EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PndGeoHandling.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PndGeoHandling.cxx
1 //
2 // C++ Implementation: PndGeoHandling
3 //
4 // Description:
5 //
6 //
7 // Author: t.stockmanns <stockman@ikp455>, (C) 2007
8 //
9 // Copyright: See COPYING file that comes with this distribution
10 //
11 //
12 #include "PndGeoHandling.h"
13 #include "PndStringSeparator.h"
14 
15 #include "PndSensorNameContFact.h"
16 #include <vector>
17 #include <string>
18 #include "TROOT.h"
19 #include "TGeoVolume.h"
20 #include "TGeoShape.h"
21 #include "TGeoBBox.h"
22 #include "TList.h"
23 #include "TTree.h"
24 #include "FairMCEventHeader.h"
25 #include "FairParRootFileIo.h"
26 #include "FairBaseParSet.h"
27 
28 #include <math.h>
29 #include "stdlib.h"
30 
31 ClassImp(PndGeoHandling);
32 
33 
34 PndGeoHandling * PndGeoHandling::fInstance= NULL;
35 
36 PndGeoHandling* PndGeoHandling::Instance(){
37  if ( !fInstance){
38  std::cout<<"Info in (PndGeoHandling::Instance): Making a new instance using the framework."<<std::endl;
39  fInstance = new PndGeoHandling();
40  }
41  return fInstance;
42 }
43 
44 PndGeoHandling::PndGeoHandling():fVerbose(0),fGeoMan(),fSensorNamePar(),fRtdb(),fLevel(0),fFullPath(true),fRunId(0),fLevelNames()
45 {
46  if(fInstance) return;
47  fInstance = this;
48  FairRun* run = FairRun::Instance();
49 // if (!run) Fatal("PndGeoHandling","No FairRun object found. If used in a macro take another constructor.");
50 // run->AddTask((FairTask*)this);
51  if (run) {
52  this->SetName("PndGeoHandling");
53  this->SetTitle("FairTask");
54  run->AddTask((FairTask*)this);
55  }
56  else {
57  std::cout << "PndGeoHandling. No FairRun object found. If used in a macro take another constructor." << std::endl;
58  }
59 
60 }
61 
62 void PndGeoHandling::SetParContainers()
63 {
64  FairRun* run = FairRun::Instance();
65  if (!run) Fatal("PndGeoHandling","No FairRun object found.");
66  fRtdb = run->GetRuntimeDb();
67  if (!fRtdb) Fatal("PndGeoHandling","No runtime database found.");
68  fRunId = run->GetRunId();
69  if (fRunId < 0) Error("PndGeoHandling", "No valid run ID? %i",fRunId);
70  //if (!gGeoManager) GetGeoManager();
71  fGeoMan = gGeoManager;
72  if (!fGeoMan) Fatal("PndGeoHandling","No gGeoManager found.");
73  fSensorNamePar = (PndSensorNamePar*) (fRtdb->getContainer("PndSensorNamePar"));
74  if ( ! fSensorNamePar) Fatal("PndGeoHandling","No PndSensorNamePar parameters found.");
75  //fRtdb->initContainers(fRunId);
76  FairTask::SetParContainers();
77 }
78 
79 PndGeoHandling::PndGeoHandling(TString mcFile, TString parFile):fVerbose(0),fGeoMan(),fSensorNamePar(),fRtdb(),fLevel(0),fFullPath(true),fRunId(0),fLevelNames()
80 {
81  if(fInstance) return;
82  fInstance = this;
83  InitRuntimeDb(parFile);
84  GetRunId(mcFile);
85  if (gGeoManager) {
86  fGeoMan = gGeoManager;
87  } else
88  {
89  //Fatal("PndGeoHandling","No gGeoManager");
90  //return;
91  GetGeoManager();
92  fGeoMan = gGeoManager;
93  }
94 
95 
96  GetSensorNamePar();
97 }
98 
99 PndGeoHandling::PndGeoHandling(Int_t runId, TString parFile):fVerbose(0),fGeoMan(),fSensorNamePar(),fRtdb(),fLevel(0),fFullPath(true),fRunId(0),fLevelNames()
100 {
101  if(fInstance) return;
102  fInstance = this;
103  InitRuntimeDb(parFile);
104  fRunId = runId;
105  if (gGeoManager) {
106  fGeoMan = gGeoManager;
107  } else
108  {
109  //Fatal("PndGeoHandling","No gGeoManager");
110  //return;
111  GetGeoManager();
112  fGeoMan = gGeoManager;
113  }
114  GetSensorNamePar();
115 
116 }
117 
118 Int_t PndGeoHandling::GetRunId(TString mcFile)
119 {
120  TFile f(mcFile.Data());
121  TTree* t = (TTree*)f.Get("cbmsim");
122  FairMCEventHeader* header= new FairMCEventHeader();
123  t->SetBranchStatus("MCEventHeader.",1);
124  t->SetBranchAddress("MCEventHeader.", &header);
125  t->GetEntry(0);
126  fRunId = header->GetRunID();
127 
128  t->SetBranchStatus("MCEventHeader.",0);
129  return fRunId;
130 }
131 
132 void PndGeoHandling::GetSensorNamePar()
133 {
134  fSensorNamePar = (PndSensorNamePar*) (fRtdb->getContainer("PndSensorNamePar"));
135 
136  fRtdb->initContainers(fRunId);
137 
138  if (fVerbose > 1){
139  std::cout << "PndGeoHandling::GetSensorNamePar()" << std::endl;
140  fRtdb->Print();
141  fSensorNamePar->Print();
142  }
143 }
144 
145 void PndGeoHandling::InitRuntimeDb(TString parFileName)
146 {
147  fRtdb = FairRuntimeDb::instance();
148  FairParRootFileIo* parInput1 = new FairParRootFileIo(kTRUE);
149  parInput1->open(parFileName.Data(),"UPDATE");
150  fRtdb->setFirstInput(parInput1);
151  fRtdb->setOutput(parInput1);
152 }
153 
154 void PndGeoHandling::GetGeoManager()
155 {
156  if (fRunId < 0)
157  return;
158  FairBaseParSet* par=(FairBaseParSet*)(fRtdb->getContainer("FairBaseParSet"));
159  fRtdb->initContainers(fRunId);
160  if(fVerbose>0) par->Print();
161 }
162 
163 /*
164  TString PndGeoHandling::GetCurrentID()
165  {
166  Int_t level;
167  Int_t copyNr[100];
168  Int_t volNr[100];
169  TString result;
170 
171  level = fGeoMan->GetLevel();
172  level++;
173 
174  fGeoMan->GetBranchNumbers(copyNr, volNr);
175  for (int i=0; i<level; i++){
176  result += volNr[i];
177  result += "_";
178  result += copyNr[i];
179  result += "/";
180  }
181  return result;
182  }
183  */
184 
185 /*
186  TString PndGeoHandling::GetID(TString path)
187  {
188  TString result;
189  TString currentPath = fGeoMan->GetPath();
190  fGeoMan->cd(path.Data());
191  result = GetCurrentID();
192  fGeoMan->cd(currentPath.Data());
193  return result;
194  }
195  */
196 
197 Int_t PndGeoHandling::GetShortID(TString path)
198 {
199  TObjString myPath(path.Data());
200  if (fSensorNamePar != 0){
201 // std::cout << "PndGeoHandling::GetShortID: " << std::endl;
202 // fSensorNamePar->Print();
203  return fSensorNamePar->SensorInList(&myPath);
204  }
205  else
206  std::cout << "-E- PndGeoHandling::GetShortID: SensorNamePar is missing!" << std::endl;
207  return -1;
208 }
209 
210 TString PndGeoHandling::GetPath(Int_t shortID)
211 {
212  if (fSensorNamePar != 0){
213  return (fSensorNamePar->GetSensorName(shortID));
214  }
215  else {
216  std::cout << "-E- PndGeoHandling::GetPath(Int_t shortID): Missing SensorNamePar" << std::endl;
217  abort();
218  }
219 }
220 
221 /*
222  TString PndGeoHandling::GetPath(TString id)
223  {
224  TString result;
225  SetGeoManager(gGeoManager);
226  //std::cout << "-I- PndGeoHandling::GetPath : " << id.Data() << std::endl;
227  std::vector<std::string> idVector;
228  PndStringSeparator pathAna(id.Data(), "/_");
229  idVector = pathAna.GetStringVector();
230 
231  for(UInt_t i = 0; i < idVector.size(); i+=2){
232  result += "/";
233  Int_t VolId = atoi(idVector[i].c_str());
234  Int_t CopyNr = atoi(idVector[i+1].c_str());
235  // if(fVerbose>3) std::cout<<" -I- PndGeoHandling::GetPath: VolId = "<<VolId<<std::endl;
236  //std::cout << "-I- PndGeoHandling::GetPath : " << VolId;
237  //std::cout << " : " << fGeoMan->GetVolume(VolId)->GetName() << std::endl;
238  result += fGeoMan->GetVolume(VolId)->GetName();
239  result += "_";
240  result += CopyNr;
241  }
242  // if(fVerbose>2) std::cout<<" -I- PndGeoHandling::GetPath: result = "<<result.Data()<<std::endl;
243  return result;
244  }
245  */
246 
247 /*
248  Bool_t PndGeoHandling::cd(TString id)
249  {
250  return fGeoMan->cd(GetPath(id).Data());
251  }
252  */
253 Bool_t PndGeoHandling::cd(Int_t id)
254 {
255  return fGeoMan->cd(GetPath(id).Data());
256 }
257 
258 
259 TString PndGeoHandling::GetVolumeID(TString name)
260 {
261  TString result;
262  TGeoVolume* vol = fGeoMan->FindVolumeFast(name);
263  if (vol == 0)
264  return result;
265  result += vol->GetNumber();
266  return result;
267 }
268 
269 std::vector<TString> PndGeoHandling::GetNamesLevel(Int_t level, TString startPath, bool fullPath)
270 {
271  TString actPath = fGeoMan->GetPath();
272  fLevelNames.clear();
273 
274  if (startPath == ""){
275  fGeoMan->CdTop();
276  fLevel = level;
277  }
278  else{
279  if (fGeoMan->cd(startPath.Data())== 0)
280  return fLevelNames;
281  else {
282  fLevel = fGeoMan->GetLevel() + level;
283  }
284  }
285  FillLevelNames();
286  return fLevelNames;
287 }
288 
289 void PndGeoHandling::FillLevelNames()
290 {
291  TGeoNode* myNode = fGeoMan->GetCurrentNode();
292  if (fLevel == fGeoMan->GetLevel()){
293  if (fFullPath)
294  fLevelNames.push_back(fGeoMan->GetPath());
295  else
296  fLevelNames.push_back(myNode->GetName());
297  }
298  else {
299  Int_t nDaughters = myNode->GetNdaughters();
300  for (Int_t i = 0; i < nDaughters; i++){
301  fGeoMan->CdDown(i);
302  FillLevelNames();
303  fGeoMan->CdUp();
304  }
305  }
306 }
307 
308 void PndGeoHandling::GetOUVPath(TString path, TVector3& o, TVector3& u, TVector3& v)
309 {
310  Double_t result[3];
311  Double_t* temp;
312  TString actPath = fGeoMan->GetPath();
313  fGeoMan->cd(path);
314 
315  TGeoHMatrix* currMatrix = fGeoMan->GetCurrentMatrix();
316  temp = currMatrix->GetTranslation();
317  o.SetXYZ(temp[0], temp[1], temp[2]);
318 
319  temp[0] = 1;
320  temp[1] = 0;
321  temp[2] = 0;
322  fGeoMan->LocalToMasterVect(temp, result);
323  u.SetXYZ(result[0], result[1], result[2]);
324 
325  temp[0] = 0;
326  temp[1] = 1;
327  temp[2] = 0;
328  fGeoMan->LocalToMasterVect(temp, result);
329  v.SetXYZ(result[0], result[1], result[2]);
330 
331  if(actPath!="" && actPath!=" ") fGeoMan->cd(actPath);
332 }
333 
334 /*
335  void PndGeoHandling::GetOUVId(TString id, TVector3& o, TVector3& u, TVector3& v)
336  {
337  GetOUVPath(GetPath(id),o,u,v);
338  }
339  */
340 
341 TVector3 PndGeoHandling::GetSensorDimensionsPath(TString path)
342 {
343  TVector3 dim;
344  TString actPath = fGeoMan->GetPath();
345  fGeoMan->cd(path);
346 
347  TGeoVolume* actVolume = gGeoManager->GetCurrentVolume();
348  TGeoBBox* actBox = (TGeoBBox*)(actVolume->GetShape());
349  dim.SetX(actBox->GetDX());
350  dim.SetY(actBox->GetDY());
351  dim.SetZ(actBox->GetDZ());
352 
353  if(actPath!="" && actPath!=" ") fGeoMan->cd(actPath);
354  return dim;
355 }
356 
357 /*
358  TVector3 PndGeoHandling::GetSensorDimensionsId(TString id)
359  {
360  return GetSensorDimensionsPath(GetPath(id));
361  }
362  */
363 
364 TGeoHMatrix* PndGeoHandling::GetMatrixPath(TString path)
365 {
366  TString actPath = fGeoMan->GetPath();
367  fGeoMan->cd(path);
368 
369  TGeoHMatrix* currMatrix = fGeoMan->GetCurrentMatrix();
370  if(actPath!="" && actPath!=" ") fGeoMan->cd(actPath);
371 
372  return currMatrix;
373 
374 }
375 
376 /*
377  TGeoHMatrix* PndGeoHandling::GetMatrixId(TString id)
378  {
379  return GetMatrixPath(GetPath(id));
380  }
381  */
382 
383 // ----- conversions of POINTS (not vectors) here -----
384 /*
385  TVector3 PndGeoHandling::MasterToLocalId(const TVector3& master, const TString& id)
386  { return MasterToLocalPath(master, GetPath(id) ); }
387  */
388 
389 TVector3 PndGeoHandling::MasterToLocalPath(const TVector3& master, const TString& path)
390 {
391  // if(fVerbose>1) std::cout<<" -I- PndGeoHandling::MasterToLocalPath"<<std::endl;
392  Double_t result[3];
393  Double_t temp[3];
394 
395  temp[0] = master.X();
396  temp[1] = master.Y();
397  temp[2] = master.Z();
398 
399  TString actPath = fGeoMan->GetPath();
400  fGeoMan->cd(path);
401  fGeoMan->MasterToLocal(temp, result);
402  if(actPath != "" && actPath != " ") fGeoMan->cd(actPath);
403  return TVector3(result[0],result[1],result[2]);
404 }
405 
406 
407 /*
408  TVector3 PndGeoHandling::LocalToMasterId(const TVector3& local, const TString& id)
409  { return LocalToMasterPath(local, GetPath(id) ); }
410  */
411 
412 TVector3 PndGeoHandling::LocalToMasterPath(const TVector3& local, const TString& path)
413 {
414  Double_t result[3];
415  Double_t temp[3];
416 
417  temp[0] = local.X();
418  temp[1] = local.Y();
419  temp[2] = local.Z();
420 
421  TString actPath = fGeoMan->GetPath();
422  fGeoMan->cd(path);
423  fGeoMan->LocalToMaster(temp, result);
424  if(actPath != "" && actPath != " ") fGeoMan->cd(actPath);
425  return TVector3(result[0],result[1],result[2]);
426 }
427 
428 
429 // ROTATION of error values, CAUTION - these are always psitive defined
430 /*
431  TVector3 PndGeoHandling::MasterToLocalErrorsId(const TVector3& master, const TString& id)
432  { return MasterToLocalErrorsPath(master, GetPath(id) ); }
433  */
434 TMatrixD PndGeoHandling::MasterToLocalErrorsPath(const TMatrixD& master, const TString& path)
435 {
436  TString actPath = fGeoMan->GetPath();
437  fGeoMan->cd(path);
438  TMatrixD rot=GetCurrentRotationMatrix();
439  TMatrixD result = rot;
440  result*=master;
441  rot.T();
442  result*=rot;
443  if(fVerbose>1){
444  std::cout<<" -I- PndGeoHandling::MasterToLocalErrorsPath: print matrices: master, rotation, result=R*M*R^T"<<std::endl;
445  master.Print();
446  rot.Print();
447  result.Print();
448  }
449  if(actPath != "" && actPath != " ") fGeoMan->cd(actPath);
450  return result;
451 }
452 
453 
454 /*
455  TVector3 PndGeoHandling::LocalToMasterErrorsId(const TVector3& local, const TString& id)
456  { return LocalToMasterErrorsPath(local, GetPath(id) ); }
457  */
458 
459 TMatrixD PndGeoHandling::LocalToMasterErrorsPath(const TMatrixD& local, const TString& path)
460 {
461  TString actPath = fGeoMan->GetPath();
462  fGeoMan->cd(path);
463  TMatrixD rot=GetCurrentRotationMatrix();
464  TMatrixD result = rot;
465  result.T();
466  result*=local;
467  result*=rot;
468  if(fVerbose>1){
469  std::cout<<" -I- PndGeoHandling::LocalToMasterErrorsPath: print matrices: master, rotation, result=R^T*M*R"<<std::endl;
470  local.Print();
471  rot.Print();
472  result.Print();
473  }
474  if(actPath != "" && actPath != " ") fGeoMan->cd(actPath);
475  return result;
476 }
477 
478 TMatrixD PndGeoHandling::GetCurrentRotationMatrix()
479 {
480  TMatrixD rot(3,3,fGeoMan->GetCurrentMatrix()->GetRotationMatrix());
481 // TMatrixD rot(3,3);
482 // rot[0][0] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[0];
483 // rot[0][1] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[1];
484 // rot[0][2] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[2];
485 // rot[1][0] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[3];
486 // rot[1][1] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[4];
487 // rot[1][2] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[5];
488 // rot[2][0] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[6];
489 // rot[2][1] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[7];
490 // rot[2][2] = fGeoMan->GetCurrentMatrix()->GetRotationMatrix()[8];
491  return rot;
492 }
493 
494 
495 TString PndGeoHandling::FindNodePath(TGeoNode* node)
496 {
497  // Find a nodes full path by going there in the gGeoManager
498  // With many volumes this becomes surely slow.
499  const char* oldpath = fGeoMan->GetPath();
500  fGeoMan->CdTop(); // dive down from top node
501  DiveDownToNode(node);
502  TString pathname = fGeoMan->GetPath();
503  fGeoMan->cd(oldpath);
504  return pathname;
505 }
506 
507 void PndGeoHandling::DiveDownToNode(TGeoNode* node)
508 {
509  // cd gGeoManager from the current node to a given node
510  TGeoNode *currentNode = fGeoMan->GetCurrentNode();
511  if (currentNode == node) return;
512  for (Int_t iNod=0; iNod<currentNode->GetNdaughters();iNod++)
513  {
514  fGeoMan->CdDown(iNod);
515  DiveDownToNode(node);
516  if (fGeoMan->GetCurrentNode() == node) return;
517  fGeoMan->CdUp();
518  }
519 }
520 
521 void PndGeoHandling::DiveDownToNodeContainingString(TString name)
522 {
523  TGeoNode *currentNode = fGeoMan->GetCurrentNode();
524  TString nodeName(currentNode->GetName());
525  if (nodeName.Contains(name)){
526  return;
527  }
528  for (Int_t iNod = 0; iNod < currentNode->GetNdaughters(); iNod++) {
529  fGeoMan->CdDown(iNod);
530  DiveDownToNodeContainingString( name);
531  nodeName = fGeoMan->GetCurrentNode()->GetName();
532  if (nodeName.Contains(name)){
533  return;
534  }
535  fGeoMan->CdUp();
536  }
537 }
538 
539 void PndGeoHandling::DiveDownToFillSensNamePar(std::vector<std::string> listOfSensitives)
540 {
541  if (fSensorNamePar != 0){
542  TGeoNode *currentNode = fGeoMan->GetCurrentNode();
543  TString nodeName(currentNode->GetName());
544  //std::cout << nodeName.Data() << std::endl;
545 
546  if (VolumeIsSensitive(nodeName, listOfSensitives)){
547  PndStringSeparator sep(nodeName.Data(), "/");
548  std::vector<std::string> sepString = sep.GetStringVector();
549  if (sepString.size() > 0 && sepString[sepString.size() - 1].find("PartAss") == std::string::npos){
550  TObjString* myName = new TObjString(fGeoMan->GetPath());
551  fSensorNamePar->AddSensorName(myName);
552  }
553  }
554  for (Int_t iNod = 0; iNod < currentNode->GetNdaughters(); iNod++) {
555  fGeoMan->CdDown(iNod);
556  DiveDownToFillSensNamePar(listOfSensitives);
557 
558  fGeoMan->CdUp();
559  }
560  }
561  else
562  std::cout << "-E- PndMvdGeoHandling::DiveDownToFillSensNamePar: fSensorNamePar does not exist!" << std::endl;
563 }
564 
565 void PndGeoHandling::cd(TGeoNode* node)
566 {
567  // go to a node in the gGeoManager without knowing the full path
568  // With many volumes this becomes surely slow.
569  fGeoMan->CdTop(); // dive down from top node
570  DiveDownToNode(node);
571  return;
572 }
573 
574 void PndGeoHandling::CreateUniqueSensorId(TString startName, std::vector<std::string> listOfSensitives)
575 {
576  fGeoMan->CdTop();
577  DiveDownToNodeContainingString(startName);
578  if (fVerbose > 0)
579  std::cout << "-I- PndMvdGeoHandling::CreateUniqueSensorId: StartNode: " << fGeoMan->GetPath() << std::endl;
580  DiveDownToFillSensNamePar(listOfSensitives);
581 }
582 
583 
584 bool PndGeoHandling::VolumeIsSensitive(TString& path, std::vector<std::string>& listOfSensitives)
585 {
586  for (unsigned int i = 0; i < listOfSensitives.size(); i++){
587  if (path.Contains(listOfSensitives[i].c_str()))
588  return true;
589  }
590  return false;
591 }
592 
593 
594 
595 
596