MpdTof.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdTof, LRectangle, LStrip
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <iostream>
  8. #include <algorithm>
  9. #include <cmath>
  10. #include <TClonesArray.h>
  11. #include <TVirtualMC.h>
  12. #include <TGeoManager.h>
  13. #include <TGeoBBox.h>
  14. #include <TGeoMatrix.h>
  15. #include "FairGeoInterface.h"
  16. #include "FairGeoLoader.h"
  17. #include "FairGeoNode.h"
  18. #include "FairGeoRootBuilder.h"
  19. #include "FairRootManager.h"
  20. #include "MpdStack.h"
  21. #include "FairRuntimeDb.h"
  22. #include "FairRunAna.h"
  23. #include "FairVolume.h"
  24. #include "MpdMCTrack.h"
  25. #include "MpdTofGeo.h"
  26. #include "MpdTofGeoPar.h"
  27. #include "MpdTofPoint.h"
  28. #include "MpdTofHit.h"
  29. #include "MpdTof.h"
  30. using namespace std;
  31. ClassImp(MpdTof)
  32. //------------------------------------------------------------------------------------------------------------------------
  33. //------------------------------------------------------------------------------------------------------------------------
  34. //------------------------------------------------------------------------------------------------------------------------
  35. MpdTof::MpdTof(const char* name, Bool_t active)
  36. : FairDetector(name, active)
  37. {
  38. aTofHits = new TClonesArray("MpdTofPoint");
  39. fVerboseLevel = 1;
  40. }
  41. //------------------------------------------------------------------------------------------------------------------------
  42. MpdTof::~MpdTof()
  43. {
  44. aTofHits->Delete();
  45. delete aTofHits;
  46. }
  47. //------------------------------------------------------------------------------------------------------------------------
  48. void MpdTof::ResetParameters()
  49. {
  50. fTrackID = fVolumeID = 0;
  51. fPos.SetXYZM(NAN, NAN, NAN, NAN);
  52. fMom = fPos;
  53. fTime = fLength = fELoss = NAN;
  54. fPosIndex = 0;
  55. }
  56. //------------------------------------------------------------------------------------------------------------------------
  57. Bool_t MpdTof::ProcessHits(FairVolume* vol)
  58. {
  59. Int_t sector, detector, strip;
  60. // Set parameters at entrance of volume. Reset ELoss.
  61. if(gMC->IsTrackEntering())
  62. {
  63. fELoss = 0.;
  64. fTime = gMC->TrackTime() * 1.0e09;
  65. fLength = gMC->TrackLength();
  66. gMC->TrackPosition(fPos);
  67. gMC->TrackMomentum(fMom);
  68. }
  69. // Sum energy loss for all steps in the active volume
  70. fELoss += gMC->Edep();
  71. // Create MpdTofPoint at exit of active volume
  72. if(fELoss > 0 && (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) )
  73. {
  74. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  75. gMC->CurrentVolOffID(0, strip); // [1,72]
  76. gMC->CurrentVolOffID(1, detector); // [1,20]
  77. gMC->CurrentVolOffID(2, sector); // [1,14]
  78. fVolumeID = MpdTofPoint::GetSuid72(sector, detector, strip);
  79. AddPoint(fTrackID, fVolumeID, fPos.Vect(), fMom.Vect(), fTime, fLength, fELoss);
  80. ((MpdStack*)gMC->GetStack())->AddPoint(kTOF);
  81. ResetParameters();
  82. }
  83. return kTRUE;
  84. }
  85. //------------------------------------------------------------------------------------------------------------------------
  86. void MpdTof::EndOfEvent()
  87. {
  88. if(fVerboseLevel) Print();
  89. aTofHits->Delete();
  90. fPosIndex = 0;
  91. }
  92. //------------------------------------------------------------------------------------------------------------------------
  93. void MpdTof::Register()
  94. {
  95. FairRootManager::Instance()->Register("TOFPoint", "Tof", aTofHits, kTRUE);
  96. }
  97. //------------------------------------------------------------------------------------------------------------------------
  98. TClonesArray* MpdTof::GetCollection(Int_t iColl) const
  99. {
  100. if(iColl == 0) return aTofHits;
  101. return nullptr;
  102. }
  103. //------------------------------------------------------------------------------------------------------------------------
  104. void MpdTof::Print() const
  105. {
  106. Int_t nHits = aTofHits->GetEntriesFast();
  107. cout << "-I- MpdTof: " << nHits << " points registered in this event." << endl;
  108. if(fVerboseLevel > 1)
  109. for(Int_t i=0; i<nHits; i++) (*aTofHits)[i]->Print();
  110. }
  111. //------------------------------------------------------------------------------------------------------------------------
  112. void MpdTof::Dump(TClonesArray *aPoints, TClonesArray *aHits, TClonesArray *aTracks, const char* comment, ostream& os)
  113. {
  114. if(comment != nullptr) os<<comment;
  115. // Fill & sort MpdTofPoint&MpdTofHit by mctid key.
  116. Int_t nHits = aHits->GetEntriesFast();
  117. Int_t nPoints = aPoints->GetEntriesFast();
  118. Int_t nTracks = aTracks->GetEntriesFast();
  119. multimap<Int_t, Int_t> mmPoints, mmHits, mmTracks; // < mctid, poinIndex or hitIndex> < mctid parent, mctid>
  120. set<Int_t> sMcIndex;
  121. vector<Int_t> vec;
  122. for(Int_t index=0; index < nHits; index++)
  123. {
  124. auto *pHit = (MpdTofHit*) aHits->At(index);
  125. vec.clear();
  126. pHit->getLinks(MpdTofUtils::mcTrackIndex, vec);
  127. for(auto it=vec.begin(), itEnd = vec.end(); it != itEnd; it++)
  128. {
  129. Int_t mctid = *it;
  130. mmHits.insert({mctid, index});
  131. sMcIndex.insert(mctid);
  132. }
  133. }
  134. for(Int_t index=0; index < nPoints; index++)
  135. {
  136. auto *pPoint = (MpdTofPoint*) aPoints->At(index);
  137. Int_t mctid = pPoint->GetTrackID();
  138. mmPoints.insert({mctid, index});
  139. sMcIndex.insert(mctid);
  140. }
  141. for(Int_t index=0; index < nTracks; index++)
  142. {
  143. auto pMCtrack = (MpdMCTrack*) aTracks->UncheckedAt(index);
  144. auto ptid = pMCtrack->GetMotherId();
  145. mmTracks.insert({ptid, index});
  146. }
  147. // Print to ostream os.
  148. os<<"\n[MpdTof::Dump] points: "<<nPoints<<", hits: "<<nHits<<" --------------------------------------------------------->>>";
  149. for(auto mctid : sMcIndex)
  150. {
  151. os<<"\ntid="<<mctid<<" points: ("<<mmPoints.count(mctid)<<") ";
  152. auto range = mmPoints.equal_range(mctid);
  153. for(auto it = range.first; it != range.second; ++it) os<<" pid="<<it->second;
  154. os<<" hits: ("<<mmHits.count(mctid)<<") ";
  155. auto range2 = mmHits.equal_range(mctid);
  156. for(auto it = range2.first; it != range2.second; ++it)
  157. {
  158. auto hid = it->second;
  159. auto pHit = (MpdTofHit*) aHits->At(hid);
  160. Int_t suid = pHit->GetDetectorID();
  161. Int_t sector, detector, gap, strip;
  162. MpdTofPoint::ParseSuid(suid, sector, detector, gap, strip);
  163. os<<" hid="<<hid<<"["<<suid<<"]{"<<sector<<","<<detector<<","<<gap<<","<<strip<<"}";
  164. }
  165. }
  166. os<<"\nTrack tree ---------------------------------------------------------------------------------------------------------";
  167. for(Int_t index=0; index < nTracks; index++)
  168. {
  169. if(mmTracks.find(index) != mmTracks.end()) continue; // pass only ended track line
  170. os<<"\n tid="<<index;
  171. auto pMCtrack = (MpdMCTrack*) aTracks->At(index);
  172. do
  173. {
  174. Int_t ptid = pMCtrack->GetMotherId();
  175. if(ptid != -1)
  176. {
  177. os<<"<tid="<<ptid;
  178. pMCtrack = (MpdMCTrack*) aTracks->At(ptid);
  179. }
  180. else pMCtrack = nullptr;
  181. }
  182. while(pMCtrack);
  183. }
  184. os<<"\n[MpdTof::Dump] --------------------------------------------------------------------------------------------------<<<";
  185. }
  186. //------------------------------------------------------------------------------------------------------------------------
  187. void MpdTof::Reset()
  188. {
  189. aTofHits->Delete();
  190. ResetParameters();
  191. }
  192. //------------------------------------------------------------------------------------------------------------------------
  193. void MpdTof::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
  194. {
  195. Int_t nEntries = cl1->GetEntriesFast();
  196. cout << "-I- MpdTof: " << nEntries << " entries to add." << endl;
  197. TClonesArray& clref = *cl2;
  198. MpdTofPoint* oldpoint = nullptr;
  199. for(Int_t i=0; i<nEntries; i++)
  200. {
  201. oldpoint = (MpdTofPoint*) cl1->At(i);
  202. Int_t index = oldpoint->GetTrackID() + offset;
  203. oldpoint->SetTrackID(index);
  204. new (clref[fPosIndex]) MpdTofPoint(*oldpoint);
  205. fPosIndex++;
  206. }
  207. cout << "-I- MpdTof: " << cl2->GetEntriesFast() << " merged entries." << endl;
  208. }
  209. //--------------------------------------------------------------------------------------------------------------------------------------
  210. void MpdTof::ConstructGeometry()
  211. {
  212. TString fileName = GetGeometryFileName();
  213. if(fileName.EndsWith(".root"))
  214. {
  215. LOG(INFO)<<"Constructing TOF geometry from ROOT file "<<fileName.Data();
  216. ConstructRootGeometry();
  217. }
  218. else if ( fileName.EndsWith(".geo") )
  219. {
  220. LOG(INFO)<<"Constructing TOF geometry from ASCII file "<<fileName.Data();
  221. ConstructAsciiGeometry();
  222. }
  223. else LOG(FATAL)<<"Geometry format of TOF file "<<fileName.Data()<<" not supported.";
  224. }
  225. //------------------------------------------------------------------------------------------------------------------------
  226. void MpdTof::ConstructAsciiGeometry()
  227. {
  228. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  229. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  230. MpdTofGeo* tofGeo = new MpdTofGeo();
  231. tofGeo->setGeomFile(GetGeometryFileName());
  232. geoFace->addGeoModule(tofGeo);
  233. Bool_t rc = geoFace->readSet(tofGeo);
  234. if(rc) tofGeo->create(geoLoad->getGeoBuilder());
  235. TList* volList = tofGeo->getListOfVolumes();
  236. // store geo parameter
  237. FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
  238. MpdTofGeoPar* par =(MpdTofGeoPar*)(rtdb->getContainer("MpdTofGeoPar"));
  239. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  240. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  241. FairGeoNode *node = nullptr;
  242. FairGeoVolume *aVol = nullptr;
  243. TListIter iter(volList);
  244. while((node = (FairGeoNode*)iter.Next()))
  245. {
  246. aVol = dynamic_cast<FairGeoVolume*> (node);
  247. if(node->isSensitive()) fSensNodes->AddLast(aVol);
  248. else fPassNodes->AddLast(aVol);
  249. }
  250. par->setChanged();
  251. par->setInputVersion(FairRun::Instance()->GetRunId(), 1);
  252. ProcessNodes(volList);
  253. }
  254. //------------------------------------------------------------------------------------------------------------------------
  255. MpdTofPoint* MpdTof::AddPoint(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eLoss)
  256. {
  257. return new((*aTofHits)[aTofHits->GetEntriesFast()]) MpdTofPoint(trackID, detID, pos, mom, time, length, eLoss);
  258. }
  259. //--------------------------------------------------------------------------------------------------------------------------------------
  260. Bool_t MpdTof::CheckIfSensitive(string name)
  261. {
  262. if(name.find("Active") != string::npos) return kTRUE;
  263. return kFALSE;
  264. }
  265. //------------------------------------------------------------------------------------------------------------------------
  266. void MpdTof::Print(const TVector3& v, const char* comment, ostream& os)
  267. {
  268. if(comment != nullptr) os<<comment;
  269. os<<"("<<v.X()<<","<<v.Y()<<","<<v.Z()<<"; "<<v.Perp()<<","<<v.Mag()<<")";
  270. }
  271. //------------------------------------------------------------------------------------------------------------------------
  272. void MpdTof::GetDelta(const TVector3& mcPos, const TVector3& estPos, double& dev, double& devZ, double& devR, double& devPhi)
  273. {
  274. dev = (mcPos - estPos).Mag();
  275. devZ = mcPos.Z() - estPos.Z();
  276. devR = mcPos.Perp() - estPos.Perp();
  277. devPhi = sqrt(dev*dev - devZ*devZ - devR*devR);
  278. devPhi = (mcPos.Phi() > estPos.Phi()) ? devPhi : -1.* devPhi;
  279. }
  280. //------------------------------------------------------------------------------------------------------------------------