MpdFsaHitProducer_old.cxx 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169
  1. /////////////////////////////////////////////////////////////
  2. //
  3. // MpdFsaHitProducer
  4. //
  5. // Filler of MpdFsaHit
  6. //
  7. //
  8. //
  9. ///////////////////////////////////////////////////////////////
  10. #include "TClonesArray.h"
  11. #include "FairRootManager.h"
  12. #include "FairDetector.h"
  13. #include "MpdFsaHitProducer.h"
  14. #include "MpdFsaHit.h"
  15. #include "MpdFsaPoint.h"
  16. #include "TGeoManager.h"
  17. #include "TGeoVolume.h"
  18. #include "TGeoNode.h"
  19. #include "TGeoMatrix.h"
  20. #include "TVector3.h"
  21. #include "FairRun.h"
  22. #include "FairRuntimeDb.h"
  23. #include "FairMCTrack.h"
  24. // ----- Default constructor -------------------------------------------
  25. MpdFsaHitProducer::MpdFsaHitProducer(const char* fileGeo) :
  26. FairTask("Ideal FSA hit Producer") {
  27. fFileGeo=fileGeo;
  28. eneThr = 0.001; // Energy threshold for FSA
  29. }
  30. // -------------------------------------------------------------------------
  31. // ----- Destructor ----------------------------------------------------
  32. MpdFsaHitProducer::~MpdFsaHitProducer() { }
  33. // -------------------------------------------------------------------------
  34. // ----- Public method Init --------------------------------------------
  35. InitStatus MpdFsaHitProducer::Init() {
  36. cout << " INITIALIZATION *********************" << endl;
  37. //FairDetector::Initialize();
  38. //FairRun* sim = FairRun::Instance();
  39. //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  40. // Get RootManager
  41. FairRootManager* ioman = FairRootManager::Instance();
  42. if ( ! ioman ) {
  43. cout << "-E- MpdFsaHitProducer::Init: "
  44. << "RootManager not instantiated!" << endl;
  45. return kFATAL;
  46. }
  47. // Get input array
  48. fPointArray = (TClonesArray*) ioman->GetObject("FSAPoint");
  49. if ( ! fPointArray ) {
  50. cout << "-W- MpdFsaHitProducer::Init: "
  51. << "No FsaPoint array!" << endl;
  52. return kERROR;
  53. }
  54. fListMCtracks = (TClonesArray *)ioman->GetObject("MCTrack");
  55. // Create and register output array
  56. fDigiArray = new TClonesArray("MpdFsaHit");
  57. ioman->Register("FsaHit","Fsa",fDigiArray,kTRUE);
  58. CreateStructure();
  59. cout << "-I- MpdFsaHitProducer: Intialization successfull" << endl;
  60. return kSUCCESS;
  61. }
  62. // -------------------------------------------------------------------------
  63. // ----- Public method Exec --------------------------------------------
  64. void MpdFsaHitProducer::Exec(Option_t* opt) {
  65. //cout << " DIGI EXECUTION *********************" << endl;
  66. // Reset output array
  67. if ( ! fDigiArray ) Fatal("Exec", "No DigiArray");
  68. fDigiArray->Clear();
  69. // Declare some variables
  70. MpdFsaPoint* point = NULL;
  71. map<Int_t, Float_t> fTrackEnergy;
  72. fTrackEnergy.clear();
  73. map<Int_t, Float_t>::const_iterator p;
  74. // Loop over FsaPoints
  75. Int_t nPoints = fPointArray->GetEntriesFast();
  76. Int_t nTracks = 0;
  77. Int_t iTrack = -1;
  78. for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
  79. point = (MpdFsaPoint*) fPointArray->At(iPoint);
  80. fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
  81. point->Print("");
  82. if(iTrack != point->GetTrackID()){
  83. nTracks++;
  84. iTrack = point->GetTrackID();
  85. }
  86. FairMCTrack *gtrack = (FairMCTrack *)fListMCtracks->At(point->GetTrackID());
  87. cout << "\n PDG " << gtrack->GetPdgCode() << "\n\n";
  88. }
  89. cout << "\n count tracks " << nTracks;
  90. cout << "\n load " << nTracks/9891330.816630875 << " tracks/cm^2";
  91. // Loop over ZdcPoint
  92. // Loop to register ZdcHit
  93. for(p=fTrackEnergy.begin(); p!=fTrackEnergy.end(); ++p) {
  94. if ((*p).second>eneThr)
  95. AddHit(1, (*p).first, (*p).second);
  96. }
  97. }
  98. // -------------------------------------------------------------------------
  99. // ----- Public method Create Structure --------------------------------
  100. void MpdFsaHitProducer::CreateStructure() {
  101. /*
  102. TString work = getenv("VMCWORKDIR");
  103. work = work + "/geometry/" + fFileGeo;
  104. cout << "-I- <MpdFsaHitProducer::CreateStructure> Zdc geometry loaded from: "
  105. << work << endl;
  106. Int_t detId = -1;
  107. MpdFsaReader read(work);
  108. for(Int_t module=1; module<=read.GetMaxModules(); module++)
  109. for(Int_t row=1; row<=read.GetMaxRows(module); row++)
  110. for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
  111. DataG4 data = read.GetData(module,row,crystal);
  112. for(Int_t copy=1; copy<=20; copy++) {
  113. detId = module*100000000 + row*1000000 + copy*10000 + crystal;
  114. emcX[detId] = data.posX; emcY[detId] = data.posY; emcZ[detId] = data.posZ;
  115. emcTheta[detId] = data.theta; emcTau[detId] = data.tau;
  116. if (module==3)
  117. emcPhi[detId] = fmod(data.phi+90.*(copy-1),360);
  118. else
  119. emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
  120. }
  121. }
  122. */
  123. }
  124. // ----- Private method AddDigi --------------------------------------------
  125. MpdFsaHit* MpdFsaHitProducer::AddHit(Int_t trackID,Int_t detID, Float_t energy){
  126. // It fills the MpdFsaHit category
  127. // cout << "MpdFsaHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
  128. TClonesArray& clref = *fDigiArray;
  129. Int_t size = clref.GetEntriesFast();
  130. return new(clref[size]) MpdFsaHit(); // FIXME: real hit info needed here
  131. }
  132. // ----
  133. ClassImp(MpdFsaHitProducer)