MpdEmcHitProducer.cxx 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207
  1. /////////////////////////////////////////////////////////////
  2. //
  3. // MpdEmcHitProducer
  4. //
  5. // Filler of MpdEmcHit
  6. //
  7. ///////////////////////////////////////////////////////////////
  8. #include "TClonesArray.h"
  9. #include "FairRootManager.h"
  10. #include "TMath.h"
  11. #include "FairRunAna.h"
  12. #include "MpdEmcHitProducer.h"
  13. #include "MpdEmcGeoPar.h"
  14. #include "MpdEmcHit.h"
  15. #include "MpdEmcPoint.h"
  16. #include "FairMCPoint.h"
  17. #include "MpdMCTrack.h"
  18. using namespace std;
  19. using namespace TMath;
  20. // ----- Default constructor -------------------------------------------
  21. MpdEmcHitProducer::MpdEmcHitProducer() :
  22. FairTask("Ideal EMC hit Producer") {
  23. }
  24. // ----- Destructor ----------------------------------------------------
  25. MpdEmcHitProducer::~MpdEmcHitProducer() {
  26. }
  27. // -------------------------------------------------------------------------
  28. // ----- Public method Init --------------------------------------------
  29. InitStatus MpdEmcHitProducer::Init() {
  30. cout << "******************* EMC INIT *********************" << endl;
  31. // Get RootManager
  32. FairRootManager* ioman = FairRootManager::Instance();
  33. if (!ioman) {
  34. cout << "-E- MpdEmcHitProducer::Init: "
  35. << "RootManager not instantiated!" << endl;
  36. return kFATAL;
  37. }
  38. fGeoPar = new MpdEmcGeoPar();
  39. // Get input array
  40. fPointArray = (TClonesArray*) ioman->GetObject("MpdEmcPoint");
  41. if (!fPointArray) {
  42. cout << "-W- MpdEmcHitProducer::Init: " << "No EmcPoint array!" << endl;
  43. return kERROR;
  44. }
  45. fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  46. if (!fMcTrackArray) {
  47. cout << "-W- MpdEmcHitProducer::Init: " << "No MCTrack array!" << endl;
  48. return kERROR;
  49. }
  50. // Create and register output array
  51. fDigiArray = new TClonesArray("MpdEmcHit", 100);
  52. ioman->Register("MpdEmcHit", "EMC", fDigiArray, kTRUE);
  53. cout << "-I- MpdEmcHitProducer: Intialization successfull" << endl;
  54. return kSUCCESS;
  55. }
  56. void MpdEmcHitProducer::Finish() {
  57. cout << "-I- MpdEmcHitProducer: Finish" << endl;
  58. }
  59. void MpdEmcHitProducer::Exec(Option_t* opt) {
  60. cout << "MpdEmcHitProducer::Exec started" << endl;
  61. // Reset output Array
  62. if (!fDigiArray) Fatal("MpdEmcHitProducer::Exec)", "No array of digits");
  63. fDigiArray->Delete();
  64. for (UInt_t iPnt = 0; iPnt < (UInt_t) fPointArray->GetEntriesFast(); ++iPnt) {
  65. FairMCPoint* pnt = (FairMCPoint*) fPointArray->At(iPnt);
  66. Int_t trId = pnt->GetTrackID();
  67. MpdMCTrack* tr = (MpdMCTrack*) fMcTrackArray->At(trId);
  68. // if (tr->GetMotherId() != -1) continue;
  69. Int_t pdg = tr->GetPdgCode();
  70. Float_t x = pnt->GetX();
  71. Float_t y = pnt->GetY();
  72. Float_t z = pnt->GetZ();
  73. UInt_t sec = GetSecId(x, y, z);
  74. UInt_t row = GetRowId(z);
  75. UInt_t supMod = GetSupModId(x, y, z, sec);
  76. UInt_t mod = GetModId(x, y, supMod, sec);
  77. Float_t e = pnt->GetEnergyLoss();
  78. MpdEmcHit* hit = SearchHit(sec, supMod, row, mod);
  79. if (hit == NULL) {
  80. hit = new((*fDigiArray)[fDigiArray->GetEntriesFast()]) MpdEmcHit(sec, row, supMod, mod, e);
  81. } else {
  82. hit->IncreaseEnergy(e);
  83. }
  84. hit->SetNumTracks(hit->GetNumTracks() + 1);
  85. hit->SetTrackId(trId);
  86. hit->SetPdg(pdg);
  87. hit->SetZCenter(CalcZCenter(sec, row, mod));
  88. hit->SetPhiCenter(CalcPhiCenter(sec, supMod, mod));
  89. }
  90. for (UInt_t iHit = 0; iHit < (UInt_t) fDigiArray->GetEntriesFast(); ++iHit) {
  91. MpdEmcHit* hit = (MpdEmcHit*) fDigiArray->At(iHit);
  92. if (hit->GetNumTracks() > 1) {
  93. hit->SetPdg(0);
  94. hit->SetTrackId(-1);
  95. }
  96. //hit->Print();
  97. }
  98. }
  99. UInt_t MpdEmcHitProducer::GetSecId(Float_t x, Float_t y, Float_t z) {
  100. Float_t ang = ATan2(y, x);
  101. if (ang < 0) ang += TwoPi();
  102. ang *= RadToDeg();
  103. Int_t nSec = fGeoPar->GetNsec();
  104. if (z > 0.0)
  105. return UInt_t(ang / 360 * nSec);
  106. else
  107. return nSec + UInt_t(ang / 360 * nSec);
  108. }
  109. UInt_t MpdEmcHitProducer::GetSupModId(Float_t x, Float_t y, Float_t z, UInt_t sec) {
  110. Float_t ang = ATan2(y, x);
  111. if (ang < 0) ang += TwoPi();
  112. ang *= RadToDeg();
  113. Float_t nDegreesInOneSector = fGeoPar->GetAngleOfSector();
  114. if (z < 0.0)
  115. sec -= fGeoPar->GetNsec();
  116. Float_t secStartAng = sec * nDegreesInOneSector;
  117. //Float_t secFinishAng = secStartAng + nDegreesInOneSector;
  118. Float_t localAng = ang - secStartAng;
  119. return UInt_t(localAng * fGeoPar->GetNsupMod() / nDegreesInOneSector);
  120. }
  121. UInt_t MpdEmcHitProducer::GetModId(Float_t x, Float_t y, UInt_t supMod, UInt_t sec) {
  122. Float_t ang = ATan2(y, x);
  123. if (ang < 0) ang += TwoPi();
  124. ang *= RadToDeg();
  125. Float_t secStartAng = (sec % fGeoPar->GetNsec()) * fGeoPar->GetAngleOfSector();
  126. Float_t supModStartAng = secStartAng + supMod * fGeoPar->GetAngleOfSuperModule();
  127. Float_t localAng = ang - supModStartAng;
  128. return UInt_t(localAng * fGeoPar->GetNmod() / fGeoPar->GetAngleOfSuperModule());
  129. }
  130. UInt_t MpdEmcHitProducer::GetRowId(Float_t z) {
  131. return UInt_t(Abs(z) / fGeoPar->GetLength() * fGeoPar->GetNrows());
  132. }
  133. Float_t MpdEmcHitProducer::CalcZCenter(UInt_t sec, UInt_t row, UInt_t mod) {
  134. Float_t lengthOfModuleByZ = fGeoPar->GetLengthOfModuleByZ();
  135. Float_t lengthOfSuperModuleByZ = fGeoPar->GetLengthOfSuperModuleByZ();
  136. Float_t halfLengthOfModuleByZ = lengthOfModuleByZ / 2.0;
  137. Float_t leftEdgeOfModuleByZ = row * lengthOfSuperModuleByZ + (mod % fGeoPar->GetNModInSuperModByZ()) * lengthOfModuleByZ;
  138. Float_t z = leftEdgeOfModuleByZ + halfLengthOfModuleByZ;
  139. return (sec < fGeoPar->GetNsec()) ? z : -z;
  140. }
  141. Float_t MpdEmcHitProducer::CalcPhiCenter(UInt_t sec, UInt_t supMod, UInt_t mod) {
  142. Float_t sectorAngle = fGeoPar->GetAngleOfSector();
  143. Float_t supModAngle = fGeoPar->GetAngleOfSuperModule();
  144. Float_t modAngle = fGeoPar->GetAngleOfModule();
  145. Int_t modIdInSupModByPhi = mod / fGeoPar->GetNModInSuperModByPhi(); // 0, 1, 2
  146. Float_t sectorAngleEdge = sectorAngle * (sec % fGeoPar->GetNsec());
  147. Float_t supModAngleEdge = supModAngle * supMod;
  148. Float_t modAngleEdge = modAngle * modIdInSupModByPhi;
  149. return sectorAngleEdge + supModAngleEdge + modAngleEdge + modAngle / 2.0;
  150. }
  151. MpdEmcHit* MpdEmcHitProducer::SearchHit(UInt_t sec, UInt_t supMod, UInt_t row, UInt_t mod) {
  152. MpdEmcHit* foundHit = NULL;
  153. for (Int_t i = 0; i < fDigiArray->GetEntriesFast(); ++i) {
  154. MpdEmcHit* hit = (MpdEmcHit*) fDigiArray->At(i);
  155. if ((UInt_t)hit->GetSec() != sec) continue;
  156. if ((UInt_t)hit->GetSupMod() != supMod) continue;
  157. if ((UInt_t)hit->GetRow() != row) continue;
  158. if ((UInt_t)hit->GetMod() != mod) continue;
  159. foundHit = hit;
  160. }
  161. return foundHit;
  162. }
  163. // -------------------------------------------------------------------------
  164. ClassImp(MpdEmcHitProducer)