MpdEmcKI.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcKI //
  4. // Cluster production for EMC //
  5. // Author List : D.Peresunko., RRCKI, 2019 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include "MpdEmcKI.h"
  9. #include "MpdEmcGeo.h"
  10. #include "MpdEmcGeoPar.h"
  11. #include "MpdEmcGeoUtils.h"
  12. #include "MpdEmcPointKI.h"
  13. #include "FairGeoInterface.h"
  14. #include "FairGeoLoader.h"
  15. #include "FairGeoNode.h"
  16. #include "FairGeoVolume.h"
  17. #include "FairRootManager.h"
  18. #include "FairRun.h"
  19. #include "FairRuntimeDb.h"
  20. #include "MpdStack.h"
  21. #include "FairVolume.h"
  22. #include "MpdDetectorList.h"
  23. #include "TClonesArray.h"
  24. #include "TGeoManager.h"
  25. #include "TGeoNode.h"
  26. #include "TGeoTube.h"
  27. #include "TParticle.h"
  28. #include "TVirtualMC.h"
  29. #include <iostream>
  30. using std::cout;
  31. using std::endl;
  32. MpdEmcKI::MpdEmcKI()
  33. : FairDetector("EMC", kTRUE, kECAL),
  34. fNhits(0),
  35. fCurrentTrackID(-1),
  36. fCurrentCellID(-1),
  37. fCurentSuperParent(-1),
  38. fCurrentHit(nullptr),
  39. fGeom(nullptr),
  40. fEcalRmin(0),
  41. fEcalRmax(0)
  42. {
  43. }
  44. MpdEmcKI::MpdEmcKI(const char* name, Bool_t active)
  45. : FairDetector(name, active, kECAL),
  46. fNhits(0),
  47. fCurrentTrackID(-1),
  48. fCurrentCellID(-1),
  49. fCurentSuperParent(-1),
  50. fCurrentHit(nullptr),
  51. fGeom(nullptr)
  52. {
  53. fMpdEmcPointCollection = new TClonesArray("MpdEmcPointKI");
  54. fVerboseLevel = 1;
  55. }
  56. MpdEmcKI::~MpdEmcKI()
  57. {
  58. if (fMpdEmcPointCollection) {
  59. fMpdEmcPointCollection->Delete();
  60. delete fMpdEmcPointCollection;
  61. fMpdEmcPointCollection = nullptr;
  62. }
  63. }
  64. void MpdEmcKI::Reset()
  65. {
  66. fSuperParents.clear();
  67. fMpdEmcPointCollection->Clear();
  68. fNhits = 0;
  69. fCurrentTrackID = -1;
  70. fCurrentCellID = -1;
  71. fCurentSuperParent = -1;
  72. fCurrentHit = nullptr;
  73. }
  74. /*
  75. void MpdEmcKI::Initialize()
  76. {
  77. //Extract inner and outer radii of ECAL from current geometry
  78. TGeoVolume * v =gGeoManager->GetVolume("emc1Chamber1_0");
  79. TGeoTube *sh = dynamic_cast<TGeoTube*>(v->GetShape()) ;
  80. // double towerHalfSizeZ = sh->GetDz() ;
  81. fEcalRmin = sh->GetRmin() ;
  82. fEcalRmax = sh->GetRmax() ;
  83. }
  84. */
  85. Bool_t MpdEmcKI::ProcessHits(FairVolume* vol)
  86. {
  87. // This method is called from the MC stepping
  88. // 1. Remember all particles first entered ECAL (active medium declared below: _sc* and _pl*)
  89. // 2. If this is _pl, do not create hit, just remember entered particle
  90. // 3. Collect all energy depositions in sc* by all secondaries from particle first entered ECAL
  91. // Check if this is first entered ECAL particle ("SuperParent")
  92. MpdStack* stack = static_cast<MpdStack*>(gMC->GetStack());
  93. const Int_t partID = stack->GetCurrentTrackNumber();
  94. Int_t superParent = -1;
  95. Bool_t isNewPartile = false; // Create Hit even if zero energy deposition
  96. if (partID != fCurrentTrackID) { // not same track as before, check: same SuperParent or new one?
  97. auto itTr = fSuperParents.find(partID);
  98. if (itTr == fSuperParents.end()) {
  99. // Search parent
  100. Int_t parentID = stack->GetCurrentTrack()->GetMother(0);
  101. itTr = fSuperParents.find(parentID);
  102. if (itTr == fSuperParents.end()) { // Neither track or its parent found: new SuperParent
  103. // Go back to parents untill found particle entered ECAL
  104. TParticle* part = stack->GetCurrentTrack();
  105. Double_t rvert = part->R(); // Transverse radius despite the name!
  106. TParticle* mctr = part;
  107. Int_t idmother = partID;
  108. while (rvert > fEcalRmin && rvert < fEcalRmax) {
  109. // Born inside EMC - find ancestor
  110. Int_t tmp = mctr->GetMother(0);
  111. if (tmp < 0) {
  112. break;
  113. }
  114. idmother = tmp;
  115. mctr = stack->GetParticle(idmother);
  116. rvert = mctr->R(); // transverse radius despite the name!!!
  117. }
  118. fSuperParents[partID] = idmother;
  119. if (partID != idmother) { // Add intermediate tracks if necessary
  120. Int_t tmp = part->GetMother(0);
  121. while (tmp >= 0 && tmp != idmother) {
  122. itTr = fSuperParents.find(tmp);
  123. if (itTr == fSuperParents.end()) {
  124. fSuperParents[tmp] = idmother;
  125. } else {
  126. break;
  127. }
  128. mctr = stack->GetParticle(tmp);
  129. tmp = mctr->GetMother(0);
  130. }
  131. }
  132. superParent = idmother;
  133. isNewPartile = true;
  134. // Mark current track to keep in stack
  135. stack->AddPoint(kECAL, superParent);
  136. } else { // parent found, this track - not
  137. superParent = itTr->second;
  138. fSuperParents[partID] = superParent;
  139. fCurrentTrackID = partID;
  140. }
  141. } else {
  142. superParent = itTr->second;
  143. fCurrentTrackID = partID;
  144. }
  145. } else { // Same track as before
  146. superParent = fCurentSuperParent;
  147. }
  148. Double_t lostenergy = gMC->Edep();
  149. if (strstr(vol->GetName(), "cl_pl") != 0) { // Do not count energy deposited in cover
  150. lostenergy = 0.;
  151. }
  152. if (lostenergy < DBL_EPSILON && !isNewPartile) {
  153. fCurentSuperParent = superParent; // If we switched to another track/superparent, remember it
  154. return false; // do not create hits with zero energy deposition
  155. }
  156. if (!fGeom) {
  157. fGeom = MpdEmcGeoUtils::GetInstance();
  158. }
  159. //printf("GEANT path: <%s>\n",gMC->CurrentVolPath()) ;
  160. Int_t chamberH, chamber, sector, crate, module, boxA, boxB, dummyA, dummyB ;
  161. char dummyC1,dummyC2;
  162. /* gMC->CurrentVolOffID(5,
  163. chamberH); // 5: number of geom. levels between emc1_cl_sc* and sector: get the sector number ;
  164. Int_t chamber;
  165. gMC->CurrentVolOffID(4, chamber); // 4: number of geom. levels between emc1_cl_sc* and sector: get the sector number ;
  166. Int_t sector;
  167. gMC->CurrentVolOffID(3, sector); // 3: number of geom. levels between emc1_cl_sc* and sector: get the sector number ;
  168. Int_t crate;
  169. gMC->CurrentVolOffID(2, crate); // 2: Crate in sector: number of geom levels between emc1_cl_sc get the crate number
  170. Int_t box;
  171. gMC->CurrentVolOffID(1, box); // 1: Tower in crate: number of geom levels between emc1_cl_sc and box
  172. */
  173. //Example path: </cave_1/emcChamber_0/emcChH_1/emcSector_20/emcCrate_1/emcModule6_0/emc_box50_99/emc_cl_sc81_161>
  174. // /cave_1/emcChamber_0/emcChH_1/emcSector_0/emcCrate_5/emcModule6_0/emc_box53_105/emc_cl_pl1_0
  175. int nRead= sscanf(gMC->CurrentVolPath(),"/cave_1/emcChamber_%d/emcChH_%d/emcSector_%d/emcCrate_%d/emcModule%d_0/emc_box%d_%d/emc_cl_%c%c%d_%d",
  176. &chamberH,&chamber,&sector,&crate,&module,&boxA,&boxB, &dummyC1,&dummyC2,&dummyA, &dummyB) ;
  177. if(nRead!=11){
  178. LOG(FATAL)<<"Can not parse Geant path:" << gMC->CurrentVolPath() << " nRead = " << nRead<<endl ;
  179. return false ;
  180. }
  181. Int_t detID = fGeom->GeantToDetId(chamberH, chamber, sector, crate, module, boxA, boxB);
  182. if (superParent == fCurentSuperParent && detID == fCurrentCellID && fCurrentHit) {
  183. // continue with current hit
  184. fCurrentHit->AddEnergyLoss(lostenergy); // TODO implement light attenuation vs. scinitllator number
  185. return true;
  186. }
  187. // try to find existing Hit and add energy to it
  188. if (!isNewPartile) {
  189. for (Int_t itr = fNhits - 1; itr >= 0; itr--) {
  190. MpdEmcPointKI* h = static_cast<MpdEmcPointKI*>(fMpdEmcPointCollection->At(itr));
  191. if (h->GetTrackID() != superParent) // switched to another SuperParent, do not search further
  192. break;
  193. if (h->GetDetectorID() == detID) { // found correct hit
  194. h->AddEnergyLoss(lostenergy);
  195. fCurentSuperParent = superParent;
  196. fCurrentTrackID = partID;
  197. fCurrentCellID = detID;
  198. fCurrentHit = h;
  199. return true;
  200. }
  201. }
  202. }
  203. // Create new Hit
  204. Float_t posX = 0., posY = 0., posZ = 0., momX = 0, momY = 0., momZ = 0., energy = 0.;
  205. gMC->TrackPosition(posX, posY, posZ);
  206. gMC->TrackMomentum(momX, momY, momZ, energy);
  207. Double_t time = gMC->TrackTime() * 1.e+9; // time in ns?? To be consistent with EMCAL
  208. Float_t length = gMC->TrackLength();
  209. fCurrentHit =
  210. AddHit(superParent, detID, TVector3(posX, posY, posZ), TVector3(momX, momY, momZ), time, length, lostenergy);
  211. fCurentSuperParent = superParent;
  212. fCurrentTrackID = partID;
  213. fCurrentCellID = detID;
  214. return true;
  215. }
  216. void MpdEmcKI::EndOfEvent()
  217. {
  218. Print();
  219. Reset();
  220. }
  221. void MpdEmcKI::Register()
  222. {
  223. /** This will create a branch in the output tree called
  224. MpdEmcPointKI, setting the last parameter to kFALSE means:
  225. this collection will not be written to the file, it will exist
  226. only during the simulation.
  227. */
  228. FairRootManager::Instance()->Register("EmcPoint", "MpdEmc", fMpdEmcPointCollection, kTRUE);
  229. }
  230. TClonesArray* MpdEmcKI::GetCollection(Int_t iColl) const
  231. {
  232. if (iColl == 0)
  233. return fMpdEmcPointCollection;
  234. else
  235. return NULL;
  236. }
  237. void MpdEmcKI::FinishEvent()
  238. {
  239. // Sort Hits
  240. // Add duplicates if any and remove them
  241. if (fNhits == 0) {
  242. return;
  243. }
  244. fMpdEmcPointCollection->Sort();
  245. for (Int_t i = 0; i < fNhits; i++) {
  246. MpdEmcPointKI* firstHit = static_cast<MpdEmcPointKI*>(fMpdEmcPointCollection->At(i));
  247. if (!firstHit) { // hit already removed
  248. continue;
  249. }
  250. for (Int_t j = i + 1; j < fNhits; j++) {
  251. MpdEmcPointKI* secondHit = static_cast<MpdEmcPointKI*>(fMpdEmcPointCollection->At(j));
  252. if ((*firstHit) == (*secondHit)) { // same volume and same parent: add and remove second
  253. (*firstHit) += (*secondHit);
  254. fMpdEmcPointCollection->RemoveAt(j);
  255. } else { // hits are sorted, no need to scan further
  256. break;
  257. }
  258. }
  259. // Clean, do not store empty points
  260. if (firstHit->GetEnergyLoss() == 0) {
  261. fMpdEmcPointCollection->RemoveAt(i);
  262. }
  263. }
  264. // Remove empty slots and shink array
  265. fMpdEmcPointCollection->Compress();
  266. fNhits = fMpdEmcPointCollection->GetEntriesFast();
  267. fMpdEmcPointCollection->Expand(fNhits);
  268. }
  269. void MpdEmcKI::Print() const
  270. {
  271. Int_t nHits = fMpdEmcPointCollection->GetEntriesFast();
  272. cout << "-I- MpdEmcKI: " << nHits << " points registered in this event." << endl;
  273. if (fVerboseLevel > 1)
  274. for (Int_t i = 0; i < nHits; i++)
  275. (*fMpdEmcPointCollection)[i]->Print();
  276. }
  277. void MpdEmcKI::ConstructGeometry()
  278. {
  279. TString fileName = GetGeometryFileName();
  280. if (fileName.EndsWith(".root")) {
  281. LOG(INFO) << "Constructing EMC geometry from ROOT file " << fileName.Data();
  282. ConstructRootGeometry();
  283. } else if (fileName.EndsWith(".geo")) {
  284. LOG(INFO) << "Constructing EMC geometry from ASCII file " << fileName.Data();
  285. ConstructAsciiGeometry();
  286. } else {
  287. LOG(FATAL) << "Geometry format of EMC file " << fileName.Data() << " not supported.";
  288. }
  289. // Extract inner and outer radii of ECAL from current geometry
  290. //TODO!!!!!!!!!
  291. /* TGeoVolume* v = gGeoManager->GetVolume("emc1Chamber1");
  292. TGeoTube* sh = dynamic_cast<TGeoTube*>(v->GetShape());
  293. // double towerHalfSizeZ = sh->GetDz() ;
  294. fEcalRmin = sh->GetRmin();
  295. fEcalRmax = sh->GetRmax();
  296. */}
  297. void MpdEmcKI::ConstructAsciiGeometry()
  298. {
  299. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  300. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  301. MpdEmcGeo* Geo = new MpdEmcGeo();
  302. Geo->setGeomFile(GetGeometryFileName());
  303. geoFace->addGeoModule(Geo);
  304. Bool_t rc = geoFace->readSet(Geo);
  305. if (rc)
  306. Geo->create(geoLoad->getGeoBuilder());
  307. TList* volList = Geo->getListOfVolumes();
  308. // store geo parameter
  309. FairRun* fRun = FairRun::Instance();
  310. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  311. MpdEmcGeoPar* par = (MpdEmcGeoPar*)(rtdb->getContainer("MpdEmcGeoPar"));
  312. TObjArray* fSensNodes = par->GetGeoSensitiveNodes();
  313. TObjArray* fPassNodes = par->GetGeoPassiveNodes();
  314. TListIter iter(volList);
  315. FairGeoNode* node = NULL;
  316. FairGeoVolume* aVol = NULL;
  317. while ((node = (FairGeoNode*)iter.Next())) {
  318. aVol = dynamic_cast<FairGeoVolume*>(node);
  319. if (node->isSensitive()) {
  320. fSensNodes->AddLast(aVol);
  321. } else {
  322. fPassNodes->AddLast(aVol);
  323. }
  324. }
  325. par->setChanged();
  326. par->setInputVersion(fRun->GetRunId(), 1);
  327. ProcessNodes(volList);
  328. }
  329. // Check sensitivity
  330. Bool_t MpdEmcKI::CheckIfSensitive(std::string name)
  331. {
  332. return (name.find("cl_sc") != std::string::npos || // Active scinillator
  333. name.find("cl_pl") != std::string::npos); // forward plate to remember primary track
  334. }
  335. MpdEmcPointKI* MpdEmcKI::AddHit(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, Double_t time, Double_t length,
  336. Double_t ELoss)
  337. {
  338. return new ((*fMpdEmcPointCollection)[fNhits++]) MpdEmcPointKI(trackID, detID, pos, mom, time, length, ELoss);
  339. }
  340. ClassImp(MpdEmcKI)