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