StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTofrGeometry.cxx
1 /*******************************************************************
2  *
3  * $Id: StTofrGeometry.cxx,v 1.12 2018/02/26 23:26:51 smirnovd Exp $
4  *
5  * Authors: Shuwei Ye, Xin Dong
6  *******************************************************************
7  *
8  * Description: Collection of geometry classes for the TOF-MRPC
9  * initializes from GEANT geometry
10  *
11  *******************************************************************/
12 #include "Stiostream.h"
13 #include <math.h>
14 #include <vector>
15 #include <stdlib.h>
16 #include <stdio.h>
17 
18 #include "StTofrGeometry.h"
19 #include "TFile.h"
20 #include "TDataSet.h"
21 #include "TDataSetIter.h"
22 #include "StMaker.h"
23 //#include "StMemStat.h"
24 #include "StMessMgr.h"
25 
26 //#include "Debugger.h"
27 
29 //
30 // group of classes for Tofr Geometry:
31 //
32 // StTofrGeometry(v1), StTofrNode(v2), StTofrGeomNode(v2)(off)
33 // StTofrGeomTray(v1), StTofrGeomSensor(v2), StTofrGeomCell(off)
34 //
35 // Usage:
36 //
37 // StTofrGeometry* geo = new StTofrGeometry("tofr","tofr geometry");
38 // // geo->Init(TVolume *starHall, const Int_t TofrConf);
39 // geo->Init(TVolume *starHall);
40 //
41 // ---------------------------------------------------------------------------
42 //
43 // StTofrNode
44 // ==============
45 //
47 
48 
49 Bool_t StTofrNode::mDebug = kFALSE;
50 Double_t const StTofrGeomSensor::mSensorDy = 10.35; // Actual module length;
51 const char* const StTofrGeometry::sectorPref = "BSEC";
52 const char* const StTofrGeometry::trayPref = "BTRA";
53 const char* const StTofrGeometry::senPref = "BRMD";
54 
55 //_____________________________________________________________________________
56 StTofrNode::StTofrNode(TVolumeView *element, TVolumeView *top)
57  : fView(element), pView(element->GetPosition()), mMasterNode(top), mTransFlag(kFALSE)
58 {
59  UpdateMatrix();
60  BuildMembers();
61 }
62 
63 //_____________________________________________________________________________
64 StTofrNode::~StTofrNode()
65 {
66  fView = 0;
67  pView = 0;
68  mMasterNode = 0;
69  mTransFlag = kFALSE;
70 }
71 
72 //_____________________________________________________________________________
73 void StTofrNode::UpdateMatrix()
74 {
75  //
76  //Update the Translation and RotateMatrix between local and master
77  // and store them as members: mTransMRS, mRotMRS
78  //while TNode stores them locally to share among all TNode objects,
79  // thus they may be changed afterward !
80  //
81 
82  mTransFlag = kFALSE;
83  CalcMatrix(this, mTransMRS, mRotMRS);
84  mTransFlag = kTRUE;
85 
86 }
87 
88 //_____________________________________________________________________________
89 void StTofrNode::CalcMatrix(StTofrNode* son,
90  Double_t* trans, Double_t* rot, StTofrNode* mother)
91 {
92  //
93  //Translation vector and Rotation matrix from TNode "mother" to "son"
94  //
95  // 1. Transformation: mother(x,y,z) ==> son(x',y',z')
96  //
97  // [x'] [ m[0] m[1] m[2] ] [ [x] [dx] ]
98  // [y'] = [ m[3] m[4] m[5] ] . [ [y] - [dy] ]
99  // [z'] [ m[6] m[7] m[8] ] [ [z] [dz] ]
100  //
101  // Transpose([x',y',z']) = M . Transpose([x-dx,y-dy,z-dz])
102  // and
103  // Transpose([x-dx,y-dy,z-dz]) = inverse(M) . Transpose([x',y',z'])
104  //
105  // where transpose(M) = inverse(M)
106  //
107  // 2. In the representation of vector:
108  //
109  // Vectors: V' = x' I' + y' J' + z' K'
110  // V = x I + y J + z K = V' + dV
111  // dV = dx I + dy J + dz K
112  //
113  // V' . I' = x' = (V-dV) . I'
114  // = (x-dx) (I.I') + (y-dy) (J.I') + (z-dz) (K.I')
115  // = (x-dx) m[0] + (y-dy) m[1] + (z-dz) m[2]
116  //
117  // so, the projections of unit vector I' on I, J, K
118  // yield m[0], m[1], m[2]
119  //
120  //********************************************************************
121 
122 
123  Double_t xl[3], xm[3];
124 
125  // son->TNode::UpdateMatrix();
126  // son->GetpView()->UpdateMatrix();
127  // LOG_INFO << "StTofrNode::CalcMatrix" << endm;
128 
129  xl[0] = 0; xl[1] = 0; xl[2] = 0;
130  ConvertPos(son,xl, mother,xm);
131  trans[0] = xm[0];
132  trans[1] = xm[1];
133  trans[2] = xm[2];
134 
135  xl[0] = 1; xl[1] = 0; xl[2] = 0;
136  ConvertPos(son,xl, mother,xm);
137  rot[0] = xm[0]-trans[0];
138  rot[1] = xm[1]-trans[1];
139  rot[2] = xm[2]-trans[2];
140 
141  xl[0] = 0; xl[1] = 1; xl[2] = 0;
142  ConvertPos(son,xl, mother,xm);
143  rot[3] = xm[0]-trans[0];
144  rot[4] = xm[1]-trans[1];
145  rot[5] = xm[2]-trans[2];
146 
147  xl[0] = 0; xl[1] = 0; xl[2] = 1;
148  ConvertPos(son,xl, mother,xm);
149  rot[6] = xm[0]-trans[0];
150  rot[7] = xm[1]-trans[1];
151  rot[8] = xm[2]-trans[2];
152 
153  return;
154 }
155 
156 //_____________________________________________________________________________
157 void StTofrNode::ConvertPos(
158  StTofrNode* from, const Double_t* pos_from,
159  StTofrNode* to, Double_t* pos_to)
160 {
161  if (to==0) from->Local2Master(pos_from,pos_to);
162  else {
163  Double_t xg[3];
164  from->Local2Master(pos_from,xg);
165  to->Master2Local(xg,pos_to);
166  }
167 
168  return;
169 }
170 
171 //_____________________________________________________________________________
172 void StTofrNode::BuildMembers()
173 {
174  //
175  //Build date members: mCenter{Rxy,Eta,Phi}, m{Eta,Phi}{Min,Max}
176  //
177 
178  Double_t xl[3];
179  Double_t xg[3];
180  TBRIK *brik = dynamic_cast<TBRIK*>(GetShape());
181  // LOG_INFO << " Get shape ready" << endm;
182  if(!brik) { LOG_INFO << " no brik " << endm; return; }
183  // Double_t dx = brik->GetDx();
184  Double_t dy = brik->GetDy();
185  Double_t dz = brik->GetDz();
186 
187  // LOG_INFO << " size = " <<dx << " " << dy << " " << dz << endm;
188  /*
189  //center point
190  Double_t xc[3];
191  xc[0] = 0; xc[1] = 0; xc[2] = 0;
192  Local2Master(xc, xg);
193  StThreeVectorD center(xg[0], xg[1], xg[2]);
194  mCenterRxy = center.perp();
195  mCenterEta = center.pseudoRapidity();
196  mCenterPhi = center.phi();
197  */
198 
199  // -dz
200  xl[0] = 0; xl[1] = 0; xl[2] = -dz;
201  Local2Master(xl, xg);
202  mEtaMin = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
203 
204  // +dz
205  xl[0] = 0; xl[1] = 0; xl[2] = dz;
206  Local2Master(xl, xg);
207  mEtaMax = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
208 
209  // -dy
210  xl[0] = 0; xl[1] = -dy; xl[2] = 0;
211  Local2Master(xl, xg);
212  mPhiMin = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
213 
214  // +dy
215  xl[0] = 0; xl[1] = dy; xl[2] = 0;
216  Local2Master(xl, xg);
217  mPhiMax = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
218 
219  if (mEtaMin>mEtaMax) {
220  Double_t tmp = mEtaMin;
221  mEtaMin = mEtaMax;
222  mEtaMax = tmp;
223  }
224 
225  if (mPhiMin>mPhiMax) {
226  Double_t tmp = mPhiMin;
227  mPhiMin = mPhiMax;
228  mPhiMax = tmp;
229  }
230 
231  return;
232 }
233 
234 //_____________________________________________________________________________
235 void StTofrNode::Local2Master(const Double_t* local, Double_t* master)
236 {
237  //
238  //Transform local coordinate into global coordinate
239  //
240  // pView->UpdateMatrix();
241  if (!mTransFlag) {
242  if (!mMasterNode) { LOG_INFO << " no Master! " << endm; return;}
243  TVolumeView *son = GetfView();
244  TVolumeView *mrs = GetTopNode();
245 
246  TVolumePosition *pos = 0;
247  pos = son->Local2Master(son, mrs);
248  pos->Local2Master(local, master);
249  delete pos;
250  return;
251  }
252 
253  //mTransFlag==kTRUE, i.e. StTofrGeomNode::UpdateMatrix() invoked already
254  Double_t x, y, z;
255  x = local[0];
256  y = local[1];
257  z = local[2];
258 
259  master[0] = mTransMRS[0] + mRotMRS[0]*x + mRotMRS[3]*y + mRotMRS[6]*z;
260  master[1] = mTransMRS[1] + mRotMRS[1]*x + mRotMRS[4]*y + mRotMRS[7]*z;
261  master[2] = mTransMRS[2] + mRotMRS[2]*x + mRotMRS[5]*y + mRotMRS[8]*z;
262 
263 }
264 
265 //_____________________________________________________________________________
266 void StTofrNode::Master2Local(const Double_t* master, Double_t* local)
267 {
268  //
269  //Transform global coordinate into local coordinate
270  //
271  // pView->UpdateMatrix();
272  // pView->Master2Local(master, local);
273  if (!mTransFlag) {
274  LOG_INFO << " and No TVolumePosition::Master2Local is wrong, so do nothing" << endm;
275  return;
276  }
277 
278  //mTransFlag==kTRUE, i.e. StTofrGeomNode::UpdateMatrix() invoked already
279  Double_t x, y, z;
280  x = master[0] - mTransMRS[0];
281  y = master[1] - mTransMRS[1];
282  z = master[2] - mTransMRS[2];
283 
284  local[0] = mRotMRS[0]*x + mRotMRS[1]*y + mRotMRS[2]*z;
285  local[1] = mRotMRS[3]*x + mRotMRS[4]*y + mRotMRS[5]*z;
286  local[2] = mRotMRS[6]*x + mRotMRS[7]*y + mRotMRS[8]*z;
287 }
288 
289 //_____________________________________________________________________________
290 StThreeVectorD StTofrNode::YZPlaneNormal()
291 {
292  //
293  //Calculate the vector unit of normal to YZ-plane
294  // i.e. the global representation of local unit vector (1,0,0)
295  //
296 
297  Double_t ux[3], nx[3];
298 
299  ux[0] = 1; ux[1] = 0; ux[2] = 0;
300  Local2Master(ux,nx);
301 
302  nx[0] -= mTransMRS[0];
303  nx[1] -= mTransMRS[1];
304  nx[2] -= mTransMRS[2];
305 
306  return StThreeVectorD(nx[0],nx[1],nx[2]);
307 }
308 
309 //_____________________________________________________________________________
310 StThreeVectorD StTofrNode::GetCenterPosition()
311 const
312 {
313  //
314  //Return the global representation of this node(TBRIk-shape) center
315  //
316 
317  Double_t xg[3];
318  xg[0] = mTransMRS[0];
319  xg[1] = mTransMRS[1];
320  xg[2] = mTransMRS[2];
321 
322  return StThreeVectorD(xg[0],xg[1],xg[2]);
323 }
324 
325 //_____________________________________________________________________________
326 Bool_t StTofrNode::IsLocalPointIn(const Double_t x, const Double_t y,
327  const Double_t z)
328 {
329  TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
330  Double_t dx = brik->GetDx();
331  Double_t dy = brik->GetDy();
332  Double_t dz = brik->GetDz();
333  Bool_t ret = -dx<x && x<dx && -dy<y && y<dy && -dz<z && z<dz;
334 
335  return ret;
336 }
337 
338 //_____________________________________________________________________________
339 Bool_t StTofrNode::IsGlobalPointIn(const StThreeVectorD &global)
340 {
341  Double_t xl[3], xg[3];
342  xg[0] = global.x();
343  xg[1] = global.y();
344  xg[2] = global.z();
345  Master2Local(xg, xl);
346  Bool_t ret = IsLocalPointIn(xl[0], xl[1], xl[2]);
347 
348  return ret;
349 }
350 
351 //_____________________________________________________________________________
352 Bool_t StTofrNode::HelixCross(const StHelixD &helix, Double_t &pathLen,
353  StThreeVectorD &cross)
354 {
355  //
356  // check if helix go through this node(TBRIK)
357  // and return the path length of helix before crossing this node
358  //
359  // static const Float_t MaxPathLength = 1000.; //Maximum path length
360  Float_t MaxPathLength = 1000.;
361 
362  Bool_t ret = kFALSE;
363  pathLen = 0;
364 
365  //
366  // Get the normal to the YZ-plane
367  //
368  StThreeVectorD planeNormal = YZPlaneNormal();
369 
370  //
371  // Get the center position
372  //
373  StThreeVectorD centerPos = GetCenterPosition();
374 
375  //
376  // Find the intersection point between the helix and the cell plane
377  //
378  pathLen = helix.pathLength(centerPos,planeNormal);
379  if ( pathLen>0 && pathLen<MaxPathLength ) {
380  cross = helix.at(pathLen);
381  //
382  // Check if the intersected point is really in the cell
383  //
384  ret = IsGlobalPointIn(cross);
385  }
386  return ret;
387 }
388 
389 //_____________________________________________________________________________
390 void StTofrNode::Print(Option_t *opt) const
391 {
392  TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
393  LOG_INFO << "Name=" << GetName() << "\t TBRIK-dimension=" << brik->GetDx()
394  << " : " << brik->GetDy() << " : " << brik->GetDz()
395  // << "\n Center Rxy:Eta:Phi=" << mCenterRxy << " : "
396  // << mCenterEta << " : " << mCenterPhi
397  << "\n EtaRange=" << mEtaMin << " : " << mEtaMax
398  << "\t PhiRange=" << mPhiMin << " : " << mPhiMax
399  << endm;
400  LOG_INFO <<"trans[0-2]=" << mTransMRS[0] <<" " <<mTransMRS[1]
401  <<" " <<mTransMRS[2]
402  <<"\nmRotMRS[0-2]=" <<mRotMRS[0] <<" " <<mRotMRS[1] <<" " <<mRotMRS[2]
403  <<"\nmRotMRS[3-5]=" <<mRotMRS[3] <<" " <<mRotMRS[4] <<" " <<mRotMRS[5]
404  <<"\nmRotMRS[6-8]=" <<mRotMRS[6] <<" " <<mRotMRS[7] <<" " <<mRotMRS[8]
405  <<endm;
406 }
407 
408 #if 0
409 // - - - - - -
410 //
411 // StTofrGeomNode
412 // ==============
413 /*
414 //_____________________________________________________________________________
415 StTofrGeomNode::StTofrGeomNode(const char* name, const char* title,
416  const TBRIK* brik, const Double_t x, const Double_t y, const Double_t z,
417  const Double_t theta)
418  : TNode(name, title, brik, x, y, z), mTransFlag(kFALSE)
419 {
420  fMatrix = new TRotMatrix("rot","rot",90.+theta,0.,90.,90.,theta,0.);
421  UpdateMatrix();
422  BuildMembers();
423 }
424 */
425 
426 
427 Bool_t StTofrGeomNode::mDebug = kFALSE;
428 //_____________________________________________________________________________
429 /*
430 StTofrGeomNode::StTofrGeomNode(const char* name, const char* title,
431  const TBRIK* brik, const Double_t x, const Double_t y, const Double_t z,
432  const TRotMatrix* matrix)*/
433 
434 StTofrGeomNode::StTofrGeomNode(const char* name, const char* title,
435  TBRIK* brik, const Double_t x, const Double_t y, const Double_t z,
436  TRotMatrix* matrix)
437  : TNode(name, title, brik, x, y, z, matrix), mTransFlag(kFALSE)
438 {
439  UpdateMatrix();
440  BuildMembers();
441 }
442 
443 //_____________________________________________________________________________
444 StTofrGeomNode::~StTofrGeomNode()
445 {
446 }
447 
448 //_____________________________________________________________________________
449 void StTofrGeomNode::UpdateMatrix()
450 {
451  //
452  //Update the Translation and RotateMatrix between local and master
453  // and store them as members: mTransMRS, mRotMRS
454  //while TNode stores them locally to share among all TNode objects,
455  // thus they may be changed afterward !
456  //
457 
458  mTransFlag = kFALSE;
459  CalcMatrix(this, mTransMRS, mRotMRS);
460  mTransFlag = kTRUE;
461 }
462 
463 //_____________________________________________________________________________
464 void StTofrGeomNode::CalcMatrix(TNode* son,
465  Double_t* trans, Double_t* rot, StTofrGeomNode* mother)
466 {
467  //
468  //Translation vector and Rotation matrix from TNode "mother" to "son"
469  //
470  // 1. Transformation: mother(x,y,z) ==> son(x',y',z')
471  //
472  // [x'] [ m[0] m[1] m[2] ] [ [x] [dx] ]
473  // [y'] = [ m[3] m[4] m[5] ] . [ [y] - [dy] ]
474  // [z'] [ m[6] m[7] m[8] ] [ [z] [dz] ]
475  //
476  // Transpose([x',y',z']) = M . Transpose([x-dx,y-dy,z-dz])
477  // and
478  // Transpose([x-dx,y-dy,z-dz]) = inverse(M) . Transpose([x',y',z'])
479  //
480  // where transpose(M) = inverse(M)
481  //
482  // 2. In the representation of vector:
483  //
484  // Vectors: V' = x' I' + y' J' + z' K'
485  // V = x I + y J + z K = V' + dV
486  // dV = dx I + dy J + dz K
487  //
488  // V' . I' = x' = (V-dV) . I'
489  // = (x-dx) (I.I') + (y-dy) (J.I') + (z-dz) (K.I')
490  // = (x-dx) m[0] + (y-dy) m[1] + (z-dz) m[2]
491  //
492  // so, the projections of unit vector I' on I, J, K
493  // yield m[0], m[1], m[2]
494  //
495  //********************************************************************
496 
497  Double_t xl[3], xm[3];
498 
499  son->TNode::UpdateMatrix();
500 
501  xl[0] = 0; xl[1] = 0; xl[2] = 0;
502  ConvertPos(son,xl, mother,xm);
503  trans[0] = xm[0];
504  trans[1] = xm[1];
505  trans[2] = xm[2];
506 
507  xl[0] = 1; xl[1] = 0; xl[2] = 0;
508  ConvertPos(son,xl, mother,xm);
509  rot[0] = xm[0]-trans[0];
510  rot[1] = xm[1]-trans[1];
511  rot[2] = xm[2]-trans[2];
512 
513  xl[0] = 0; xl[1] = 1; xl[2] = 0;
514  ConvertPos(son,xl, mother,xm);
515  rot[3] = xm[0]-trans[0];
516  rot[4] = xm[1]-trans[1];
517  rot[5] = xm[2]-trans[2];
518 
519  xl[0] = 0; xl[1] = 0; xl[2] = 1;
520  ConvertPos(son,xl, mother,xm);
521  rot[6] = xm[0]-trans[0];
522  rot[7] = xm[1]-trans[1];
523  rot[8] = xm[2]-trans[2];
524 
525  return;
526 }
527 
528 //_____________________________________________________________________________
529 void StTofrGeomNode::ConvertPos(
530  TNode* from, const Double_t* pos_from,
531  StTofrGeomNode* to, Double_t* pos_to)
532 {
533  //
534  //Convert the coordinate "Double_t* pos_from" in TNode* from
535  // into the coordinate "Doulbe_t *pos_to" in TNode* to
536  //
537 
538  if (to==0) from->Local2Master(pos_from,pos_to);
539  else {
540  Double_t xg[3];
541  from->Local2Master(pos_from,xg);
542  to->Master2Local(xg,pos_to);
543  }
544 
545  return;
546 }
547 
548 //_____________________________________________________________________________
549 void StTofrGeomNode::BuildMembers()
550 {
551  //
552  //Build date members: mCenter{Rxy,Eta,Phi}, m{Eta,Phi}{Min,Max}
553  //
554 
555  Double_t xl[3];
556  Double_t xg[3];
557  TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
558  //Double_t dx = brik->GetDx();
559  Double_t dy = brik->GetDy();
560  Double_t dz = brik->GetDz();
561 
562  /*
563  //center point
564  Double_t xc[3];
565  xc[0] = 0; xc[1] = 0; xc[2] = 0;
566  Local2Master(xc, xg);
567  StThreeVectorD center(xg[0], xg[1], xg[2]);
568  mCenterRxy = center.perp();
569  mCenterEta = center.pseudoRapidity();
570  mCenterPhi = center.phi();
571  */
572 
573  // -dz
574  xl[0] = 0; xl[1] = 0; xl[2] = -dz;
575  Local2Master(xl, xg);
576  mEtaMin = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
577 
578  // +dz
579  xl[0] = 0; xl[1] = 0; xl[2] = dz;
580  Local2Master(xl, xg);
581  mEtaMax = StThreeVectorD(xg[0],xg[1],xg[2]).pseudoRapidity();
582 
583  // -dy
584  xl[0] = 0; xl[1] = -dy; xl[2] = 0;
585  Local2Master(xl, xg);
586  mPhiMin = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
587 
588  // +dy
589  xl[0] = 0; xl[1] = dy; xl[2] = 0;
590  Local2Master(xl, xg);
591  mPhiMax = StThreeVectorD(xg[0],xg[1],xg[2]).phi();
592 
593  if (mEtaMin>mEtaMax) {
594  Double_t tmp = mEtaMin;
595  mEtaMin = mEtaMax;
596  mEtaMax = tmp;
597  }
598 
599  if (mPhiMin>mPhiMax) {
600  Double_t tmp = mPhiMin;
601  mPhiMin = mPhiMax;
602  mPhiMax = tmp;
603  }
604 
605 }
606 
607 //_____________________________________________________________________________
608 void StTofrGeomNode::Local2Master(const Double_t* local, Double_t* master)
609 {
610  //
611  //Transform local coordinate into global coordinate
612  //
613 
614  if (!mTransFlag) {
615  // TNode::UpdateMatrix();
616  TNode::Local2Master(local,master);
617  return;
618  }
619 
620  //mTransFlag==kTRUE, i.e. StTofrGeomNode::UpdateMatrix() invoked already
621  Double_t x, y, z;
622  x = local[0];
623  y = local[1];
624  z = local[2];
625 
626  master[0] = mTransMRS[0] + mRotMRS[0]*x + mRotMRS[3]*y + mRotMRS[6]*z;
627  master[1] = mTransMRS[1] + mRotMRS[1]*x + mRotMRS[4]*y + mRotMRS[7]*z;
628  master[2] = mTransMRS[2] + mRotMRS[2]*x + mRotMRS[5]*y + mRotMRS[8]*z;
629 
630 }
631 
632 //_____________________________________________________________________________
633 void StTofrGeomNode::Master2Local(const Double_t* master, Double_t* local)
634 {
635  //
636  //Transform global coordinate into local coordinate
637  // Please notice that TNode::Master2Local is INCORRECT !!!
638  //
639 
640  if (!mTransFlag) {
641  LOG_INFO << "Warning in StTofrGeomNode::Master2Local\n"
642  << " StTofrGeomNode::UpdateMatrix not yet invoked\n"
643  << " and TNode::Master2Local is wrong, so do nothing" << endm;
644  return;
645  }
646 
647  //mTransFlag==kTRUE, i.e. StTofrGeomNode::UpdateMatrix() invoked already
648  Double_t x, y, z;
649  x = master[0] - mTransMRS[0];
650  y = master[1] - mTransMRS[1];
651  z = master[2] - mTransMRS[2];
652 
653  local[0] = mRotMRS[0]*x + mRotMRS[1]*y + mRotMRS[2]*z;
654  local[1] = mRotMRS[3]*x + mRotMRS[4]*y + mRotMRS[5]*z;
655  local[2] = mRotMRS[6]*x + mRotMRS[7]*y + mRotMRS[8]*z;
656 
657 }
658 
659 //_____________________________________________________________________________
660 StThreeVectorD StTofrGeomNode::YZPlaneNormal()
661 {
662  //
663  //Calculate the vector unit of normal to YZ-plane
664  // i.e. the global representation of local unit vector (1,0,0)
665  //
666 
667  Double_t ux[3], nx[3];
668 
669  ux[0] = 1; ux[1] = 0; ux[2] = 0;
670  Local2Master(ux,nx);
671 
672  nx[0] -= mTransMRS[0];
673  nx[1] -= mTransMRS[1];
674  nx[2] -= mTransMRS[2];
675 
676  return StThreeVectorD(nx[0],nx[1],nx[2]);
677 }
678 
679 //_____________________________________________________________________________
680 StThreeVectorD StTofrGeomNode::GetCenterPosition()
681 const
682 {
683  //
684  //Return the global representation of this node(TBRIk-shape) center
685  //
686 
687  Double_t xg[3];
688  xg[0] = mTransMRS[0];
689  xg[1] = mTransMRS[1];
690  xg[2] = mTransMRS[2];
691 
692  return StThreeVectorD(xg[0],xg[1],xg[2]);
693 }
694 
695 //_____________________________________________________________________________
696 Bool_t StTofrGeomNode::IsLocalPointIn(const Double_t x, const Double_t y,
697  const Double_t z)
698 const
699 {
700  TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
701  Double_t dx = brik->GetDx();
702  Double_t dy = brik->GetDy();
703  Double_t dz = brik->GetDz();
704  Bool_t ret = -dx<x && x<dx && -dy<y && y<dy && -dz<z && z<dz;
705 
706  return ret;
707 }
708 
709 //_____________________________________________________________________________
710 Bool_t StTofrGeomNode::IsGlobalPointIn(const StThreeVectorD &global)
711 {
712  Double_t xl[3], xg[3];
713  xg[0] = global.x();
714  xg[1] = global.y();
715  xg[2] = global.z();
716  Master2Local(xg, xl);
717  Bool_t ret = IsLocalPointIn(xl[0], xl[1], xl[2]);
718 
719  return ret;
720 }
721 
722 //_____________________________________________________________________________
723 Bool_t StTofrGeomNode::HelixCross(const StHelixD &helix, Double_t &pathLen,
724  StThreeVectorD &cross)
725 {
726  //
727  // check if helix go through this node(TBRIK)
728  // and return the path length of helix before crossing this node
729  //
730 
731  static const Float_t MaxPathLength = 1000.; //Maximum path length
732 
733  Bool_t ret = kFALSE;
734  pathLen = 0;
735 
736  //
737  // Get the normal to the YZ-plane
738  //
739  StThreeVectorD planeNormal = YZPlaneNormal();
740 
741  //
742  // Get the center position
743  //
744  StThreeVectorD centerPos = GetCenterPosition();
745 
746  //
747  // Find the intersection point between the helix and the cell plane
748  //
749  pathLen = helix.pathLength(centerPos,planeNormal);
750  if ( pathLen>0 && pathLen<MaxPathLength ) {
751  cross = helix.at(pathLen);
752  //
753  // Check if the intersected point is really in the cell
754  //
755  ret = IsGlobalPointIn(cross);
756  }
757 
758  return ret;
759 }
760 
761 //_____________________________________________________________________________
762 void StTofrGeomNode::Print(Option_t *opt) const
763 const
764 {
765  TBRIK *brik = dynamic_cast<TBRIK*> (GetShape());
766  LOG_INFO << "Name=" << GetName() << "\t TBRIK-dimension=" << brik->GetDx()
767  << " : " << brik->GetDy() << " : " << brik->GetDz()
768  // << "\n Center Rxy:Eta:Phi=" << mCenterRxy << " : "
769  // << mCenterEta << " : " << mCenterPhi
770  << "\n EtaRange=" << mEtaMin << " : " << mEtaMax
771  << "\t PhiRange=" << mPhiMin << " : " << mPhiMax
772  << endm;
773  LOG_INFO <<"trans[0-2]=" << mTransMRS[0] <<" " <<mTransMRS[1]
774  <<" " <<mTransMRS[2]
775  <<"\nmRotMRS[0-2]=" <<mRotMRS[0] <<" " <<mRotMRS[1] <<" " <<mRotMRS[2]
776  <<"\nmRotMRS[3-5]=" <<mRotMRS[3] <<" " <<mRotMRS[4] <<" " <<mRotMRS[5]
777  <<"\nmRotMRS[6-8]=" <<mRotMRS[6] <<" " <<mRotMRS[7] <<" " <<mRotMRS[8]
778  <<endm;
779 }
780 #endif
781 
782 
784 //
785 // StTofrGeomTray
786 // ==============
787 //
789 
790 
791 Bool_t StTofrGeomTray::mDebug = kFALSE;
792 
793 //_____________________________________________________________________________
794 StTofrGeomTray::StTofrGeomTray(const Int_t ibtoh, TVolumeView *sector, TVolumeView *top)
795  : StTofrNode((TVolumeView *)sector->First(), top)
796 {
797  mSectorsInBTOH = top->GetListSize()/2;
798  mBTOHIndex = ibtoh + 1;
799  mTrayIndex = ibtoh * mSectorsInBTOH + sector->GetPosition()->GetId();
800 }
801 
802 //_____________________________________________________________________________
803 StTofrGeomTray::~StTofrGeomTray()
804 {
805  mBTOHIndex = 0;
806  mTrayIndex = 0;
807 }
808 //_____________________________________________________________________________
809 /*StTofrGeomTray::StTofrGeomTray(const char* name, const char* title,
810  TBRIK *brik, const Double_t x, const Double_t y, const Double_t z,
811  TRotMatrix *matrix, const Int_t itray)
812  : StTofrGeomNode(name,title,brik,x,y,z,matrix)
813 {
814  mTrayIndex = itray;
815 }
816 */
817 
818 //_____________________________________________________________________________
819 /*void StTofrGeomTray::PrepareCopyNode(TNode* node, StTofrGeomNode* top,
820  TShape* &shape, Double_t* pos, TRotMatrix* &newrot)
821 {
822  //
823  //Re-create/Clone components of a node to prepare for CopyNode or AddNode
824  //
825 
826  // static const char* where = "StTofrGeomTray::PrepareCopyNode";
827  // gDebugger.Print(mDebug,"==>Come to %s(...)\n",where);
828 
829  shape = 0;
830  newrot = 0;
831  if (node==0) return;
832 
833  THashList* shapeList = gGeometry->GetListOfShapes();
834  shape = dynamic_cast<TShape*>
835  ( shapeList->FindObject(node->GetShape()->GetName()) );
836 
837  if (shape == 0) {
838  //Clone the shape of the node
839  shape = dynamic_cast<TShape*> (node->GetShape()->Clone());
840  shapeList->Add(shape);
841  THashList* materialList = gGeometry->GetListOfMaterials();
842  materialList->Add(shape->GetMaterial());
843  }
844 
845  const char* name;
846  const char* title;
847  TRotMatrix *oldrot= node->GetMatrix();
848  name = oldrot->GetName();
849  title = oldrot->GetTitle();
850 
851  Double_t trans[3];
852  Double_t matrix[9];
853  CalcMatrix(node,trans,matrix,top);
854  newrot = new TRotMatrix(name,title,matrix);
855 
856  pos[0] = trans[0];
857  pos[1] = trans[1];
858  pos[2] = trans[2];
859 
860  // gDebugger.Print(mDebug," <== Leaving %s(...)\n",where);
861  return;
862 }
863 
864 //_____________________________________________________________________________
865 StTofrGeomTray* StTofrGeomTray::CopyNode(TNode* node, const Int_t itray)
866 {
867  //
868  //Clone-Convert "TNode* node" into "StTofrGeomTray*" under top TNode(MRS)
869  //
870 
871  // static const char* where = "StTofrGeomTray::CopyNode";
872  // gDebugger.Print(mDebug,"==>Come to %s(.,itray=%d)\n",where, itray);
873 
874  if (node==0) return NULL;
875 
876  TShape *shape;
877  Double_t trans[3];
878  TRotMatrix *newrot;
879 
880  PrepareCopyNode(node,0,shape,trans,newrot);
881 
882  const char* name;
883  const char* title;
884  TBRIK* brik = dynamic_cast<TBRIK*> (shape);
885  name = node->GetName();
886  title = node->GetTitle();
887 
888  //create StTofrGeomTray directly under gGeometry
889  // gGeometry->SetCurrentNode(0);
890 
891  StTofrGeomTray* tray = new StTofrGeomTray(name,title,brik,
892  trans[0],trans[1],trans[2],newrot, itray);
893 
894  // gDebugger.Print(mDebug," <== Leaving %s(.,itray=%d)\n",where, itray);
895  return tray;
896 }
897 
898 //_____________________________________________________________________________
899 StTofrGeomSensor* StTofrGeomTray::AddNode(TNode* node, const Int_t imodule)
900 {
901  //
902  //Clone-onvert "TNode* node" into "StTofrGeomSensor*" and put under this
903  // the translation/rotation will be built between "node" and this
904  //
905 
906  // static const char* where = "StTofrGeomTray::AddNode";
907  // gDebugger.Print(mDebug,"==>Come to %s(.,imodule=%d)\n",
908  // where, imodule);
909 
910  if (node==0) return NULL;
911 
912  const char* name;
913  const char* title;
914  TShape *shape;
915  Double_t trans[3];
916  TRotMatrix *newrot;
917 
918  PrepareCopyNode(node,this,shape,trans,newrot);
919  TBRIK* brik = dynamic_cast<TBRIK*> (shape);
920  name = node->GetName();
921  title = node->GetTitle();
922 
923  this->cd();
924  StTofrGeomSensor* sensor = new StTofrGeomSensor(name,title,brik,
925  trans[0],trans[1],trans[2],newrot, imodule);
926 
927 // gDebugger.Print(mDebug," <== Leaving %s(.,,imodule=%d)\n",where, imodule);
928  return sensor;
929 }12
930 */
931 //_____________________________________________________________________________
932 StTofrGeomSensor* StTofrGeomTray::GetSensor(const Int_t imodule) const
933 {
934 
935  StTofrGeomSensor* sensor = 0;
936 
937  TVolumeView *volume = GetfView();
938  // TList *list = fView->Nodes();
939  if ( !(volume->GetListSize()) ) {
940  LOG_INFO << " No Modules in this tray " << endm;
941  return sensor;
942  }
943  // list->Delete();
944 
945  TDataSetIter nextSensor(volume);
946  TVolumeView *sensorVolume = 0;
947  TVolumeView *top = GetTopNode();
948 
949  while ( (sensorVolume = (TVolumeView *)nextSensor()) )
950  {
951  if ( sensorVolume && (int)(sensorVolume->GetPosition()->GetId()) == imodule ) {
952  sensor = new StTofrGeomSensor(sensorVolume, top);
953  }
954  }
955 
956  return sensor;
957 }
958 
959 //_____________________________________________________________________________
960 void StTofrGeomTray::Print(const Option_t *opt) const
961 {
962  LOG_INFO << "StTofrGeomTray, tray#=" << mTrayIndex << endm;
963  StTofrNode::Print(opt);
964 }
965 
966 
968 //
969 // StTofrGeomSensor
970 // ================
971 //
973 
974 
975 Bool_t StTofrGeomSensor::mDebug = kFALSE;
976 
977 //_____________________________________________________________________________
978 StTofrGeomSensor::StTofrGeomSensor(TVolumeView *element, TVolumeView *top)
979  : StTofrNode(element, top)
980 {
981  mModuleIndex = element->GetPosition()->GetId();
982  CreateGeomCells();
983 }
984 
985 //_____________________________________________________________________________
986 StTofrGeomSensor::~StTofrGeomSensor()
987 {
988  mModuleIndex = 0;
989 }
990 //_____________________________________________________________________________
991 /*
992 StTofrGeomSensor::StTofrGeomSensor(const char* name, const char* title,
993  const TBRIK *brik, const Double_t x, const Double_t y, const Double_t z,
994  const TRotMatrix* matrix, const Int_t imodule)
995 StTofrGeomSensor::StTofrGeomSensor(const char* name, const char* title,
996  TBRIK *brik, const Double_t x, const Double_t y, const Double_t z,
997  TRotMatrix* matrix, const Int_t imodule)
998  : StTofrGeomNode(name,title,brik,x,y,z,matrix)
999 {
1000  mModuleIndex = imodule;
1001  CreateGeomCells();
1002 }
1003 */
1004 //_____________________________________________________________________________
1005 void StTofrGeomSensor::CreateGeomCells()
1006 {
1007  //
1008  //Divide this sensor to creat cells
1009  //
1010  /*
1011  TBRIK *thisBrik = dynamic_cast<TBRIK*> (GetShape());
1012  Double_t sensor_dy = thisBrik->GetDy();
1013  */
1014 
1015  //
1016  // change the size according to the real cells
1017  // mSensorDy -- defined by myself
1018  //
1019  Double_t sensor_dy = mSensorDy;
1020  Double_t cell_width = 2*sensor_dy/mCells;
1021 
1022  for (int i=0; i<=mCells; i++) mCellY[i] = cell_width*i - sensor_dy;
1023 
1024 }
1025 
1026 //_____________________________________________________________________________
1027 Double_t StTofrGeomSensor::GetCellYMin(const Int_t icell)
1028 const
1029 {
1030  Double_t ret = 0;
1031  if (icell<=0 || icell>mCells) {
1032  LOG_INFO << "cell#=" << icell <<" is out range=[0," << mCells << "]"
1033  << endm;
1034  } else ret = mCellY[icell-1];
1035 
1036  return ret;
1037 }
1038 
1039 //_____________________________________________________________________________
1040 Double_t StTofrGeomSensor::GetCellYMax(const Int_t icell)
1041 const
1042 {
1043  Double_t ret = 0;
1044  if (icell<=0 || icell>mCells) {
1045  LOG_INFO << "cell#=" << icell <<" is out range=[0," << mCells << "]"
1046  << endm;
1047  } else ret = mCellY[icell];
1048 
1049  return ret;
1050 }
1051 
1052 //_____________________________________________________________________________
1053 StThreeVectorD StTofrGeomSensor::GetCellPosition(const Int_t icell)
1054 {
1055  //
1056  // Get the center position of cell in this sensor
1057  //
1058 
1059  static const char* where = "StTofrGeomSensor::GetCellPosition";
1060  /*
1061  TBRIK *thisBrik = dynamic_cast<TBRIK*> (GetShape());
1062  Double_t sensor_dy = thisBrik->GetDy();
1063  */
1064 
1065  // change the size according to the real cells
1066  Double_t sensor_dy = mSensorDy;
1067  Double_t cell_dy = 2*sensor_dy/mCells;
1068 
1069  Double_t xl[3], xg[3];
1070  if (icell>=1 && icell<=mCells) {
1071  xl[0] = 0;
1072  xl[1] = (icell*2-1)*cell_dy - sensor_dy;
1073  xl[2] = 0;
1074  Local2Master(xl,xg);
1075  } else { //invalid cell, return (0.,0.,0.)
1076  LOG_INFO << "Warning in " << where <<" Invalid cell# =" << icell << endm;
1077  xg[0] = 0.;
1078  xg[1] = 0.;
1079  xg[2] = 0.;
1080  }
1081 
1082  return StThreeVectorD(xg[0],xg[1],xg[2]);
1083 }
1084 
1085 
1086 //_____________________________________________________________________________
1087 Int_t StTofrGeomSensor::FindCellIndex(const Double_t* local)
1088 {
1089  //
1090  //Look up the cell the local point in
1091  //
1092 
1093  Int_t icell=-1;
1094  if ( IsLocalPointIn(local[0],local[1],local[2]) ) {
1095 
1096  for (Int_t i=0; i<mCells; i++) {
1097  if (mCellY[i]<= local[1] && local[1]<=mCellY[i+1]) {
1098  icell = i+1;
1099  break;
1100 
1101  }
1102  }
1103  }
1104 
1105  return icell;
1106 }
1107 
1108 //_____________________________________________________________________________
1109 void StTofrGeomSensor::Print(const Option_t *opt) const
1110 {
1111  LOG_INFO << "StTofrGeomSensor, module#=" << mModuleIndex << endm;
1112  StTofrNode::Print(opt);
1113 
1114  LOG_INFO << " Cells=" << mCells << "\t Y range for cells=\n";
1115  for (Int_t i=0; i<=mCells; i++) LOG_INFO << " : " << mCellY[i];
1116  LOG_INFO << endm;
1117 }
1118 
1119 
1121 //
1122 // StTofrGeometry
1123 // ==============
1124 //
1126 
1127 StTofrGeometry *gTofrGeometry = 0;
1128 static const Int_t CELLSINMODULE = 6;
1129 
1130 
1131 Bool_t StTofrGeometry::mDebug = kFALSE;
1132 
1133 //_____________________________________________________________________________
1134 StTofrGeometry::StTofrGeometry(const char* name, const char* title)
1135  : TNamed(name,title)
1136 {
1137  mCellsInModule = StTofrGeomSensor::GetCells();
1138  mModulesInTray = 0;
1139  mTrays = 0;
1140  mRootFile = 0;
1141  mInitFlag = kFALSE;
1142  mTopNode = 0;
1143  mStarHall = 0;
1144 
1145  for(int i=0;i<mNTrays;i++) {
1146  mTofrTray[i] = 0;
1147  for(int j=0;j<mNModules;j++) {
1148  mTofrSensor[i][j] = 0;
1149  }
1150  }
1151 
1152  //
1153  //We only need one instance of StTofrGeometry
1154  //
1155  if (gTofrGeometry) {
1156  LOG_INFO << "Warning !! There is already StTofrGeometry at pointer="
1157  << (void*)gTofrGeometry << ", so it is deleted"
1158  << endm;
1159  delete gTofrGeometry;
1160  }
1161  gTofrGeometry = this;
1162 }
1163 
1164 //_____________________________________________________________________________
1165 StTofrGeometry::~StTofrGeometry()
1166 {
1167  LOG_INFO << "Warning !! StTofrGeometry at pointer =" << (void*)gTofrGeometry
1168  << " is deleted" << endm;
1169  gTofrGeometry = 0;
1170 
1171  for(int i=0;i<mNTrays;i++) {
1172  if(mTofrTray[i]) delete mTofrTray[i];
1173  mTofrTray[i] = 0;
1174  for(int j=0;j<mNModules;j++) {
1175  if(mTofrSensor[i][j]) delete mTofrSensor[i][j];
1176  mTofrSensor[i][j] = 0;
1177  }
1178  }
1179 
1180 }
1181 
1182 //_____________________________________________________________________________
1183 //void StTofrGeometry::Init(const char *file, Option_t *option)
1184 //void StTofrGeometry::Init(TVolume *starHall, const Int_t TofrConf)
1185 void StTofrGeometry::Init(TVolume *starHall)
1186 {
1187  //
1188  //Define geometry parameters and establish the geometry
1189  // current available options: "root", and "xdf"
1190  //
1191 
1192  // InitFromStar(starHall, TofrConf);
1193  InitFromStar(starHall);
1194  mStarHall = starHall;
1195 
1196  //save current gGeometry to resume later
1197  //
1198  //TGeometry *currentGeo = gGeometry;
1199 
1200  //static TString optRoot("root");
1201  //static TString optXdf("xdf");
1202 
1203  //TString opt(option);
1204  //opt.ToLower();
1205 
1206  //if (opt==optRoot) {
1207  // //
1208  // //Initialization from ROOT file of geometry (converted from RZ file)
1209  // //
1210  // InitFromRoot(file);
1211  // mRootFile = file;
1212  //} else if (opt==optXdf) {
1213  // //
1214  // //Initialize from XDF file, not ready
1215  // //
1216  // //InitFromXdf(file);
1217  // //Xdf2Geometry();
1218  // LOG_INFO << "StTofrGeom::Init, sorry! option \"xdf\" not yet ready" << endm;
1219  //} else {
1220  // LOG_INFO << "StTofrGeom::Init, Warning!! not yet implemented option="
1221  // << option << endm;
1222  //}
1223 
1226  //gGeometry = currentGeo;
1227 
1230  //mInitFlag = kTRUE;
1231 
1232 }
1233 //_____________________________________________________________________________
1234 //void StTofrGeometry::InitFromStar(TVolume *starHall, const Int_t TofrConf=0)
1235 void StTofrGeometry::InitFromStar(TVolume *starHall)
1236 {
1237  // Initialize TOFr geometry from STAR geometry
1238  // TofrConf -- 0 tray_Tofr (default)
1239  // 1 full_Tofr
1240 
1241  // TVolume *starHall = (TVolume *) GetDataSet("HALL");
1242  // mStarHall = starHall;
1243 
1244  // mTofrConf = TofrConf;
1245 // static const char* sectorPref="BSEC";
1246 // static const char* trayPref="BTRA";
1247 // static const char* senPref ="BRMD";
1248 // static const char* topName ="HALL1";
1249 
1250 // static const char* tofElements[] = { sectorPref, trayPref, senPref, topName };
1251 
1252  // Loop over the STAR geometry and mark the volume needed
1253  TDataSetIter volume(starHall,0);
1254 
1255  TVolume *starDetectorElement = 0;
1256  while ( (starDetectorElement = ( TVolume *)volume()) )
1257  {
1258  // const char *elementName = starDetectorElement->GetName();
1259  // Bool_t found = ! ( strcmp(elementName,tofElements[0]) && strcmp(elementName,tofElements[1]) && strcmp(elementName,tofElements[2]) );
1260  Bool_t found = ( IsBSEC(starDetectorElement) || IsBTRA(starDetectorElement) || IsBRMD(starDetectorElement) );
1261  if (found) {
1262 
1263  starDetectorElement->SetVisibility(TVolume::kBothVisible);
1264  starDetectorElement->Mark();
1265  if (starDetectorElement->GetLineColor()==1 || starDetectorElement->GetLineColor()==7)
1266  starDetectorElement->SetLineColor(14);
1267 
1268  } else {
1269 
1270  starDetectorElement->UnMark();
1271  starDetectorElement->SetVisibility(TVolume::kThisUnvisible);
1272 
1273  }
1274  }
1275 
1276  starHall->SetVisibility(TVolume::kBothVisible);
1277  mTopNode = new TVolumeView(*starHall,10);
1278 
1279  // mTrays = 120;
1280  // mModulesInTray = 33;
1281 
1282  mSectorsInBTOH = mTopNode->GetListSize()/2; // # of sectors in one half
1283 
1284  // check tray-tofr or full-tofr
1285  TDataSetIter nextSector(mTopNode);
1286  TVolumeView *secVolume = 0;
1287  mTrays = 0; // non-emtry tray number
1288  while ( (secVolume = (TVolumeView *)nextSector()) ) {
1289  TVolumeView *trayVolume = (TVolumeView *)secVolume->First();
1290  if ( trayVolume->GetListSize() ) {
1291  mTrays++;
1292  mModulesInTray = trayVolume->GetListSize();
1293  }
1294  }
1295  mTofrConf = 0;
1296  if (mTrays==120) mTofrConf = 1;
1297  //
1298 
1300  // save the sensors and trays
1302  gMessMgr->Info("","OS") << " # of trays = " << mTopNode->GetListSize() << endm;
1303  TList *list = mTopNode->Nodes();
1304  Int_t ibtoh =0;
1305  TVolumeView *sectorVolume = 0;
1306  mNValidTrays = 0;
1307  mNValidModules = 0;
1308  for(Int_t i=0;i<list->GetSize();i++) {
1309  sectorVolume = dynamic_cast<TVolumeView*> (list->At(i));
1310  TVolumeView *trayVolume = (TVolumeView *)sectorVolume->First();
1311  if( !trayVolume->GetListSize() ) continue;
1312  if ( i>=60 ) ibtoh = 1;
1313  // gMessMgr->Info("","OS") << " test sector size = " << trayVolume->GetListSize() << endm;
1314  mTofrTray[mNValidTrays] = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1315  TList *list1 = trayVolume->Nodes();
1316  // gMessMgr->Info("","OS") << " # of modules in tray " << mTofrTray[mNValidTrays]->Index() << " = " << trayVolume->GetListSize() << endm;
1317  if (!list1 ) continue;
1318  TVolumeView *sensorVolume = 0;
1319  if(list1->GetSize()>mNValidModules) mNValidModules=list1->GetSize();
1320  for(Int_t j=0;j<list1->GetSize();j++) {
1321  sensorVolume = dynamic_cast<TVolumeView*> (list1->At(j));
1322  mTofrSensor[mNValidTrays][j] = new StTofrGeomSensor(sensorVolume, mTopNode);
1323  }
1324  mNValidTrays++;
1325  }
1326  gMessMgr->Info("","OS") << "\n-------------------------------------------\n"
1327  << " Summary of initialization: "
1328  << " NValidTrays = " << mNValidTrays << " NValidModules = " << mNValidModules << endm;
1329 
1330 
1331 // if ( !TofrConf ) { // tray Tofr --- delete garbages
1332 // // We should not delete the volumes , just keep them
1333 // /*
1334 // TList garbage;
1335 // TDataSetIter nextSector(mTopNode);
1336 // TVolumeView *sector = 0;
1337 // Int_t ibtoh = 0, i = 0;
1338 // while ( (sector = (TVolumeView *)nextSector()) ) {
1339 // i++;
1340 // if ( i>mSectorsInBTOH ) ibtoh = 1;
1341 // Int_t itray = mSectorsInBTOH * ibtoh + sector->GetPosition()->GetId();
1342 // if ( itray!=83 ) garbage.Add(sector);
1343 // }
1344 // garbage.Delete();
1345 // */
1346 // mTrays = 1;
1347 // mModulesInTray = 33;
1348 // }
1349 
1350  return;
1351 
1352 }
1353 
1354 //_____________________________________________________________________________
1355 Bool_t StTofrGeometry::ContainOthers(TVolume *element)
1356 {
1357  TVolumeView *elementView = new TVolumeView(*element,1);
1358  TList *list = elementView->GetListOfShapes();
1359  if ( list ) {
1360  LOG_INFO << " yes list in " << element->GetName() << endm;
1361  return kTRUE;
1362  } else {
1363  LOG_INFO << " no list in " << element->GetName() << endm;
1364  return kFALSE;
1365  }
1366 }
1367 
1368 /*
1369 //_____________________________________________________________________________
1370 Bool_t StTofrGeometry::InitFromRoot(const char* geofile)
1371 {
1372  //
1373  //Read geometry root file and setup own simple StTofrGeometry frame
1374  // // The geometry root file can be created in the following procedure:
1375  //
1376  // 1. staf> make cavegeo
1377  // > make btofgeo
1378  // > rz/file 80 btofr_geo.rz on
1379  //
1380  // 2. shell> g2root btofr_geo.rz btofr_geo.root
1381  //
1382  //-------------------------------------------------------------------
1383 
1384  static const char* trayPref="BTR";
1385  static const char* senPref ="BRMD";
1386  static const char* topName ="HALL1";
1387 
1388  TObject *obj;
1389  TFile file(geofile);
1390  if (! file.IsOpen() ) return kFALSE;
1391 
1392  // Get the geometry in the root file
1393  //
1394  TList *list = file.GetListOfKeys();
1395  TObject *key;
1396  TGeometry *btofr_geo;
1397  TIter next(list);
1398  while ( (key=next()) != 0 ) {
1399  obj = file.Get(key->GetName());
1400  if ( strcmp(obj->ClassName(),"TGeometry") == 0) {
1401  btofr_geo = dynamic_cast<TGeometry*> (obj);
1402  break;
1403  }
1404  }
1405  if (!btofr_geo) return kFALSE;
1406  file.Close();
1407 
1408  TNode *topNode = btofr_geo->GetNode(topName);
1409  if (!topNode) return kFALSE;
1410 
1411  //Create the simplified geometry "top(MRS)->Tray->Module(Sensor)->Cell"
1412  // under this Geometry
1413  //
1414  this->cd();
1415  CopyTopNode(topNode);
1416 
1417  //Create StTofrGeomTrays and its descedenets(StTofrGeomSensors)
1418  //
1419  TList trayList;
1420  TList senList;
1421  GetPrefixNodes(topNode,trayPref,trayList);
1422  Int_t itray = 0 ;
1423  Int_t ntray = 0;
1424  Int_t imodule;
1425  TIter trayNext(&trayList);
1426  TNode* node;
1427  while ( (obj=trayNext()) != 0 ) {
1428  itray++;
1429 
1430  // if (itray>=10 && itray<=20) break; //Skip tray=[10,20] for yesw-Test
1431 
1432  node = dynamic_cast<TNode*> (obj);
1433  senList.Clear();
1434  GetPrefixNodes(node,senPref,senList);
1435  //gDebugger.Print(mDebug,"%s, in itray=%d found module=%d\n",where,
1436  // itray, senList.GetSize() );
1437 
1438  if (senList.GetSize() > 0) {
1439  ntray++;
1440  //
1441  //Create StTofrGeomTray under TopNode
1442  //
1443  mTopNode->cd();
1444  StTofrGeomTray *tray = StTofrGeomTray::CopyNode(node,itray);
1445  //gDebugger.Print(mDebug,"%s, itray=%d copied as %d tray\n",where,
1446  // itray, ntray);
1447  if (!tray) {
1448  LOG_INFO << "StTofrGeometry::InitFromRoot, Warning!!"
1449  << "failed in CopyNode=" <<node->GetName() << endm;
1450  break;
1451  }
1452  tray->cd();
1453  TIter senNext(&senList);
1454  TObject* obj2;
1455  TNode* node2;
1456  imodule = 0;
1457  while ( (obj2=senNext()) != 0 ) {
1458  imodule++;
1459  node2 = dynamic_cast<TNode*> (obj2);
1460  //gDebugger.Print(mDebug,"%s, adding sensor %d name=%s\n",where,
1461  //imodule, node2->GetName() );
1462  //
1463  //Create StTofrGeomSensor under current GeomTray
1464  //
1465  tray->AddNode(node2,imodule);
1466  }
1467  if (mModulesInTray==0 && imodule!=0) mModulesInTray=imodule;
1468  }
1469 
1470  }
1471  if (mTrays==0 && ntray!=0) mTrays=ntray;
1472 
1473  delete btofr_geo;
1474 
1475  return kTRUE;
1476 }
1477 */
1478 /*
1479 //_____________________________________________________________________________
1480 Bool_t StTofrGeometry::CopyTopNode(TNode* top)
1481 {
1482  //
1483  //Create the top TNode under StTofrGeometry
1484  //
1485 
1486  //check if it is a top TNode
1487  //
1488  if (top==0) return kFALSE;
1489  if (top->GetParent() != 0) return kFALSE;
1490 
1491  TShape *shape;
1492  Double_t trans[3];
1493  TRotMatrix *newrot;
1494  StTofrGeomTray::PrepareCopyNode(top,0,shape,trans,newrot);
1495 
1496  const char* name;
1497  const char* title;
1498  name = top->GetName();
1499  title = top->GetTitle();
1500 
1501  //
1502  // right under StTofrGeometry
1503  //
1504  gGeometry->SetCurrentNode(0);
1505  mTopNode = new TNode(name,title,shape,trans[0],trans[1],trans[2]);
1506 
1507  return kTRUE;
1508 }
1509 
1510 //_____________________________________________________________________________
1511 void StTofrGeometry::GetPrefixNodes(const TNode* topNode,
1512  const char* key, TList &list)
1513 {
1514  //
1515  // Recusively loop over all descendent nodes and find out those
1516  // name prefixed with key
1517  //
1518 
1519  static Int_t counts = 0;
1520 
1521  if (!topNode || !key) return;
1522  if (strlen(key)==0) return;
1523 
1524  TList *sons = topNode->GetListOfNodes();
1525  if (sons==0) return;
1526 
1527  TIter next(sons); //create class TIter for iteration over sons
1528  TObject* obj;
1529  TNode* node;
1530  const char *name;
1531 
1532  while ( (obj=next()) !=0 ) {
1533  node = dynamic_cast<TNode*> (obj);
1534  name = node->GetName();
1535  if (strstr(name,key)==name) {
1536  counts++;
1537  list.Add(node);
1538  } else {
1539  GetPrefixNodes(node,key,list);
1540  }
1541  }
1542 
1543 
1544 }
1545 */
1546 
1547 //_____________________________________________________________________________
1548 Bool_t StTofrGeometry::LackThis(const char* fromWhere)
1549 {
1550  if (gTofrGeometry == 0) {
1551  LOG_INFO << " !! Warning from " << fromWhere
1552  << "\n no StTofrGeometry existing, create one instance first"
1553  << endm;
1554  return kTRUE;
1555  } else return kFALSE;
1556 }
1557 
1558 /*
1559 //_____________________________________________________________________________
1560 TRotMatrix* StTofrGeometry::CreateMatrix(const Double_t theta)
1561 {
1562  return new TRotMatrix("rot","rot",90.+theta,0.,90.,90.,theta,0.);
1563 }
1564 */
1565 
1566 //_____________________________________________________________________________
1567 Int_t StTofrGeometry::CalcCellId(const Int_t volumeId, const Float_t* local)
1568 const
1569 {
1570  Double_t dlocal[3];
1571  dlocal[0] = local[0];
1572  dlocal[1] = local[1];
1573  dlocal[2] = local[2];
1574  return CalcCellId(volumeId,dlocal);
1575 }
1576 
1577 //_____________________________________________________________________________
1578 Int_t StTofrGeometry::CalcCellId(const Int_t volumeId, const Double_t* local)
1579 const
1580 {
1581  //
1582  // Calculate cellID based on volumeId and local position in a sensor
1583  //
1584 
1585  Int_t icell, imodule, itray;
1586  DecodeVolumeId(volumeId, imodule, itray);
1587  StTofrGeomSensor *sensor = GetGeomSensor(imodule, itray);
1588  if (!sensor) return -1;
1589  icell = sensor->FindCellIndex(local);
1590  Int_t ret = CalcCellId(icell, imodule, itray);
1591 
1592  return ret;
1593 }
1594 
1595 //_____________________________________________________________________________
1596 Bool_t StTofrGeometry::IsCellValid(const Int_t icell)
1597 const
1598 {
1599  //
1600  //Check the validity of cell# = [1,mCellsInModule]
1601  //
1602  return (icell>=1 && icell<=mCellsInModule);
1603 }
1604 
1605 //_____________________________________________________________________________
1606 Bool_t StTofrGeometry::IsSensorValid(const Int_t imodule)
1607 const
1608 {
1609  //
1610  //Check the validity of module# = [1,mModulesInTray]
1611  //
1612  return (imodule>=1 && imodule<=mModulesInTray);
1613 }
1614 
1615 //_____________________________________________________________________________
1616 Bool_t StTofrGeometry::IsTrayValid(const Int_t itray)
1617 const
1618 {
1619  //
1620  //Check the validity of tray#
1621  // return kTRUE if found in the tray list
1622  //
1623 
1624  Bool_t ret = kFALSE;
1625  // TList *list = mTopNode->GetListOfNodes();
1626  TDataSetIter nextSector(mTopNode);
1627  TVolumeView* sectorVolume = 0;
1628  Int_t ibtoh = 0, i = 0;
1629  // StTofrGeomTray *tray = 0;
1630  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
1631  i++;
1632  if ( i>mSectorsInBTOH ) ibtoh = 1;
1633  int trayIndex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
1634  // tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1635  // if (tray->Index() == itray) {
1636  if (trayIndex == itray) {
1637  ret = kTRUE;
1638  break;
1639  }
1640  // delete tray;
1641  }
1642 
1643  /* TList *list = mTopNode->Nodes();
1644  if (!list || mTrays==0) return ret;
1645 
1646  StTofrGeomTray* tray;
1647  for (Int_t i=0; i<list->GetSize(); i++) {
1648  tray = dynamic_cast<StTofrGeomTray*> (list->At(i));
1649  if (tray->Index() == itray) {
1650  ret = kTRUE;
1651  break;
1652  }
1653  }
1654  */
1655 
1656  return ret;
1657 }
1658 
1659 //_____________________________________________________________________________
1660 Int_t StTofrGeometry::CalcSensorId(const Int_t imodule, const Int_t itray)
1661 const
1662 {
1663  //
1664  //Encode imodule and itray into SensorId
1665  //
1666 
1667  Int_t sensorId = -1;
1668  if (!IsSensorValid(imodule)) return sensorId;
1669 
1670  Int_t idx = GetAtOfTray(itray);
1671 
1672  if (idx<0) return sensorId;
1673 
1674  sensorId = imodule-1 + mModulesInTray*idx;
1675 
1676  return sensorId;
1677 }
1678 
1679 //_____________________________________________________________________________
1680 Int_t StTofrGeometry::PrevCellId(const Int_t cellId)
1681 const
1682 {
1683  //
1684  //Look up the cell prior to cellId in the same module
1685  // and return the cellId of cell found
1686  //
1687 
1688  Int_t found = -1;
1689  Int_t icell, imodule, itray;
1690  DecodeCellId(cellId,icell,imodule,itray);
1691  StTofrGeomSensor* sensor = GetGeomSensor(imodule,itray);
1692  Int_t icell_p = sensor->PrevCellIndex(icell);
1693 
1694  found = CalcCellId(icell_p,imodule,itray);
1695 
1696  return found;
1697 }
1698 
1699 //_____________________________________________________________________________
1700 Int_t StTofrGeometry::NextCellId(const Int_t cellId)
1701 const
1702 {
1703  //
1704  //Look up the cell prior to cellId in the same module
1705  // and return the cellId of cell found
1706  //
1707 
1708  Int_t found = -1;
1709  Int_t icell, imodule, itray;
1710  DecodeCellId(cellId,icell,imodule,itray);
1711  StTofrGeomSensor* sensor = GetGeomSensor(imodule,itray);
1712  Int_t icell_p = sensor->NextCellIndex(icell);
1713 
1714  found = CalcCellId(icell_p,imodule,itray);
1715 
1716  return found;
1717 }
1718 
1719 //_____________________________________________________________________________
1720 Int_t StTofrGeometry::CalcCellId(const Int_t icell, const Int_t imodule,
1721  const Int_t itray)
1722 const
1723 {
1724  //
1725  //Encode icell, imodule and itray into CellId
1726  //
1727 
1728  Int_t cellId = -1;
1729  if (!IsCellValid(icell)) return cellId;
1730 
1731  Int_t sensorId = CalcSensorId(imodule,itray);
1732  if (sensorId<0) return cellId; //Invalid sensorId
1733 
1734  cellId = icell-1 + mCellsInModule*sensorId;
1735 
1736  return cellId;
1737 }
1738 
1739 //_____________________________________________________________________________
1740 Bool_t StTofrGeometry::DecodeSensorId(const Int_t sensorId,
1741  Int_t &imodule, Int_t &itray)
1742 const
1743 {
1744  //
1745  //decode sensorId into tray#, module# in itray, imodule
1746  //
1747 
1748  imodule = sensorId%mModulesInTray + 1;
1749  if (!IsSensorValid(imodule)) return kFALSE;
1750 
1751  Int_t idx = sensorId/mModulesInTray;
1752  StTofrGeomTray *tray = GetGeomTrayAt(idx);
1753  if (!tray) return kFALSE;
1754  itray = tray->Index();
1755  delete tray;
1756  return kTRUE;
1757 }
1758 
1759 //_____________________________________________________________________________
1760 Bool_t StTofrGeometry::DecodeCellId(const Int_t cellId, Int_t &icell,
1761  Int_t &imodule, Int_t &itray)
1762 const
1763 {
1764  //
1765  //decode cellId into tray#, module#, cell# in itray, imodule, icell
1766  //
1767 
1768  Int_t sensorId = cellId/mCellsInModule;
1769  if (!DecodeSensorId(sensorId,imodule,itray)) return kFALSE;
1770 
1771  icell = cellId%mCellsInModule + 1;
1772 
1773  return IsCellValid(icell);
1774 }
1775 
1776 //_____________________________________________________________________________
1777 Int_t StTofrGeometry::GetCellIndex(const Int_t cellId)
1778 const
1779 {
1780  //
1781  //decode cellId and return cell#
1782  //
1783 
1784  Int_t icell = cellId%mCellsInModule + 1;
1785 
1786  return icell;
1787 }
1788 
1789 //_____________________________________________________________________________
1790 void StTofrGeometry::DecodeVolumeId(const Int_t volumeId,
1791  Int_t &imodule, Int_t &itray)
1792 const
1793 {
1794  //
1795  //decode the volumeId into tray# and module# in itray, imodule
1796  // see the definition of TOFr's volumeId in "pams/sim/g2t/g2t_volume_id.g"
1797  //
1798 
1799  Int_t ires = volumeId;
1800 
1801  Int_t rileft = int(ires/10/100/100);
1802  ires = ires-rileft*100*100*10;
1803 
1804  itray = int(ires/10/100);
1805  ires = ires-itray*100*10;
1806 
1807  imodule = int(ires/10);
1808  // Int_t icell = ires%10; //always 0 since not defined in geant geometry
1809 
1810  itray = itray + (rileft-1)*mSectorsInBTOH;
1811 
1812  return;
1813 }
1814 
1815 //_____________________________________________________________________________
1816 StTofrGeomSensor* StTofrGeometry::GetGeomCell(const Int_t cellId)
1817 const
1818 {
1819  //
1820  //Return StTofrGeomSensor* where the cell of cellId is in
1821  //
1822 
1823  Int_t icell, imodule, itray;
1824  DecodeCellId(cellId, icell, imodule, itray);
1825  StTofrGeomSensor* sensor = GetGeomSensor(imodule, itray);
1826 
1827  return sensor;
1828 }
1829 
1830 //_____________________________________________________________________________
1831 StTofrGeomSensor* StTofrGeometry::GetGeomSensor(const Int_t imodule,
1832  const Int_t itray)
1833 const
1834 {
1835  //
1836  //itray is dummy if itray==0 and it is the current single tray
1837  //
1838 
1839  StTofrGeomTray *tray = GetGeomTray(itray);
1840  StTofrGeomSensor* sensor = NULL;
1841  if (tray) sensor = tray->GetSensor(imodule);
1842  delete tray;
1843  return sensor;
1844 }
1845 
1846 //_____________________________________________________________________________
1847 StTofrGeomTray* StTofrGeometry::GetGeomTray(const Int_t itray)
1848 const
1849 {
1850  //
1851  //itray is dummy if itray==0 and it is the current single tray
1852  //
1853 
1854  StTofrGeomTray* found = 0;
1855  // StTofrGeomTray* tray = 0;
1856  TDataSetIter nextSector(mTopNode);
1857  TVolumeView *sectorVolume = 0;
1858  Int_t ibtoh = 0, i = 0;
1859  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
1860  i++;
1861  if (i>mSectorsInBTOH ) ibtoh = 1;
1862  int trayIndex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
1863  // tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1864  // if (tray->Index() == itray) {
1865  if (trayIndex == itray) {
1866  // found = tray;
1867  found = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1868  break;
1869  }
1870  // delete tray;
1871  }
1872 
1873  /*
1874  // TList *list = mTopNode->GetListOfNodes();
1875  TList *list = mTopNode->Nodes();
1876  if (!list || mTrays==0) return NULL;
1877 
1878  //case of single tray, itray is dummy
1879  //
1880  if (mTrays==1 && itray==0) {
1881  found = dynamic_cast<StTofrGeomTray*> (list->At(0));
1882  return found;
1883  }
1884 
1885  //case of multiple trays
1886  //
1887  TIter next(list);
1888  TObject *obj;
1889  StTofrGeomTray* tray;
1890  while ( (obj=next()) != 0 ) {
1891  tray = dynamic_cast<StTofrGeomTray*> (obj);
1892  if (tray->Index() == itray) {
1893  found = tray;
1894  break;
1895  }
1896  }
1897  */
1898 
1899  return found;
1900 }
1901 
1902 //_____________________________________________________________________________
1903 StTofrGeomTray* StTofrGeometry::GetGeomTrayAt(const Int_t idx)
1904 const
1905 {
1906  //
1907  //Get the StTofrGeomTray at index of the list
1908  //
1909 
1910  StTofrGeomTray* found = 0;
1911 
1912  TDataSetIter nextSector(mTopNode);
1913  TVolumeView *sectorVolume = 0;
1914  Int_t ibtoh = 0, i = 0;
1915  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
1916  i++;
1917  if ( i>mSectorsInBTOH ) ibtoh = 1;
1918  if (i==idx) {
1919  found = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1920  }
1921  }
1922 
1923  // TList *list = mTopNode->GetListOfNodes();
1924  /*
1925  TList *list = mTopNode->Nodes();
1926  if (!list || mTrays==0) return NULL;
1927 
1928  if (idx<0) return found; //bug in list->At(i) if i<0
1929 
1930  TVolumeView *trayVolume = dynamic_cast<StTofrGeomTray*> (list->At(idx));
1931  found = new StTofrGeomTray(trayVolume);
1932  */
1933  // found = dynamic_cast<StTofrGeomTray*> (list->At(idx));
1934 
1935  return found;
1936 }
1937 
1938 //_____________________________________________________________________________
1939 Int_t StTofrGeometry::GetAtOfTray(const Int_t itray)
1940 const
1941 {
1942  //
1943  //Find out the list-index of StTofrGeomTray with TrayIndex=itray
1944  // itray is dummy if itray==0 and it is the current single tray
1945  //
1946 
1947  Int_t at = -1;
1948 
1949  // TList *list = mTopNode->GetListOfNodes();
1950 
1951  TDataSetIter nextSector(mTopNode);
1952  TVolumeView *sectorVolume = 0;
1953  // StTofrGeomTray *tray = 0;
1954  Int_t ibtoh = 0, i = 0;
1955  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
1956  i++;
1957  if ( i>mSectorsInBTOH ) ibtoh = 1;
1958  int trayIndex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
1959  // tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
1960  // if (tray->Index() == itray) {
1961  if (trayIndex == itray) {
1962  at = i;
1963  break;
1964  }
1965  // delete tray;
1966  }
1967 
1968  /*
1969  TList *list = mTopNode->Nodes();
1970  if (!list || mTrays==0) return at;
1971 
1972  //case of single tray, itray is dummy
1973  //
1974  if (mTrays==1 && itray==0) {
1975  if (list->At(0) != 0) return 0;
1976  else return at;
1977  }
1978 
1979  //case of multiple trays
1980  //
1981  StTofrGeomTray* tray;
1982  for (Int_t i=0; i<list->GetSize(); i++) {
1983  tray = dynamic_cast<StTofrGeomTray*> (list->At(i));
1984  if (tray->Index() == itray) {
1985  at = i;
1986  break;
1987  }
1988  }
1989  */
1990 
1991  return at;
1992 }
1993 
1994 //_____________________________________________________________________________
1995 void StTofrGeometry::Print(Option_t *opt) const
1996 {
1997  LOG_INFO << "Trays=" << mTrays <<"\t ModulesInTray=" << mModulesInTray
1998  << "\t CellsInModule=" << mCellsInModule << endm;
1999 }
2000 
2001 //_____________________________________________________________________________
2002 Int_t StTofrGeometry::CellIdPointIn(const StThreeVectorD& point)
2003 const
2004 {
2005  //
2006  //Return the ID of cell which containing the global point
2007  //
2008 
2009  Int_t cellId = -1;
2010  Double_t xl[3], xg[3];
2011  xg[0] = point.x();
2012  xg[1] = point.y();
2013  xg[2] = point.z();
2014 
2015  //Loop over trays
2016  //
2017 
2018  Int_t itray = -1, imodule = -1, icell = -1;
2019  for(int i=0;i<mNValidTrays;i++) {
2020  if(!mTofrTray[i]) continue;
2021  if ( mTofrTray[i]->IsGlobalPointIn(point) ) {
2022  itray = mTofrTray[i]->Index();
2023  if ( !(mTofrTray[i]->GetfView()->GetListSize()) ) {
2024  LOG_INFO << " No sensors in tray " << itray << endm;
2025  return cellId;
2026  }
2027 
2028  for( int j=0;j<mNValidModules;j++) {
2029  if(!mTofrSensor[i][j]) continue;
2030  if ( mTofrSensor[i][j]->IsGlobalPointIn(point) ) {
2031  imodule = mTofrSensor[i][j]->Index();
2032  mTofrSensor[i][j]->Master2Local(xg,xl);
2033  icell = mTofrSensor[i][j]->FindCellIndex(xl);
2034  }
2035  } // end for (j)
2036  } // end if
2037  } // end for (i)
2038 
2039  /*
2040  TDataSetIter nextSector(mTopNode);
2041  TVolumeView *sectorVolume = 0;
2042  // StTofrGeomTray *tray = 0;
2043 
2044  Int_t itray = -1, imodule = -1, icell = -1;
2045  Int_t ibtoh = 0, i = 0;
2046  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
2047  i++;
2048  if ( i>mSectorsInBTOH ) ibtoh = 1;
2049  StTofrGeomTray *tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
2050  if ( tray->IsGlobalPointIn(point) ) {
2051  itray = tray->Index();
2052 
2053  if ( !(tray->GetfView()->GetListSize()) ) {
2054  LOG_INFO << " No sensors in tray " << tray->Index() << endm;
2055  delete tray;
2056  return cellId;
2057  }
2058 
2059  TDataSetIter nextSensor(tray->GetfView());
2060  TVolumeView *sensorVolume = 0;
2061  // StTofrGeomSensor *sensor = 0;
2062  while ( (sensorVolume = (TVolumeView *)nextSensor()) ) {
2063  StTofrGeomSensor *sensor = new StTofrGeomSensor(sensorVolume, mTopNode);
2064  if ( sensor->IsGlobalPointIn(point) ) {
2065  imodule = sensor->Index();
2066  sensor->Master2Local(xg,xl);
2067  icell = sensor->FindCellIndex(xl);
2068  }
2069  delete sensor;
2070  }
2071  }
2072  delete tray;
2073  }
2074  */
2075 
2076  if ( itray <= 0 || imodule <= 0 ) return cellId;
2077  cellId = CalcCellId(icell, imodule, itray);
2078  return cellId;
2079  /*
2080  StTofrGeomTray* tray;
2081  Int_t itray = -1;
2082  // TIter next( mTopNode->GetListOfNodes() );
2083  TIter next( mTopNode->Nodes() );
2084  while ( (obj=next()) != 0 ) {
2085  tray = dynamic_cast<StTofrGeomTray*> (obj);
2086  if ( tray->IsGlobalPointIn(point) ) {
2087  itray = tray->Index();
2088  break;
2089  }
2090  }
2091  if (itray <= 0) return cellId;
2092 
2093  //Loop over sensors in a tray
2094  //
2095  // TIter next2( tray->GetListOfNodes() );
2096  TIter next2( tray->GetfView()->Nodes() );
2097  StTofrGeomSensor *sensor;
2098  Int_t imodule = -1;
2099  while ( (obj=next2()) != 0 ) {
2100  sensor = dynamic_cast<StTofrGeomSensor*> (obj);
2101  if ( sensor->IsGlobalPointIn(point) ) {
2102  imodule = sensor->Index();
2103  break;
2104  }
2105  }
2106  if (imodule <= 0) return cellId;
2107  */
2108  // Find cell# in a sensor
2109  /* sensor->Master2Local(xg,xl);
2110  Int_t icell = sensor->FindCellIndex(xl);
2111 
2112  cellId = CalcCellId(icell, imodule, itray);
2113 
2114  return cellId;
2115  */
2116 }
2117 
2118 //_____________________________________________________________________________
2119 Bool_t StTofrGeometry::HelixCross(const StHelixD &helix)
2120 const
2121 {
2122  //
2123  // return "kTRUE" if any cell is crossed by this helix
2124  //
2125 
2126  IntVec idVec;
2127  DoubleVec pathVec;
2128  PointVec crossVec;
2129 
2130  Bool_t crossed = HelixCrossCellIds(helix,idVec,pathVec,crossVec);
2131 
2132  return crossed;
2133 }
2134 
2135 //_____________________________________________________________________________
2136 Bool_t StTofrGeometry::HelixCrossCellIds(const StHelixD &helix,
2137  IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec)
2138 const
2139 {
2140  //
2141  // return "kTRUE" if any cell is crossed by this helix
2142  // and also fill the cellIds which are crossed by the helix
2143  // and the path length of the helix before crossing the cell
2144  //
2145 
2146  Double_t pathLen;
2147  Int_t cellId;
2148  StThreeVectorD cross;
2149  idVec.clear();
2150  pathVec.clear();
2151  crossVec.clear();
2152 
2153  for(int i=0;i<mNValidTrays;i++) {
2154  if(!mTofrTray[i]) continue;
2155  int trayId = mTofrTray[i]->Index();
2156 
2157  for(int j=0;j<mNValidModules;j++) {
2158  if(!mTofrSensor[i][j]) continue;
2159  int moduleId = mTofrSensor[i][j]->Index();
2160  if ( mTofrSensor[i][j]->HelixCross(helix,pathLen,cross) ) {
2161  Double_t global[3], local[3];
2162  global[0] = cross.x();
2163  global[1] = cross.y();
2164  global[2] = cross.z();
2165  mTofrSensor[i][j]->Master2Local(global,local);
2166  Int_t icell = mTofrSensor[i][j]->FindCellIndex(local);
2167  cellId = CalcCellId(icell, moduleId, trayId);
2168  if (cellId>=0) { // a bug before // reject hit in the edge of module;
2169  pathVec.push_back(pathLen);
2170  idVec.push_back(cellId);
2171  crossVec.push_back(cross);
2172  }
2173  }
2174  } // end for (j)
2175  } // end for (i)
2176 
2177 
2178  /*
2179  TDataSetIter nextSector(mTopNode);
2180  TVolumeView *sectorVolume = 0;
2181  // StTofrGeomTray *tray = 0;
2182 
2183  Int_t ibtoh = 0, i = 0;
2184  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
2185  i++;
2186  if ( i>mSectorsInBTOH ) ibtoh = 1;
2187  //
2188  // For tray-tofr, select tray #83 directly to save time
2189  //
2190  if ( !mTofrConf ) {
2191  int trayindex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
2192  if ( trayindex!=mY03TrayIndex ) {
2193  // cout << " skip tray " << trayindex << endl;
2194  continue;
2195  }
2196  }
2197  //
2198  StTofrGeomTray *tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
2199  // TVolumeView *trayVolume = tray->GetfView();
2200  if ( tray->HelixCross(helix,pathLen,cross) ) {
2201 
2202  cout << " Helix cross tray # = " << tray->Index() << endl;
2203  if (!tray->GetfView()->GetListSize()) {
2204  delete tray;
2205  return kFALSE;
2206  }
2207 
2208  TDataSetIter nextSensor(tray->GetfView());
2209 
2210  TVolumeView *sensorVolume = 0;
2211 
2212  while ( (sensorVolume = (TVolumeView *)nextSensor()) ) {
2213  StTofrGeomSensor *sensor = new StTofrGeomSensor(sensorVolume, mTopNode);
2214  if ( sensor->HelixCross(helix,pathLen,cross) ) {
2215 
2216  // cout << " hit the sensor volume " << endl;
2217  Double_t global[3], local[3];
2218  global[0] = cross.x();
2219  global[1] = cross.y();
2220  global[2] = cross.z();
2221  sensor->Master2Local(global,local);
2222  Int_t icell = sensor->FindCellIndex(local);
2223  cellId = CalcCellId(icell, sensor->Index(), tray->Index());
2224  if (cellId>0) { // reject hit in the edge of module;
2225  pathVec.push_back(pathLen);
2226  idVec.push_back(cellId);
2227  crossVec.push_back(cross);
2228  }
2229  }
2230  delete sensor;
2231  }
2232 
2233  } // end if
2234  delete tray;
2235 
2236  }
2237  */
2238  /*
2239  TObject *obj;
2240  StTofrNode *node;
2241 
2242  // find out the tray crossed by the helix, and one tray can be crossed
2243  // but there is totally only one tray of TOFr at this moment
2244  StTofrGeomTray *tray=NULL;
2245  // TList *list = mTopNode->GetListOfNodes();
2246  TList *list = mTopNode->Nodes();
2247  TIter next(list);
2248  while ( (obj=next()) != 0) {
2249  node = dynamic_cast<StTofrNode*> (obj);
2250  if ( node->HelixCross(helix,pathLen,cross) ) {
2251  tray = dynamic_cast<StTofrGeomTray*> (obj);
2252  break;
2253  }
2254  }
2255 
2256  if (!tray) return kFALSE;
2257 
2258  // find out modules crossed by the helix
2259  // and next cells further
2260  //
2261  StTofrGeomSensor *sensor;
2262  // list = tray->GetListOfNodes();
2263  list = tray->GetfView()->Nodes();
2264  TIter next2(list);
2265  while ( (obj=next2()) != 0 ) {
2266  node = dynamic_cast<StTofrNode*> (obj);
2267  if ( node->HelixCross(helix,pathLen,cross) ) { //module crossed
2268  sensor = dynamic_cast<StTofrGeomSensor*> (obj);
2269 
2270  // Method-1: just Look up the cell the point of cross is in
2271  //
2272  Double_t global[3], local[3];
2273  global[0] = cross.x();
2274  global[1] = cross.y();
2275  global[2] = cross.z();
2276  sensor->Master2Local(global,local);
2277  Int_t icell = sensor->FindCellIndex(local);
2278  cellId = CalcCellId(icell, sensor->Index(), tray->Index());
2279  pathVec.push_back(pathLen);
2280  idVec.push_back(cellId);
2281  crossVec.push_back(cross);
2282  }
2283  }
2284  */
2285 
2286  if (idVec.size()>0) {
2287  // mleak1.PrintMem(" normal return true");
2288  return kTRUE;
2289  }
2290  else {
2291  // mleak1.PrintMem(" normal return false");
2292  return kFALSE;
2293  }
2294 }
2295 
2296 
2297 //_____________________________________________________________________________
2298 Bool_t StTofrGeometry::HelixCross(const StHelixD &helix, IntVec validModuleVec, IntVec projTrayVec)
2299 const
2300 {
2301  //
2302  // return "kTRUE" if any cell is crossed by this helix
2303  //
2304 
2305  IntVec idVec;
2306  DoubleVec pathVec;
2307  PointVec crossVec;
2308 
2309  Bool_t crossed = HelixCrossCellIds(helix,validModuleVec, projTrayVec, idVec,pathVec,crossVec);
2310 
2311  return crossed;
2312 }
2313 
2314 
2315 //_____________________________________________________________________________
2316 Bool_t StTofrGeometry::HelixCrossCellIds(const StHelixD &helix, IntVec validModuleVec, IntVec projTrayVec, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec)
2317 const
2318 {
2320  // optimized :
2321  // input : helix, validModuleVec from DAQ,
2322  // projTrayVec possible hit tray from this helix
2323  // Xin Dong
2325  //
2326  // return "kTRUE" if any cell is crossed by this helix
2327  // and also fill the cellIds which are crossed by the helix
2328  // and the path length of the helix before crossing the cell
2329  //
2330 
2331  if(validModuleVec.size()==0) return kFALSE;
2332  if(projTrayVec.size()==0) return kFALSE;
2333 
2334  Double_t pathLen;
2335  Int_t cellId;
2336  StThreeVectorD cross;
2337  idVec.clear();
2338  pathVec.clear();
2339  crossVec.clear();
2340 
2341  for(int i=0;i<mNValidTrays;i++) {
2342  if(!mTofrTray[i]) continue;
2343  int trayId = mTofrTray[i]->Index();
2344  bool itrayFind = kFALSE;
2345 
2346  for(size_t it=0;it<projTrayVec.size();it++) {
2347  int validtrayId = projTrayVec[it];
2348  if(validtrayId==trayId) {
2349  itrayFind = kTRUE;
2350  break;
2351  }
2352  }
2353  if(!itrayFind) continue;
2354 
2355  // LOG_INFO << " Helix cross sensitive tray " << trayId << endm;
2356 
2357  for(int j=0;j<mNValidModules;j++) {
2358  if(!mTofrSensor[i][j]) continue;
2359  int moduleId = mTofrSensor[i][j]->Index();
2360  for(size_t iv=0;iv<validModuleVec.size();iv++) {
2361  int validtrayId = validModuleVec[iv]/100;
2362  int validmoduleId = validModuleVec[iv]%100;
2363  if(validtrayId==trayId&&validmoduleId==moduleId) {
2364  if ( mTofrSensor[i][j]->HelixCross(helix,pathLen,cross) ) {
2365  Double_t global[3], local[3];
2366  global[0] = cross.x();
2367  global[1] = cross.y();
2368  global[2] = cross.z();
2369  mTofrSensor[i][j]->Master2Local(global,local);
2370  Int_t icell = mTofrSensor[i][j]->FindCellIndex(local);
2371  cellId = CalcCellId(icell, moduleId, trayId);
2372  if (cellId>=0) { // a bug before // reject hit in the edge of module;
2373  pathVec.push_back(pathLen);
2374  idVec.push_back(cellId);
2375  crossVec.push_back(cross);
2376  }
2377  }
2378  } // endif (tray && module)
2379  } // end for (iv)
2380  } // end for (j)
2381  } // end for (i)
2382 
2383  /*
2384  TDataSetIter nextSector(mTopNode);
2385  TVolumeView *sectorVolume = 0;
2386  // StTofrGeomTray *tray = 0;
2387 
2388  Int_t ibtoh = 0, i = 0;
2389  while ( (sectorVolume = (TVolumeView *)nextSector()) ) {
2390  i++;
2391  if ( i>mSectorsInBTOH ) ibtoh = 1;
2392  //
2393  // For tray-tofr, select tray #83 directly to save time
2394  //
2395  if ( !mTofrConf ) {
2396  int trayindex = ibtoh * mSectorsInBTOH + sectorVolume->GetPosition()->GetId();
2397  if ( trayindex!=mY03TrayIndex ) {
2398  // LOG_INFO << " skip tray " << trayindex << endm;
2399  continue;
2400  }
2401  }
2402  //
2403  StTofrGeomTray *tray = new StTofrGeomTray(ibtoh, sectorVolume, mTopNode);
2404  // TVolumeView *trayVolume = tray->GetfView();
2405  if ( tray->HelixCross(helix,pathLen,cross) ) {
2406 
2407  LOG_INFO << " Helix cross tray # = " << tray->Index() << endm;
2408  if (!tray->GetfView()->GetListSize()) {
2409  delete tray;
2410  return kFALSE;
2411  }
2412 
2413  TDataSetIter nextSensor(tray->GetfView());
2414 
2415  TVolumeView *sensorVolume = 0;
2416 
2417  while ( (sensorVolume = (TVolumeView *)nextSensor()) ) {
2418  StTofrGeomSensor *sensor = new StTofrGeomSensor(sensorVolume, mTopNode);
2419  if ( sensor->HelixCross(helix,pathLen,cross) ) {
2420 
2421  // LOG_INFO << " hit the sensor volume " << endm;
2422  Double_t global[3], local[3];
2423  global[0] = cross.x();
2424  global[1] = cross.y();
2425  global[2] = cross.z();
2426  sensor->Master2Local(global,local);
2427  Int_t icell = sensor->FindCellIndex(local);
2428  cellId = CalcCellId(icell, sensor->Index(), tray->Index());
2429  if (cellId>0) { // reject hit in the edge of module;
2430  pathVec.push_back(pathLen);
2431  idVec.push_back(cellId);
2432  crossVec.push_back(cross);
2433  }
2434  }
2435  delete sensor;
2436  }
2437 
2438  } // end if
2439  delete tray;
2440 
2441  }
2442  */
2443  /*
2444  TObject *obj;
2445  StTofrNode *node;
2446 
2447  // find out the tray crossed by the helix, and one tray can be crossed
2448  // but there is totally only one tray of TOFr at this moment
2449  StTofrGeomTray *tray=NULL;
2450  // TList *list = mTopNode->GetListOfNodes();
2451  TList *list = mTopNode->Nodes();
2452  TIter next(list);
2453  while ( (obj=next()) != 0) {
2454  node = dynamic_cast<StTofrNode*> (obj);
2455  if ( node->HelixCross(helix,pathLen,cross) ) {
2456  tray = dynamic_cast<StTofrGeomTray*> (obj);
2457  break;
2458  }
2459  }
2460 
2461  if (!tray) return kFALSE;
2462 
2463  // find out modules crossed by the helix
2464  // and next cells further
2465  //
2466  StTofrGeomSensor *sensor;
2467  // list = tray->GetListOfNodes();
2468  list = tray->GetfView()->Nodes();
2469  TIter next2(list);
2470  while ( (obj=next2()) != 0 ) {
2471  node = dynamic_cast<StTofrNode*> (obj);
2472  if ( node->HelixCross(helix,pathLen,cross) ) { //module crossed
2473  sensor = dynamic_cast<StTofrGeomSensor*> (obj);
2474 
2475  // Method-1: just Look up the cell the point of cross is in
2476  //
2477  Double_t global[3], local[3];
2478  global[0] = cross.x();
2479  global[1] = cross.y();
2480  global[2] = cross.z();
2481  sensor->Master2Local(global,local);
2482  Int_t icell = sensor->FindCellIndex(local);
2483  cellId = CalcCellId(icell, sensor->Index(), tray->Index());
2484  pathVec.push_back(pathLen);
2485  idVec.push_back(cellId);
2486  crossVec.push_back(cross);
2487  }
2488  }
2489  */
2490 
2491  if (idVec.size()>0) {
2492  // mleak1.PrintMem(" normal return true");
2493  return kTRUE;
2494  }
2495  else {
2496  // mleak1.PrintMem(" normal return false");
2497  return kFALSE;
2498  }
2499 }
2500 
2501 //---------------------------------------------------------------------------
2502 // estimate the possible projection on the TOF tray
2503 Bool_t StTofrGeometry::projTrayVector(const StHelixD &helix, IntVec &trayVec) const {
2504 
2505  trayVec.clear();
2506  double R_tof = 215.;
2507  double res = 5.0;
2508  double s1 = helix.pathLength(R_tof).first;
2509  if(s1<0.) s1 = helix.pathLength(R_tof).second;
2510  StThreeVectorD point = helix.at(s1);
2511  double phi = point.phi()*180/3.14159;
2512 
2513  // east ring, start from 108 deg (id=61) , clock-wise from east facing west
2514  int itray_east = (255+(int)phi)%360/6+61;
2515  trayVec.push_back(itray_east);
2516 
2517  int itray_east1 = (255+(int)(phi+res))%360/6+61;
2518  int itray_east2 = (255+(int)(phi-res))%360/6+61;
2519  if(itray_east1!=itray_east) {
2520  trayVec.push_back(itray_east1);
2521  }
2522  if(itray_east2!=itray_east&&itray_east2!=itray_east1) {
2523  trayVec.push_back(itray_east2);
2524  }
2525 
2526  // west ring, start from 72 deg (id=1) , clock-wise from west facing east
2527  int itray_west = (435-(int)phi)%360/6+1;
2528  trayVec.push_back(itray_west);
2529 
2530  int itray_west1 = (435-(int)(phi+res))%360/6+1;
2531  int itray_west2 = (435-(int)(phi-res))%360/6+1;
2532  if(itray_west1!=itray_west) {
2533  trayVec.push_back(itray_west1);
2534  }
2535  if(itray_west2!=itray_west&&itray_west2!=itray_west1) {
2536  trayVec.push_back(itray_west2);
2537  }
2538 
2539 // LOG_INFO << " proj tray id = ";
2540 // for(size_t it=0;it<trayVec.size();it++) {
2541 // LOG_INFO << trayVec[it] << " ";
2542 // }
2543 // LOG_INFO << endm;
2544 
2545  if(trayVec.size()>0) return kTRUE;
2546  else return kFALSE;
2547 }
2548 
2549 /*
2550 //_____________________________________________________________________________
2551 //
2552 // non-member function
2553 // ==========
2554 //_____________________________________________________________________________
2555 Bool_t lackTofrgeometry(const char* fromWhere)
2556 {
2557  if (gTofrGeometry == 0) {
2558  LOG_INFO << " !! Warning from " << fromWhere
2559  << "\n no StTofrGeometry existing, create one instance first"
2560  << endm;
2561  return kTRUE;
2562  } else return kFALSE;
2563 }
2564 */
2565 
2566 /*******************************************************************
2567  * $Log: StTofrGeometry.cxx,v $
2568  * Revision 1.12 2018/02/26 23:26:51 smirnovd
2569  * StTof: Remove outdated ClassImp macro
2570  *
2571  * Revision 1.11 2018/02/26 23:13:21 smirnovd
2572  * Move embedded CVS log messages to the end of file
2573  *
2574  * Revision 1.10 2009/01/26 15:05:33 fisyak
2575  * rename TMemStat => StMemStat due to clash with ROOT class
2576  *
2577  * Revision 1.9 2008/03/27 00:15:39 dongx
2578  * Update for Run8 finished.
2579  *
2580  * Revision 1.8 2007/04/17 23:01:52 dongx
2581  * replaced with standard STAR Loggers
2582  *
2583  * Revision 1.7 2007/02/28 23:28:17 dongx
2584  * R_tof used for pre-match calculation updated to ~215cm
2585  *
2586  * Revision 1.6 2004/05/03 23:07:49 dongx
2587  * -Introduce data members to save the Tray and Sensor geometries in the initialization.
2588  * -Optimize the HelixCrossCellIds() function to save CPU time
2589  * -Introduce a new function projTrayVector()
2590  * -Update the ClassDef number 1->2
2591  *
2592  *
2593  * Revision 1.4 2004/03/09 16:45:16 dongx
2594  * Remove InitDaqMap() since a StTofrDaqMap is introduced
2595  *
2596  * Revision 1.3 2003/09/11 05:49:23 perev
2597  * ansi corrs
2598  *
2599  * Revision 1.2 2003/09/07 03:49:06 perev
2600  * gcc 3.2 + WarnOff
2601  *
2602  * Revision 1.1 2003/08/06 23:00:53 geurts
2603  * First Release
2604  */
TDataSet::First
virtual TDataSet * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TDataSet.cxx:403
StTofrGeometry::sectorPref
static const char *const sectorPref
Control message printing of this class.
Definition: StTofrGeometry.h:365
StTofrGeomTray::StTofrGeomTray
StTofrGeomTray(const Int_t ibtoh, TVolumeView *sector, TVolumeView *top)
Control message printing of this class.
Definition: StTofrGeometry.cxx:794
StTofrGeomSensor::CreateGeomCells
void CreateGeomCells()
Control message printing of this class.
Definition: StTofrGeometry.cxx:1005
TVolume::SetVisibility
virtual void SetVisibility(ENodeSEEN vis=TVolume::kBothVisible)
Definition: TVolume.cxx:753
TVolumeView::Local2Master
virtual TVolumePosition * Local2Master(const TVolumeView *localNode, const TVolumeView *masterNode=0)
to be documented
Definition: TVolumeView.cxx:675
TVolumePosition::Local2Master
virtual Double_t * Local2Master(const Double_t *local, Double_t *master, Int_t nPoints=1) const
Definition: TVolumePosition.cxx:352
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StTofrGeometry::mSectorsInBTOH
Int_t mSectorsInBTOH
the root file of geometry
Definition: StTofrGeometry.h:348
StTofrNode::StTofrNode
StTofrNode(TVolumeView *element, TVolumeView *top)
Control message printing of this class.
Definition: StTofrGeometry.cxx:56
Generated by   doxygen 1.8.5