TpcDetector.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class TpcDetector
  7. // see TpcDetector.hh for details
  8. //
  9. // Environment:
  10. // Software developed for the PANDA Detector at FAIR.
  11. //
  12. // Author List:
  13. // Sebastian Neubert TUM (original author)
  14. //
  15. //
  16. //-----------------------------------------------------------
  17. // Panda Headers ----------------------
  18. // This Class' Header ------------------
  19. #include "TpcDetector.h"
  20. #include "MpdTpcEDepParams.h"
  21. // C/C++ Headers ----------------------
  22. #include <iostream>
  23. using namespace std;
  24. // Collaborating Class Headers --------
  25. #include "TClonesArray.h"
  26. #include "FairRootManager.h"
  27. #include "TpcPoint.h"
  28. #include "TVirtualMC.h"
  29. #include "TLorentzVector.h"
  30. #include "TpcGeo.h"
  31. #include "FairGeoLoader.h"
  32. #include "FairGeoInterface.h"
  33. #include "FairGeoNode.h"
  34. #include "FairGeoRootBuilder.h"
  35. #include "TList.h"
  36. #include "FairRun.h"
  37. #include "FairRuntimeDb.h"
  38. #include "TpcGeoPar.h"
  39. #include "TObjArray.h"
  40. #include <TSystem.h>
  41. #include "FairGeoNode.h"
  42. #include "FairGeoVolume.h"
  43. #include "FairVolume.h"
  44. #include "MpdStack.h"
  45. #include "TFile.h"
  46. #include "TGeoManager.h"
  47. #include "TGDMLParse.h"
  48. #include "FairGeoMedia.h"
  49. // Class Member definitions -----------
  50. TpcDetector::TpcDetector(const char * Name, Bool_t Active)
  51. : FairDetector(Name, Active)
  52. {
  53. fTpcPointCollection= new TClonesArray("TpcPoint");
  54. fVerboseLevel = 1;
  55. }
  56. TpcDetector::TpcDetector()
  57. {
  58. fTpcPointCollection= new TClonesArray("TpcPoint");
  59. fVerboseLevel = 1;
  60. }
  61. TpcDetector::~TpcDetector()
  62. {
  63. if (fTpcPointCollection) {
  64. fTpcPointCollection->Delete();
  65. delete fTpcPointCollection;
  66. }
  67. }
  68. void TpcDetector::EndOfEvent()
  69. {
  70. if(fVerboseLevel) Print();
  71. fTpcPointCollection->Delete();
  72. //AZ - memory monitor
  73. ProcInfo_t proc;
  74. gSystem->GetProcInfo(&proc);
  75. cout << " User CPU time: " << proc.fCpuUser << " Memory: resident " << proc.fMemResident << ", virtual " << proc.fMemVirtual << endl;
  76. }
  77. void TpcDetector::Register() {
  78. /** This will create a branch in the output tree called TpcDetectorPoint, setting the last parameter to kFALSE means:
  79. this collection will not be written to the file, it will exist only during the simulation. */
  80. FairRootManager::Instance()->Register("TpcPoint", "Tpc", fTpcPointCollection, kTRUE);
  81. }
  82. //---------------------------------------------------------------------------
  83. Bool_t TpcDetector::ProcessHits( FairVolume *v)
  84. {
  85. // create Hit for every MC step
  86. Double_t time = gMC->TrackTime() * 1.0e09;
  87. Double_t length = gMC->TrackLength();
  88. TLorentzVector pos;
  89. gMC->TrackPosition(pos);
  90. // pos.Print();
  91. TLorentzVector mom;
  92. gMC->TrackMomentum(mom);
  93. // mom.Print();
  94. Double_t eLoss = gMC->Edep();
  95. Int_t trackID = gMC->GetStack()->GetCurrentTrackNumber();
  96. Int_t volumeID = v->getMCid();
  97. Float_t pX=0., pY=0., pZ=0., eTot=0.;
  98. gMC->TrackMomentum( pX, pY, pZ, eTot);
  99. if (eLoss > 0 && eTot>0.001) { //AZ - hard-coded 1 MeV cut
  100. //if (eLoss > 0) { //AZ
  101. //AZ - dE/dx simulation from the new prescription
  102. //cout << " xxxxxxxxxxxxx " << gMC->TrackMass() << endl;
  103. Double_t betgam = mom.Vect().Mag() / gMC->TrackMass();
  104. //cout << betgam << " " << mom.Beta() * mom.Gamma() << endl;
  105. if (betgam > 0.01) {
  106. betgam = TMath::Log10(betgam);
  107. MpdTpcEDepParams *deParams = MpdTpcEDepParams::Instance();
  108. //IR 05-APR-2021 if (betgam >= deParams->GetMinLimit() && betgam < deParams->GetMaxLimit()) {
  109. betgam = TMath::Min( betgam, deParams->GetMaxLimit());
  110. if (betgam >= deParams->GetMinLimit()) {
  111. // New value
  112. eLoss = deParams->GetEloss (betgam, gMC->TrackCharge(), gMC->TrackStep());
  113. //cout << " xxxxxx " << gMC->Edep() << " " << eLoss << endl;
  114. //AZ
  115. TpcPoint* p = AddHit(trackID, volumeID, pos.Vect(), mom.Vect(), time, length, eLoss);
  116. p->SetStep(gMC->TrackStep());
  117. ((MpdStack*)gMC->GetStack())->AddPoint(kTPC);
  118. }
  119. }
  120. }
  121. // p->Print("");
  122. return kTRUE;
  123. }
  124. //---------------------------------------------------------------------------
  125. TClonesArray* TpcDetector::GetCollection(Int_t iColl) const {
  126. if (iColl == 0) return fTpcPointCollection;
  127. else return nullptr;
  128. }
  129. TpcPoint* TpcDetector::AddHit(Int_t trackID, Int_t detID, TVector3 pos,
  130. TVector3 mom, Double_t time, Double_t length,
  131. Double_t eLoss) {
  132. TClonesArray& clref = *fTpcPointCollection;
  133. Int_t size = clref.GetEntriesFast();
  134. return new(clref[size]) TpcPoint(trackID, detID, pos, mom,
  135. time, length, eLoss);
  136. }
  137. void TpcDetector::Reset() {
  138. fTpcPointCollection->Delete();
  139. //ResetParameters();
  140. }
  141. void TpcDetector::Print() const
  142. {
  143. Int_t nHits = fTpcPointCollection->GetEntriesFast();
  144. cout<<"-I- MpdTpc: " << nHits << " points registered in this event." << endl;
  145. if(fVerboseLevel > 1)
  146. for(Int_t i=0; i<nHits; i++) (*fTpcPointCollection)[i]->Print();
  147. }
  148. void TpcDetector::ConstructGeometry() {
  149. TString fileName = GetGeometryFileName();
  150. if ( fileName.EndsWith(".root") ) {
  151. LOG(INFO) << "Constructing TPC geometry from ROOT file " << fileName.Data();
  152. ConstructRootGeometry();
  153. }
  154. else if ( fileName.EndsWith(".geo") ) {
  155. LOG(INFO) << "Constructing TPC geometry from ASCII file " << fileName.Data();
  156. ConstructAsciiGeometry();
  157. }
  158. else if ( fileName.EndsWith(".gdml") )
  159. {
  160. LOG(INFO) << "Constructing TPC geometry from GDML file " << fileName.Data();
  161. //ConstructGDMLGeometry();
  162. }
  163. else
  164. {
  165. LOG(FATAL) << "Geometry format of TPC file " << fileName.Data() << " not supported.";
  166. }
  167. }
  168. // ----- ConstructAsciiGeometry -------------------------------------------
  169. void TpcDetector::ConstructAsciiGeometry() {
  170. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  171. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  172. TpcGeo* TPCGeo = new TpcGeo();
  173. TPCGeo->setGeomFile(GetGeometryFileName());
  174. geoFace->addGeoModule(TPCGeo);
  175. Bool_t rc = geoFace->readSet(TPCGeo);
  176. if (rc) TPCGeo->create(geoLoad->getGeoBuilder());
  177. TList* volList = TPCGeo->getListOfVolumes();
  178. // store geo parameter
  179. FairRun *fRun = FairRun::Instance();
  180. FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb();
  181. TpcGeoPar* par=(TpcGeoPar*)(rtdb->getContainer("TpcGeoPar"));
  182. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  183. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  184. TListIter iter(volList);
  185. FairGeoNode* node = nullptr;
  186. FairGeoVolume *aVol=nullptr;
  187. while( (node = (FairGeoNode*)iter.Next()) ) {
  188. aVol = dynamic_cast<FairGeoVolume*> ( node );
  189. if ( node->isSensitive() ) {
  190. fSensNodes->AddLast( aVol );
  191. }else{
  192. fPassNodes->AddLast( aVol );
  193. }
  194. }
  195. par->setChanged();
  196. par->setInputVersion(fRun->GetRunId(),1);
  197. ProcessNodes( volList );
  198. }
  199. // ----------------------------------------------------------------------------
  200. // ----- ConstructGDMLGeometry -------------------------------------------
  201. void TpcDetector::ConstructGDMLGeometry()
  202. {
  203. TFile *old = gFile;
  204. TGDMLParse parser;
  205. TGeoVolume* gdmlTop;
  206. // Before importing GDML
  207. Int_t maxInd = gGeoManager->GetListOfMedia()->GetEntries() - 1;
  208. gdmlTop = parser.GDMLReadFile(GetGeometryFileName());
  209. // Cheating - reassigning media indices after GDML import (need to fix this in TGDMLParse class!!!)
  210. // for (Int_t i=0; i<gGeoManager->GetListOfMedia()->GetEntries(); i++)
  211. // gGeoManager->GetListOfMedia()->At(i)->Dump();
  212. // After importing GDML
  213. Int_t j = gGeoManager->GetListOfMedia()->GetEntries() - 1;
  214. Int_t curId;
  215. TGeoMedium* m;
  216. do {
  217. m = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(j);
  218. curId = m->GetId();
  219. m->SetId(curId+maxInd);
  220. j--;
  221. } while (curId > 1);
  222. // LOG(DEBUG) << "====================================================================";
  223. // for (Int_t i=0; i<gGeoManager->GetListOfMedia()->GetEntries(); i++)
  224. // gGeoManager->GetListOfMedia()->At(i)->Dump();
  225. Int_t newMaxInd = gGeoManager->GetListOfMedia()->GetEntries() - 1;
  226. gGeoManager->GetTopVolume()->AddNode(gdmlTop, 1, 0);
  227. ExpandNodeForGdml(gGeoManager->GetTopVolume()->GetNode(gGeoManager->GetTopVolume()->GetNdaughters()-1));
  228. for (Int_t k = maxInd+1; k < newMaxInd+1; k++) {
  229. TGeoMedium* medToDel = (TGeoMedium*)(gGeoManager->GetListOfMedia()->At(maxInd+1));
  230. LOG(DEBUG) << " removing media " << medToDel->GetName() << " with id " << medToDel->GetId() << " (k=" << k << ")";
  231. gGeoManager->GetListOfMedia()->Remove(medToDel);
  232. }
  233. gGeoManager->SetAllIndex();
  234. gFile = old;
  235. }
  236. void TpcDetector::ExpandNodeForGdml(TGeoNode* node)
  237. {
  238. LOG(DEBUG) << "----------------------------------------- ExpandNodeForGdml for node " << node->GetName();
  239. TGeoVolume* curVol = node->GetVolume();
  240. LOG(DEBUG) << " volume: " << curVol->GetName();
  241. if (curVol->IsAssembly()) {
  242. LOG(DEBUG) << " skipping volume-assembly";
  243. }
  244. else
  245. {
  246. TGeoMedium* curMed = curVol->GetMedium();
  247. TGeoMaterial* curMat = curVol->GetMaterial();
  248. TGeoMedium* curMedInGeoManager = gGeoManager->GetMedium(curMed->GetName());
  249. TGeoMaterial* curMatOfMedInGeoManager = curMedInGeoManager->GetMaterial();
  250. TGeoMaterial* curMatInGeoManager = gGeoManager->GetMaterial(curMat->GetName());
  251. // Current medium and material assigned to the volume from GDML
  252. LOG(DEBUG2) << " curMed\t\t\t\t" << curMed << "\t" << curMed->GetName() << "\t" << curMed->GetId();
  253. LOG(DEBUG2) << " curMat\t\t\t\t" << curMat << "\t" << curMat->GetName() << "\t" << curMat->GetIndex();
  254. // Medium and material found in the gGeoManager - either the pre-loaded one or one from GDML
  255. LOG(DEBUG2) << " curMedInGeoManager\t\t" << curMedInGeoManager
  256. << "\t" << curMedInGeoManager->GetName() << "\t" << curMedInGeoManager->GetId();
  257. LOG(DEBUG2) << " curMatOfMedInGeoManager\t\t" << curMatOfMedInGeoManager
  258. << "\t" << curMatOfMedInGeoManager->GetName() << "\t" << curMatOfMedInGeoManager->GetIndex();
  259. LOG(DEBUG2) << " curMatInGeoManager\t\t" << curMatInGeoManager
  260. << "\t" << curMatInGeoManager->GetName() << "\t" << curMatInGeoManager->GetIndex();
  261. TString matName = curMat->GetName();
  262. TString medName = curMed->GetName();
  263. if (curMed->GetId() != curMedInGeoManager->GetId()) {
  264. if (fFixedMedia.find(medName) == fFixedMedia.end()) {
  265. LOG(DEBUG) << " Medium needs to be fixed";
  266. fFixedMedia[medName] = curMedInGeoManager;
  267. Int_t ind = curMat->GetIndex();
  268. gGeoManager->RemoveMaterial(ind);
  269. LOG(DEBUG) << " removing material " << curMat->GetName()
  270. << " with index " << ind;
  271. for (Int_t i=ind; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++) {
  272. TGeoMaterial* m = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(i);
  273. m->SetIndex(m->GetIndex()-1);
  274. }
  275. LOG(DEBUG) << " Medium fixed";
  276. }
  277. else
  278. {
  279. LOG(DEBUG) << " Already fixed medium found in the list ";
  280. }
  281. }
  282. else
  283. {
  284. if (fFixedMedia.find(medName) == fFixedMedia.end()) {
  285. LOG(DEBUG) << " There is no correct medium in the memory yet";
  286. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  287. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  288. FairGeoMedia* geoMediaBase = geoFace->getMedia();
  289. FairGeoBuilder* geobuild = geoLoad->getGeoBuilder();
  290. FairGeoMedium* curMedInGeo = geoMediaBase->getMedium(medName);
  291. if (curMedInGeo == 0)
  292. {
  293. LOG(FATAL) << " Media not found in Geo file: " << medName;
  294. //! This should not happen.
  295. //! This means that somebody uses material in GDML that is not in the media.geo file.
  296. //! Most probably this is the sign to the user to check materials' names in the CATIA model.
  297. }
  298. else
  299. {
  300. LOG(DEBUG) << " Found media in Geo file" << medName;
  301. Int_t nmed = geobuild->createMedium(curMedInGeo);
  302. LOG(DEBUG) << " nmed " << nmed;
  303. fFixedMedia[medName] = (TGeoMedium*)gGeoManager->GetListOfMedia()->Last();
  304. gGeoManager->RemoveMaterial(curMatOfMedInGeoManager->GetIndex());
  305. LOG(DEBUG) << " removing material " << curMatOfMedInGeoManager->GetName()
  306. << " with index " << curMatOfMedInGeoManager->GetIndex();
  307. for (Int_t i=curMatOfMedInGeoManager->GetIndex(); i<gGeoManager->GetListOfMaterials()->GetEntries(); i++) {
  308. TGeoMaterial* m = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(i);
  309. m->SetIndex(m->GetIndex()-1);
  310. }
  311. }
  312. if (curMedInGeo->getSensitivityFlag()) {
  313. LOG(DEBUG) << " Adding sensitive " << curVol->GetName();
  314. AddSensitiveVolume(curVol);
  315. }
  316. }
  317. else
  318. {
  319. LOG(DEBUG) << " Already fixed medium found in the list";
  320. LOG(DEBUG) << "!!! Sensitivity: " << fFixedMedia[medName]->GetParam(0);
  321. if (fFixedMedia[medName]->GetParam(0) == 1) {
  322. LOG(DEBUG) << " Adding sensitive " << curVol->GetName();
  323. AddSensitiveVolume(curVol);
  324. }
  325. }
  326. }
  327. curVol->SetMedium(fFixedMedia[medName]);
  328. gGeoManager->SetAllIndex();
  329. // gGeoManager->GetListOfMaterials()->Print();
  330. // gGeoManager->GetListOfMedia()->Print();
  331. }
  332. //! Recursevly go down the tree of nodes
  333. if (curVol->GetNdaughters() != 0)
  334. {
  335. TObjArray* NodeChildList = curVol->GetNodes();
  336. TGeoNode* curNodeChild;
  337. for (Int_t j=0; j<NodeChildList->GetEntriesFast(); j++)
  338. {
  339. curNodeChild = (TGeoNode*)NodeChildList->At(j);
  340. ExpandNodeForGdml(curNodeChild);
  341. }
  342. }
  343. }
  344. //-----------------------------------------------------------------------------
  345. //Check if Sensitive-----------------------------------------------------------
  346. Bool_t TpcDetector::CheckIfSensitive(std::string name) {
  347. TString tsname = name;
  348. if (tsname.Contains("Active") || tsname.Contains("tpc01sv")) {
  349. return kTRUE;
  350. }
  351. return kFALSE;
  352. // if(0 == TString(name).CompareTo("DCH1DetV")) {
  353. // return kTRUE;
  354. // }
  355. // return kFALSE;
  356. }
  357. //---------------------------------------------------------
  358. ClassImp(TpcDetector)