MpdEmc.cxx 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229
  1. #include "MpdEmc.h"
  2. #include "MpdEmcPoint.h"
  3. #include "MpdEmcGeo.h"
  4. #include "MpdEmcGeoPar.h"
  5. #include "FairVolume.h"
  6. #include "FairGeoVolume.h"
  7. #include "FairGeoNode.h"
  8. #include "FairRootManager.h"
  9. #include "FairGeoLoader.h"
  10. #include "FairGeoInterface.h"
  11. #include "FairRun.h"
  12. #include "FairRuntimeDb.h"
  13. #include "MpdDetectorList.h"
  14. #include "MpdStack.h"
  15. #include "TClonesArray.h"
  16. #include "TVirtualMC.h"
  17. #include <iostream>
  18. using std::cout;
  19. using std::endl;
  20. MpdEmc::MpdEmc() :
  21. FairDetector("EMC", kTRUE, kECAL) {
  22. /** create your collection for data points */
  23. fMpdEmcPointCollection = new TClonesArray("MpdEmcPoint");
  24. }
  25. MpdEmc::MpdEmc(const char* name, Bool_t active)
  26. : FairDetector(name, active, kECAL) {
  27. fMpdEmcPointCollection = new TClonesArray("MpdEmcPoint");
  28. fVerboseLevel=1;
  29. }
  30. MpdEmc::~MpdEmc() {
  31. if (fMpdEmcPointCollection) {
  32. fMpdEmcPointCollection->Delete();
  33. delete fMpdEmcPointCollection;
  34. }
  35. }
  36. /*
  37. void MpdEmc::Initialize()
  38. {
  39. FairDetector::Initialize();
  40. FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb();
  41. //MpdEmcGeoPar* par=(MpdEmcGeoPar*)(rtdb->getContainer("MpdEmcGeoPar"));
  42. }
  43. */
  44. Bool_t MpdEmc::ProcessHits(FairVolume* vol)
  45. {
  46. /** This method is called from the MC stepping */
  47. //#define EDEBUG
  48. #ifdef EDEBUG
  49. static int my_counter=0;
  50. if (my_counter<10) {
  51. cout << "-I- MpdEmc::ProcessHits IsTrackEntering=" << ( gMC->IsTrackEntering() ) << endl;
  52. my_counter++;
  53. }
  54. #endif
  55. #undef EDEBUG
  56. //Set parameters at entrance of volume. Reset ELoss.
  57. if ( gMC->IsTrackEntering() ) {
  58. fELoss = 0.;
  59. fTime = gMC->TrackTime() * 1.0e09;
  60. fLength = gMC->TrackLength();
  61. gMC->TrackPosition(fPos);
  62. gMC->TrackMomentum(fMom);
  63. }
  64. // Sum energy loss for all steps in the active volume
  65. fELoss += gMC->Edep();
  66. // Create MpdEmcPoint at exit of active volume
  67. if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared() )
  68. {
  69. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  70. fVolumeID = vol->getMCid();
  71. if (fELoss == 0. ) return kFALSE;
  72. AddHit(fTrackID, fVolumeID, TVector3(fPos.X(), fPos.Y(), fPos.Z()),
  73. TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, fLength,
  74. fELoss);
  75. // Increment number of tutorial det points in TParticle
  76. // printf("Name x, y, z %s %f %f %f \n", vol->GetName(), fPos.X(), fPos.Y(), fPos.Z());
  77. MpdStack* stack = (MpdStack*) gMC->GetStack();
  78. stack->AddPoint(kECAL);
  79. }
  80. return kTRUE;
  81. }
  82. void MpdEmc::EndOfEvent() {
  83. Print();
  84. fMpdEmcPointCollection->Clear();
  85. }
  86. void MpdEmc::Register() {
  87. /** This will create a branch in the output tree called
  88. MpdEmcPoint, setting the last parameter to kFALSE means:
  89. this collection will not be written to the file, it will exist
  90. only during the simulation.
  91. */
  92. FairRootManager::Instance()->Register("EmcPoint", "MpdEmc",
  93. fMpdEmcPointCollection, kTRUE);
  94. }
  95. TClonesArray* MpdEmc::GetCollection(Int_t iColl) const {
  96. if (iColl == 0) return fMpdEmcPointCollection;
  97. else return nullptr;
  98. }
  99. void MpdEmc::Reset() {
  100. fMpdEmcPointCollection->Clear();
  101. }
  102. void MpdEmc::Print() const {
  103. Int_t nHits = fMpdEmcPointCollection->GetEntriesFast();
  104. cout << "-I- MpdEmc: " << nHits << " points registered in this event."
  105. << endl;
  106. if (fVerboseLevel>1)
  107. for (Int_t i=0; i<nHits; i++) (*fMpdEmcPointCollection)[i]->Print();
  108. }
  109. void MpdEmc::ConstructGeometry() {
  110. TString fileName = GetGeometryFileName();
  111. if ( fileName.EndsWith(".root") ) {
  112. LOG(info)<<"Constructing EMC geometry from ROOT file "<<fileName.Data();
  113. ConstructRootGeometry();
  114. }
  115. else if ( fileName.EndsWith(".geo") ) {
  116. LOG(info)<<"Constructing EMC geometry from ASCII file "<<fileName.Data();
  117. ConstructAsciiGeometry();
  118. }
  119. else {
  120. LOG(fatal)<<"Geometry format of EMC file %s not supported."<<fileName.Data();
  121. }
  122. }
  123. void MpdEmc::ConstructAsciiGeometry() {
  124. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  125. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  126. MpdEmcGeo* Geo = new MpdEmcGeo();
  127. Geo->setGeomFile(GetGeometryFileName());
  128. geoFace->addGeoModule(Geo);
  129. Bool_t rc = geoFace->readSet(Geo);
  130. if (rc) Geo->create(geoLoad->getGeoBuilder());
  131. TList* volList = Geo->getListOfVolumes();
  132. // store geo parameter
  133. FairRun *fRun = FairRun::Instance();
  134. FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb();
  135. MpdEmcGeoPar* par=(MpdEmcGeoPar*)(rtdb->getContainer("MpdEmcGeoPar"));
  136. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  137. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  138. TListIter iter(volList);
  139. FairGeoNode* node = nullptr;
  140. FairGeoVolume *aVol=nullptr;
  141. while( (node = (FairGeoNode*)iter.Next()) ) {
  142. aVol = dynamic_cast<FairGeoVolume*> ( node );
  143. if ( node->isSensitive() ) {
  144. fSensNodes->AddLast( aVol );
  145. }else{
  146. fPassNodes->AddLast( aVol );
  147. }
  148. }
  149. par->setChanged();
  150. par->setInputVersion(fRun->GetRunId(),1);
  151. ProcessNodes ( volList );
  152. }
  153. // Check sensitivity
  154. Bool_t MpdEmc::CheckIfSensitive(std::string name) {
  155. TString tsname = name;
  156. if (tsname.Contains("cl_sc")) {
  157. return kTRUE;
  158. }
  159. return kFALSE;
  160. }
  161. MpdEmcPoint* MpdEmc::AddHit(Int_t trackID, Int_t detID,
  162. TVector3 pos, TVector3 mom,
  163. Double_t time, Double_t length,
  164. Double_t ELoss) {
  165. TClonesArray& clref = *fMpdEmcPointCollection;
  166. Int_t size = clref.GetEntriesFast();
  167. return new(clref[size]) MpdEmcPoint(trackID, detID, pos, mom,
  168. time, length, ELoss);
  169. }
  170. ClassImp(MpdEmc)