BmdDetector.cxx 16 KB

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