MpdBbcHitProducer.cxx 6.2 KB

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