Class Reference for E1039 Core & Analysis Software
GeomSvc.cxx
Go to the documentation of this file.
1 /*
2 GeomSvc.cxx
3 
4 Implementation of class GeomSvc.
5 
6 Author: Kun Liu, liuk@fnal.gov
7 Created: 10-19-2011
8 */
9 
10 #include <iostream>
11 #include <fstream>
12 #include <sstream>
13 #include <cmath>
14 #include <algorithm>
15 #include <iomanip>
16 
17 #include <TROOT.h>
18 #include <TTree.h>
19 #include <TSQLServer.h>
20 #include <TSQLResult.h>
21 #include <TSQLRow.h>
22 #include <TMath.h>
23 #include <TString.h>
24 #include <TPRegexp.h>
25 #include <TRotation.h>
26 #include <TMatrixD.h>
27 
28 #include <phool/recoConsts.h>
29 
30 #include "GeomParamPlane.h"
31 #include "GeomSvc.h"
32 
33 Plane::Plane()
34 {
35  detectorID = -1;
36  planeType = -1;
37 
38  nElements = 0;
39  x0 = 0.;
40  y0 = 0.;
41  z0 = 0.;
42 
43  x1 = 1.E6;
44  y1 = 1.E6;
45  z1 = 1.E6;
46  x2 = -1.E6;
47  y2 = -1.E6;
48  z2 = -1.E6;
49 
50  thetaX = 0.;
51  thetaY = 0.;
52  thetaZ = 0.;
53  rotX = 0.;
54  rotY = 0.;
55  rotZ = 0.;
56 
57  deltaX = 0.;
58  deltaY = 0.;
59  deltaZ = 0.;
60  deltaW = 0.;
61  for(int i = 0; i < 9; ++i) deltaW_module[i] = 0.;
62 
63  xc = 0.;
64  yc = 0.;
65  zc = 0.;
66  wc = 0.;
67  rX = 0.;
68  rY = 0.;
69  rZ = 0.;
70 
71  tmin = -1.E6;
72  tmax = 1.E6;
73 
74  rtprofile = NULL;
75  elementPos.clear();
76 }
77 
78 void Plane::update()
79 {
80  xc = x0 + deltaX;
81  yc = y0 + deltaY;
82  zc = z0 + deltaZ;
83 
84  rX = thetaX + rotX;
85  rY = thetaY + rotY;
86  rZ = thetaZ + rotZ;
87 
88  sintheta = sin(angleFromVert + rZ);
89  costheta = cos(angleFromVert + rZ);
90  tantheta = tan(angleFromVert + rZ);
91 
92  wc = getW(xc, yc);
93 
94  rotM[0][0] = cos(rZ)*cos(rY);
95  rotM[0][1] = cos(rZ)*sin(rX)*sin(rY) - cos(rX)*sin(rZ);
96  rotM[0][2] = cos(rX)*cos(rZ)*sin(rY) + sin(rX)*sin(rZ);
97  rotM[1][0] = sin(rZ)*cos(rY);
98  rotM[1][1] = sin(rZ)*sin(rX)*sin(rY) + cos(rZ)*cos(rX);
99  rotM[1][2] = sin(rZ)*sin(rY)*cos(rX) - cos(rZ)*sin(rX);
100  rotM[2][0] = -sin(rY);
101  rotM[2][1] = cos(rY)*sin(rX);
102  rotM[2][2] = cos(rY)*cos(rX);
103 
104  uVec[0] = cos(angleFromVert);
105  uVec[1] = sin(angleFromVert);
106  uVec[2] = 0.;
107  uVec = rotM*uVec;
108 
109  vVec[0] = -sin(angleFromVert);
110  vVec[1] = cos(angleFromVert);
111  vVec[2] = 0.;
112  vVec = rotM*vVec;
113 
114  xVec[0] = 1.;
115  xVec[1] = 0.;
116  xVec[2] = 0.;
117  xVec = rotM*xVec;
118 
119  yVec[0] = 0.;
120  yVec[1] = 1.;
121  yVec[2] = 0.;
122  yVec = rotM*yVec;
123 
124  nVec[0] = uVec[1]*vVec[2] - vVec[1]*uVec[2];
125  nVec[1] = uVec[2]*vVec[0] - vVec[2]*uVec[0];
126  nVec[2] = uVec[0]*vVec[1] - vVec[0]*uVec[1];
127 }
128 
129 double Plane::intercept(double tx, double ty, double x0_track, double y0_track) const
130 {
131  //double mom[3] = {tx, ty, 1.};
132  //double pos[3] = {x0_track, y0_track, 0.};
133 
134  double det = -(tx*nVec[0] + ty*nVec[1] + nVec[2]);
135  double dpos[3] = {x0_track - xc, y0_track - yc, -zc};
136 
137  double vcp[3];
138  vcp[0] = vVec[1] - vVec[2]*ty;
139  vcp[1] = vVec[2]*tx - vVec[0];
140  vcp[2] = vVec[0]*ty - vVec[1]*tx;
141 
142  //LogInfo(detectorID << " " << detectorName << " " << tx << " " << ty << " " << x0_track << " " << y0_track << " " << -(vcp[0]*dpos[0] + vcp[1]*dpos[1] + vcp[2]*dpos[2])/det);
143  return -(vcp[0]*dpos[0] + vcp[1]*dpos[1] + vcp[2]*dpos[2])/det + wc;
144 }
145 
146 double Plane::getWirePosition(int elementID) const
147 {
148  return (elementID - (nElements+1.)/2.)*spacing + xoffset + x0*costheta + y0*sintheta + deltaW;
149 }
150 
151 TVectorD Plane::getEndPoint(int elementID, int sign) const
152 {
153  //sign = -1 means start point, sign = 1 means end point, wire vector points from start to the end
154  TVectorD detCenter(3);
155  detCenter[0] = xc;
156  detCenter[1] = yc;
157  detCenter[2] = zc;
158 
159  TVectorD ep(3);
160  if(planeType != 4) //special treatment for Y planes
161  {
162  double cellLength = fabs(y2 - y1)/cos(angleFromVert);
163  double hspacing = spacing/cos(angleFromVert);
164  double hoffset = xoffset/cos(angleFromVert);
165 
166  ep = detCenter + ((elementID - (nElements+1.)/2.)*hspacing + hoffset)*xVec + 0.5*sign*cellLength*vVec;
167  }
168  else
169  {
170  double cellLength = fabs(x2 - x1)/sin(angleFromVert);
171  double vspacing = spacing/sin(angleFromVert);
172  double voffset = xoffset/sin(angleFromVert);
173 
174  ep = detCenter + ((elementID - (nElements+1.)/2.)*vspacing + voffset)*yVec + 0.5*sign*cellLength*vVec;
175  }
176 
177  return ep;
178 }
179 
180 std::ostream& operator << (std::ostream& os, const Plane& plane)
181 {
182  os << std::setw(6) << std::setiosflags(std::ios::right) << plane.detectorID
183  << std::setw(6) << std::setiosflags(std::ios::right) << plane.detectorName
184  << std::setw(6) << std::setiosflags(std::ios::right) << plane.planeType
185  << std::setw(10) << std::setiosflags(std::ios::right) << plane.spacing
186  << std::setw(10) << std::setiosflags(std::ios::right) << plane.cellWidth
187  << std::setw(10) << std::setiosflags(std::ios::right) << plane.xoffset
188  << std::setw(10) << std::setiosflags(std::ios::right) << plane.overlap
189  << std::setw(10) << std::setiosflags(std::ios::right) << plane.x2 - plane.x1
190  << std::setw(10) << std::setiosflags(std::ios::right) << plane.y2 - plane.y1
191  << std::setw(10) << std::setiosflags(std::ios::right) << plane.z2 - plane.z1
192  << std::setw(10) << std::setiosflags(std::ios::right) << plane.nElements
193  << std::setw(10) << std::setiosflags(std::ios::right) << plane.angleFromVert
194  << std::setw(10) << std::setiosflags(std::ios::right) << plane.xc
195  << std::setw(10) << std::setiosflags(std::ios::right) << plane.yc
196  << std::setw(10) << std::setiosflags(std::ios::right) << plane.zc
197  << std::setw(10) << std::setiosflags(std::ios::right) << plane.deltaW << "\n"
198  << ", nVec: {" << plane.nVec[0] << ", " << plane.nVec[1] << ", " << plane.nVec[2] << "} "
199  << ", uVec: {" << plane.uVec[0] << ", " << plane.uVec[1] << ", " << plane.uVec[2] << "} "
200  << ", vVec: {" << plane.vVec[0] << ", " << plane.vVec[1] << ", " << plane.vVec[2] << "} ";
201  if(plane.detectorID>=nChamberPlanes+nHodoPlanes+1 && plane.detectorID<=nChamberPlanes+nHodoPlanes+nPropPlanes) {
202  os << "\n";
203  for(int i=0; i<9; ++i) os << std::setw(10) << std::setiosflags(std::ios::right) << plane.deltaW_module[i];
204  }
205 
206  return os;
207 }
208 
209 GeomSvc* GeomSvc::p_geometrySvc = NULL;
210 bool GeomSvc::use_dbsvc = true;
211 
212 GeomSvc* GeomSvc::instance()
213 {
214  if(p_geometrySvc == NULL)
215  {
216  p_geometrySvc = new GeomSvc;
217  p_geometrySvc->init();
218  }
219 
220  return p_geometrySvc;
221 }
222 
223 void GeomSvc::close()
224 {
225  if(!p_geometrySvc)
226  {
227  for(int i = 0; i < nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes; ++i)
228  {
229  if(planes[i].rtprofile != NULL) delete planes[i].rtprofile;
230  }
231 
232  delete p_geometrySvc;
233  }
234  else
235  {
236  std::cout << "Error: no instance of geometry service found! " << std::endl;
237  }
238 }
239 
240 void GeomSvc::init()
241 {
242  using namespace std;
243  rc = recoConsts::instance();
244 
245  //Initialize the detectorID --- detectorName convention
246  typedef std::map<std::string, int>::value_type nameToID;
247 
248  int idx = 0;
249  map_detectorID.insert(nameToID("D0U" , ++idx));
250  map_detectorID.insert(nameToID("D0Up" , ++idx));
251  map_detectorID.insert(nameToID("D0X" , ++idx));
252  map_detectorID.insert(nameToID("D0Xp" , ++idx));
253  map_detectorID.insert(nameToID("D0V" , ++idx));
254  map_detectorID.insert(nameToID("D0Vp" , ++idx));
255  map_detectorID.insert(nameToID("D1V" , ++idx));
256  map_detectorID.insert(nameToID("D1Vp" , ++idx));
257  map_detectorID.insert(nameToID("D1X" , ++idx));
258  map_detectorID.insert(nameToID("D1Xp" , ++idx));
259  map_detectorID.insert(nameToID("D1U" , ++idx));
260  map_detectorID.insert(nameToID("D1Up" , ++idx));
261  map_detectorID.insert(nameToID("D2V" , ++idx));
262  map_detectorID.insert(nameToID("D2Vp" , ++idx));
263  map_detectorID.insert(nameToID("D2Xp" , ++idx));
264  map_detectorID.insert(nameToID("D2X" , ++idx));
265  map_detectorID.insert(nameToID("D2U" , ++idx));
266  map_detectorID.insert(nameToID("D2Up" , ++idx));
267  map_detectorID.insert(nameToID("D3pVp", ++idx));
268  map_detectorID.insert(nameToID("D3pV" , ++idx));
269  map_detectorID.insert(nameToID("D3pXp", ++idx));
270  map_detectorID.insert(nameToID("D3pX" , ++idx));
271  map_detectorID.insert(nameToID("D3pUp", ++idx));
272  map_detectorID.insert(nameToID("D3pU" , ++idx));
273  map_detectorID.insert(nameToID("D3mVp", ++idx));
274  map_detectorID.insert(nameToID("D3mV" , ++idx));
275  map_detectorID.insert(nameToID("D3mXp", ++idx));
276  map_detectorID.insert(nameToID("D3mX" , ++idx));
277  map_detectorID.insert(nameToID("D3mUp", ++idx));
278  map_detectorID.insert(nameToID("D3mU" , ++idx));
279  // Please make sure of "idx = nChamberPlanes" here
280 
281  map_detectorID.insert(nameToID("H1B" , ++idx));
282  map_detectorID.insert(nameToID("H1T" , ++idx));
283  map_detectorID.insert(nameToID("H1L" , ++idx));
284  map_detectorID.insert(nameToID("H1R" , ++idx));
285  map_detectorID.insert(nameToID("H2L" , ++idx));
286  map_detectorID.insert(nameToID("H2R" , ++idx));
287  map_detectorID.insert(nameToID("H2B" , ++idx));
288  map_detectorID.insert(nameToID("H2T" , ++idx));
289  map_detectorID.insert(nameToID("H3B" , ++idx));
290  map_detectorID.insert(nameToID("H3T" , ++idx));
291  map_detectorID.insert(nameToID("H4Y1L", ++idx));
292  map_detectorID.insert(nameToID("H4Y1R", ++idx));
293  map_detectorID.insert(nameToID("H4Y2L", ++idx));
294  map_detectorID.insert(nameToID("H4Y2R", ++idx));
295  map_detectorID.insert(nameToID("H4B" , ++idx));
296  map_detectorID.insert(nameToID("H4T" , ++idx));
297  // Please make sure of "idx = nChamberPlanes + nHodoPlanes" here
298 
299  map_detectorID.insert(nameToID("P1Y1", ++idx)); // 47
300  map_detectorID.insert(nameToID("P1Y2", ++idx));
301  map_detectorID.insert(nameToID("P1X1", ++idx));
302  map_detectorID.insert(nameToID("P1X2", ++idx));
303  map_detectorID.insert(nameToID("P2X1", ++idx));
304  map_detectorID.insert(nameToID("P2X2", ++idx));
305  map_detectorID.insert(nameToID("P2Y1", ++idx));
306  map_detectorID.insert(nameToID("P2Y2", ++idx));
307  // Please make sure of "idx = nChamberPlanes + nHodoPlanes + nPropPlanes" here
308 
309  map_detectorID.insert(nameToID("DP1TL", ++idx)); // 55
310  map_detectorID.insert(nameToID("DP1TR", ++idx));
311  map_detectorID.insert(nameToID("DP1BL", ++idx));
312  map_detectorID.insert(nameToID("DP1BR", ++idx));
313  map_detectorID.insert(nameToID("DP2TL", ++idx));
314  map_detectorID.insert(nameToID("DP2TR", ++idx));
315  map_detectorID.insert(nameToID("DP2BL", ++idx));
316  map_detectorID.insert(nameToID("DP2BR", ++idx));
317  // Please make sure of "idx = nChamberPlanes + nHodoPlanes + nPropPlanes + nDarkPhotonPlanes" here
318 
319  map_detectorID.insert(nameToID("BeforeInhNIM" , ++idx)); // 63
320  map_detectorID.insert(nameToID("BeforeInhMatrix", ++idx));
321  map_detectorID.insert(nameToID("AfterInhNIM" , ++idx));
322  map_detectorID.insert(nameToID("AfterInhMatrix" , ++idx));
323  map_detectorID.insert(nameToID("BOS" , ++idx));
324  map_detectorID.insert(nameToID("EOS" , ++idx));
325  map_detectorID.insert(nameToID("L1" , ++idx));
326  map_detectorID.insert(nameToID("RF" , ++idx));
327  map_detectorID.insert(nameToID("STOP" , ++idx));
328  map_detectorID.insert(nameToID("L1PXtp" , ++idx));
329  map_detectorID.insert(nameToID("L1PXtn" , ++idx));
330  map_detectorID.insert(nameToID("L1PXbp" , ++idx));
331  map_detectorID.insert(nameToID("L1PXbn" , ++idx));
332  map_detectorID.insert(nameToID("L1NIMxt" , ++idx));
333  map_detectorID.insert(nameToID("L1NIMxb" , ++idx));
334  map_detectorID.insert(nameToID("L1NIMyt" , ++idx));
335  map_detectorID.insert(nameToID("L1NIMyb" , ++idx));
336  // No constant for this name group is defined in GlobalConsts.h.
337 
338  map_detectorID.insert(nameToID("H4Y1Ll", ++idx)); // 80
339  map_detectorID.insert(nameToID("H4Y1Lr", ++idx));
340  map_detectorID.insert(nameToID("H4Y1Rl", ++idx));
341  map_detectorID.insert(nameToID("H4Y1Rr", ++idx));
342  map_detectorID.insert(nameToID("H4Y2Ll", ++idx));
343  map_detectorID.insert(nameToID("H4Y2Lr", ++idx));
344  map_detectorID.insert(nameToID("H4Y2Rl", ++idx));
345  map_detectorID.insert(nameToID("H4Y2Rr", ++idx));
346  map_detectorID.insert(nameToID("H4Tu" , ++idx));
347  map_detectorID.insert(nameToID("H4Td" , ++idx));
348  map_detectorID.insert(nameToID("H4Bu" , ++idx));
349  map_detectorID.insert(nameToID("H4Bd" , ++idx));
350  // No constant for this name group is defined in GlobalConsts.h.
351 
352  // TODO temp solution
353  for(int i=1; i<=nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes; ++i) {
354  map_detid_triggerlv[i] = -1;
355  }
356  for(int i=nChamberPlanes+1; i<=nChamberPlanes+4; ++i) {
357  map_detid_triggerlv[i] = 0;
358  }
359  for(int i=nChamberPlanes+5; i<=nChamberPlanes+8; ++i) {
360  map_detid_triggerlv[i] = 1;
361  }
362  for(int i=nChamberPlanes+9; i<=nChamberPlanes+10; ++i) {
363  map_detid_triggerlv[i] = 2;
364  }
365  for(int i=nChamberPlanes+11; i<=nChamberPlanes+16; ++i) {
366  map_detid_triggerlv[i] = 3;
367  }
368 
369  typedef std::map<int, std::string>::value_type idToName;
370  for(std::map<std::string, int>::iterator iter = map_detectorID.begin(); iter != map_detectorID.end(); ++iter)
371  {
372  map_detectorName.insert(idToName(iter->second, iter->first));
373  }
374 
375  //----------------------- hard-coded part is over-----------------------
376  if (use_dbsvc) initPlaneDbSvc();
377  else initPlaneDirect();
378 
379  // Alignment and calibration
380  calibration_loaded = false;
381  if(rc->get_BoolFlag("OnlineAlignment"))
382  {
383  loadOnlineAlignment(rc->get_CharFlag("AlignmentMille"));
384  }
385  else if(!rc->get_BoolFlag("IdealGeom"))
386  {
387  loadAlignment("NULL", rc->get_CharFlag("AlignmentHodo"), rc->get_CharFlag("AlignmentProp"));
388  loadMilleAlignment(rc->get_CharFlag("AlignmentMille"));
389  loadCalibration(rc->get_CharFlag("Calibration"));
390  }
391 
392  initWireLUT();
393 
395  xmin_kmag = -57.*2.54;
396  xmax_kmag = 57.*2.54;
397  ymin_kmag = -40.*2.54;
398  ymax_kmag = 40.*2.54;
399 
400  zmin_kmag = 1064.26 - 120.*2.54;
401  zmax_kmag = 1064.26 + 120.*2.54;
402 
403 #ifdef _DEBUG_ON
404  for(int i = 1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes; ++i)
405  {
406  cout << planes[i] << endl;
407  }
408 #endif
409 }
410 
411 void GeomSvc::initPlaneDirect() {
412  std::cerr << "!! ERROR !!\n"
413  << "!! GeomSvc::initPlaneDirect() is called but no longer available.\n"
414  << "!! Probably you have 'GeomSvc::UseDbSvc(false)' in your macro.\n"
415  << "!! Please delete it. Abort for now." << std::endl;
416  exit(1);
417 }
418 
419 void GeomSvc::initPlaneDbSvc() {
420  using namespace std;
421  recoConsts* rc = recoConsts::instance();
422  GeomParamPlane* geom = new GeomParamPlane();
423 
424  if (rc->FlagExist("GeomPlane") && rc->get_CharFlag("GeomPlane") != "") {
425  string map_id = rc->get_CharFlag("GeomPlane");
426  cout << "GeomSvc: Load the plane geometry info for map_id = " << map_id << "." << endl;
427  geom->SetMapID(map_id);
428  } else {
429  int run_id = rc->get_IntFlag("RUNNUMBER");
430  cout << "GeomSvc: Load the plane geometry info for run_id = " << run_id << "." << endl;
431  geom->SetMapIDbyDB(run_id);
432  rc->set_CharFlag("GeomPlane", geom->GetMapID());
433  }
434  geom->ReadFromDB();
435 
436  int dummy = 0;
437  for (int ii = 0; ii < geom->GetNumPlanes(); ii++) {
438  GeomParamPlane::Plane* pl = geom->GetPlane(ii);
439 
440  string detectorName(pl->det_name);
441  toLocalDetectorName(detectorName, dummy);
442 
443  int detectorID = map_detectorID[detectorName];
444  //cout << "ii: " << ii << " " << detectorID << " " << detectorName << endl;
445  planes[detectorID].detectorID = detectorID;
446  planes[detectorID].detectorName = detectorName;
447  planes[detectorID].spacing = pl->cell_spacing;
448  planes[detectorID].cellWidth = pl->cell_width;
449  planes[detectorID].overlap = pl->cell_width - pl->cell_spacing;
450  planes[detectorID].angleFromVert = pl->angle_from_vert;
451  planes[detectorID].xoffset = pl->xoffset;
452  planes[detectorID].thetaX = pl->theta_x;
453  planes[detectorID].thetaY = pl->theta_y;
454  planes[detectorID].thetaZ = pl->theta_z;
455 
456  //Following items need to be sumed or averaged over all modules
457  planes[detectorID].nElements += pl->n_ele; // still needed in E1039??
458  double x0_i = pl->x0;
459  double y0_i = pl->y0;
460  double z0_i = pl->z0;
461  double width_i = pl->width;
462  double height_i = pl->height;
463 
464  double x1_i = x0_i - 0.5*width_i;
465  double x2_i = x0_i + 0.5*width_i;
466  double y1_i = y0_i - 0.5*height_i;
467  double y2_i = y0_i + 0.5*height_i;
468 
469  planes[detectorID].x0 += x0_i;
470  planes[detectorID].y0 += y0_i;
471  planes[detectorID].z0 += z0_i;
472  if(planes[detectorID].x1 > x1_i) planes[detectorID].x1 = x1_i;
473  if(planes[detectorID].x2 < x2_i) planes[detectorID].x2 = x2_i;
474  if(planes[detectorID].y1 > y1_i) planes[detectorID].y1 = y1_i;
475  if(planes[detectorID].y2 < y2_i) planes[detectorID].y2 = y2_i;
476 
477  //Calculated value
478  planes[detectorID].resolution = planes[detectorID].cellWidth/sqrt(12.);
479  planes[detectorID].update();
480 
481  //Set the plane type
482  if(detectorName.find("X") != string::npos || detectorName.find("T") != string::npos || detectorName.find("B") != string::npos)
483  {
484  planes[detectorID].planeType = 1;
485  }
486  else if((detectorName.find("U") != string::npos || detectorName.find("V") != string::npos) && planes[detectorID].angleFromVert > 0)
487  {
488  planes[detectorID].planeType = 2;
489  }
490  else if((detectorName.find("U") != string::npos || detectorName.find("V") != string::npos) && planes[detectorID].angleFromVert < 0)
491  {
492  planes[detectorID].planeType = 3;
493  }
494  else if(detectorName.find("Y") != string::npos || detectorName.find("L") != string::npos || detectorName.find("R") != string::npos)
495  {
496  planes[detectorID].planeType = 4;
497  }
498  }
499 
500  //For prop. tube only, average over 9 modules
501  for(int i = nChamberPlanes+nHodoPlanes+1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes; i++)
502  {
503  planes[i].x0 = planes[i].x0/9.;
504  planes[i].y0 = planes[i].y0/9.;
505  planes[i].z0 = planes[i].z0/9.;
506 
507  planes[i].update();
508  }
509 
510  //Set the depth of the plane for Geant4 simulation simplification
511  for(int i = 1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes; ++i)
512  {
513  if(i <= nChamberPlanes) //For drift chambers the depth equals the wire spacing
514  {
515  planes[i].z1 = planes[i].z0 - planes[i].cellWidth/2.;
516  planes[i].z2 = planes[i].z0 + planes[i].cellWidth/2.;
517  }
518  else if(i <= nChamberPlanes+nHodoPlanes) //For hodos the depth is hard-coded
519  {
520  planes[i].z1 = planes[i].z0 - 0.635/2.;
521  planes[i].z2 = planes[i].z0 + 0.635/2.;
522  }
523  else if(i <= nChamberPlanes+nHodoPlanes+nPropPlanes) //For prop. tubes the depth is defined by the equivalent gas volume
524  {
525  planes[i].z1 = planes[i].z0 - planes[i].cellWidth*TMath::Pi()/8.;
526  planes[i].z2 = planes[i].z0 + planes[i].cellWidth*TMath::Pi()/8.;
527  }
528  else //For dark photon hodos the depth equals the cell width
529  {
530  planes[i].z1 = planes[i].z0 - planes[i].cellWidth/2.;
531  planes[i].z2 = planes[i].z0 + planes[i].cellWidth/2.;
532  }
533  }
534 }
535 
536 
537 void GeomSvc::initWireLUT() {
538 
540  typedef std::unordered_map<int, double>::value_type posType;
541  typedef std::unordered_map<int, TVectorD>::value_type epType;
542  for(int i = 1; i <= nChamberPlanes; ++i)
543  {
544  map_wirePosition[i].clear();
545  map_endPoint1 [i].clear();
546  map_endPoint2 [i].clear();
547  planes [i].elementPos.clear();
548  for(int j = 1; j <= planes[i].nElements; ++j)
549  {
550  double pos = (j - (planes[i].nElements+1.)/2.)*planes[i].spacing + planes[i].xoffset + planes[i].x0*planes[i].costheta + planes[i].y0*planes[i].sintheta + planes[i].deltaW;
551  map_wirePosition[i].insert(posType(j, pos));
552 
553  map_endPoint1[i].insert(epType(j, planes[i].getEndPoint(j, -1)));
554  map_endPoint2[i].insert(epType(j, planes[i].getEndPoint(j, 1)));
555 
556  planes[i].elementPos.push_back(pos);
557  }
558  std::sort(planes[i].elementPos.begin(), planes[i].elementPos.end());
559  }
560 
561  // 2. for hodoscopes and prop. tubes
562  for(int i = nChamberPlanes + 1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes; ++i)
563  {
564  map_wirePosition[i].clear();
565  map_endPoint1 [i].clear();
566  map_endPoint2 [i].clear();
567  planes [i].elementPos.clear();
568  for(int j = 1; j <= planes[i].nElements; ++j)
569  {
570  double pos;
571  if(i <= nChamberPlanes+nHodoPlanes)
572  {
573  pos = planes[i].x0*planes[i].costheta + planes[i].y0*planes[i].sintheta + planes[i].xoffset + (j - (planes[i].nElements+1)/2.)*planes[i].spacing + planes[i].deltaW;
574  }
575  else if(i <= nChamberPlanes+nHodoPlanes+nPropPlanes)
576  {
577  int moduleID = 8 - int((j - 1)/8); //Need to re-define moduleID for run2, note it's reversed compared to elementID
578  pos = planes[i].x0*planes[i].costheta + planes[i].y0*planes[i].sintheta + planes[i].xoffset + (j - (planes[i].nElements+1)/2.)*planes[i].spacing + planes[i].deltaW_module[moduleID];
579  }
580  else
581  {
582  int elementID = ((i-55) & 2) > 0 ? planes[i].nElements + 1 - j : j;
583  pos = planes[i].x0*planes[i].costheta + planes[i].y0*planes[i].sintheta + planes[i].xoffset + (elementID - (planes[i].nElements+1)/2.)*planes[i].spacing + planes[i].deltaW;
584  }
585  map_wirePosition[i].insert(posType(j, pos));
586  map_endPoint1[i].insert(epType(j, planes[i].getEndPoint(j, -1)));
587  map_endPoint2[i].insert(epType(j, planes[i].getEndPoint(j, 1)));
588  planes[i].elementPos.push_back(pos);
589  }
590  std::sort(planes[i].elementPos.begin(), planes[i].elementPos.end());
591  }
592 }
593 
594 std::vector<int> GeomSvc::getDetectorIDs(std::string pattern)
595 {
596  TPRegexp pattern_re(pattern.c_str());
597 
598  std::vector<int> detectorIDs;
599  detectorIDs.clear();
600 
601  for(std::map<std::string, int>::iterator iter = map_detectorID.begin(); iter != map_detectorID.end(); ++iter)
602  {
603  TString detectorName((*iter).first.c_str());
604  if(detectorName(pattern_re) != "")
605  {
606  detectorIDs.push_back((*iter).second);
607  }
608  }
609  std::sort(detectorIDs.begin(), detectorIDs.end());
610 
611  return detectorIDs;
612 }
613 
614 bool GeomSvc::findPatternInDetector(int detectorID, std::string pattern)
615 {
616  TPRegexp pattern_re(pattern.c_str());
617  TString detectorName(map_detectorName[detectorID]);
618 
619  return detectorName(pattern_re) != "";
620 }
621 
622 bool GeomSvc::isInPlane(int detectorID, double x, double y)
623 {
624  if(x < planes[detectorID].x1 || x > planes[detectorID].x2) return false;
625  if(y < planes[detectorID].y1 || y > planes[detectorID].y2) return false;
626 
627  return true;
628 }
629 
630 bool GeomSvc::isInElement(int detectorID, int elementID, double x, double y, double tolr)
631 {
632  double x_min, x_max, y_min, y_max;
633  get2DBoxSize(detectorID, elementID, x_min, x_max, y_min, y_max);
634 
635  x_min -= (tolr*(x_max - x_min));
636  x_max += (tolr*(x_max - x_min));
637  y_min -= (tolr*(y_max - y_min));
638  y_max += (tolr*(y_max - y_min));
639 
640  return x > x_min && x < x_max && y > y_min && y < y_max;
641 }
642 
643 bool GeomSvc::isInKMAG(double x, double y)
644 {
645  if(x < xmin_kmag || x > xmax_kmag) return false;
646  if(y < ymin_kmag || y > ymax_kmag) return false;
647 
648  return true;
649 }
650 
651 void GeomSvc::getMeasurement(int detectorID, int elementID, double& measurement, double& dmeasurement)
652 {
653  measurement = map_wirePosition[detectorID][elementID];
654  dmeasurement = planes[detectorID].resolution;
655 }
656 
657 double GeomSvc::getMeasurement(int detectorID, int elementID)
658 {
659  return map_wirePosition[detectorID][elementID];
660 }
661 
662 void GeomSvc::getEndPoints(int detectorID, int elementID, TVectorD& ep1, TVectorD& ep2)
663 {
664  ep1 = map_endPoint1[detectorID][elementID];
665  ep2 = map_endPoint2[detectorID][elementID];
666 }
667 
668 void GeomSvc::getEndPoints(int detectorID, int elementID, TVector3& ep1, TVector3& ep2)
669 {
670  TVectorD vp1(map_endPoint1[detectorID][elementID]);
671  TVectorD vp2(map_endPoint2[detectorID][elementID]);
672 
673  ep1.SetXYZ(vp1[0], vp1[1], vp1[2]);
674  ep2.SetXYZ(vp2[0], vp2[1], vp2[2]);
675 }
676 
677 int GeomSvc::getExpElementID(int detectorID, double pos_exp)
678 {
679  //check if pos_exp is within the correct range
680  double pos_min = planes[detectorID].elementPos.front();
681  double pos_max = planes[detectorID].elementPos.back();
682  pos_min -= (0.5*planes[detectorID].cellWidth);
683  pos_max += (0.5*planes[detectorID].cellWidth);
684 
685  if(pos_exp > pos_max) return planes[detectorID].nElements+1;
686  if(pos_exp < pos_min) return 0;
687 
688  int index = std::lower_bound(planes[detectorID].elementPos.begin(), planes[detectorID].elementPos.end(), pos_exp) - planes[detectorID].elementPos.begin();
689  int elementID = index + 1 + (planes[detectorID].elementPos[index] - pos_exp > 0.5*planes[detectorID].spacing ? -1 : 0);
690 
691  if(detectorID > nChamberPlanes+nHodoPlanes+nPropPlanes)
692  {
693  bool bottom = ((detectorID-55) & 2) > 0;
694  if(bottom) elementID = planes[detectorID].nElements+1-elementID;
695  }
696 
697  return elementID;
698 }
699 
700 void GeomSvc::get2DBoxSize(int detectorID, int elementID, double& x_min, double& x_max, double& y_min, double& y_max)
701 {
702  std::string detectorName = getDetectorName(detectorID);
703  if(planes[detectorID].planeType == 1)
704  {
705  double x_center = map_wirePosition[detectorID][elementID];
706  double x_width = 0.5*planes[detectorID].cellWidth;
707  x_min = x_center - x_width;
708  x_max = x_center + x_width;
709 
710  y_min = planes[detectorID].y1;
711  y_max = planes[detectorID].y2;
712  }
713  else
714  {
715  double y_center = map_wirePosition[detectorID][elementID];
716  double y_width = 0.5*planes[detectorID].cellWidth;
717  y_min = y_center - y_width;
718  y_max = y_center + y_width;
719 
720  x_min = planes[detectorID].x1;
721  x_max = planes[detectorID].x2;
722  }
723 }
724 
725 void GeomSvc::getWireEndPoints(int detectorID, int elementID, double& x_min, double& x_max, double& y_min, double& y_max)
726 {
727  y_min = planes[detectorID].y1;
728  y_max = planes[detectorID].y2;
729 
730  double dw = planes[detectorID].xoffset + planes[detectorID].spacing*(elementID - (planes[detectorID].nElements + 1.)/2.);
731  x_min = planes[detectorID].xc + dw*cos(planes[detectorID].rZ) - fabs(0.5*planes[detectorID].tantheta*(y_max - y_min));
732  x_max = planes[detectorID].xc + dw*cos(planes[detectorID].rZ) + fabs(0.5*planes[detectorID].tantheta*(y_max - y_min));
733 }
734 
740 void GeomSvc::toLocalDetectorName(std::string& detectorName, int& eID)
741 {
742  using namespace std;
743 
744  if(detectorName[0] == 'P' && (detectorName[2] == 'H' || detectorName[2] == 'V'))
745  {
746  string XY = detectorName[2] == 'H' ? "Y" : "X";
747  string FB = (detectorName[3] == 'f' || detectorName[4] == 'f') ? "1" : "2"; //temporary solution
748  int moduleID = std::isdigit(detectorName[3]) == 1 ? atoi(&detectorName[3]) : atoi(&detectorName[4]);
749 
750  detectorName.replace(2, detectorName.length(), "");
751  detectorName += XY;
752  detectorName += FB;
753 
754  if(eID <= 8) //Means either it's at lowest module or elementID is only from 1 to 8, repeatedly
755  {
756  eID = (9 - moduleID)*8 + eID;
757  }
758  }
759  //else if(detectorName.substr(0, 2) == "H4")
760  //{
761  // if(detectorName.find("T") != string::npos || detectorName.find("B") != string::npos)
762  // {
763  // detectorName.replace(3, detectorName.length(), "");
764  // }
765  // else
766  // {
767  // detectorName.replace(5, detectorName.length(), "");
768  // }
769  //}
770 }
771 
772 bool GeomSvc::isChamber(const int detectorID) const
773 {
774  return isChamber(getDetectorName(detectorID));
775 }
776 
777 bool GeomSvc::isChamber(const std::string detectorName) const
778 {
779  return detectorName[0] == 'D' && '0' <= detectorName[1] && detectorName[1] <= '3';
780 }
781 
782 bool GeomSvc::isHodo(const int detectorID) const
783 {
784  return isHodo(getDetectorName(detectorID));
785 }
786 
787 bool GeomSvc::isHodo(const std::string detectorName) const
788 {
789  return detectorName[0] == 'H' && '1' <= detectorName[1] && detectorName[1] <= '4';
790 }
791 
792 bool GeomSvc::isPropTube(const int detectorID) const
793 {
794  return isPropTube(getDetectorName(detectorID));
795 }
796 
797 bool GeomSvc::isPropTube(const std::string detectorName) const
798 {
799  return detectorName[0] == 'P' && '1' <= detectorName[1] && detectorName[1] <= '2';
800 }
801 
802 bool GeomSvc::isDPHodo(const int detectorID) const
803 {
804  return isDPHodo(getDetectorName(detectorID));
805 }
806 
807 bool GeomSvc::isDPHodo(const std::string detectorName) const
808 {
809  return detectorName.substr(0, 2) == "DP";
810 }
811 
812 bool GeomSvc::isInh(const int detectorID) const
813 {
814  return isInh(getDetectorName(detectorID));
815 }
816 
817 bool GeomSvc::isInh(const std::string detectorName) const
818 {
819  return detectorName.substr(0, 6) == "Before" || detectorName.substr(0, 5) == "After";
820 }
821 
822 int GeomSvc::getHodoStation(const int detectorID) const
823 {
824  return getHodoStation(getDetectorName(detectorID));
825 }
826 
827 int GeomSvc::getHodoStation(const std::string detectorName) const
828 {
829  if (detectorName.size() == 0 || detectorName[0] != 'H') return 0;
830  int num = detectorName[1] - '0';
831  return 1 <= num && num <= 4 ? num : 0;
832 }
833 
834 double GeomSvc::getDriftDistance(int detectorID, double tdcTime)
835 {
836  if(!calibration_loaded)
837  {
838  return 0.;
839  }
840  else if(planes[detectorID].rtprofile == NULL)
841  {
842  return 0.;
843  }
844  else if(tdcTime < planes[detectorID].tmin)
845  {
846  return 0.5*planes[detectorID].cellWidth;
847  }
848  else if(tdcTime > planes[detectorID].tmax)
849  {
850  return 0.;
851  }
852  else
853  {
854  return planes[detectorID].rtprofile->Eval(tdcTime);
855  }
856 
857  return 0.;
858 }
859 
860 double GeomSvc::getInterceptionFast(int detectorID, double tx, double ty, double x0, double y0) const
861 {
862  return (tx*planes[detectorID].zc + x0)*planes[detectorID].costheta + (ty*planes[detectorID].zc + y0)*planes[detectorID].sintheta;
863 }
864 
865 double GeomSvc::getDCA(int detectorID, int elementID, double tx, double ty, double x0, double y0)
866 {
867  TVector3 trkp0(x0, y0, 0.);
868  TVector3 trkdir(tx, ty, 1.);
869 
870  TVector3 ep1, ep2;
871  getEndPoints(detectorID, elementID, ep1, ep2);
872  TVector3 wiredir = ep2 - ep1;
873 
874  TVector3 norm = trkdir.Cross(wiredir);
875  norm.SetMag(1.);
876  return (ep1 - trkp0).Dot(norm);
877 }
878 
880 
885 void GeomSvc::loadOnlineAlignment(const std::string& alignmentFile_mille)
886 {
887  loadMilleAlignment(alignmentFile_mille); //final chance of overwrite resolution numbers in online mode
888 
889  //TSQLServer* con = TSQLServer::Connect(rc->get_CharFlag("MySQLURL").c_str(), "seaguest","qqbar2mu+mu-");
890  //char query[300];
891  //const char* buf_offsets = "SELECT detectorName,deltaX,deltaY,deltaZ,rotateAboutZ FROM %s.PlaneOffsets WHERE"
892  // " detectorName LIKE 'D%%' OR detectorName LIKE 'H__' OR detectorName LIKE 'H____' OR detectorName LIKE 'P____'";
893  //sprintf(query, buf_offsets, rc->get_CharFlag("Geometry").c_str());
894  //TSQLResult* res = con->Query(query);
895  //unsigned int nRows = res->GetRowCount();
896  //if(nRows >= nChamberPlanes) cout << "GeomSvc: loaded chamber alignment parameters from database: " << rc->get_CharFlag("Geometry") << endl;
897  //for(unsigned int i = 0; i < nRows; ++i)
898  //{
899  // TSQLRow* row = res->Next();
900  // string detectorName(row->GetField(0));
901  // toLocalDetectorName(detectorName, dummy);
902  //
903  // int detectorID = map_detectorID[detectorName];
904  // if(detectorID > nChamberPlanes) // ?? WTF?
905  // {
906  // delete row;
907  // continue;
908  // }
909  //
910  // planes[detectorID].deltaX = atof(row->GetField(1));
911  // planes[detectorID].deltaY = atof(row->GetField(2));
912  // planes[detectorID].deltaZ = atof(row->GetField(3));
913  // planes[detectorID].rotZ = atof(row->GetField(4));
914  //
915  // planes[detectorID].update();
916  // planes[detectorID].deltaW = planes[detectorID].getW(planes[detectorID].deltaX, planes[detectorID].deltaY);
917  //
918  // delete row;
919  //}
920  //delete res;
921  //delete con;
922 }
923 
924 void GeomSvc::loadAlignment(const std::string& alignmentFile_chamber, const std::string& alignmentFile_hodo, const std::string& alignmentFile_prop)
925 {
926  using namespace std;
927 
928  //load alignment numbers for chambers
929  fstream _align_chamber;
930  _align_chamber.open(alignmentFile_chamber.c_str(), ios::in);
931 
932  char buf[300];
933  if(_align_chamber)
934  {
935  for(int i = 1; i <= nChamberPlanes; i++)
936  {
937  _align_chamber.getline(buf, 100);
938  istringstream stringBuf(buf);
939 
940  stringBuf >> planes[i].deltaW >> planes[i].resolution;
941  //if(planes[i].resolution < RESOLUTION_DC) planes[i].resolution = RESOLUTION_DC;
942 
943  planes[i].deltaX = planes[i].deltaW*planes[i].costheta;
944  planes[i].deltaY = planes[i].deltaW*planes[i].sintheta;
945  planes[i].update();
946  }
947 
948  for(int i = 1; i <= nChamberPlanes; i += 2)
949  {
950  double resol = planes[i].resolution > planes[i+1].resolution ? planes[i].resolution : planes[i+1].resolution;
951  planes[i ].resolution = resol;
952  planes[i+1].resolution = resol;
953  }
954 
955  cout << "GeomSvc: loaded chamber alignment parameters from " << alignmentFile_chamber << endl;
956  }
957  _align_chamber.close();
958 
959  if (! alignmentFile_hodo.empty()) {
960  //load alignment numbers for hodos
961  fstream _align_hodo;
962  _align_hodo.open(alignmentFile_hodo.c_str(), ios::in);
963 
964  if(_align_hodo)
965  {
966  for(int i = nChamberPlanes+1; i <= nChamberPlanes+nHodoPlanes; i++)
967  {
968  _align_hodo.getline(buf, 100);
969  istringstream stringBuf(buf);
970 
971  stringBuf >> planes[i].deltaW;
972  planes[i].deltaX = planes[i].deltaW*planes[i].costheta;
973  planes[i].deltaY = planes[i].deltaW*planes[i].sintheta;
974  planes[i].update();
975  }
976  cout << "GeomSvc: loaded hodoscope alignment parameters from " << alignmentFile_hodo << endl;
977  }
978  else
979  {
980  cout << "GeomSvc: failed to load hodoscope alignment parameters from " << alignmentFile_hodo << endl;
981  }
982  _align_hodo.close();
983  }
984 
985  if (! alignmentFile_prop.empty()) {
986  //load alignment numbers for prop. tubes
987  fstream _align_prop;
988  _align_prop.open(alignmentFile_prop.c_str(), ios::in);
989 
990  if(_align_prop)
991  {
992  for(int i = nChamberPlanes+nHodoPlanes+1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes; i += 2)
993  {
994  planes[i].deltaW = 0.;
995  planes[i+1].deltaW = 0.;
996  for(int j = 0; j < 9; j++)
997  {
998  _align_prop.getline(buf, 100);
999  istringstream stringBuf(buf);
1000 
1001  stringBuf >> planes[i].deltaW_module[j];
1002  planes[i+1].deltaW_module[j] = planes[i].deltaW_module[j];
1003 
1004  planes[i].deltaW += planes[i].deltaW_module[j];
1005  planes[i+1].deltaW += planes[i+1].deltaW_module[j];
1006  }
1007 
1008  planes[i].deltaW /= 9.;
1009  planes[i].deltaX = planes[i].deltaW*planes[i].costheta;
1010  planes[i].deltaY = planes[i].deltaW*planes[i].sintheta;
1011 
1012  planes[i+1].deltaW /= 9.;
1013  planes[i+1].deltaX = planes[i+1].deltaW*planes[i+1].costheta;
1014  planes[i+1].deltaY = planes[i+1].deltaW*planes[i+1].sintheta;
1015  }
1016  cout << "GeomSvc: loaded prop. tube alignment parameters from " << alignmentFile_prop << endl;
1017  }
1018  else
1019  {
1020  cout << "GeomSvc: failed to load prop. tube alignment parameters from " << alignmentFile_prop << endl;
1021  }
1022  }
1023 }
1024 
1025 void GeomSvc::loadMilleAlignment(const std::string& alignmentFile_mille)
1026 {
1027  using namespace std;
1028  if (alignmentFile_mille.empty()) return;
1029 
1030  //load alignment numbers for chambers
1031  fstream _align_mille;
1032  _align_mille.open(alignmentFile_mille.c_str(), ios::in);
1033 
1034  char buf[300];
1035  if(_align_mille)
1036  {
1037  for(int i = 1; i <= nChamberPlanes; i++)
1038  {
1039  _align_mille.getline(buf, 100);
1040  istringstream stringBuf(buf);
1041 
1042  stringBuf >> planes[i].deltaZ >> planes[i].rotZ >> planes[i].deltaW >> planes[i].resolution;
1043  planes[i].deltaX = planes[i].deltaW*planes[i].costheta;
1044  planes[i].deltaY = planes[i].deltaW*planes[i].sintheta;
1045  planes[i].update();
1046 
1047  //if(planes[i].resolution < RESOLUTION_DC) planes[i].resolution = RESOLUTION_DC;
1048  }
1049  cout << "GeomSvc: loaded millepede-based alignment parameters from " << alignmentFile_mille << endl;
1050 
1051  for(int i = 1; i <= nChamberPlanes; i+=2)
1052  {
1053  planes[i].resolution = rc->get_DoubleFlag("RESOLUTION_FACTOR")*0.5*(planes[i].resolution + planes[i+1].resolution);
1054  planes[i+1].resolution = planes[i].resolution;
1055  }
1056  }
1057  else
1058  {
1059  cout << "GeomSvc: failed to load mullepede-based alignment parameters from " << alignmentFile_mille << endl;
1060  }
1061 
1062  _align_mille.close();
1063 }
1064 
1065 void GeomSvc::loadCalibration(const std::string& calibrationFile)
1066 {
1067  using namespace std;
1068  if (calibrationFile.empty()) return;
1069 
1070  fstream _cali_file;
1071  _cali_file.open(calibrationFile.c_str(), ios::in);
1072 
1073  char buf[300];
1074  int iBin, nBin, detectorID;
1075  double tmin_temp, tmax_temp;
1076  double R[500], T[500];
1077  if(_cali_file)
1078  {
1079  calibration_loaded = true;
1080 
1081  while(_cali_file.getline(buf, 100))
1082  {
1083  istringstream detector_info(buf);
1084  detector_info >> detectorID >> nBin >> tmin_temp >> tmax_temp;
1085  planes[detectorID].tmax = tmax_temp;
1086  planes[detectorID].tmin = tmin_temp;
1087 
1088  for(int i = 0; i < nBin; i++)
1089  {
1090  _cali_file.getline(buf, 100);
1091  istringstream cali_line(buf);
1092 
1093  cali_line >> iBin >> T[i] >> R[i];
1094  }
1095 
1096  if(planes[detectorID].rtprofile != NULL) delete planes[detectorID].rtprofile;
1097  if(nBin > 0) planes[detectorID].rtprofile = new TSpline3(getDetectorName(detectorID).c_str(), T, R, nBin, "b1e1");
1098  }
1099  cout << "GeomSvc: loaded calibration parameters from " << calibrationFile << endl;
1100  }
1101  _cali_file.close();
1102 }
1103 
1104 bool GeomSvc::isInTime(int detectorID, double tdcTime)
1105 {
1106  return tdcTime > planes[detectorID].tmin && tdcTime < planes[detectorID].tmax;
1107 }
1108 
1109 void GeomSvc::printWirePosition()
1110 {
1111  for(std::map<std::string, int>::iterator iter = map_detectorID.begin(); iter != map_detectorID.end(); ++iter)
1112  {
1113  if(iter->second>nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes) continue;
1114  int detectorID = (*iter).second;
1115  std::cout << " ====================== " << (*iter).first << " ==================== " << std::endl;
1116  for(int i = 1; i <= planes[detectorID].nElements; ++i)
1117  {
1118  std::cout << std::setw(6) << std::setiosflags(std::ios::right) << detectorID;
1119  std::cout << std::setw(6) << std::setiosflags(std::ios::right) << (*iter).first;
1120  std::cout << std::setw(6) << std::setiosflags(std::ios::right) << i;
1121  std::cout << std::setw(10) << std::setiosflags(std::ios::right) << map_wirePosition[detectorID][i];
1122  std::cout << std::setw(10) << std::setiosflags(std::ios::right) << map_wirePosition[detectorID][i] - 0.5*planes[detectorID].cellWidth;
1123  std::cout << std::setw(10) << std::setiosflags(std::ios::right) << map_wirePosition[detectorID][i] + 0.5*planes[detectorID].cellWidth;
1124  std::cout << std::endl;
1125  }
1126  }
1127 }
1128 
1129 void GeomSvc::printAlignPar()
1130 {
1131  std::cout << "detectorID DetectorName offset_pos offset_z offset_phi" << std::endl;
1132  for(std::map<std::string, int>::iterator iter = map_detectorID.begin(); iter != map_detectorID.end(); ++iter)
1133  {
1134  if(iter->second>nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes) continue;
1135  std::cout << iter->second << " " << iter->first << " " << planes[(*iter).second].deltaW << " " << planes[(*iter).second].deltaZ << " " << planes[(*iter).second].rotZ << std::endl;
1136  }
1137 }
1138 
1139 void GeomSvc::printTable()
1140 {
1141  std::cout << "detectorID detectorName planeType Spacing Xoffset overlap width height nElement angleFromVertical z0 x0 y0 deltaW" << std::endl;
1142  for(std::map<std::string, int>::iterator iter = map_detectorID.begin(); iter != map_detectorID.end(); ++iter)
1143  {
1144  if(iter->second>nChamberPlanes+nHodoPlanes+nPropPlanes+nDarkPhotonPlanes) continue;
1145  std::cout << planes[iter->second] << std::endl;
1146  }
1147 }
operator<<
std::ostream & operator<<(std::ostream &os, const Plane &plane)
Definition: GeomSvc.cxx:180
nPropPlanes
#define nPropPlanes
Definition: GlobalConsts.h:8
nChamberPlanes
#define nChamberPlanes
Definition: GlobalConsts.h:6
nDarkPhotonPlanes
#define nDarkPhotonPlanes
Definition: GlobalConsts.h:9
nHodoPlanes
#define nHodoPlanes
Definition: GlobalConsts.h:7
NULL
#define NULL
Definition: Pdb.h:9
GeomParamPlane::GetPlane
Plane * GetPlane(const int idx)
Definition: GeomParamPlane.h:38
GeomSvc
User interface class about the geometry of detector planes.
Definition: GeomSvc.h:164
GeomSvc::printWirePosition
void printWirePosition()
Definition: GeomSvc.cxx:1109
GeomSvc::findPatternInDetector
bool findPatternInDetector(int detectorID, std::string pattern)
Definition: GeomSvc.cxx:614
GeomSvc::printTable
void printTable()
Definition: GeomSvc.cxx:1139
GeomSvc::initWireLUT
void initWireLUT()
Definition: GeomSvc.cxx:537
GeomSvc::initPlaneDirect
void initPlaneDirect()
Definition: GeomSvc.cxx:411
GeomSvc::getDriftDistance
double getDriftDistance(int detectorID, double tdcTime)
Definition: GeomSvc.cxx:834
GeomSvc::isInElement
bool isInElement(int detectorID, int elementID, double x, double y, double tolr=0.)
Definition: GeomSvc.cxx:630
GeomSvc::isPropTube
bool isPropTube(const int detectorID) const
Return "true" for prop tube planes.
Definition: GeomSvc.cxx:792
GeomSvc::isHodo
bool isHodo(const int detectorID) const
Return "true" for hodo planes.
Definition: GeomSvc.cxx:782
GeomSvc::getWireEndPoints
void getWireEndPoints(int detectorID, int elementID, double &x_min, double &x_max, double &y_min, double &y_max)
Definition: GeomSvc.cxx:725
GeomSvc::printAlignPar
void printAlignPar()
Debugging print of the content.
Definition: GeomSvc.cxx:1129
GeomSvc::instance
static GeomSvc * instance()
singlton instance
Definition: GeomSvc.cxx:212
GeomSvc::toLocalDetectorName
void toLocalDetectorName(std::string &detectorName, int &eID)
Convert the official detectorName to local detectorName.
Definition: GeomSvc.cxx:740
GeomSvc::loadMilleAlignment
void loadMilleAlignment(const std::string &alignmentFile_mille)
Definition: GeomSvc.cxx:1025
GeomSvc::getMeasurement
void getMeasurement(int detectorID, int elementID, double &measurement, double &dmeasurement)
Convert the detectorID and elementID to the actual hit position.
Definition: GeomSvc.cxx:651
GeomSvc::loadOnlineAlignment
void loadOnlineAlignment(const std::string &alignmentFile_mille)
Load parameters used in the online-alignment mode.
Definition: GeomSvc.cxx:885
GeomSvc::getEndPoints
void getEndPoints(int detectorID, int elementID, TVectorD &ep1, TVectorD &ep2)
Definition: GeomSvc.cxx:662
GeomSvc::isInKMAG
bool isInKMAG(double x, double y)
Definition: GeomSvc.cxx:643
GeomSvc::isInh
bool isInh(const int detectorID) const
Return "true" for BeforeInh and AfterInh signals.
Definition: GeomSvc.cxx:812
GeomSvc::init
void init()
Initialization, either from MySQL or from ascii file.
Definition: GeomSvc.cxx:240
GeomSvc::getHodoStation
int getHodoStation(const int detectorID) const
Return a station number (1-4) for hodo planes or "0" for others.
Definition: GeomSvc.cxx:822
GeomSvc::getExpElementID
int getExpElementID(int detectorID, double pos_exp)
Definition: GeomSvc.cxx:677
GeomSvc::isInTime
bool isInTime(int detectorID, double tdcTime)
Definition: GeomSvc.cxx:1104
GeomSvc::getInterceptionFast
double getInterceptionFast(int detectorID, double tx, double ty, double x0, double y0) const
Definition: GeomSvc.cxx:860
GeomSvc::loadAlignment
void loadAlignment(const std::string &alignmentFile_chamber, const std::string &alignmentFile_hodo, const std::string &alignmentFile_prop)
Definition: GeomSvc.cxx:924
GeomSvc::loadCalibration
void loadCalibration(const std::string &calibrateFile)
Definition: GeomSvc.cxx:1065
GeomSvc::getDetectorName
std::string getDetectorName(const int &detectorID) const
Definition: GeomSvc.h:223
GeomSvc::isChamber
bool isChamber(const int detectorID) const
Return "true" for chamber planes.
Definition: GeomSvc.cxx:772
GeomSvc::getDCA
double getDCA(int detectorID, int elementID, double tx, double ty, double x0, double y0)
Definition: GeomSvc.cxx:865
GeomSvc::initPlaneDbSvc
void initPlaneDbSvc()
Definition: GeomSvc.cxx:419
GeomSvc::isInPlane
bool isInPlane(int detectorID, double x, double y)
See if a point is in a plane.
Definition: GeomSvc.cxx:622
GeomSvc::close
void close()
Close the geometry service before exit or starting a new one.
Definition: GeomSvc.cxx:223
GeomSvc::get2DBoxSize
void get2DBoxSize(int detectorID, int elementID, double &x_min, double &x_max, double &y_min, double &y_max)
Definition: GeomSvc.cxx:700
GeomSvc::isDPHodo
bool isDPHodo(const int detectorID) const
Return "true" for DP hodo planes.
Definition: GeomSvc.cxx:802
GeomSvc::getDetectorIDs
std::vector< int > getDetectorIDs(std::string pattern)
Definition: GeomSvc.cxx:594
PHFlag::FlagExist
virtual int FlagExist(const std::string &name) const
Definition: PHFlag.cc:255
PHFlag::get_IntFlag
virtual int get_IntFlag(const std::string &name) const
Definition: PHFlag.cc:117
PHFlag::get_CharFlag
virtual const std::string get_CharFlag(const std::string &flag) const
Definition: PHFlag.cc:13
Plane
Definition: GeomSvc.h:37
Plane::nElements
int nElements
Definition: GeomSvc.h:69
Plane::resolution
double resolution
Definition: GeomSvc.h:102
Plane::spacing
double spacing
Definition: GeomSvc.h:70
Plane::rotY
double rotY
Definition: GeomSvc.h:100
Plane::angleFromVert
double angleFromVert
Definition: GeomSvc.h:74
Plane::nVec
TVectorD nVec
Definition: GeomSvc.h:114
Plane::deltaZ
double deltaZ
Definition: GeomSvc.h:96
Plane::rotX
double rotX
Definition: GeomSvc.h:99
Plane::yVec
TVectorD yVec
Definition: GeomSvc.h:118
Plane::overlap
double overlap
Definition: GeomSvc.h:73
Plane::x2
double x2
Definition: GeomSvc.h:86
Plane::elementPos
std::vector< double > elementPos
Definition: GeomSvc.h:127
Plane::rotM
TMatrixD rotM
Definition: GeomSvc.h:119
Plane::sintheta
double sintheta
Definition: GeomSvc.h:75
Plane::deltaW_module
double deltaW_module[9]
Definition: GeomSvc.h:98
Plane::zc
double zc
Definition: GeomSvc.h:107
Plane::rotZ
double rotZ
Definition: GeomSvc.h:101
Plane::rtprofile
TSpline3 * rtprofile
Definition: GeomSvc.h:124
Plane::rZ
double rZ
Definition: GeomSvc.h:111
Plane::costheta
double costheta
Definition: GeomSvc.h:76
Plane::tantheta
double tantheta
Definition: GeomSvc.h:77
Plane::z1
double z1
Definition: GeomSvc.h:85
Plane::cellWidth
double cellWidth
Definition: GeomSvc.h:71
Plane::xc
double xc
Definition: GeomSvc.h:105
Plane::getWirePosition
double getWirePosition(int elementID) const
Definition: GeomSvc.cxx:146
Plane::update
void update()
Definition: GeomSvc.cxx:78
Plane::planeType
int planeType
Definition: GeomSvc.h:66
Plane::tmax
double tmax
Definition: GeomSvc.h:123
Plane::y2
double y2
Definition: GeomSvc.h:87
Plane::xoffset
double xoffset
Definition: GeomSvc.h:72
Plane::getEndPoint
TVectorD getEndPoint(int elementID, int sign=-1) const
Definition: GeomSvc.cxx:151
Plane::deltaW
double deltaW
Definition: GeomSvc.h:97
Plane::thetaX
double thetaX
Definition: GeomSvc.h:89
Plane::xVec
TVectorD xVec
Definition: GeomSvc.h:117
Plane::x0
double x0
Definition: GeomSvc.h:80
Plane::detectorID
int detectorID
Definition: GeomSvc.h:64
Plane::z0
double z0
Definition: GeomSvc.h:82
Plane::x1
double x1
Definition: GeomSvc.h:83
Plane::uVec
TVectorD uVec
Definition: GeomSvc.h:115
Plane::thetaZ
double thetaZ
Definition: GeomSvc.h:91
Plane::Plane
Plane()
Definition: GeomSvc.cxx:33
Plane::rY
double rY
Definition: GeomSvc.h:110
Plane::yc
double yc
Definition: GeomSvc.h:106
Plane::detectorName
std::string detectorName
Definition: GeomSvc.h:65
Plane::deltaX
double deltaX
Definition: GeomSvc.h:94
Plane::intercept
double intercept(double tx, double ty, double x0_track, double y0_track) const
Definition: GeomSvc.cxx:129
Plane::deltaY
double deltaY
Definition: GeomSvc.h:95
Plane::tmin
double tmin
Definition: GeomSvc.h:122
Plane::getW
double getW(double x, double y) const
Definition: GeomSvc.h:48
Plane::y1
double y1
Definition: GeomSvc.h:84
Plane::y0
double y0
Definition: GeomSvc.h:81
Plane::vVec
TVectorD vVec
Definition: GeomSvc.h:116
Plane::thetaY
double thetaY
Definition: GeomSvc.h:90
Plane::z2
double z2
Definition: GeomSvc.h:88
Plane::wc
double wc
Definition: GeomSvc.h:108
Plane::rX
double rX
Definition: GeomSvc.h:109
RunParamBase::SetMapID
void SetMapID(const std::string map_id)
Definition: RunParamBase.h:28
RunParamBase::ReadFromDB
void ReadFromDB()
Definition: RunParamBase.cc:102
RunParamBase::SetMapIDbyDB
void SetMapIDbyDB(const std::string map_id)
Definition: RunParamBase.cc:31
RunParamBase::GetMapID
std::string GetMapID()
Definition: RunParamBase.h:27
recoConsts::instance
static recoConsts * instance()
Definition: recoConsts.cc:7
recoConsts::set_CharFlag
virtual void set_CharFlag(const std::string &name, const std::string &flag)
overide the virtual function to expand the environmental variables
Definition: recoConsts.cc:21