StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMtdGeometry.cxx
1 /********************************************************************
2  * $Id: StMtdGeometry.cxx,v 1.20 2018/04/25 15:58:35 marr Exp $
3  ********************************************************************
4  *
5  * $Log: StMtdGeometry.cxx,v $
6  * Revision 1.20 2018/04/25 15:58:35 marr
7  * Fix the bug of converting projected localy to cellId
8  *
9  * Revision 1.19 2017/03/08 20:40:38 marr
10  * Add back the old implementation of GetCellLocalYCenter() function to make
11  * the class backward compatible.
12  *
13  * Revision 1.18 2017/02/13 02:56:11 marr
14  * From 2017, do not move BL 8&24 along y direction by hand since this is already
15  * done in the geometry file. Calibration, production and analysis should use
16  * the new version consistently.
17  *
18  * Revision 1.17 2016/08/05 16:12:34 marr
19  * Add MTD hit IdTruth to avoid applying dy shift for BL 8 and 24 for MC hits
20  *
21  * Revision 1.16 2015/07/29 14:52:47 smirnovd
22  * Added class scope to static members (overlooked in previous commit)
23  *
24  * Revision 1.15 2015/07/29 01:11:25 smirnovd
25  * Initialize static constants outside of class definition
26  *
27  * C++ forbids initialization of non-integral static const members within the class
28  * definition. The syntax is allowed only for integral type variables.
29  *
30  * Revision 1.14 2015/07/24 15:56:05 marr
31  * 1. Remove calling a macro in Init() to create geometry. It should be done within
32  * the maker that uses this utility class.
33  * 2. Add the TGeoManager parameter to the default constructor to force the existance
34  * of the gometry when using this utility class.
35  * 3. Simplify the code for getting the pointer to the magnetic field
36  *
37  * Revision 1.13 2015/05/01 01:55:34 marr
38  * Fix the geometry of shifted backleg 8 and 24
39  *
40  * Revision 1.12 2015/04/07 16:23:33 marr
41  * 1. Make use the constants defined in StMtdConstants.h
42  * 2. Cleaning up
43  *
44  * Revision 1.11 2015/02/09 21:29:23 marr
45  * Fix an overlook in extracting geometry for 2015 during bfc chain running
46  *
47  * Revision 1.10 2014/12/23 16:32:33 marr
48  * In 2015 geometry file, MTD modules are placed under directory MagRefSys_1.
49  * Modification is maded to reflect this change.
50  *
51  * Revision 1.9 2014/09/24 18:14:47 marr
52  * Flip localY for cells in modules 4 and 5
53  *
54  * Revision 1.8 2014/09/18 21:59:08 marr
55  * 1. Disable retrieving table geant2backlegID for year 2012 ana before
56  * 2. The default geometry tag is set to yYYYYa
57  *
58  * Revision 1.7 2014/07/25 19:44:18 marr
59  * Fix a minor inconsistency in using the fNExtraCells
60  *
61  * Revision 1.6 2014/07/24 17:02:30 huangbc
62  * Add protection for reading magnetic field in case of track projection position is (nan,nan,nan).
63  *
64  * Revision 1.5 2014/07/16 15:31:01 huangbc
65  * Add an option to lock bfield to FF.
66  *
67  * Revision 1.4 2014/07/10 20:45:13 huangbc
68  * New geometry class for MTD, load geometry from geant geometry. Need gGeoManager.
69  *
70  * Revision 1.3 2013/08/07 18:25:30 geurts
71  * include CVS Id and Log tags
72  *
73  *
74  *******************************************************************/
75 #include <vector>
76 #include <assert.h>
77 #include <stdlib.h>
78 #include "StMtdGeometry.h"
79 #include "StMaker.h"
80 #include "StMessMgr.h"
81 #include "TGeoManager.h"
82 #include "TGeoVolume.h"
83 #include "TGeoBBox.h"
84 #include "TGeoNode.h"
85 #include "TROOT.h"
86 #include "StThreeVectorD.hh"
87 #include "StPhysicalHelixD.hh"
88 #include "PhysicalConstants.h"
89 #include "SystemOfUnits.h" // has "tesla" in it
90 #include "StDetectorDbMaker/St_MagFactorC.h"
91 #include "tables/St_mtdGeant2BacklegIDMap_Table.h"
92 
93 
94 const Double_t StMtdGeometry::mMtdMinR = 392.802;
95 const Double_t StMtdGeometry::mMtdMaxR = 418.865;
96 const Double_t StMtdGeometry::mMagInR = 303.290;
97 const Double_t StMtdGeometry::mMagOutR = 364.290;
98 const Double_t StMtdGeometry::mEmcInR = 223.505;
99 const Double_t StMtdGeometry::mEmcOutR = 248.742;
100 
101 const char* StMtdGeometry::backlegPref[4] = {"MTMT","MTMF","MTTG","MTT1"};
102 const char* StMtdGeometry::modulePref = "MTRA";
103 //const char* StMtdGeometry::sensorPref = "MIGG";
104 
105 //----------------------------------------------------//
106 // //
107 // StMtdGeoNode //
108 // //
109 //----------------------------------------------------//
110 
111 //#ifdef __ROOT__
112 //ClassImp(StMtdGeoNode)
113 //#endif
114 
115 // ___________________________________________________________________________
116 StMtdGeoNode::StMtdGeoNode(TGeoVolume *vol, TGeoHMatrix *mat, StThreeVectorD point, Int_t nExtraCells)
117 : fPoint(point)
118 {
119 
120  fMatrix = new TGeoHMatrix();
121  fMatrix->CopyFrom(mat);
122  fVolume = vol->CloneVolume();
123  fTransFlag = kFALSE;
124  fNExtraCells = nExtraCells;
125  UpdateMatrix();
126  fNormal = YZPlaneNormal();
127 }
128 
129 StMtdGeoNode::~StMtdGeoNode()
130 {
131 
132  delete fVolume;
133  delete fMatrix;
134  delete fTransMRS;
135  delete fRotMRS;
136 }
137 
138 
139 //_____________________________________________________________________________
140 void StMtdGeoNode::UpdateMatrix()
141 {
142  //
143  //Update the Translation and RotateMatrix between local and master
144  // and store them as members: mTransMRS, mRotMRS
145  //while TNode stores them locally to share among all TNode objects,
146  // thus they may be changed afterward !
147  //
148 
149  fTransFlag = kFALSE;
150  if(fMatrix){
151  fTransMRS = fMatrix->GetTranslation();
152  fRotMRS = fMatrix->GetRotationMatrix();
153  fTransFlag = kTRUE;
154  }else{
155  LOG_WARN<<"No TGeoHMatrix!"<<endm;
156  }
157 
158 }
159 //_____________________________________________________________________________
160 void StMtdGeoNode::LocalToMaster(const Double_t* local, Double_t* master)
161 {
162  //
163  //Transform local coordinate into global coordinate
164  //
165  // UpdateMatrix();
166 
167  if (!fMatrix) {
168  if(fTransFlag){
169  //mTransFlag==kTRUE, i.e. StMtdGeomNode::UpdateMatrix() invoked already
170  Double_t x, y, z;
171  x = local[0];
172  y = local[1];
173  z = local[2];
174 
175  master[0] = fTransMRS[0] + fRotMRS[0]*x + fRotMRS[1]*y + fRotMRS[2]*z;
176  master[1] = fTransMRS[1] + fRotMRS[3]*x + fRotMRS[4]*y + fRotMRS[5]*z;
177  master[2] = fTransMRS[2] + fRotMRS[6]*x + fRotMRS[7]*y + fRotMRS[8]*z;
178  }else{
179  LOG_WARN << " No TGeoHMatrix::LocalToMaster is wrong, so do nothing" << endm;
180  }
181  } else {
182  fMatrix->LocalToMaster(local, master);
183  }
184 
185  //definition is not same as StBTofGeometry, belows are wrong
186  //master[0] = fTransMRS[0] + fRotMRS[0]*x + fRotMRS[3]*y + fRotMRS[6]*z;
187  //master[1] = fTransMRS[1] + fRotMRS[1]*x + fRotMRS[4]*y + fRotMRS[7]*z;
188  //master[2] = fTransMRS[2] + fRotMRS[2]*x + fRotMRS[5]*y + fRotMRS[8]*z;
189  //LOG_INFO<<"master ("<<master[0]<<","<<master[1]<<","<<master[2]<<")"<<endm;
190  //LOG_INFO<<"test ("<<test[0]<<","<<test[1]<<","<<test[2]<<")"<<endm;
191 
192  return;
193 
194 }
195 
196 //_____________________________________________________________________________
197 void StMtdGeoNode::MasterToLocal(const Double_t* master, Double_t* local)
198 {
199  //
200  //Transform global coordinate into local coordinate
201  //
202  if(!fMatrix){
203  if (fTransFlag) {
204  //fTransFlag==kTRUE, i.e. StMtdGeoNode::UpdateMatrix() invoked already
205  Double_t x, y, z;
206  x = master[0] - fTransMRS[0];
207  y = master[1] - fTransMRS[1];
208  z = master[2] - fTransMRS[2];
209 
210  local[0] = fRotMRS[0]*x + fRotMRS[3]*y + fRotMRS[6]*z;
211  local[1] = fRotMRS[1]*x + fRotMRS[4]*y + fRotMRS[7]*z;
212  local[2] = fRotMRS[2]*x + fRotMRS[5]*y + fRotMRS[8]*z;
213 
214  }else{
215  LOG_WARN<< " No TGeoHMatrix::MasterToLocal is wrong, so do nothing" << endm;
216  }
217  }else{
218  fMatrix->MasterToLocal(master,local);
219  }
220 
221  Double_t x, y, z;
222  x = master[0] - fTransMRS[0];
223  y = master[1] - fTransMRS[1];
224  z = master[2] - fTransMRS[2];
225 
226  Double_t test[3];
227  test[0] = fRotMRS[0]*x + fRotMRS[3]*y + fRotMRS[6]*z;
228  test[1] = fRotMRS[1]*x + fRotMRS[4]*y + fRotMRS[7]*z;
229  test[2] = fRotMRS[2]*x + fRotMRS[5]*y + fRotMRS[8]*z;
230 
231  return;
232 
233 }
234 
235 //_____________________________________________________________________________
236 StThreeVectorD StMtdGeoNode::YZPlaneNormal()
237 {
238  //
239  //Calculate the vector unit of normal to YZ-plane
240  // i.e. the global representation of local unit vector (1,0,0)
241  //
242 
243  Double_t ux[3], nx[3];
244 
245  ux[0] = 1; ux[1] = 0; ux[2] = 0;
246  LocalToMaster(ux,nx);
247 
248  nx[0] -= fTransMRS[0];
249  nx[1] -= fTransMRS[1];
250  nx[2] -= fTransMRS[2];
251 
252  return StThreeVectorD(nx[0],nx[1],nx[2]);
253 }
254 //_____________________________________________________________________________
255 Bool_t StMtdGeoNode::HelixCross(const StPhysicalHelixD helix, const Double_t pathToMagOutR, const Double_t tofToMagOutR, Double_t &pathL, Double_t &tof, StThreeVectorD &cross){
256 
257  Float_t MaxPathLength = 1000.;
258 
259  Bool_t ret = kFALSE;
260  pathL = -9999.;
261  tof = -9999.;
262 
263  StThreeVectorD oh = helix.origin();
264  if(TMath::IsNaN(oh.x())||TMath::IsNaN(oh.y())||TMath::IsNaN(oh.z())||TMath::Abs(oh.perp())>450.||TMath::Abs(oh.z())>300.) return ret;
265  pathL = helix.pathLength(fPoint, fNormal);
266  double betaGam = helix.momentum(gMtdGeometry->GetFieldZ(oh)).mag()/muonMass;
267  double velocity = sqrt(betaGam*betaGam/(1+betaGam*betaGam))*c_light*1e-9;
268  tof = pathL/velocity;
269  LOG_DEBUG<<"StMtdGeoNode::HelixCross() pathL = "<<pathL<<" point ("<<fPoint.x()<<","<<fPoint.y()<<","<<fPoint.z()<<") normal("<<fNormal.x()<<","<<fNormal.y()<<","<<fNormal.z()<<") radius = "<<fPoint.perp()<<endm;
270  if(pathL>0 && pathL<MaxPathLength){
271  cross = helix.at(pathL);
272  ret = IsGlobalPointIn(cross);
273  LOG_DEBUG<<"track cross point "<<cross.x()<<","<<cross.y()<<","<<cross.z()<<endm;
274  pathL += pathToMagOutR;
275  tof += tofToMagOutR;
276  }
277  return ret;
278 }
279 
280 //_____________________________________________________________________________
281 Bool_t StMtdGeoNode::IsGlobalPointIn(StThreeVectorD &global){
282  Double_t xl[3], xg[3];
283  xg[0] = global.x();
284  xg[1] = global.y();
285  xg[2] = global.z();
286  MasterToLocal(xg, xl);
287  Bool_t ret = IsLocalPointIn(xl[0], xl[1], xl[2]);
288  return ret;
289 }
290 
291 //_____________________________________________________________________________
292 Bool_t StMtdGeoNode::IsLocalPointIn(const Double_t x, const Double_t y, const Double_t z){
293  TGeoBBox *brik = (TGeoBBox*)fVolume->GetShape();
294  Double_t dx = brik->GetDX();
295  Float_t nExtraCells = fNExtraCells>1.66?fNExtraCells-1.66:0;
296  Double_t dy = brik->GetDY()+nExtraCells*(gMtdCellWidth+gMtdCellGap);
297  Double_t dz = brik->GetDZ();
298  Bool_t ret = -dx<x && x<dx && -dy<y && y<dy && -dz<z && z<dz;
299 
300  return ret;
301 }
302 
303 //_____________________________________________________________________________
304 StThreeVectorD StMtdGeoNode::GetNodePoint(){
305  return fPoint;
306 }
307 
308 //----------------------------------------------------//
309 // //
310 // StMtdGeoBackleg //
311 // //
312 //----------------------------------------------------//
313 
314 //#ifdef __ROOT__
315 //ClassImp(StMtdGeoBackleg)
316 //#endif
317 
318 // ___________________________________________________________________________
319  StMtdGeoBackleg::StMtdGeoBackleg (Int_t iMTTG, Int_t iBL, TGeoVolume *vol, TGeoHMatrix *mat, StThreeVectorD point, Int_t nExtraCells)
320 : StMtdGeoNode(vol, mat, point, nExtraCells), mMTTGIndex(iMTTG), mBacklegIndex(iBL)
321 {
322 
323 }
324 
325 StMtdGeoBackleg::~StMtdGeoBackleg()
326 {
327 
328 }
329 
330 //----------------------------------------------------//
331 // //
332 // StMtdGeoModule //
333 // //
334 //----------------------------------------------------//
335 
336 //#ifdef __ROOT__
337 //ClassImp(StMtdGeoModule)
338 //#endif
339 
340 // ___________________________________________________________________________
341  StMtdGeoModule::StMtdGeoModule (Int_t iMTRA, Int_t iMod, TGeoVolume *vol, TGeoHMatrix *mat, StThreeVectorD point, Int_t nExtraCells)
342 : StMtdGeoNode(vol, mat, point, nExtraCells), mMTRAIndex(iMTRA), mModuleIndex(iMod)
343 {
344 
345 }
346 
347 StMtdGeoModule::~StMtdGeoModule()
348 {
349 
350 }
351 
352 //_____________________________________________________________________________
353 Int_t StMtdGeoModule::FindCellId(const Double_t *local){
354  Int_t cellId = -99;
355  if ( IsLocalPointIn(local[0],local[1],local[2]) ) {
356  cellId = TMath::FloorNint((local[1]+(gMtdCellWidth+gMtdCellGap)*gMtdNCells/2.)/(gMtdCellWidth+gMtdCellGap));
357  }
358  return cellId;
359 }
360 
361 //_____________________________________________________________________________
362 Float_t StMtdGeoModule::GetCellPhiCenter(Int_t iCell){
363  Float_t phi = fPoint.phi();
364  Float_t r = fPoint.perp();
365  Float_t stripPhiCen = 0.;
366  if(mModuleIndex>0&&mModuleIndex<4){
367  stripPhiCen = phi-(gMtdNCells/2.-0.5-iCell)*(gMtdCellWidth+gMtdCellGap)/r; // approximation
368  }else{
369  stripPhiCen = phi+(gMtdNCells/2.-0.5-iCell)*(gMtdCellWidth+gMtdCellGap)/r;
370  }
371  if(stripPhiCen>2.*TMath::Pi()) stripPhiCen -= 2.*TMath::Pi();
372  if(stripPhiCen<0.) stripPhiCen += 2.*TMath::Pi();
373  return stripPhiCen;
374 }
375 
376 //_____________________________________________________________________________
377 Float_t StMtdGeoModule::GetCellZCenter(Int_t iCell){
378  return fPoint.z();
379 }
380 
381 //_____________________________________________________________________________
382 Float_t StMtdGeoModule::GetCellLocalYCenter(Int_t iCell, Int_t iBL, Int_t idTruth){
383  // Suitable before Run17
384 
385  Float_t cell_width = gMtdCellWidth+gMtdCellGap;
386  Float_t y_center = (mModuleIndex<4? 1 : -1) * (iCell-gMtdNCells/2+0.5) * cell_width;
387 
388  if(idTruth==0)
389  {
390  // only needed for real hits
391  // MC hits are produced with ideal geometry
392  if(iBL==8) y_center -= 3 * cell_width;
393  if(iBL==24) y_center += 2 * cell_width;
394  }
395 
396  return y_center;
397 }
398 
399 //_____________________________________________________________________________
400 Float_t StMtdGeoModule::GetCellLocalYCenter(Int_t iCell){
401  // Suitable starting from Run17
402 
403  Float_t cell_width = gMtdCellWidth+gMtdCellGap;
404  Float_t y_center = (mModuleIndex<4? 1 : -1) * (iCell-gMtdNCells/2+0.5) * cell_width;
405  return y_center;
406 }
407 
408 //----------------------------------------------------//
409 // //
410 // StMtdGeometry //
411 // //
412 //----------------------------------------------------//
413 
414 StMtdGeometry* gMtdGeometry = 0;
415 
416 #ifdef __ROOT__
417 ClassImp(StMtdGeometry)
418 #endif
419 
420 // ___________________________________________________________________________
421 StMtdGeometry::StMtdGeometry(const char* name, const char* title, TGeoManager *manager)
422 : TNamed(name,title)
423 {
424  //
425  //We only need one instance of StMtdGeometry
426  //
427 
428  mDebug = kFALSE;
429  mCosmicFlag = kFALSE;
430  mELossFlag = kTRUE;
431  mNExtraCells = 0;
432  mNValidBLs = 0;
433  mStarBField = 0;
434  mBFactor = -1.;
435  mLockBField = 0;
436  assert(manager);
437  mGeoManager = manager;
438 
439  fMagEloss = new TF1("f2","[0]*exp(-pow([1]/x,[2]))",0.,100);
440  fMagEloss->SetParameters(1.38147e+00,6.08655e-02,5.03337e-01);
441 
442  for(int i=0;i<gMtdNBacklegs;i++) {
443  mMtdGeoBackleg[i] = 0;
444  for(int j=0;j<gMtdNModules;j++) {
445  mMtdGeoModule[i][j] = 0;
446  }
447  }
448 
449  if (gMtdGeometry) {
450  LOG_INFO << "Warning !! There is already StMtdGeometry at pointer="
451  << (void*)gMtdGeometry << ", so it is deleted"
452  << endm;
453  delete gMtdGeometry;
454  }
455  gMtdGeometry = this;
456 }
457 
458 StMtdGeometry::~StMtdGeometry()
459 {
460  if(IsDebugOn()){
461  LOG_INFO << "StMtdGeometry at pointer =" << (void*)gMtdGeometry
462  << " will be deleted" << endm;
463  }
464  gMtdGeometry = 0;
465 }
466 
467 void StMtdGeometry::Init(StMaker *maker){
468 
469  if(maker->Debug()) DebugOn();
470 
471  Int_t mYear = maker->GetDateTime().GetYear();
472  if(IsDebugOn()){
473  LOG_INFO<<"Input data from year "<<mYear<<endm;
474  }
475 
476  // Get magnetic field map
477  if(!StarMagField::Instance() && mLockBField)
478  {
479  LOG_INFO<<" Initializing locked mag.field for simulation! "<<endm;
480  new StarMagField ( StarMagField::kMapped, mBFactor);
481  }
482  assert(StarMagField::Instance());
483  mStarBField = StarMagField::Instance();
484  LOG_INFO<<" Initializing mag.field from StarMagField! fScale = "<<mStarBField->GetFactor()<<endm;
485 
486  // Load geant2backlegID Map
487  // Extract MTD maps from database
488  if(mYear<=2012)
489  {
490  LOG_WARN<<"No geant2backlegID table for 2012 and before"<<endm;
491  }
492  else
493  {
494  LOG_INFO << "Retrieving geant2backlegID table from database ..." << endm;
495 
496  TDataSet *dataset = maker->GetDataBase("Geometry/mtd/mtdGeant2BacklegIDMap");
497  St_mtdGeant2BacklegIDMap *mtdGeant2BacklegIDMap = static_cast<St_mtdGeant2BacklegIDMap*>(dataset->Find("mtdGeant2BacklegIDMap"));
498  if ( !mtdGeant2BacklegIDMap )
499  {
500  LOG_ERROR << "No mtdGeant2BacklegIDMap found in database" << endm;
501  return;
502  }
503 
504  mtdGeant2BacklegIDMap_st *mGeant2BLTable = static_cast<mtdGeant2BacklegIDMap_st*>(mtdGeant2BacklegIDMap->GetTable());
505  if(mGeant2BLTable)
506  {
507  for ( Int_t i = 0; i < 30; i++ )
508  {
509  mMTTG2BL[i] = (Int_t)mGeant2BLTable->geant2backlegID[i];
510  if(mMTTG2BL[i]>30 || mMTTG2BL[i]<0) LOG_ERROR<<" Wrong map! check database! "<<endm;
511  }
512  }
513  else
514  {
515  LOG_ERROR << "No geant2backlegIDMap table found in database" << endm;
516  }
517  }
518 
519  for(int i=0;i<30;++i)
520  {
521  for(int j=0;j<5;++j)
522  {
523  if(i<11||i>19) mMTRA2Mod[i][j] = j+1;
524  else
525  {
526  if(j==3||j==4) mMTRA2Mod[i][j] = 0;
527  else mMTRA2Mod[i][j] = j+2;
528  }
529  }
530  }
531 
532  Int_t mGeoYear = 0;
533  if(mGeoManager->CheckPath("/HALL_1/CAVE_1/MUTD_1/MTMT_1")){
534  LOG_INFO<<"found y2012 geometry"<<endm;
535  mGeoYear=2012;
536  }
537  else if(mGeoManager->CheckPath("/HALL_1/CAVE_1/MagRefSys_1/MUTD_1")){
538  LOG_INFO<<"found y2015 geometry"<<endm;
539  mGeoYear=2015;
540  }
541 
542  // intialize backleg/module geometry
543  TGeoVolume *mMtdGeom = mGeoManager->FindVolumeFast("MUTD");
544  const char *elementName = mMtdGeom->GetName();
545  if(elementName){
546  if(IsDebugOn()) LOG_INFO <<" found detector:"<<elementName<<endm;
547  TGeoIterator next(mMtdGeom);
548  if(mGeoYear==2015) next.SetTopName("/HALL_1/CAVE_1/MagRefSys_1/MUTD_1");
549  else next.SetTopName("/HALL_1/CAVE_1/MUTD_1");
550  TGeoNode *node = 0;
551  TGeoVolume *blVol = 0;
552  TGeoVolume *modVol = 0;
553  Int_t ibackleg = 0;
554  Int_t imodule = 0;
555  Double_t minR = 999.;
556  Double_t maxR = 0.;
557  while ( (node=(TGeoNode*)next()) )
558  {
559  TString name = node->GetName();
560  TString path;
561  next.GetPath(path);
562  if(!mGeoManager->CheckPath(path.Data())){
563  LOG_WARN<<"Path "<<path.Data()<<" is not found"<<endm;
564  continue;
565  }
566  mGeoManager->cd(path.Data());
567  TGeoVolume *detVol = mGeoManager->GetCurrentVolume();
568  Bool_t found = ( IsMTTG(detVol) || IsMTRA(detVol) );
569  if (found) {
570  detVol->SetVisibility(kTRUE);
571  //if (detVol->GetLineColor()==1 || detVol->GetLineColor()==7)
572  // detVol->SetLineColor(14);
573  } else {
574  detVol->SetVisibility(kFALSE);
575  continue;
576  }
577  if(IsDebugOn()) LOG_INFO<<"currentpath = "<<mGeoManager->GetPath()<<" node name="<<mGeoManager->GetCurrentNode()->GetName()<<endm;
578 
579  //fill GeoBLs and GeoModules
580  if(IsMTTG(detVol)){
581  blVol = (TGeoVolume *)detVol;
582  Int_t iMTTG = node->GetNumber();
583  if(IsDebugOn()) LOG_INFO<<"Node name = "<<node->GetName()<<" iMTTG="<<iMTTG<<endm;
584  if(iMTTG>0){
585  Int_t iBL = 0;
586  if(mGeoYear==2012)
587  {
588  if(!strcmp(blVol->GetName(), backlegPref[0]))
589  {
590  iBL = 26;
591  }
592  else if(!strcmp(blVol->GetName(), backlegPref[1]))
593  {
594  if(iMTTG==1) iBL = 27;
595  else if(iMTTG==2) iBL = 28;
596  else
597  {
598  LOG_ERROR<<"Wrong BL id, this is not Y2012 geometry!"<<endm;
599  iBL = 0;
600  continue;
601  }
602  }
603  else
604  {
605  LOG_ERROR<<"Wrong BL id, this is not Y2012 geometry!"<<endm;
606  iBL = 0;
607  continue;
608  }
609  }
610  else
611  {
612  iBL = mMTTG2BL[iMTTG-1];
613  }
614  Double_t op[3];
615  Double_t local[3] = {0,0,0};
616  mGeoManager->LocalToMaster(local,op);
617  ++mNValidBLs;
618  TGeoHMatrix *mat = mGeoManager->GetCurrentMatrix();
619  StThreeVectorD point(op[0],op[1],op[2]);
620  if(mGeoYear==2012)
621  {
622  Double_t *trans = mat->GetTranslation();
623  Double_t *rot = mat->GetRotationMatrix();
624  trans[0] = -1*trans[0];
625  rot[0] = -1*rot[0];
626  rot[4] = -1*rot[4];
627  mat->SetTranslation(trans);
628  mat->SetRotation(rot);
629  point.setX(-1*op[0]);
630  }
631  mMtdGeoBackleg[iBL-1] = new StMtdGeoBackleg(iMTTG, iBL, blVol, mat, point, mNExtraCells);
632  if(IsDebugOn()) LOG_INFO<<"iMTTG="<<iMTTG<<" iBL="<<iBL<<endm;
633  imodule = 0; // clear for this backleg
634  ibackleg = iBL;
635  }
636  }
637 
638  if(IsMTRA(detVol)){
639  modVol = (TGeoVolume *)detVol;
640  Int_t iMTRA = node->GetNumber();
641  if(iMTRA>0){
642  Int_t iMod = mMTRA2Mod[ibackleg-1][iMTRA-1];
643  if(mGeoYear==2012&&ibackleg==26){
644  iMod = iMTRA+1;
645  }
646  Double_t op[3];
647  Double_t local[3] = {0,0,0};
648  mGeoManager->LocalToMaster(local,op);
649  ++imodule;
650  double R = sqrt(op[0]*op[0]+op[1]*op[1]);
651  if(R<minR) minR = R;
652  if(R>maxR) maxR = R;
653 
654  TGeoHMatrix *mat = mGeoManager->GetCurrentMatrix();
655  StThreeVectorD point(op[0],op[1],op[2]);
656  if(mGeoYear==2012)
657  {
658  Double_t *trans = mat->GetTranslation();
659  Double_t *rot = mat->GetRotationMatrix();
660  trans[0] = -1*trans[0];
661  rot[0] = -1*rot[0];
662  rot[4] = -1*rot[4];
663  mat->SetTranslation(trans);
664  mat->SetRotation(rot);
665  point.setX(-1*op[0]);
666  }
667 
668  mMtdGeoModule[ibackleg-1][iMod-1] = new StMtdGeoModule(iMTRA, iMod, modVol, mat, point, mNExtraCells);
669  if(IsDebugOn()) LOG_INFO<<"iMTRA="<<iMTRA<<" iMod="<<iMod<<endm;
670  }
671  }
672  }
673  }
674 
675  if(IsDebugOn()){
676  for(int i=0;i<gMtdNBacklegs;i++){
677  for(int j=0;j<gMtdNModules;j++){
678  if(mMtdGeoModule[i][j])
679  LOG_INFO<<"valid (backleg,module) = "<<i+1<<","<<j+1<<endm;
680  }
681  }
682  }
683 }
684 
685 Bool_t StMtdGeometry::ProjToMagOutR(const StPhysicalHelixD helix, const StThreeVectorD vertex, StPhysicalHelixD &outHelix, Double_t &pathL, Double_t &tof, StThreeVectorD &pos){
686 
687  // 1. project to vertex
688  pathL = 0.;
689  tof = -9999.;
690  pos.set(0,0,0);
691  StThreeVectorD dcaPos(0,0,0);
692  Double_t pathL2Vtx = 0.;
693  ProjToVertex(helix,vertex,pathL2Vtx,tof,dcaPos);
694 
695  if(IsDebugOn()){
696  LOG_INFO<<"------------ start of projection to magOutR ------------"<<endm;
697  LOG_INFO<<" to vertex("<< vertex.x()<<","<<vertex.y()<<","<<vertex.z()<<") dca = "<<dcaPos.mag()<<endm;
698  }
699  //if(!mCosmicFlag && dcaPos.mag()>10) return kFALSE;
700 
701  //2. project to Emc
702  pairD sInnerEmc = helix.pathLength(mEmcInR);
703  if(sInnerEmc.first<=0 && sInnerEmc.second<=0) return kFALSE;
704  double rInnerEmc = (sInnerEmc.first < 0 || sInnerEmc.second < 0)
705  ? max(sInnerEmc.first, sInnerEmc.second) : min(sInnerEmc.first, sInnerEmc.second);
706  StThreeVectorD innerEmcPos = helix.at(rInnerEmc);
707 
708  if(TMath::IsNaN(innerEmcPos.x())||TMath::IsNaN(innerEmcPos.y())||TMath::IsNaN(innerEmcPos.z())||TMath::Abs(innerEmcPos.perp())>450.||TMath::Abs(innerEmcPos.z())>300.) return kFALSE;
709  Double_t bField = GetFieldZ(innerEmcPos);
710  Int_t charge = helix.charge(bField);
711  if(IsDebugOn()) LOG_INFO<<" charge = "<<charge<<" bField = "<<bField<<endm;
712  StThreeVectorD innerEmcMom = helix.momentumAt(rInnerEmc,bField*tesla);
713  StPhysicalHelixD helixInEmc(innerEmcMom,innerEmcPos,bField*tesla,charge);
714 
715  double betaGam = innerEmcMom.mag()/muonMass;
716  double vInner = sqrt(betaGam*betaGam/(1+betaGam*betaGam))*c_light*1e-9;
717  double tof2InnerEmc = -9999.;
718  tof2InnerEmc = (pathL2Vtx+rInnerEmc)/vInner;
719 
720  if(IsDebugOn()){
721  LOG_INFO<<" to Emcinner: pos x,y,z:"<<innerEmcPos.x()<<","<<innerEmcPos.y()<<","<<innerEmcPos.z()<<endm;
722  LOG_INFO<<" to Emcinner: mom p,pt,eta,phi:"<<innerEmcMom.mag()<<","<<innerEmcMom.perp()<<","<<innerEmcMom.pseudoRapidity()<<","<<innerEmcMom.phi()<<endm;
723  LOG_INFO<<" to Emcinner: tof:"<<tof2InnerEmc<<endm;
724  }
725 
726  //3. outer Emc
727  const int nEmcStep = 4;
728  double rEmcStep = (mEmcOutR-mEmcInR)/nEmcStep; //cm
729  double pathLEmcLayer[nEmcStep];
730  double tofEmcLayer[nEmcStep];
731  double vEmc = -9999.;
732  double betaGamEmc= -9999.;
733  StThreeVectorD EmcLayerPos = innerEmcPos;
734  StThreeVectorD EmcLayerMom = innerEmcMom;
735  Double_t elossEmc = 0.;
736  if(mELossFlag){
737  if(mCosmicFlag&&EmcLayerPos.phi()>0&&EmcLayerPos.phi()<TMath::Pi()) elossEmc = -1.*gMtdMuonELossInEmc/nEmcStep;
738  else elossEmc = gMtdMuonELossInEmc/nEmcStep;
739  }
740  for( int i=0; i<nEmcStep; i++){
741  double EmcLayerRadius = mEmcInR+rEmcStep*(i+1);
742 
743  if(TMath::IsNaN(EmcLayerPos.x())||TMath::IsNaN(EmcLayerPos.y())||TMath::IsNaN(EmcLayerPos.z())||TMath::Abs(EmcLayerPos.perp())>450.||TMath::Abs(EmcLayerPos.z())>300.) return kFALSE;
744  bField = GetFieldZ(EmcLayerPos);
745  StPhysicalHelixD helixInEmc(EmcLayerMom,EmcLayerPos,bField*tesla,charge);
746 
747  pairD sEmcLayer = helixInEmc.pathLength(EmcLayerRadius);
748  if(sEmcLayer.first<=0 && sEmcLayer.second<=0) return kFALSE;
749 
750  double rEmcLayer = (sEmcLayer.first < 0 || sEmcLayer.second < 0)
751  ? max(sEmcLayer.first, sEmcLayer.second) : min(sEmcLayer.first, sEmcLayer.second);
752 
753  betaGamEmc = EmcLayerMom.mag()/muonMass;
754  vEmc = sqrt(betaGamEmc*betaGamEmc/(1.+betaGamEmc*betaGamEmc))*c_light*1e-9;
755  pathLEmcLayer[i] = rEmcLayer;
756  tofEmcLayer[i] = rEmcLayer/vEmc;
757 
758  EmcLayerPos = helixInEmc.at(rEmcLayer);
759  EmcLayerMom = helixInEmc.momentumAt(rEmcLayer,bField*tesla);
760  //EmcLayerMom = (sqrt(pow((sqrt(pow((abs(EmcLayerMom)),2)+pow(muonMass,2))-rEmcLayer*elossEmc),2)-pow(muonMass,2)))/(abs(EmcLayerMom))*EmcLayerMom;
761  EmcLayerMom = (EmcLayerMom.mag()-elossEmc)/EmcLayerMom.mag()*EmcLayerMom;
762  }
763 
764  if(IsDebugOn()){
765  double tof2OuterEmc = 0.;
766  for(int i=0;i<nEmcStep;i++) tof2OuterEmc += tofEmcLayer[i];
767  LOG_INFO<<" to Emcouter: pos x,y,z:"<<EmcLayerPos.x()<<","<<EmcLayerPos.y()<<","<<EmcLayerPos.z()<<endm;
768  LOG_INFO<<" to Emcouter: mom p,pt,eta,phi:"<<EmcLayerMom.mag()<<","<<EmcLayerMom.perp()<<","<<EmcLayerMom.pseudoRapidity()<<","<<EmcLayerMom.phi()<<endm;
769  LOG_INFO<<" to Emcouter: tof:"<<tof2OuterEmc<<endm;
770  }
771 
772  //4: coil
773  if(TMath::IsNaN(EmcLayerPos.x())||TMath::IsNaN(EmcLayerPos.y())||TMath::IsNaN(EmcLayerPos.z())||TMath::Abs(EmcLayerPos.perp())>450.||TMath::Abs(EmcLayerPos.z())>300.) return kFALSE;
774  bField = GetFieldZ(EmcLayerPos);
775  const int nCoilStep = 5;
776  double rCoilStep = (mMagInR-mEmcOutR)/nCoilStep; //cm
777  double pathLCoilLayer[nCoilStep];
778  double tofCoilLayer[nCoilStep];
779  double vCoil=-9999.;
780  double betaGamCoil=-9999.;
781  StThreeVectorD CoilLayerPos = EmcLayerPos;
782  StThreeVectorD CoilLayerMom = EmcLayerMom;
783  Double_t elossCoil = 0.;
784  if(mELossFlag){
785  if(mCosmicFlag&&CoilLayerPos.phi()>0&&CoilLayerPos.phi()<TMath::Pi()) elossCoil = -1.*gMtdMuonELossInCoil/nCoilStep;
786  else elossCoil = gMtdMuonELossInCoil/nCoilStep;
787  }
788  for( int i=0; i<nCoilStep; i++){
789  double CoilLayerRadius = mEmcOutR+rCoilStep*(i+1);
790 
791  if(TMath::IsNaN(CoilLayerPos.x())||TMath::IsNaN(CoilLayerPos.y())||TMath::IsNaN(CoilLayerPos.z())||TMath::Abs(CoilLayerPos.perp())>450.||TMath::Abs(CoilLayerPos.z())>300.) return kFALSE;
792  bField = GetFieldZ(CoilLayerPos);
793  StPhysicalHelixD helixInCoil(CoilLayerMom,CoilLayerPos,bField*tesla,charge);
794 
795  pairD sCoilLayer = helixInCoil.pathLength(CoilLayerRadius);
796  if(sCoilLayer.first<=0 && sCoilLayer.second<=0) return kFALSE;
797 
798  double rCoilLayer = (sCoilLayer.first < 0 || sCoilLayer.second < 0)
799  ? max(sCoilLayer.first, sCoilLayer.second) : min(sCoilLayer.first, sCoilLayer.second);
800 
801  betaGamCoil = CoilLayerMom.mag()/muonMass;
802  vCoil = sqrt(betaGamCoil*betaGamCoil/(1.+betaGamCoil*betaGamCoil))*c_light*1e-9;
803  pathLCoilLayer[i] = rCoilLayer;
804  tofCoilLayer[i] = rCoilLayer/vCoil;
805 
806  CoilLayerPos = helixInCoil.at(rCoilLayer);
807  CoilLayerMom = helixInCoil.momentumAt(rCoilLayer,bField*tesla);
808  if(IsDebugOn()){
809  LOG_INFO<<" to Coil Layer "<<i<<" bField = "<<bField<<endm;
810  LOG_INFO<<" to Coil Layer"<< i <<" : pos x,y,z:"<<CoilLayerPos.x()<<","<<CoilLayerPos.y()<<","<<CoilLayerPos.z()<<endm;
811  LOG_INFO<<" to Coil Layer"<< i <<" : mom p,pt,eta,phi:"<<CoilLayerMom.mag()<<","<<CoilLayerMom.perp()<<","<<CoilLayerMom.pseudoRapidity()<<","<<CoilLayerMom.phi()<<endm;
812  }
813  //CoilLayerMom=(sqrt(pow((sqrt(pow((abs(CoilLayerMom)),2)+pow(muonMass,2))-rCoilLayer*elossCoil),2)-pow(muonMass,2)))/(abs(CoilLayerMom))*CoilLayerMom;
814  CoilLayerMom = (CoilLayerMom.mag()-elossCoil)/CoilLayerMom.mag()*CoilLayerMom;
815  }
816 
817  if(IsDebugOn()){
818  double tof2OuterCoil = 0.;
819  for(int i=0;i<nCoilStep;i++) tof2OuterCoil += tofCoilLayer[i];
820  LOG_INFO<<" to Coilouter: pos x,y,z:"<<CoilLayerPos.x()<<","<<CoilLayerPos.y()<<","<<CoilLayerPos.z()<<endm;
821  LOG_INFO<<" to Coilouter: mom p,pt,eta,phi:"<<CoilLayerMom.mag()<<","<<CoilLayerMom.perp()<<","<<CoilLayerMom.pseudoRapidity()<<","<<CoilLayerMom.phi()<<endm;
822  LOG_INFO<<" to Coilouter: tof:"<<tof2OuterCoil<<endm;
823  }
824 
825  //5. Mag
826  const Int_t nMagStep = 10;
827  double pathLMagLayer[nMagStep];
828  double tofMagLayer[nMagStep];
829  double vMag;
830  double betaGamMag;
831  double rStep = (mMagOutR-mMagInR)/nMagStep; //cm
832  StThreeVector<double> MagLayerPos = CoilLayerPos;
833  StThreeVector<double> MagLayerMom = CoilLayerMom;
834  Double_t elossMag = 0.;
835  double mMagELoss = 0.;
836  if(mELossFlag){
837  mMagELoss = fMagEloss->Eval(innerEmcMom.mag())-gMtdMuonELossInEmc-gMtdMuonELossInCoil;
838  if(mCosmicFlag&&MagLayerPos.phi()>0&&MagLayerPos.phi()<TMath::Pi()) elossMag = -1.*mMagELoss/nMagStep;
839  else elossMag = mMagELoss/nMagStep;
840  }
841  if(IsDebugOn()){
842  LOG_INFO<<" mMagELoss = "<<mMagELoss<<" p = "<<innerEmcMom.mag()<<endm;
843  }
844  for( int i=0; i<nMagStep; i++){
845 
846  double MagLayerRadius = mMagInR+rStep*(i+1);
847  if(TMath::IsNaN(MagLayerPos.x())||TMath::IsNaN(MagLayerPos.y())||TMath::IsNaN(MagLayerPos.z())||TMath::Abs(MagLayerPos.perp())>450.||TMath::Abs(MagLayerPos.z())>300.) return kFALSE;
848  bField = GetFieldZ(MagLayerPos);
849  StPhysicalHelixD helixInMag(MagLayerMom,MagLayerPos,bField*tesla,charge);
850 
851  if(IsDebugOn()){
852  LOG_INFO<<" to Magnet layer "<<i<<" bField = "<<bField<<endm;
853  LOG_INFO<<" to Magnet Layer "<< i <<" : pos x,y,z:"<<MagLayerPos.x()<<","<<MagLayerPos.y()<<","<<MagLayerPos.z()<<endm;
854  LOG_INFO<<" to Magnet Layer "<< i <<" : mom p,pt,eta,phi:"<<MagLayerMom.mag()<<","<<MagLayerMom.perp()<<","<<MagLayerMom.pseudoRapidity()<<","<<MagLayerMom.phi()<<endm;
855  }
856  pairD sMagLayer = helixInMag.pathLength(MagLayerRadius);
857  if(sMagLayer.first<=0 && sMagLayer.second<=0) return kFALSE;
858 
859  double rMagLayer = (sMagLayer.first < 0 || sMagLayer.second < 0)
860  ? max(sMagLayer.first, sMagLayer.second) : min(sMagLayer.first, sMagLayer.second);
861 
862  betaGamMag = MagLayerMom.mag()/muonMass;
863  vMag = sqrt(betaGamMag*betaGamMag/(1+betaGamMag*betaGamMag))*c_light*1e-9;
864  pathLMagLayer[i] = rMagLayer;
865  tofMagLayer[i] = rMagLayer/vMag;
866 
867  MagLayerPos = helixInMag.at(rMagLayer);
868  MagLayerMom = helixInMag.momentumAt(rMagLayer,bField*tesla);
869 
870  double momMag = MagLayerMom.mag();
871  if(momMag<elossMag) return kFALSE;
872  MagLayerMom = (momMag-elossMag)/momMag*MagLayerMom;
873  }
874 
875  if(IsDebugOn()){
876  double tof2OuterMag = 0.;
877  for(int i=0;i<nMagStep;i++) tof2OuterMag += tofMagLayer[i];
878  LOG_INFO<<" to OuterMag: pos x,y,z:"<<MagLayerPos.x()<<","<<MagLayerPos.y()<<","<<MagLayerPos.z()<<endm;
879  LOG_INFO<<" to OuterMag: mom p,pt,eta,phi:"<<MagLayerMom.mag()<<","<<MagLayerMom.perp()<<","<<MagLayerMom.pseudoRapidity()<<","<<MagLayerMom.phi()<<endm;
880  LOG_INFO<<" to OuterMag: tof:"<<tof2OuterMag<<endm;
881  }
882 
883  double tof2MagOuter = -9999.;
884  tof2MagOuter = tof2InnerEmc;
885  for(int i=0;i<nEmcStep;i++) tof2MagOuter += tofEmcLayer[i];
886  for(int i=0;i<nCoilStep;i++) tof2MagOuter += tofCoilLayer[i];
887  for(int i=0;i<nMagStep;i++) tof2MagOuter += tofMagLayer[i];
888 
889  double pathL2MagOuter = -9999.;
890  if(IsDebugOn()) LOG_INFO<<" path to vertex("<< vertex.x()<<","<<vertex.y()<<","<<vertex.z()<<") = "<<pathL2Vtx<<endm;
891  pathL2MagOuter = pathL2Vtx+rInnerEmc;
892  if(IsDebugOn()) LOG_INFO<<" path to inner Emc = "<<pathL2MagOuter<<endm;
893  for(int i=0;i<nEmcStep;i++) pathL2MagOuter += pathLEmcLayer[i];
894  if(IsDebugOn()) LOG_INFO<<" path to outer Emc = "<<pathL2MagOuter<<endm;
895  for(int i=0;i<nCoilStep;i++) pathL2MagOuter += pathLCoilLayer[i];
896  if(IsDebugOn()) LOG_INFO<<" path to outer Coil = "<<pathL2MagOuter<<endm;
897  for(int i=0;i<nMagStep;i++) pathL2MagOuter += pathLMagLayer[i];
898  if(IsDebugOn()) LOG_INFO<<" path to outer mag = "<<pathL2MagOuter<<endm;
899 
900  pathL = pathL2MagOuter;
901  tof = tof2MagOuter;
902  pos = MagLayerPos;
903  if(TMath::IsNaN(MagLayerPos.x())||TMath::IsNaN(MagLayerPos.y())||TMath::IsNaN(MagLayerPos.z())||TMath::Abs(MagLayerPos.perp())>450.||TMath::Abs(MagLayerPos.z())>300.) return kFALSE;
904  bField = GetFieldZ(MagLayerPos);
905  //bField = 0.;
906  outHelix = StPhysicalHelixD(MagLayerMom,MagLayerPos,bField*tesla,charge);
907  if(IsDebugOn()){
908  LOG_INFO<<"bField of magnet outside is "<<bField<<" from magMap = "<<GetFieldZ(MagLayerPos)<<" mom = "<<MagLayerMom.x()<<","<<MagLayerMom.y()<<","<<MagLayerMom.z()<<" mag = "<<MagLayerMom.mag()<<endm;
909  LOG_INFO<<"------------ end of projection to magOutR ------------"<<endm;
910  }
911  return kTRUE;
912 }
913 
914 void StMtdGeometry::ProjToVertex(const StPhysicalHelixD helix, const StThreeVectorD vertex, Double_t &pathL, Double_t &tof, StThreeVectorD &dcaPos){
915  pathL = -9999.;
916  tof = -9999.;
917  pathL = TMath::Abs(helix.pathLength(vertex));
918  dcaPos = helix.at(helix.pathLength(vertex));
919 
920  StThreeVectorD oh = helix.origin();
921  if(TMath::IsNaN(oh.x())||TMath::IsNaN(oh.y())||TMath::IsNaN(oh.z())||TMath::Abs(oh.perp())>450.||TMath::Abs(oh.z())>300.) return;
922  Double_t betaGam = helix.momentum(GetFieldZ(oh)).mag()/muonMass;
923  Double_t v = sqrt(betaGam*betaGam/(1+betaGam*betaGam))*c_light*1e-9;
924  if(v!=0){
925  tof = pathL/v;
926  }else{
927  tof = -9999.;
928  }
929 }
930 
931 Bool_t StMtdGeometry::ProjToBLModVect(const StPhysicalHelixD helix, IntVec &blVect, IntVec &modVect){
932 
933  blVect.clear();
934  modVect.clear();
935 
936  Double_t R_mtd[2] = {mMtdMinR,mMtdMaxR};
937  Double_t iBL[2] = {0,0};
938  Double_t iMod[2] = {0,0};
939  for(int i=0;i<2;i++) {
940  Double_t s = helix.pathLength(R_mtd[i]).first;
941  if(s<0.) s = helix.pathLength(R_mtd[i]).second;
942  StThreeVectorD point = helix.at(s);
943  Double_t phi = point.phi();
944  Double_t z = point.z();
945 
946  iBL[i] = FindBLId(phi);
947  iMod[i] = FindModId(z);
948  }
949 
950  for (Int_t i = 0; i < 2; i++) {
951  for(Int_t j=0;j<3;j++){
952  Int_t idx = iBL[i]-1+j;
953  if(idx>gMtdNBacklegs) idx -= gMtdNBacklegs;
954  if(idx<1) idx += gMtdNBacklegs;
955  if(idx>0&&idx<=gMtdNBacklegs) blVect.push_back(idx);
956  }
957  }
958  RemoveDuplicate(blVect);
959  for (Int_t i = 0; i < 2; i++) {
960  for(Int_t j=0;j<3;j++){
961  Int_t idx = iMod[i]-1+j;
962  if(idx>0&&idx<=gMtdNModules) modVect.push_back(idx);
963  }
964  }
965  RemoveDuplicate(modVect);
966  return kTRUE;
967 }
968 
969 void StMtdGeometry::RemoveDuplicate(IntVec &vec){
970 
971  sort(vec.begin(),vec.end());
972  Int_t tmpid = 0;
973  for(IntVec::iterator tmpIter = vec.begin();tmpIter<vec.end();++tmpIter){
974  if(tmpIter==vec.begin()){
975  tmpid = *tmpIter;
976  continue;
977  }
978  if(tmpid==*tmpIter){
979  vec.erase(tmpIter);
980  --tmpIter;
981  }else{
982  tmpid=*tmpIter;
983  }
984  }
985 }
986 
987 Int_t StMtdGeometry::FindBLId(Double_t phi){
988 
989  Int_t iBL = -1;
990  if(phi<0) phi += 2.*(TMath::Pi()); // -pi,pi --> 0,2*pi
991  double backLegPhiWidth = 8.*(TMath::Pi())/180.; //rad, 8 degree per backLeg
992  double backLegPhiGap = 4.*(TMath::Pi())/180.; //rad, 4 degree
993  double dphi = backLegPhiWidth+backLegPhiGap;
994  iBL = (int)(phi/dphi);
995  iBL+= 24;
996  if(iBL>30) iBL-= 30;
997  if(iBL<1||iBL>30){
998  if(IsDebugOn()) LOG_WARN<<"Invalid BL id:"<<iBL<<" input phi = "<<phi<<endm;
999  return -1;
1000  }
1001 
1002  return iBL;
1003 }
1004 
1005 Int_t StMtdGeometry::FindModId(Double_t z){
1006 
1007  Int_t iMod = -1;
1008 
1009  iMod = (int)((z+2.5*gMtdCellLength)/gMtdCellLength+1);
1010 
1011  if(iMod<0||iMod>6){
1012  if(IsDebugOn()) LOG_WARN<<"Invalid Module id:"<<iMod<<" input z = "<<z<<endm;
1013  return -1;
1014  }
1015 
1016  return iMod;
1017 }
1018 
1019 Bool_t StMtdGeometry::HelixCrossCellIds(const StPhysicalHelixD helix, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec, DoubleVec &tofVec){
1020  const StThreeVectorD vertex(0,0,0);
1021  return HelixCrossCellIds(helix,vertex,idVec,pathVec,crossVec,tofVec);
1022 }
1023 
1024 Bool_t StMtdGeometry::HelixCrossCellIds(const StPhysicalHelixD helix, const StThreeVectorD vertex, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec, DoubleVec &tofVec){
1025 
1026  // input : helix
1027  // output: hit cell id, pathlength, cross point, tof
1028 
1029  // A. project to magnet outer radius
1030  StPhysicalHelixD outHelix;
1031  Double_t pathLToMagOutR;
1032  Double_t tofToMagOutR;
1033  StThreeVectorD posAtMagOutR;
1034  if(!ProjToMagOutR(helix,vertex,outHelix,pathLToMagOutR,tofToMagOutR,posAtMagOutR)) return kFALSE;
1035 
1036  // B. project to mtd inner and outer radius, get backleg and module candidates
1037  IntVec projBLVec;
1038  IntVec projModVec;
1039  if(!ProjToBLModVect(outHelix,projBLVec,projModVec)) return kFALSE;
1040 
1041  Int_t cellId=0;
1042  idVec.clear();
1043  pathVec.clear();
1044  crossVec.clear();
1045  tofVec.clear();
1046 
1047  // C. loop over all the candidate modules, find the hit cells and positions.
1048  for (UInt_t i = 0; i < projBLVec.size(); i++) {
1049  Int_t iBL = projBLVec[i];
1050  if(!mMtdGeoBackleg[iBL-1]) continue;
1051  for (UInt_t j = 0; j < projModVec.size(); j++) {
1052  Int_t iMod = projModVec[j];
1053  if(!mMtdGeoModule[iBL-1][iMod-1]) continue;
1054  Double_t pathToMtd = 0.;
1055  Double_t tofToMtd = 0.;
1056  StThreeVectorD crossToMtd(0,0,0);
1057  if(IsDebugOn()){
1058  LOG_INFO<<" checking track cross BL = "<<iBL<<" Module = "<<iMod<<" or not ... "<<endm;
1059  }
1060  if(mMtdGeoModule[iBL-1][iMod-1]->HelixCross(outHelix,pathLToMagOutR,tofToMagOutR,pathToMtd,tofToMtd,crossToMtd)){
1061  Double_t global[3] = {crossToMtd.x(),crossToMtd.y(),crossToMtd.z()};
1062  Double_t local[3] = {0,0,0};
1063  mMtdGeoModule[iBL-1][iMod-1]->MasterToLocal(global,local);
1064  Int_t iCel = mMtdGeoModule[iBL-1][iMod-1]->FindCellId(local);
1065  if(iCel<-50) continue;
1066  if(iMod>3) iCel = 11-iCel;
1067  cellId = CalcCellId(iBL , iMod, iCel);
1068  if(IsDebugOn()){
1069  LOG_INFO<<" track cross BL = "<<iBL<<" Module = "<<iMod<< " iCel = "<<iCel<<endm;
1070  LOG_INFO<<" cellId = "<<cellId<<" pathToMtd = "<<pathToMtd<<" tofToMtd = "<<tofToMtd<<endm;
1071  }
1072  if(IsIdValid(cellId)){
1073  idVec.push_back(cellId);
1074  pathVec.push_back(pathToMtd);
1075  crossVec.push_back(crossToMtd);
1076  tofVec.push_back(tofToMtd);
1077  }
1078  }
1079  }
1080  }
1081  return kTRUE;
1082 }
1083 
1084 Int_t StMtdGeometry::CalcCellId(Int_t iBL, Int_t iMod, Int_t iCel){
1085  Int_t cellId = -999;
1086  if(iBL<1|| iBL>gMtdNBacklegs) return cellId;
1087  if(iMod<1|| iMod>gMtdNModules) return cellId;
1088  if(iCel<-mNExtraCells || iCel>11+mNExtraCells) return cellId;
1089  cellId = iBL*1000+iMod*100+(iCel+50);
1090 
1091  return cellId;
1092 }
1093 
1094 
1095 Bool_t StMtdGeometry::IsIdValid(Int_t id){
1096 
1097  Int_t iBL = id/1000;
1098  Int_t iMod = (id%1000)/100;
1099  Int_t iCell = id%100-50;
1100 
1101  if(iBL<1 || iBL>gMtdNBacklegs) return kFALSE;
1102  if(iMod<1 || iMod>gMtdNModules) return kFALSE;
1103  if(iCell<-mNExtraCells || iCell>11+mNExtraCells) return kFALSE;
1104  return kTRUE;
1105 }
1106 
1107 void StMtdGeometry::DecodeCellId(Int_t id, Int_t &iBL, Int_t &iMod, Int_t &iCell){
1108 
1109  iBL = id/1000;
1110  iMod = (id%1000)/100;
1111  iCell = id%100-50;
1112 
1113 }
1114 
1115 StThreeVectorD StMtdGeometry::GetField(StThreeVectorD pos) const{
1116  return GetField(pos.x(),pos.y(),pos.z());
1117 }
1118 
1119 StThreeVectorD StMtdGeometry::GetField(Double_t x, Double_t y, Double_t z) const{
1120  Double_t B[3] = {0,0,0};
1121  Double_t X[3] = {x,y,z};
1122  mStarBField->BField(X,B);
1123  for(int i=0;i<3;++i) B[i] /= 10.;
1124  return StThreeVectorD(B[0],B[1],B[2]);
1125 }
1126 
1127 Double_t StMtdGeometry::GetFieldZ(StThreeVectorD pos) const{
1128  StThreeVectorD bField = GetField(pos.x(),pos.y(),pos.z());
1129  return bField.z();
1130 }
1131 
1132 Double_t StMtdGeometry::GetFieldZ(Double_t x, Double_t y, Double_t z) const{
1133  StThreeVectorD bField = GetField(x,y,z);
1134  return bField.z();
1135 }
1136 
1137 StMtdGeoModule *StMtdGeometry::GetGeoModule(Int_t iBL, Int_t iMod) const{
1138  if(iBL>0&&iBL<=gMtdNBacklegs&&iMod>0&&iMod<=gMtdNModules){
1139  return mMtdGeoModule[iBL-1][iMod-1];
1140  }else{
1141  return 0;
1142  }
1143 }
1144 
1145 Bool_t StMtdGeometry::IsMTTG(const TGeoVolume * vol) const {
1146  Bool_t found = false;
1147  for(int i=0;i<4;i++){
1148  if(!strcmp(vol->GetName(), backlegPref[i])){
1149  found = true;
1150  break;
1151  }
1152  }
1153  return found;
1154 }
1155 
1156 void StMtdGeometry::SetLockBField(Bool_t val){
1157  mLockBField=val;
1158  SetBFactor(1);
1159 }
StMtdGeometry::mCosmicFlag
Bool_t mCosmicFlag
Control message printing of this class.
Definition: StMtdGeometry.h:238
StMtdGeometry::mNExtraCells
Int_t mNExtraCells
Control mag field to FF.
Definition: StMtdGeometry.h:241
StMtdGeometry::mELossFlag
Bool_t mELossFlag
Cosmic event flag.
Definition: StMtdGeometry.h:239
tof
Definition: tof.h:15
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StMtdGeometry::IsIdValid
Bool_t IsIdValid(Int_t id)
return BL*1000+Module*100+(Cell+50). BL: 1-30, Module: 1-5, Cell: 0-11.
Definition: StMtdGeometry.cxx:1095
StHelix::origin
const StThreeVector< double > & origin() const
-sign(q*B);
Definition: StHelix.hh:224
StMtdGeometry::mNValidBLs
Int_t mNValidBLs
Control matching range in the module.
Definition: StMtdGeometry.h:242
StMtdGeometry::StMtdGeometry
StMtdGeometry(const char *name="mtdGeo", const char *title="Simplified Mtd Geometry", TGeoManager *manager=0)
EMC system outer radius.
Definition: StMtdGeometry.cxx:421
StMtdGeometry::backlegPref
static const char * backlegPref[4]
EMC system outer radius.
Definition: StMtdGeometry.h:252
StMtdGeometry::mLockBField
Bool_t mLockBField
Control energy loss flag.
Definition: StMtdGeometry.h:240
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
Generated by   doxygen 1.8.5