MpdCpcHitProducer.cxx 7.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260
  1. /////////////////////////////////////////////////////////////
  2. //
  3. // MpdCpcHitProducer
  4. //
  5. // Filler of MpdCpcHit
  6. //
  7. ///////////////////////////////////////////////////////////////
  8. #include "TClonesArray.h"
  9. #include "FairRootManager.h"
  10. #include "FairDetector.h"
  11. #include "FairBaseParSet.h"
  12. #include "FairGeoNode.h"
  13. #include "TGeoManager.h"
  14. #include "TGeoVolume.h"
  15. #include "TGeoNode.h"
  16. #include "TGeoMatrix.h"
  17. #include "TVector3.h"
  18. #include "FairRun.h"
  19. #include "FairRuntimeDb.h"
  20. #include "TMath.h"
  21. #include "FairRunAna.h"
  22. #include "MpdCpcHitProducer.h"
  23. #include "MpdCpcGeoPar.h"
  24. #include "MpdCpcHit.h"
  25. #include "MpdCpcPoint.h"
  26. // ----- Default constructor -------------------------------------------
  27. MpdCpcHitProducer::MpdCpcHitProducer(const char* fileGeo) :
  28. FairTask("Ideal CPC hit Producer") {
  29. fFileGeo=fileGeo;
  30. eneThr = 0.001; // Energy threshold for CPC
  31. }
  32. // ----- Destructor ----------------------------------------------------
  33. MpdCpcHitProducer::~MpdCpcHitProducer() { }
  34. // -------------------------------------------------------------------------
  35. // ----- Public method Init --------------------------------------------
  36. InitStatus MpdCpcHitProducer::Init() {
  37. cout << "******************* INITIALIZATION *********************" << endl;
  38. //FairDetector::Initialize();
  39. //FairRun* sim = FairRun::Instance();
  40. //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  41. // Get RootManager
  42. FairRootManager* ioman = FairRootManager::Instance();
  43. if ( ! ioman ) {
  44. cout << "-E- MpdCpcHitProducer::Init: "
  45. << "RootManager not instantiated!" << endl;
  46. return kFATAL;
  47. }
  48. // Get input array
  49. fPointArray = (TClonesArray*) ioman->GetObject("CPCPoint");
  50. if ( ! fPointArray ) {
  51. cout << "-W- MpdCpcHitProducer::Init: "
  52. << "No CpcPoint array!" << endl;
  53. return kERROR;
  54. }
  55. // Create and register output array
  56. fDigiArray = new TClonesArray("MpdCpcHit");
  57. ioman->Register("CpcHit","Cpc",fDigiArray,kTRUE);
  58. CreateStructure();
  59. hlist = new TList();
  60. MakeHists();
  61. cout << "-I- MpdCpcHitProducer: Intialization successfull" << endl;
  62. return kSUCCESS;
  63. }
  64. //__________________________________________________________________
  65. void MpdCpcHitProducer::Finish() {
  66. //---
  67. cout << "-I- MpdCpcHitProducer: FinishTask" << endl;
  68. }
  69. //__________________________________________________________________
  70. void MpdCpcHitProducer::FinishTask() {
  71. //---
  72. cout << "-I- MpdCpcHitProducer: Finish" << endl;
  73. // save histograms of corresponding list
  74. if (hlist!=0) {
  75. TObject *obj;
  76. TIter next(hlist);
  77. while((obj = (TObject*)next())) obj->Write();
  78. }
  79. }
  80. //__________________________________________________________________
  81. void MpdCpcHitProducer::MakeHists() {
  82. //---
  83. fZ = new TH1F("z","", 500, 150., 250.);
  84. hlist->Add(fZ);
  85. //Float_t xbins[5] = {0., 15., 20., 30., 50.0};
  86. fR = new TH1F("r","", 100, 0., 55.);
  87. //fR = new TH1F("r","", 4, xbins);
  88. hlist->Add(fR);
  89. fXY = new TH2F("xy","", 100, -50., 50., 100, -50., 50.);
  90. hlist->Add(fXY);
  91. fRphi = new TH2F("rphi","", 100, 0., TMath::TwoPi(), 100, 0., 50.);
  92. hlist->Add(fRphi);
  93. }
  94. // ----- Private method SetParContainers -------------------------------
  95. void MpdCpcHitProducer::SetParContainers() {
  96. // Get run and runtime database
  97. FairRunAna* run = FairRunAna::Instance();
  98. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  99. FairRuntimeDb* db = run->GetRuntimeDb();
  100. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  101. // Get Cpc geometry parameter container
  102. db->getContainer("MpdCpcGeoPar");
  103. }
  104. // ----- Public method Exec --------------------------------------------
  105. void MpdCpcHitProducer::Exec(Option_t* opt) {
  106. //cout << " DIGI EXECUTION *********************" << endl;
  107. // Reset output array
  108. if ( ! fDigiArray ) Fatal("Exec", "No DigiArray");
  109. ///////////////////////////////////////
  110. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  111. //rtdb->printParamContext();
  112. MpdCpcGeoPar *geoPar = (MpdCpcGeoPar*) rtdb->getContainer("MpdCpcGeoPar");
  113. //TString volName = "cpc01l";
  114. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  115. Double_t rMin, rMax;
  116. for (Int_t i=0; i<sensNodes->GetEntriesFast(); i++){
  117. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->At(i));
  118. TArrayD* params = sensVol->getParameters();
  119. rMin = params->At(0);
  120. rMax = params->At(1);
  121. cout << " *** CPC sensitive volume " <<(i+1) <<": " << sensVol->GetName() << " " << rMin << " " << rMax << endl;
  122. }
  123. /////////////////////////////////////
  124. fDigiArray->Clear();
  125. // Declare some variables
  126. MpdCpcPoint* point = NULL;
  127. map<Int_t, Float_t> fTrackEnergy;
  128. fTrackEnergy.clear();
  129. map<Int_t, Float_t>::const_iterator p;
  130. int oldtrackID = -6;
  131. // Loop over CpcPoints
  132. Int_t nPoints = fPointArray->GetEntriesFast();
  133. for (Int_t iPoint=0; iPoint < nPoints; iPoint++) {
  134. point = (MpdCpcPoint*) fPointArray->At(iPoint);
  135. // fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
  136. int trackID = point->GetTrackID();
  137. if (oldtrackID != trackID) {
  138. fZ->Fill(point->GetZ());
  139. fXY->Fill(point->GetX(), point->GetY());
  140. TVector2 v = TVector2(point->GetX(), point->GetY());
  141. fR->Fill(v.Mod());
  142. fRphi->Fill(v.Phi(), v.Mod());
  143. }
  144. }
  145. #if 0
  146. // Loop to register CpcHit
  147. for(p=fTrackEnergy.begin(); p!=fTrackEnergy.end(); ++p) {
  148. if ((*p).second>eneThr)
  149. AddHit(1, (*p).first, (*p).second);
  150. }
  151. #endif
  152. }
  153. // -------------------------------------------------------------------------
  154. // ----- Public method Create Structure --------------------------------
  155. void MpdCpcHitProducer::CreateStructure() {
  156. /*
  157. TString work = getenv("VMCWORKDIR");
  158. work = work + "/geometry/" + fFileGeo;
  159. cout << "-I- <MpdCpcHitProducer::CreateStructure> Cpc geometry loaded from: "
  160. << work << endl;
  161. Int_t detId = -1;
  162. MpdCpcReader read(work);
  163. for(Int_t module=1; module<=read.GetMaxModules(); module++)
  164. for(Int_t row=1; row<=read.GetMaxRows(module); row++)
  165. for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
  166. DataG4 data = read.GetData(module,row,crystal);
  167. for(Int_t copy=1; copy<=20; copy++) {
  168. detId = module*100000000 + row*1000000 + copy*10000 + crystal;
  169. emcX[detId] = data.posX; emcY[detId] = data.posY; emcZ[detId] = data.posZ;
  170. emcTheta[detId] = data.theta; emcTau[detId] = data.tau;
  171. if (module==3)
  172. emcPhi[detId] = fmod(data.phi+90.*(copy-1),360);
  173. else
  174. emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
  175. }
  176. }
  177. */
  178. }
  179. // ----- Private method AddDigi --------------------------------------------
  180. MpdCpcHit* MpdCpcHitProducer::AddHit(Int_t trackID,Int_t detID, Float_t energy){
  181. // It fills the MpdCpcHit category
  182. // cout << "MpdCpcHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
  183. TClonesArray& clref = *fDigiArray;
  184. Int_t size = clref.GetEntriesFast();
  185. return new(clref[size]) MpdCpcHit(); // FIXME: real hit info needed here
  186. }
  187. // ----
  188. ClassImp(MpdCpcHitProducer)