MpdFfd.cxx 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdFfd
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include "iostream"
  8. #include "MpdMCTrack.h"
  9. #include "TClonesArray.h"
  10. #include "TLorentzVector.h"
  11. #include "TParticle.h"
  12. #include "TVirtualMC.h"
  13. #include "FairGeoInterface.h"
  14. #include "FairGeoLoader.h"
  15. #include "FairGeoNode.h"
  16. #include "FairGeoRootBuilder.h"
  17. #include "FairRootManager.h"
  18. #include "FairRun.h"
  19. #include "FairRuntimeDb.h"
  20. #include "FairVolume.h"
  21. #include "MpdStack.h"
  22. #include "TObjArray.h"
  23. #include "MpdFfd.h"
  24. #include "MpdFfdGeo.h"
  25. #include "MpdFfdGeoPar.h"
  26. #include "MpdFfdPoint.h"
  27. #include "MpdFfdHitProducer.h"
  28. class FairVolume;
  29. using namespace std;
  30. ClassImp(MpdFfd)
  31. //------------------------------------------------------------------------------------------------------------------------
  32. MpdFfd::MpdFfd()
  33. : FairDetector("FFD", kTRUE)
  34. {
  35. aFfdPoints = new TClonesArray("MpdFfdPoint");
  36. fPosIndex = 0;
  37. fVerboseLevel = 1;
  38. }
  39. //------------------------------------------------------------------------------------------------------------------------
  40. MpdFfd::MpdFfd(const char *name, Bool_t active, Int_t verbose)
  41. : FairDetector(name, active)
  42. {
  43. aFfdPoints = new TClonesArray("MpdFfdPoint");
  44. fPosIndex = 0;
  45. fVerboseLevel = verbose;
  46. }
  47. //------------------------------------------------------------------------------------------------------------------------
  48. MpdFfd::~MpdFfd()
  49. {
  50. aFfdPoints->Delete();
  51. delete aFfdPoints;
  52. }
  53. //------------------------------------------------------------------------------------------------------------------------
  54. MpdFfdPoint* MpdFfd::FindPoint(Int_t ptid, Int_t suid)
  55. {
  56. if(cPtid == ptid && cSuid == suid) return cPoint; // cached result
  57. // looking for new pair <tid, suid>
  58. for(int i = 0, N = aFfdPoints->GetEntriesFast(); i < N; i++) // cycle already created MpdFfdPoint
  59. {
  60. auto entry = (MpdFfdPoint*) aFfdPoints->UncheckedAt(i);
  61. if(entry->IsSame(ptid, suid))
  62. {
  63. cPoint = entry; // update current point
  64. cPtid = ptid;
  65. cSuid = suid;
  66. return cPoint;
  67. }
  68. }
  69. return nullptr;
  70. }
  71. //------------------------------------------------------------------------------------------------------------------------
  72. MpdFfdPoint* MpdFfd::CreatePoint(Int_t tid, Int_t suid)
  73. {
  74. return new((*aFfdPoints)[aFfdPoints->GetEntriesFast()]) MpdFfdPoint(tid, suid);
  75. }
  76. //------------------------------------------------------------------------------------------------------------------------
  77. Bool_t MpdFfd::ProcessHits(FairVolume *vol)
  78. {
  79. Int_t suid;
  80. gMC->CurrentVolID(suid);
  81. Int_t tid = gMC->GetStack()->GetCurrentTrackNumber();
  82. Int_t pid = gMC->TrackPid(); // Cherenkov 50000050 FeedbackPhoton 50000051
  83. if(gMC->IsTrackEntering())
  84. {
  85. if(pid != 50000050) // event = non-op track ENTER to quartz volume
  86. {
  87. gMC->TrackPosition(fPos);
  88. gMC->TrackMomentum(fMom);
  89. auto ret = FindTrackParams(tid, suid);
  90. if(ret.second == false) // don't found
  91. {
  92. TrackParam param(suid, fPos.Vect(), fMom, gMC->TrackTime() * 1.e+9, gMC->TrackLength());
  93. params.insert(make_pair(tid, param)); // insert new parent track parameters
  94. }
  95. }
  96. else // event = op create at quartz volume
  97. {
  98. TLorentzVector mom;
  99. gMC->TrackMomentum(mom);
  100. double energy = mom.E()* 1.e+9; // op energy [eV]
  101. bool addEntry = true;
  102. if(MpdFfdPoint::fCurrentMode == MpdFfdPoint::kPhotoElectron) // current mode = save pe only
  103. {
  104. if(! MpdFfdHitProducer::IsPeCreated(energy)) addEntry = false; // failed create pe from op
  105. }
  106. if(addEntry)
  107. {
  108. Int_t ptid = gMC->GetStack()->GetCurrentParentTrackNumber();
  109. auto point = FindPoint(ptid, suid);
  110. if(point == nullptr) point = CreatePoint(ptid, suid);
  111. if(! point->AddOp(energy, gMC->TrackTime() * 1.e+9)) // [ns]; if MpdFfdPoint not closed, need to call SaveParentTrackParams ONCE!
  112. {
  113. auto ret = FindTrackParams(ptid, suid);
  114. if(ret.second)
  115. {
  116. auto param = ret.first->second;
  117. point->SaveParentTrackParams(param.posIn, param.posOut, param.mom, param.time, param.length);
  118. ((MpdStack *)gMC->GetStack())->AddPoint(kFFD, ptid); // add marker for parent track
  119. }
  120. }
  121. } // addEntry
  122. } // op
  123. } // track entering
  124. // Update MpdFfdPoint at non-op track EXIT from quartz
  125. if(pid != 50000050 && (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) )
  126. {
  127. TLorentzVector posOut;
  128. gMC->TrackPosition(posOut);
  129. auto ret = FindTrackParams(tid, suid);
  130. if(ret.second) ret.first->second.posOut = posOut.Vect(); // update parent track output position
  131. ResetParameters();
  132. } // track exiting
  133. return kTRUE;
  134. }
  135. //------------------------------------------------------------------------------------------------------------------------
  136. pair <MpdFfd::Tparams::iterator, bool> MpdFfd::FindTrackParams(Int_t tid, Int_t suid)
  137. {
  138. auto ret = params.equal_range(tid);
  139. for(auto it = ret.first; it != ret.second; ++it)
  140. {
  141. if(it->second.suid == suid)
  142. {
  143. return make_pair(it, true);
  144. }
  145. }
  146. return make_pair(params.end(), false);
  147. }
  148. //------------------------------------------------------------------------------------------------------------------------
  149. void MpdFfd::EndOfEvent()
  150. {
  151. if(fVerboseLevel) Print();
  152. aFfdPoints->Delete();
  153. fPosIndex = 0;
  154. // event data cleanup
  155. params.clear();
  156. cPoint = nullptr;
  157. cPtid = cSuid = -1;
  158. }
  159. //------------------------------------------------------------------------------------------------------------------------
  160. void MpdFfd::Register()
  161. {
  162. FairRootManager::Instance()->Register("FFDPoint", "Ffd", aFfdPoints, kTRUE);
  163. }
  164. //------------------------------------------------------------------------------------------------------------------------
  165. TClonesArray *MpdFfd::GetCollection(Int_t iColl) const
  166. {
  167. if(iColl == 0) return aFfdPoints;
  168. return nullptr;
  169. }
  170. //------------------------------------------------------------------------------------------------------------------------
  171. void MpdFfd::Print() const
  172. {
  173. Int_t nHits = aFfdPoints->GetEntriesFast();
  174. cout<<"-I- MpdFfd: "<<nHits<<" points registered in this event.\n";
  175. if(fVerboseLevel > 1) for (Int_t i = 0; i < nHits; i++) (*aFfdPoints)[i]->Print();
  176. }
  177. //------------------------------------------------------------------------------------------------------------------------
  178. void MpdFfd::Reset()
  179. {
  180. aFfdPoints->Delete();
  181. ResetParameters();
  182. }
  183. //------------------------------------------------------------------------------------------------------------------------
  184. void MpdFfd::CopyClones(TClonesArray *cl1, TClonesArray *cl2, Int_t offset)
  185. {
  186. Int_t nEntries = cl1->GetEntriesFast();
  187. cout << "-I- MpdFfd: " << nEntries << " entries to add." << endl;
  188. TClonesArray &clref = *cl2;
  189. MpdFfdPoint *oldpoint = nullptr;
  190. for (Int_t i = 0; i < nEntries; i++) {
  191. oldpoint = (MpdFfdPoint *)cl1->At(i);
  192. Int_t index = oldpoint->GetTrackID() + offset;
  193. oldpoint->SetTrackID(index);
  194. new (clref[fPosIndex]) MpdFfdPoint(*oldpoint);
  195. fPosIndex++;
  196. }
  197. cout << "-I- MpdFfd: " << cl2->GetEntriesFast() << " merged entries." << endl;
  198. }
  199. //------------------------------------------------------------------------------------------------------------------------
  200. void MpdFfd::ConstructGeometry()
  201. {
  202. TString fileName = GetGeometryFileName();
  203. if (fileName.EndsWith(".root")) {
  204. LOG(INFO) << "Constructing FFD geometry from ROOT file " << fileName.Data()
  205. << endl;
  206. ConstructRootGeometry();
  207. } else if (fileName.EndsWith(".geo")) {
  208. LOG(INFO) << "Constructing FFD geometry from ASCII file " << fileName.Data()
  209. << endl;
  210. ConstructAsciiGeometry();
  211. }
  212. /*else if ( fileName.EndsWith(".gdml") )
  213. {
  214. LOG(INFO) << "Constructing CPC geometry from GDML file " <<
  215. fileName.Data() << endl; ConstructGDMLGeometry();
  216. }*/
  217. else {
  218. LOG(FATAL) << "Geometry format of FFD file " << fileName.Data()
  219. << " not supported." << endl;
  220. }
  221. }
  222. //------------------------------------------------------------------------------------------------------------------------
  223. void MpdFfd::ConstructAsciiGeometry()
  224. {
  225. int count = 0;
  226. int count_tot = 0;
  227. FairGeoLoader *geoLoad = FairGeoLoader::Instance();
  228. FairGeoInterface *geoFace = geoLoad->getGeoInterface();
  229. MpdFfdGeo *Geo = new MpdFfdGeo();
  230. Geo->setGeomFile(GetGeometryFileName());
  231. geoFace->addGeoModule(Geo);
  232. Bool_t rc = geoFace->readSet(Geo);
  233. if (rc)
  234. Geo->create(geoLoad->getGeoBuilder());
  235. else
  236. std::cerr << "FfdDetector:: geometry could not be read!" << std::endl;
  237. TList *volList = Geo->getListOfVolumes();
  238. // store geo parameter
  239. FairRun *fRun = FairRun::Instance();
  240. FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
  241. MpdFfdGeoPar *par = (MpdFfdGeoPar *)(rtdb->getContainer("MpdFfdGeoPar"));
  242. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  243. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  244. TListIter iter(volList);
  245. FairGeoNode *node = nullptr;
  246. FairGeoVolume *aVol = nullptr;
  247. while ((node = (FairGeoNode *)iter.Next())) {
  248. aVol = dynamic_cast<FairGeoVolume *>(node);
  249. if (node->isSensitive()) {
  250. fSensNodes->AddLast(aVol);
  251. count++;
  252. } else
  253. fPassNodes->AddLast(aVol);
  254. count_tot++;
  255. }
  256. par->setChanged();
  257. par->setInputVersion(fRun->GetRunId(), 1);
  258. ProcessNodes(volList);
  259. }
  260. //------------------------------------------------------------------------------------------------------------------------
  261. Bool_t MpdFfd::CheckIfSensitive(std::string name)
  262. {
  263. if(name.find("Active") != string::npos) return kTRUE;
  264. return kFALSE;
  265. }
  266. //------------------------------------------------------------------------------------------------------------------------