MpdEtofGeoUtils.cxx 7.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #include <assert.h>
  3. #include <iostream>
  4. #include <TGeoManager.h>
  5. #include "TClonesArray.h"
  6. #include "TH1D.h"
  7. #include "TH2D.h"
  8. #include "TGeoMatrix.h"
  9. #include "TGeoBBox.h"
  10. #include "TGeoNode.h"
  11. #include "FairRootManager.h"
  12. #include "FairLogger.h"
  13. #include "MpdTofPoint.h"
  14. #include "MpdEtofGeoUtils.h"
  15. using namespace std;
  16. MpdEtofGeoUtils* MpdEtofGeoUtils::instance = nullptr;
  17. //------------------------------------------------------------------------------------------------------------------------
  18. MpdEtofGeoUtils::MpdEtofGeoUtils()
  19. {
  20. }
  21. //------------------------------------------------------------------------------------------------------------------------
  22. void MpdEtofGeoUtils::FindNeighborStrips(Double_t thresh, bool forced)
  23. {
  24. /* static bool founded = false;
  25. if( !forced && founded) return; // neighbor strips already founded and filled
  26. size_t NR = 0, NL= 0;
  27. const LStrip *strip2; double distance;
  28. for(auto it1 = mStrips.begin(), itEnd1 = mStrips.end(); it1 != itEnd1 ; it1++) // cycle1 by strips
  29. {
  30. LStrip *strip1 = &(it1->second);
  31. for(auto it2 = mStrips.begin(), itEnd2 = mStrips.end(); it2 != itEnd2 ; it2++) // cycle2 by strips
  32. {
  33. strip2 = &(it2->second);
  34. // CATION: Ckeck only left and right sides(one row strips NOW)
  35. distance = strip1->Distance(LStrip::kRight, *strip2);
  36. if(h1) h1->Fill(distance);
  37. if(distance < thresh) // CAUTION: constant depends on the geometry layout(see h1TestDistance histo)
  38. {
  39. strip1->neighboring[LStrip::kRight] = strip2->volumeUID;
  40. if(h2) h2->Fill(strip1->stripID, strip2->stripID);
  41. }
  42. distance = strip1->Distance(LStrip::kLeft, *strip2);
  43. if(h1) h1->Fill(distance);
  44. if(distance < thresh) // CAUTION: constant depends on the geometry layout(see h1TestDistance histo)
  45. {
  46. strip1->neighboring[LStrip::kLeft] = strip2->volumeUID;
  47. if(h2) h2->Fill( strip2->stripID, strip1->stripID);
  48. }
  49. }// cycle2 by strips
  50. }// cycle1 by strips
  51. founded = true;
  52. cout<<" [MpdEtofGeoUtils::FindNeighborStrips] Neighbor strips: left = "<<NL<<", right = "<<NR<<endl;
  53. */
  54. }
  55. //------------------------------------------------------------------------------------------------------------------------
  56. void MpdEtofGeoUtils::ParseTGeoManager(bool forced)
  57. {
  58. assert(gGeoManager);
  59. //static const double degree180 = 3.14159265359; // 180 degree to rad
  60. if( !forced && !mStrips.empty()) return; // already parsed and filled
  61. /*
  62. mStrips.clear();
  63. int nBoxes = 0, nStrips = 0, nGases = 0;
  64. const char* pathToTOF[] = {"/cave_1/etof1_1", "/cave_1/etof1_2"};
  65. for(size_t side = 0, sectorID = 1; side <= 1 ; side++, sectorID++)
  66. {
  67. gGeoManager->cd(pathToTOF[side]);
  68. TGeoNode *boxNode, *stripNode, *gasNode;
  69. Int_t volumeUID, boxID, stripID, gasID; // module[1,...,30], pad [1,...,131]
  70. vector<intervalType> vDetectorsR, vDetectorsPhi;
  71. Double_t *X0Y0Z0 = new Double_t[3]; X0Y0Z0[0] = X0Y0Z0[1] = X0Y0Z0[2] = 0.; // center of sensetive detector
  72. Double_t *local = new Double_t[3], master[3], dX, dY, dZ;
  73. TObjArray *array = gGeoManager->GetCurrentVolume()->GetNodes();
  74. TIterator *it1 = array->MakeIterator();
  75. while( (boxNode = (TGeoNode*) it1->Next()) ) // BOXES
  76. {
  77. TString PATH1 = TString(pathToTOF[side]) + "/" + boxNode->GetName(); boxID = boxNode->GetNumber(); nBoxes++;
  78. ///cout<<"\n MODULE: "<<boxNode->GetName()<<", copy# "<<boxID<<" path= "<<PATH1.Data();
  79. TIterator *it2 = boxNode->GetNodes()->MakeIterator();
  80. while( (stripNode = (TGeoNode*) it2->Next()) ) // STRIPS
  81. {
  82. TString stripName = stripNode->GetName();
  83. TString PATH2 = PATH1 + "/" + stripName; nStrips++;
  84. TString stripid = stripName; stripid.Remove(0, 7); stripid.Remove(stripid.Last('_')); stripID = stripid.Atoi();
  85. ///cout<<"\n\t PAD: "<<stripName.Data()<<", copy# "<<stripID<<" path= "<<PATH2.Data();
  86. TIterator *it3 = stripNode->GetNodes()->MakeIterator();
  87. while( (gasNode = (TGeoNode*) it3->Next()) ) // GASES
  88. {
  89. TString detName = gasNode->GetName();
  90. if(!detName.Contains("gas")) continue;
  91. TString PATH3 = PATH2 + "/" + detName; nGases++;
  92. TString gasid = stripName; gasid.Remove(0, 7); gasid.Remove(gasid.Last('_')); gasID = gasid.Atoi();
  93. ///cout<<"\n\t\t GAS: "<<detName.Data()<<", copy# "<<gasID<<" path= "<<PATH3.Data();
  94. gGeoManager->cd(PATH3);
  95. TGeoMatrix *matrix = gGeoManager->GetCurrentMatrix(); // calculate global TGeoHMatrix for current node
  96. matrix->LocalToMaster(X0Y0Z0, master); // 0.0.0. --> MRS
  97. TGeoBBox *box = (TGeoBBox*) gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
  98. dX = box->GetDX(); dY = box->GetDY(); dZ = box->GetDZ();
  99. volumeUID = MpdTofPoint::GetSuid(sectorID, boxID, 1, stripID); // FIXME: now ONE detector per box(detector=1)
  100. ///cout<<"\n uid="<<volumeUID<<", center= "<<master[0]<<", "<<master[1]<<", "<<master[2]<<", dXYZ= "<<dX<<", "<<dY<<", "<<dZ;
  101. LStrip stripData(volumeUID, sectorID, boxID, 1, stripID);
  102. stripData.center.SetXYZ(master[0], master[1], master[2]);
  103. // point A
  104. local[0] = -dX; local[1] = -dY; local[2] = +dZ;
  105. matrix->LocalToMaster(local, master);
  106. stripData.A.SetXYZ(master[0], master[1], master[2]);
  107. // point B
  108. local[0] = +dX; local[1] = -dY; local[2] = +dZ;
  109. matrix->LocalToMaster(local, master);
  110. stripData.B.SetXYZ(master[0], master[1], master[2]);
  111. // point C
  112. local[0] = +dX; local[1] = -dY; local[2] = -dZ;
  113. matrix->LocalToMaster(local, master);
  114. stripData.C.SetXYZ(master[0], master[1], master[2]);
  115. // point D
  116. local[0] = -dX; local[1] = -dY; local[2] = -dZ;
  117. matrix->LocalToMaster(local, master);
  118. stripData.D.SetXYZ(master[0], master[1], master[2]);
  119. stripData.InitCenterPerp();
  120. if(h2)
  121. {
  122. h2->Fill(stripData.A.Perp(), stripData.A.Phi());
  123. h2->Fill(stripData.B.Perp(), stripData.B.Phi());
  124. h2->Fill(stripData.C.Perp(), stripData.C.Phi());
  125. h2->Fill(stripData.D.Perp(), stripData.D.Phi());
  126. }
  127. ///stripData.Dump("\n ETOFstrip:");
  128. bool IsUniqueUID = mStrips.insert(make_pair(volumeUID, stripData)).second;
  129. assert(IsUniqueUID);
  130. // Add strip location interval to intervalTree
  131. LRectangle * pStrip = new LRectangle(stripData);
  132. pStrip->InitCenterPerp();
  133. double Rmin, Rmax, Phimin, Phimax;
  134. pStrip->GetRPhiRanges(Rmin, Rmax, Phimin, Phimax);
  135. vDetectorsR.push_back(intervalType(Rmin, Rmax, pStrip));
  136. if(Phimax - Phimin > degree180) // segment overlap the boundary [-180.,+180[
  137. {
  138. vDetectorsPhi.push_back(intervalType(-degree180, Phimin, pStrip));
  139. vDetectorsPhi.push_back(intervalType(Phimax, degree180, pStrip));
  140. ///cout<<"\n OVERLAPED phimin="<<Phimin<<", phimax="<<Phimax;
  141. }
  142. else vDetectorsPhi.push_back(intervalType(Phimin, Phimax, pStrip));
  143. } // GASES
  144. } // STRIPS
  145. } // MODULES
  146. mDetectorsR[side] = intervalTreeType(vDetectorsR);
  147. mDetectorsPhi[side] = intervalTreeType(vDetectorsPhi);
  148. }
  149. FairLogger::GetLogger()->Info(MESSAGE_ORIGIN, "[MpdEtofGeoUtils::ParseTGeoManager] Regions= %d, modules= %d, strips= %d. ", 2, nBoxes, nStrips);
  150. */
  151. }
  152. //------------------------------------------------------------------------------------------------------------------------
  153. const LStrip* MpdEtofGeoUtils::FindStrip(Int_t UID)
  154. {
  155. auto cit = mStrips.find(UID);
  156. assert(cit != mStrips.end());
  157. return &(cit->second);
  158. }
  159. //----------------------------------------------------------------------------------------------------------------------