StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EEmcSmdGeom.cxx
1 
35 #include "Stiostream.h"
36 #include "EEmcSmdGeom.h"
37 #include "EEmcStripGeom.h"
38 #include <assert.h>
39 #include <TMath.h>
40 
41 // decouple from StarClassLibrary
42 //#include "PhysicalConstants.h"
43 #ifndef HEP_SYSTEM_OF_UNITS_H
44 static const double radian = 1.;
45 static const double pi = M_PI; // from <math.h>
46 static const double degree = (M_PI/180.0)*radian;
47 #endif
48 
49 
50 #include "StMessMgr.h"
51 
52 
54 ClassImp(EEmcSmdGeom)
55 
56 EEmcSmdGeom::EEmcSmdGeom()
57  : TObject()
58 {
59  for(int iSec = 0; iSec < kEEmcNumSectors; iSec++) mIsSectorIn[iSec] = true;
60 };
62 EEmcSmdGeom::~EEmcSmdGeom(){
63  delete sInstance;
64  sInstance = 0;
65 }
66 
68 void EEmcSmdGeom::init(){
69  buildSmdGeom();
70 }
71 
72 EEmcSmdGeom* EEmcSmdGeom::sInstance = 0;
73 // all setctors
74 EEmcSmdGeom* EEmcSmdGeom::instance() {
75  if(!sInstance){
76  sInstance = new EEmcSmdGeom();
77  sInstance->init();
78  }
79  return sInstance;
80 }
81 
82 // selected sectors
83 EEmcSmdGeom* EEmcSmdGeom::instance(intVec sectorIdVec) {
84  if(!sInstance){
85  sInstance = new EEmcSmdGeom();
86  sInstance->setSectors(sectorIdVec);
87  sInstance->init();
88  }
89  return sInstance;
90 }
91 
92 // build a glabal geometry database from local coordinates
93 void EEmcSmdGeom::buildSmdGeom(){
94  mEEmcSmdParam.stripWidth = 0.5;
95  for(int iPlane=0; iPlane<kEEmcNumSmdPlanes; iPlane++)
96  mEEmcSmdParam.rOffset[iPlane] = kEEmcSmdROffset[iPlane];
97 
98  float x0[kEEmcNumStrips];
99  float y1[kEEmcNumStrips];
100  float y2[kEEmcNumStrips];
101  float length[kEEmcNumStrips];
102  float x0Edge[kEEmcNumEdgeStrips];
103  float y1Edge[kEEmcNumEdgeStrips];
104  float y2Edge[kEEmcNumEdgeStrips];
105  float lengthEdge[kEEmcNumEdgeStrips];
106 
107 // fill variable arrays with data in EmcStripGeom.h
108  for (int i = 0; i < kEEmcNumStrips; i++) {
109  x0[i] = EEmcStripGeomData[i].x0;
110  y1[i] = EEmcStripGeomData[i].y1;
111  y2[i] = EEmcStripGeomData[i].y2;
112  length[i] = EEmcStripGeomData[i].length;
113  }
114  for (int i = 0; i < kEEmcNumEdgeStrips; i++) {
115  x0Edge[i] = EEmcEdgeStripGeomData[i].x0;
116  y1Edge[i] = EEmcEdgeStripGeomData[i].y1;
117  y2Edge[i] = EEmcEdgeStripGeomData[i].y2;
118  lengthEdge[i] = EEmcEdgeStripGeomData[i].length;
119  }
120 /*
121 // write x, y, and length to a file for a check
122  FILE *fp;
123  fp = fopen("arrCheck.lis", "w");
124  fprintf(fp, " Reqular\n");
125  for (int i = 0; i < kEEmcNumStrips; i++) {
126  fprintf(fp, " strip %d: x0, y1, y2, length = %f %f %f %f\n",
127  i, x0[i], y1[i], y2[i], length[i]);
128  }
129  fprintf(fp, " Edge\n");
130  for (int i = 0; i < kEEmcNumEdgeStrips; i++) {
131  fprintf(fp, " strip %d: x0, y1, y2, length = %f %f %f %f\n",
132  i, x0Edge[i], y1Edge[i], y2Edge[i], lengthEdge[i]);
133  }
134  fclose(fp);
135 */
136  float delPhi = 2*pi/degree/kEEmcNumSectors;
137  float PhiRotation[kEEmcNumSmdUVs][kEEmcNumSectors];
138 
139 // calculate rotation angles to prepare for local-global transformation
140  for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
141  for(int iSec=0; iSec<kEEmcNumSectors; iSec++){
142  PhiRotation[iUV][iSec]=(-15.0 + iSec*delPhi)*degree;
143  if(iUV == 1) PhiRotation[iUV][iSec] = -1.0*PhiRotation[iUV][iSec];
144  }
145  }
146 
147 // loop over planes
148  for (int iPlane = 0; iPlane < kEEmcNumSmdPlanes; iPlane++) {
149  float globalX1, globalY1, globalX2, globalY2;
150  float x0Corr, y1Corr, y2Corr, lengthCorr;
151  float phi1, phi2, phiMin, phiMax;
152  float r, rMin, rMax;
153 
154  mEEmcSmdParam.zPlane[iPlane] = kEEmcZSMD +
155  (iPlane - kEEmcNumSmdPlanes + 2) * kEEmcSmdZPlaneShift ;
156 // loop over UV
157  for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
158 
159 // loop over sectors
160  for(int iUVSec=iPlane+1-iUV; iUVSec<kEEmcNumSectors+1-iUV;
161  iUVSec=iUVSec+kEEmcNumSmdPlanes) {
162  int iSec;
163  if(iUVSec == 12) iSec = 0;
164  else iSec = iUVSec;
165 
166  if(IsSectorIn(iSec)) {
167 
168  rMin = 1000.0;
169  rMax = 0.0;
170  phiMin = pi;
171  phiMax = -pi;
172 
173 // loop over strips
174  for (int iStrip = 0; iStrip < kEEmcNumStrips; iStrip++) {
175  if(kEEmcSmdMapEdge[iPlane][iSec] && iStrip > kEEmcNumEdgeStrips-1)
176  break;
177 // Id = index + 1 in all cases
178  StructEEmcStripId stripStructId;
179  stripStructId.sectorId = iSec+1;
180  stripStructId.UVId = iUV+1;
181  stripStructId.stripId = iStrip + 1;
182  stripStructId.planeId = iPlane + 1;
183 
184  StructEEmcStrip* stripPtr = new StructEEmcStrip;
185  stripPtr->stripStructId = stripStructId;
186 
187 // correct for radius offset
188  x0Corr = x0[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
189  y2Corr = y2[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
190  if(kEEmcSmdMapEdge[iPlane][iSec]) {
191  y1Corr = y1Edge[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
192  lengthCorr = lengthEdge[iStrip];
193  }
194  else {
195  y1Corr = y1[iStrip] - kEEmcSmdROffset[iPlane]*::sqrt(0.5);
196  lengthCorr = length[iStrip];
197  }
198 // Transform local to gloabal by rotation. Local X axis is perpendicular to
199 // strip, Y is parallel. After rotation, x & y should be swapped for V sector.
200 
201  if(iUV == 0) {
202  globalX1 = x0Corr*cos(PhiRotation[iUV][iSec])+
203  y1Corr*sin(PhiRotation[iUV][iSec]);
204  globalY1 = y1Corr*cos(PhiRotation[iUV][iSec]) -
205  x0Corr*sin(PhiRotation[iUV][iSec]);
206  globalX2 = x0Corr*cos(PhiRotation[iUV][iSec]) +
207  y2Corr*sin(PhiRotation[iUV][iSec]);
208  globalY2 = y2Corr*cos(PhiRotation[iUV][iSec]) -
209  x0Corr*sin(PhiRotation[iUV][iSec]);
210  }
211  else {
212  globalX1 = y1Corr*cos(PhiRotation[iUV][iSec]) -
213  x0Corr*sin(PhiRotation[iUV][iSec]);
214  globalY1 = x0Corr*cos(PhiRotation[iUV][iSec]) +
215  y1Corr*sin(PhiRotation[iUV][iSec]);
216  globalX2 = y2Corr*cos(PhiRotation[iUV][iSec]) -
217  x0Corr*sin(PhiRotation[iUV][iSec]);
218  globalY2 = x0Corr*cos(PhiRotation[iUV][iSec]) +
219  y2Corr*sin(PhiRotation[iUV][iSec]);
220  }
221  r = ::sqrt(globalX1*globalX1 + globalY1*globalY1);
222  if(r < rMin) rMin = r;
223  r = ::sqrt(globalX2*globalX2 + globalY2*globalY2);
224  if(r > rMax) rMax = r;
225 
226 //Fill StripPtrVec
227  stripPtr->end1.SetX(globalX1) ;
228  stripPtr->end1.SetY(globalY1) ;
229  stripPtr->end1.SetZ(mEEmcSmdParam.zPlane[iPlane]);
230  stripPtr->end2.SetX(globalX2) ;
231  stripPtr->end2.SetY(globalY2) ;
232  stripPtr->end2.SetZ(mEEmcSmdParam.zPlane[iPlane]);
233  stripPtr->length = lengthCorr;
234 
235  phi1 = stripPtr->end1.Phi();
236  phi2 = stripPtr->end2.Phi();
237 
238  if(iSec != kEEmcSmdSectorIdPhiCrossPi - 1) {
239  if(phi1 < phiMin) phiMin = phi1;
240  if(phi1 > phiMax) phiMax = phi1;
241  if(phi2 < phiMin) phiMin = phi2;
242  if(phi2 > phiMax) phiMax = phi2;
243  }
244  else {
245  if(phi1 > 0) if(phi1 < phiMin) phiMin = phi1;
246  if(phi1 < 0 ) if(phi1 > phiMax) phiMax = phi1;
247  if(phi2 > 0) if(phi2 < phiMin) phiMin = phi2;
248  if(phi2 < 0) if(phi2 > phiMax) phiMax = phi2;
249  }
250  mEEmcSector[iUV][iSec].stripPtrVec.push_back(stripPtr);
251  } // loop over iStrip
252 
253 // Fill mEEmcSectors
254  mEEmcSector[iUV][iSec].sectorId = iSec+1;
255  mEEmcSector[iUV][iSec].planeId = iPlane+1;
256  mEEmcSector[iUV][iSec].phiMin = phiMin;
257  mEEmcSector[iUV][iSec].phiMax = phiMax;
258  mEEmcSector[iUV][iSec].rMin = rMin;
259  mEEmcSector[iUV][iSec].rMax = rMax;
260 
261  } // if sector selected
262  } // end of iUVSec loop
263  } // end of iUV loop
264  } // end of iPlane loop
265 
266 
267  // build revers mapping (iUV,iSec) --> iPlane
268  memset(kEEmcSmdMap_iPlane,0,sizeof(kEEmcSmdMap_iPlane));
269  for (int iPlane = 0; iPlane < kEEmcNumSmdPlanes; iPlane++)
270  for(int iSec=0; iSec<kEEmcNumSectors; iSec++) {
271  int iuv=kEEmcSmdMapUV[iPlane][iSec];
272  if(iuv<0 ) continue;
273  assert(iuv<kEEmcNumSmdUVs);
274  kEEmcSmdMap_iPlane[iuv][iSec]=iPlane;
275  }
276 
277  buildStripPtrVector();
278 
279 } // end of buildSmdGeom
280 //===========================================================
281 //===========================================================
282 
283 // build mStripPtrVector with getEEmcSector()
284 void EEmcSmdGeom::buildStripPtrVector() {
285  StructEEmcStrip *dummyStripPtr = new StructEEmcStrip;
286  *dummyStripPtr = initStrip();
287  EEmcStripPtrVec stripPtrVec;
288  EEmcStripPtrVecIter p;
289  for(int iSec = 0; iSec < kEEmcNumSectors; iSec++) {
290  if(mIsSectorIn[iSec]) {
291  for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
292  stripPtrVec = getEEmcSector(iUV,iSec).stripPtrVec;
293  p = stripPtrVec.begin();
294  int PlaneId = getEEmcSector(iUV,iSec).planeId;
295  while(p !=stripPtrVec.end()) {
296  mStripPtrVector.push_back(*p);
297  p++;
298 
299  }
300  if(kEEmcSmdMapEdge[PlaneId-1][iSec]) {
301  for(int i=0; i < kEEmcNumStrips - kEEmcNumEdgeStrips; i++)
302  mStripPtrVector.push_back(dummyStripPtr);
303  }
304  }
305  }
306  else {
307  for(int iUV = 0; iUV < kEEmcNumSmdUVs; iUV++) {
308  for(int iStrip = 0; iStrip < kEEmcNumStrips; iStrip++)
309  mStripPtrVector.push_back(dummyStripPtr);
310  }
311  }
312  }
313 }
314 
315 // set status of selected sectors
316 void EEmcSmdGeom::setSectors(const intVec sectorIdVec) {
317  for(int iSec = 0; iSec< kEEmcNumSectors; iSec++)
318  mIsSectorIn[iSec] = false;
319  for(unsigned int i = 0; i < sectorIdVec.size(); i++) {
320  for(int iSec = 0; iSec< kEEmcNumSectors; iSec++) {
321  if (sectorIdVec[i] == iSec+1) mIsSectorIn[iSec] = true;
322  }
323  }
324 }
325 
326 // instance and initialize a strip
327 StructEEmcStrip EEmcSmdGeom::initStrip() const {
328  TVector3 zero(0,0,0);
329  StructEEmcStrip strip;
330  strip.stripStructId.stripId = 0;
331  strip.stripStructId.UVId = 0;
332  strip.stripStructId.sectorId = 0;
333  strip.stripStructId.planeId = 0;
334  strip.end1 = zero;
335  strip.end2 = zero;
336  strip.length = 0.0;
337  return strip;
338 }
339 
340 // return index of a sector from a global point in a plane
341 Int_t EEmcSmdGeom::getEEmcISec(const Int_t iPlane,
342  const TVector3& point) const {
343  Int_t indexSec = -1;
344  float phiMin, phiMax, rMin, rMax;
345  float phi = point.Phi();
346  float r = ::sqrt(point.x()*point.x() + point.y()*point.y());
347 
348  for (int iSec = 0; iSec < kEEmcNumSectors; iSec++) {
349  int iUV = kEEmcSmdMapUV[iPlane][iSec];
350  if(iUV >= 0 && IsSectorIn(iSec)) {
351  phiMin = mEEmcSector[iUV][iSec].phiMin;
352  phiMax = mEEmcSector[iUV][iSec].phiMax;
353  rMin = mEEmcSector[iUV][iSec].rMin;
354  rMax = mEEmcSector[iUV][iSec].rMax;
355  if(iSec != kEEmcSmdSectorIdPhiCrossPi - 1) {
356  if (phi >= phiMin && phi < phiMax && r > rMin && r < rMax) {
357  indexSec = iSec;
358  break;
359  }
360  }
361 // sector9 between 165 deg (Min) and -165 deg (Max)
362  else {
363  if(((phi > 0.0 && phi >= phiMin) || (phi < 0.0 && phi < phiMax))
364  && r > rMin && r < rMax){
365  indexSec = iSec;
366  break;
367  }
368  }
369  }
370  }
371  return indexSec;
372 }
373 
374 // return a strip pointer from indices
375 StructEEmcStrip* EEmcSmdGeom::getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec) {
376  int i = iStrip + iUV*kEEmcNumStrips + iSec*kEEmcNumStrips*kEEmcNumSmdUVs;
377  return mStripPtrVector[i];
378 }
379 
380 // return a strip pointer from indices
381 const StructEEmcStrip* EEmcSmdGeom::getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec) const {
382  int i = iStrip + iUV*kEEmcNumStrips + iSec*kEEmcNumStrips*kEEmcNumSmdUVs;
383  return mStripPtrVector[i];
384 }
385 
386 
387 //==================================================================
388 // get DCA strip pointer from a point
389 // iPlane=[0,1,2] - experts only, changes meaning form sector to sector
390 const StructEEmcStrip* EEmcSmdGeom::getDcaStripPtr(const Int_t iPlane,
391  const Int_t iSec, const TVector3& point, Float_t* dca) const
392 {
393  // StructEEmcStrip* stripPtr;
394  // stripPtr = new StructEEmcStrip;
395  *dca = 1000.0;
396  int iStrip = -1;
397  //$$$ int iUV;
398  float x1,y1,x2,y2,mu,d;
399  EEmcStripPtrVec stripPtrVec;
400  EEmcStripPtrVecIter p;
401 
402 
403  Int_t iUV = kEEmcSmdMapUV[iPlane][iSec]; // jcw 2/6/04 moved here from $$$
404  // iUV may be used uninitialized
405  // below, otherwise...
406 // int iSec = getEEmcISec(iPlane, point);
407  if(iSec >= 0 && IsSectorIn(iSec)) {
408  //$$$iUV = kEEmcSmdMapUV[iPlane][iSec];
409  stripPtrVec = getEEmcSector(iUV,iSec).stripPtrVec;
410  p = stripPtrVec.begin();
411  while(p != stripPtrVec.end()) {
412  x1 = (*p)->end1.x();
413  y1 = (*p)->end1.y();
414  x2 = (*p)->end2.x();
415  y2 = (*p)->end2.y();
416  mu = -1.0/::sqrt((y2-y1)*(y2-y1) + (x2-x1)*(x2-x1)) *
417  ((x2*y1-x1*y2)/fabs(x2*y1-x1*y2));
418 // distance d carries a sign
419  d = ((y2-y1)*point.x() + (x1-x2)*point.y() + (x2*y1-x1*y2))*mu;
420 
421  if(fabs(d) < fabs(*dca)) {
422  *dca = d;
423  iStrip = (*p)->stripStructId.stripId - 1;
424  }
425  if(d < 0) break;
426  p++;
427  }
428  }
429  if(iStrip >=0) {
430  //stripPtr = getStripPtr(iStrip,iUV,iSec);
431  //return stripPtr;
432  return getStripPtr(iStrip,iUV,iSec);
433  }
434  else {
435  StructEEmcStrip *stripPtr = new StructEEmcStrip;
436  (*stripPtr) = initStrip();
437  // std::cout << "NO dca strip found in plane (sector empty or not in)" << std::endl; //silentium
438 
439  return stripPtr;
440  }
441 }
442 
443 //==================================================================
444 // get DCA strip pointer from a point
445 // iPlane=[0,1,2] - experts only, changes meaning form sector to sector
446 const StructEEmcStrip* EEmcSmdGeom::getDcaStripPtr(const Int_t iPlane,
447  const TVector3& point, Float_t* dca) const {
448  int iSec = getEEmcISec(iPlane, point);
449  return getDcaStripPtr(iPlane, iSec, point, dca);
450 }
451 
452 
453 //==================================================================
454 // get DCA strip pointer from a point
455 // iUV=[0,1] , maps [U,V]
456 // Warn, this code does not handle well sector boundaries, it ignores possible tripple overlaps and uses only nominal U or V planes in a given sector
457 const StructEEmcStrip* EEmcSmdGeom::getDca2Strip(const Int_t iUV,
458  const TVector3& point, Float_t* dca) const {
459  assert(iUV>=0 || iUV<kEEmcNumSmdUVs);
460  float phiDeg=atan2(point.y(),point.x())/3.1316*180.;
461  //printf("phiDeg=%.1f \n",phiDeg);
462  int iSec= ((int) ( 12.-(phiDeg-75.)/30.) )%12;
463  assert(iSec>=0);
464  assert( iSec<kEEmcNumSectors);
465  int iPlane= kEEmcSmdMap_iPlane[iUV][iSec];// now find mapping iUV --> iPlane
466  assert(iPlane>=0 && iPlane<kEEmcNumSmdPlanes);
467  return getDcaStripPtr(iPlane, iSec, point, dca);
468 }
469 
470 
471 //==================================================================
472 // match two strips
473  bool EEmcSmdGeom::matchStrips(const StructEEmcStripId &stripStructId1,
474  const StructEEmcStripId &stripStructId2,
475  Int_t nTolerance) const {
476  bool match = false;
477  if(stripStructId1.UVId == stripStructId2.UVId &&
478  stripStructId1.sectorId == stripStructId2.sectorId) {
479  if((TMath::Abs(stripStructId1.stripId - stripStructId2.stripId) <= nTolerance))
480  match = true;
481  }
482  return match;
483 }
484 
485 
486 
487 TVector3 EEmcSmdGeom::getstripEnd(const StructEEmcStrip &strip,
488  const Int_t endId) const {
489  TVector3 end;
490  if(endId == 1) end = strip.end1;
491  else end = strip.end2;
492 
493  return end;
494 }
495 
496 // methods of printout
498 void EEmcSmdGeom::printGeom(ostream& os) const {
499  os << "------EEmcSmdGeom::printGeom()------" << endl;
500  os << " " << "z[3] = "
501  << " " << getEEmcSmdParam().zPlane[0]
502  << " " << getEEmcSmdParam().zPlane[1]
503  << " " << getEEmcSmdParam().zPlane[2] << endl;
504  os << " " << "rOffset[3] = "
505  << " " << getEEmcSmdParam().rOffset[0]
506  << " " << getEEmcSmdParam().rOffset[1]
507  << " " << getEEmcSmdParam().rOffset[2] << endl;
508  os << " " << "stripWidth = "
509  << " " << getEEmcSmdParam().stripWidth << endl;
510  os << "---------------------------------------" << endl;
511 }
512 
514 void EEmcSmdGeom::printSector(const StructEEmcSmdSector sector, ostream& os) const {
515  float delPhi;
516  int iUV = kEEmcSmdMapUV[sector.planeId-1][sector.sectorId-1];
517  delPhi = (sector.phiMax - sector.phiMin)/degree;
518  if(sector.sectorId == kEEmcSmdSectorIdPhiCrossPi)
519  delPhi = 2*pi/degree + delPhi;
520 
521  os << "------EEmcSmdGeom::printSector()------" << endl;
522  os << kEEmcSmdUVChar[iUV] << " Sector: sectorId, planeId, nStrips = "
523  << " " << sector.sectorId
524  << " " << sector.planeId
525  << " " << sector.stripPtrVec.size() << endl;
526  os << " phiMin, phiMax, delPhi = "
527  << " " << sector.phiMin/degree
528  << " " << sector.phiMax/degree
529  << " " << delPhi << endl;
530  os << " rMin, rMax delR = "
531  << " " << sector.rMin
532  << " " << sector.rMax
533  << " " << sector.rMax - sector.rMin << endl;
534  os << "------------------------------------" << endl;
535 }
536 
538 void EEmcSmdGeom::printStrip(const StructEEmcStrip strip, ostream& os) const {
539  char UVChar;
540  if(strip.stripStructId.sectorId == 0) UVChar = 'X';
541  else
542  UVChar = kEEmcSmdUVChar[strip.stripStructId.UVId - 1];
543 
544  os << "------EEmcSmdGeom::printStrip()------" << endl;
545 
546  os << "Strip: sectorId, planeUV, stripId, planeId = "
547  << " " << strip.stripStructId.sectorId
548  << " " << UVChar
549  << " " << strip.stripStructId.stripId
550  << " " << strip.stripStructId.planeId << endl;
551  os << " x1, y1, x2, y2, z = "
552  << " " << strip.end1.x()
553  << " " << strip.end1.y()
554  << " " << strip.end2.x()
555  << " " << strip.end2.y()
556  << " " << strip.end2.z() << endl;
557  os << " phi1, phi2, length = "
558  << " " << strip.end1.Phi()/degree
559  << " " << strip.end2.Phi()/degree
560  << " " << strip.length << endl;
561  os << "------------------------------------" << endl;
562 }
563 
565 void EEmcSmdGeom::printStripId(const StructEEmcStripId stripStructId, ostream& os) const {
566  char UVChar;
567  if(stripStructId.sectorId == 0) UVChar = 'X';
568  else
569  UVChar = kEEmcSmdUVChar[stripStructId.UVId - 1];
570 
571  os << "------EEmcSmdGeom::printStripId()------" << endl;
572  os << "Strip: sectorId, stripId, planeId = "
573  << " " << stripStructId.sectorId
574  << UVChar
575  << " " << stripStructId.stripId
576  << " " << stripStructId.planeId << endl;
577  os << "------------------------------------" << endl;
578 }
579 
580 
581 
582 // Move to StEEmcSmdGeom
583 #if 0
584 // printout delPhi and centerPhi used in ITTF
585 void EEmcSmdGeom::printSectorPhis(const Int_t iPlane, const Int_t iSec,
586  ostream& os ) {
587  int iUV;
588  iUV = kEEmcSmdMapUV[iPlane][iSec];
589 
590  os << "------EEmcSmdGeom::printPhis()------" << endl;
591  os << " planeId = " << iPlane + 1 << " sectorId = " << iSec + 1 << endl;
592  if(iUV >= 0)
593  os << " " << kEEmcSmdUVChar[iUV] << " Sector" << endl;
594  else
595  os << " Empty" << endl;
596  os << " delPhi = " << getEEmcSmdDelPhi(iPlane, iSec)/degree <<
597  " " << "centerPhi = " << getEEmcSmdCenterPhi(iPlane, iSec)/degree
598  << endl;
599 
600 }
601 #endif
602 
603 
605 
606 //
607 // Function(s) to return the position of the crossing of two strips, given
608 // either the sectorID and strip ID's, or the StructEEmcStrips.
609 //
610 
611 // Wrapper function when we don't have the actual strips, just sector and
612 // strip ID's.
613 TVector3 EEmcSmdGeom::getIntersection ( Int_t sector, Int_t uId, Int_t vId ) const {
614 
615  Int_t nU = getNStrips( sector, 0 );
616  Int_t nV = getNStrips( sector, 1 );
617  if ( uId >= nU || vId >= nV ) {
618  LOG_DEBUG<<"::getIntersection(...) passed invalid strip ID sector="<<sector<<" uId="<<uId<<" vId="<<vId<<std::endl;
619  return TVector3(1.,1.,-999.0);
620  }
621 
622  return getIntersection ( getStripPtr( uId, 0, sector ),
623  getStripPtr( vId, 1, sector ) );
624 
625 }
626 
627 
628 TVector3 EEmcSmdGeom::getIntersection ( Int_t sector, Int_t uId, Int_t vId, const TVector3 &vert ) const {
629 
630  Int_t nU = getNStrips( sector, 0 );
631  Int_t nV = getNStrips( sector, 1 );
632  if ( uId >= nU || vId >= nV ) {
633  LOG_DEBUG<<"::getIntersection(...) passed invalid strip ID sector="<<sector<<" uId="<<uId<<" vId="<<vId<<std::endl;
634  return TVector3(1.,1.,-999.0);
635  }
636 
637  return getIntersection ( getStripPtr( uId, 0, sector ),
638  getStripPtr( vId, 1, sector ), vert );
639 
640 }
641 
642 
643 /*********************************************************/
644 TVector3 EEmcSmdGeom::getIntersection ( const StructEEmcStrip *u,
645  const StructEEmcStrip *v,
646  const TVector3 &vertex ) const
647 {
648  // The strips are arranged sensibly, so that the ID's and the
649  // widths of the strips basically tell us the position of the
650  // crossing point. However, we need to know a few things:
651 
652  Int_t uSectorId = u -> stripStructId.sectorId; // This would be easier if
653  Int_t vSectorId = v -> stripStructId.sectorId; // we were dealing with a
654  //Int_t uPlaneId = u -> stripStructId.planeId; // class instead of a struct
655  //Int_t vPlaneId = v -> stripStructId.planeId; //
656 
657  Int_t uId = u -> stripStructId.stripId;
658  Int_t vId = v -> stripStructId.stripId;
659 
660  // Get vectors pointing to the start of the u and v strips in question,
661  // as well as the ends. NOTE: We pass uId - 1 to getStripPtr, because
662  // that routine expects the c++ _index_... (the convention in this
663  // class is that quantities beginning with an "i" correspond to a
664  // c++ index numbered from 0, rather than a fortran index numbered
665  // from 1...)
666  TVector3 u0 = getStripPtr ( uId-1, 0, uSectorId-1 ) -> end1;
667  TVector3 v0 = getStripPtr ( vId-1, 1, vSectorId-1 ) -> end1;
668 
669 
670  TVector3 uF = getStripPtr ( uId-1, 0, uSectorId-1 ) -> end2;
671  TVector3 vF = getStripPtr ( vId-1, 1, vSectorId-1 ) -> end2;
672 
673 
674  TVector3 u0p=(u0-vertex); // uOp is vector from vertex to beginning of first strip from vertex
675  TVector3 uFp=(uF-vertex); // uFp is vector from vertex to end of first strip from vertex
676  TVector3 v0p=(v0-vertex); // vOp is vector from vertex to beginning of first strip from vertex
677  TVector3 vFp=(vF-vertex); // vFp is vector from vertex to end of first strip from vertex
678 
680  // U then V! return point on V strip
682 
683  // N is normal vector to plane
684  TVector3 Nv = u0p.Cross(uFp);
685 
686  // use components of N to establish D in plane of form A(x-x0)+B(y-y0)+C(z-z0)+D=0, use vertex for point
687  Float_t Dv = (Nv.X()*vertex.X() + Nv.Y()*vertex.Y() + Nv.Z()*vertex.Z());
688 
689  // scale determines how far along vector from v0 to vF to go to get point
690  // of intersection with plane. If it's <0 or >1, then the point of
691  // intersection is off the physical size of the strip. Return an error
692  // value
693  Float_t scale_numerv = (Nv.X()*v0.X() + Nv.Y()*v0.Y() + Nv.Z()*v0.Z() - Dv);
694  Float_t scale_denomv = (Nv.X()*(v0.X()-vF.X()) + Nv.Y()*(v0.Y()-vF.Y()) + Nv.Z()*(v0.Z()-vF.Z()));
695  Float_t scalev=scale_numerv/scale_denomv;
696  if (scalev < 0 || scalev > 1) {
697  TVector3 ErrorVector(1.,1.,-999.);
698  LOG_DEBUG<<GetName()<<"::getIntersection( ) passed non-intersecting SMD strips " << *u << " " << *v << endm;
699  return ErrorVector;
700  }
701 
702 
703  // Rv is the vector to the final point from the origin ON V STRIP
704  TVector3 Rv = (v0 + scalev*(vF-v0));
705 
706 
707 
708  // S is the vector that would go from vertex to end point
709  // TVector3 S = (R-z0);
710 
712  // V THEN U! returns point on u strip
714  // N is normal vector to plane
715  TVector3 Nu = v0p.Cross(vFp);
716 
717  // use components of N to establish D in plane of form A(x-x0)+B(y-y0)+C(z-z0)+D=0, use vertex for point
718  Float_t Du = (Nu.X()*vertex.X() + Nu.Y()*vertex.Y() + Nu.Z()*vertex.Z());
719 
720  // scale determines how far along vector from v0 to vF to go to get point of intersection with plane
721  Float_t scale_numeru = (Nu.X()*u0.X() + Nu.Y()*u0.Y() + Nu.Z()*u0.Z() - Du);
722  Float_t scale_denomu = (Nu.X()*(u0.X()-uF.X()) + Nu.Y()*(u0.Y()-uF.Y()) + Nu.Z()*(u0.Z()-uF.Z()));
723  Float_t scaleu=scale_numeru/scale_denomu;
724  if (scaleu < 0 || scaleu > 1) {
725  TVector3 ErrorVector(1.,1.,-999.);
726  return ErrorVector;
727  }
728 
729 
730  // Ru is the vector to the final point from the origin ON U STRIP
731  TVector3 Ru = (u0 + scaleu*(uF-u0));
732 
733 
734  // gives us vector direction FROM THE VERTEX
735  TVector3 R = ((Rv+Ru)*0.5);
736 
737  return R;
738 
739 
740 
741 }
742 
743 /*********************************************************/
744 
745 
746 TVector3 EEmcSmdGeom::getIntersection ( const StructEEmcStrip *u,
747  const StructEEmcStrip *v
748  ) const {
749 
750  return getIntersection( u, v, TVector3(0,0,0) );
751 
752 }
753 
754 
755 /****************************************************************************/
756 
757 // output strip for smd strips
758 ostream& operator<<(ostream &os, const StructEEmcStrip &strip)
759 {
760  const Char_t *nameUV[]={"U","V"};
761  TString stripname="";
762  if ( strip.stripStructId.sectorId < 10 ) stripname+="0"; stripname+=strip.stripStructId.sectorId;
763  stripname += nameUV[ strip.stripStructId.UVId-1 ];
764  if ( strip.stripStructId.stripId < 10 ) stripname+="0";
765  if ( strip.stripStructId.stripId < 100 ) stripname+="0";
766  stripname+=strip.stripStructId.stripId;
767  os << stripname << " depth=" << strip.stripStructId.planeId;
768  return os;
769 }
770 
771 
773 /*
774  * $Log: EEmcSmdGeom.cxx,v $
775  * Revision 1.17 2015/07/21 17:10:02 jeromel
776  * Wrong init of TVector3 zero corrected
777  *
778  * Revision 1.16 2010/08/26 22:48:55 ogrebeny
779  * Improved constness
780  *
781  * Revision 1.15 2009/08/25 18:33:11 fine
782  * fix the compilation issues under SL5_64_bits gcc 4.3.2
783  *
784  * Revision 1.14 2008/02/17 17:37:32 balewski
785  * demote warning about strips intersecting beyond modle to DEBUG
786  *
787  * Revision 1.13 2007/07/12 19:30:15 fisyak
788  * Add includes for ROOT 5.16
789  *
790  * Revision 1.12 2007/05/01 20:22:25 jwebb
791  * Added error handling to EEmcSmdGeom::getIntersection().
792  *
793  * Revision 1.11 2007/02/17 01:29:05 balewski
794  * less printout
795  *
796  * Revision 1.10 2007/02/02 02:11:11 balewski
797  * simplification of EEmcSmdGeom::getDca2Strip(..)
798  *
799  * Revision 1.9 2007/02/01 20:33:46 wzhang
800  * dca initialized in getDcaStripPtr()
801  *
802  * Revision 1.8 2007/02/01 13:47:39 balewski
803  * bug fix in getDca2Strip(), more methodhs are public
804  *
805  * Revision 1.7 2007/01/26 00:51:08 balewski
806  * too strong protection
807  *
808  * Revision 1.6 2007/01/25 22:33:21 balewski
809  * add:
810  * - better writeup
811  * - new simpler to use method calculating dca fo track to strip, it is just a wrapper, some approximations were used, may fail at the sector boundary
812  *
813  * Revision 1.5 2007/01/12 23:53:14 jwebb
814  * Fix applied to eliminate parralax error in the EEmcSmdGeom::getIntersection()
815  * method.
816  *
817  * Revision 1.4 2004/06/03 23:01:06 jwebb
818  * Fixed memory leak reported by Bob. Fixed another memory leak.
819  *
820  * Revision 1.3 2004/02/06 22:33:08 jwebb
821  * Moved statement to fix warning.
822  *
823  * Revision 1.2 2004/01/29 16:37:25 jwebb
824  * Removed dependence on StMaker.h and PhysicalConstants.h. Should be fully
825  * decoupled from Star environment now.
826  *
827  * Revision 1.1 2004/01/29 15:26:10 jwebb
828  * The StEEmcSmdGeom class was split into two classes. All StRoot-independent
829  * code has been moved to EEmcSmdGeom. TVector3 replaces StThreeVectorD in
830  * all function calls in EEmcSmdGeom. StThreeVectorD wrappers are provided
831  * in StEEmcSmdGeom, for integration into Star framework.
832  *
833  *
834  *****************************************************************************
835  * Log: StEEmcSmdGeom.cxx,v
836  * Revision 1.8 2004/01/12 14:34:09 wzhang
837  * Corrected the usage of EEmcStripPtrVecIter
838  *
839  * Revision 1.7 2003/12/05 00:06:10 jwebb
840  * Member function added to return a vector pointing to the intersection of
841  * two strips.
842  *
843  * Revision 1.6 2003/10/15 15:26:08 wzhang
844  * improved and reorganized
845  *
846  * Revision 1.5 2003/09/02 17:57:56 perev
847  * gcc 3.2 updates + WarnOff
848  *
849  * Revision 1.4 2003/08/22 15:14:26 wzhang
850  * Added ClassImp and method stripEnd
851  *
852  * Revision 1.3 2003/06/11 18:58:19 wzhang
853  * added geometry methods for StiEEmc
854  *
855  * Revision 1.2 2003/04/04 15:33:32 wzhang
856  * included EEmcGeomDefs.h & improved codes
857  *
858  * Revision 1.1 2003/03/28 15:49:59 balewski
859  * first
860  *
861  *
862  *******************************************************************/
EEmcSmdGeom::IsSectorIn
bool IsSectorIn(const Int_t iSec) const
return sector status
Definition: EEmcSmdGeom.h:130
EEmcSmdGeom::mEEmcSector
StructEEmcSmdSector mEEmcSector[kEEmcNumSmdUVs][kEEmcNumSectors]
general geometry variables
Definition: EEmcSmdGeom.h:101
EEmcSmdGeom::getEEmcSmdParam
StructEEmcSmdParam & getEEmcSmdParam()
return SMD geometry parameters
Definition: EEmcSmdGeom.h:136
EEmcSmdGeom::kEEmcSmdMap_iPlane
int kEEmcSmdMap_iPlane[kEEmcNumSmdUVs][kEEmcNumSectors]
sector status.
Definition: EEmcSmdGeom.h:105
EEmcSmdGeom::getNStrips
Int_t getNStrips(Int_t iSec, Int_t iUV) const
Definition: EEmcSmdGeom.h:192
EEmcSmdGeom::getEEmcISec
Int_t getEEmcISec(const Int_t iPlane, const TVector3 &point) const
return index of a sector from a point in a plane
Definition: EEmcSmdGeom.cxx:341
EEmcSmdGeom::getEEmcSector
StructEEmcSmdSector & getEEmcSector(const Int_t iUV, const Int_t iSec)
return structure-sector from iUV and iSec
Definition: EEmcSmdGeom.h:140
EEmcSmdGeom::getDcaStripPtr
const StructEEmcStrip * getDcaStripPtr(const Int_t iPlane, const TVector3 &point, Float_t *dca) const
Definition: EEmcSmdGeom.cxx:446
EEmcSmdGeom::~EEmcSmdGeom
virtual ~EEmcSmdGeom()
default empty destructor
Definition: EEmcSmdGeom.cxx:62
EEmcSmdGeom::getstripEnd
TVector3 getstripEnd(const StructEEmcStrip &strip, const Int_t endId) const
return strip-end of 3D-vector
Definition: EEmcSmdGeom.cxx:487
EEmcSmdGeom::mIsSectorIn
bool mIsSectorIn[kEEmcNumSectors]
storage for all strip pointers
Definition: EEmcSmdGeom.h:104
EEmcSmdGeom::printGeom
void printGeom(ostream &os=cout) const
printout global geometry parameters
Definition: EEmcSmdGeom.cxx:498
EEmcSmdGeom::init
void init()
Initialize geometry class.
Definition: EEmcSmdGeom.cxx:68
EEmcSmdGeom::initStrip
StructEEmcStrip initStrip() const
instance and initialize a strip
Definition: EEmcSmdGeom.cxx:327
EEmcSmdGeom::printSector
void printSector(const StructEEmcSmdSector Sector, ostream &os=cout) const
printout sector-specific geometry parameters
Definition: EEmcSmdGeom.cxx:514
EEmcSmdGeom::mStripPtrVector
EEmcStripPtrVec mStripPtrVector
storage for all strip pointers
Definition: EEmcSmdGeom.h:103
EEmcSmdGeom::EEmcSmdGeom
EEmcSmdGeom()
defaulty constructor
Definition: EEmcSmdGeom.cxx:56
EEmcSmdGeom::matchStrips
bool matchStrips(const StructEEmcStripId &stripStructId1, const StructEEmcStripId &stripStructId2, Int_t nTolerance) const
match two strips
Definition: EEmcSmdGeom.cxx:473
EEmcSmdGeom::buildStripPtrVector
void buildStripPtrVector()
build mStripPtrVector
Definition: EEmcSmdGeom.cxx:284
EEmcSmdGeom::getIntersection
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
Definition: EEmcSmdGeom.cxx:628
EEmcSmdGeom::setSectors
void setSectors(const intVec sectorIdVec)
set sectors for partial EEMC
Definition: EEmcSmdGeom.cxx:316
EEmcSmdGeom::printStrip
void printStrip(const StructEEmcStrip Strip, ostream &os=cout) const
printout strip-specific geometry parameters
Definition: EEmcSmdGeom.cxx:538
EEmcSmdGeom::getStripPtr
StructEEmcStrip * getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec)
return a strip pointer from indices
Definition: EEmcSmdGeom.cxx:375
EEmcSmdGeom::getDca2Strip
const StructEEmcStrip * getDca2Strip(const Int_t iUV, const TVector3 &point, Float_t *dca) const
Definition: EEmcSmdGeom.cxx:457
EEmcSmdGeom::printStripId
void printStripId(const StructEEmcStripId StripId, ostream &os=cout) const
printout stripStructId
Definition: EEmcSmdGeom.cxx:565
Generated by   doxygen 1.8.5