MpdEmcDigitizer.cxx 9.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295
  1. //--------------------------------------------------------------------
  2. //
  3. // Description:
  4. // MPD EMC Digitizer - takes EmcPoints and makes digits
  5. //
  6. //
  7. // Author List:
  8. // Alexander Zinchenko LHEP, JINR, Dubna - 8-May-2016
  9. // Alexander Zinchenko LHEP, JINR, Dubna - 24-June-2018 - adapted to projective geometry
  10. //
  11. //--------------------------------------------------------------------
  12. #include "MpdEmcDigitizer.h"
  13. #include "MpdEmcGeoParams.h"
  14. #include "MpdEmcDigit.h"
  15. #include "MpdEmcPoint.h"
  16. #include "FairMCPoint.h"
  17. #include "MpdMCTrack.h"
  18. #include "FairRootManager.h"
  19. #include "FairRunAna.h"
  20. #include <TClonesArray.h>
  21. #include <TGeoBBox.h>
  22. #include <TGeoManager.h>
  23. #include <TMath.h>
  24. using namespace std;
  25. using namespace TMath;
  26. FILE *lun = fopen ("file.txt","w"); //AZ
  27. // ----- Default constructor -------------------------------------------
  28. MpdEmcDigitizer::MpdEmcDigitizer() :
  29. FairTask("EMC digitizer") {
  30. }
  31. // ----- Destructor ----------------------------------------------------
  32. MpdEmcDigitizer::~MpdEmcDigitizer() {
  33. }
  34. // -------------------------------------------------------------------------
  35. // ----- Public method Init --------------------------------------------
  36. InitStatus MpdEmcDigitizer::Init() {
  37. cout << "\n******************* EMC Digitizer INIT *********************" << endl;
  38. // Get RootManager
  39. FairRootManager* ioman = FairRootManager::Instance();
  40. if (!ioman) {
  41. cout << "-E- MpdEmcDigitizer::Init: "
  42. << "RootManager not instantiated!" << endl;
  43. return kFATAL;
  44. }
  45. fGeoPar = new MpdEmcGeoParams();
  46. // Get input array
  47. fPointArray = (TClonesArray*) ioman->GetObject("EmcPoint");
  48. if (!fPointArray) {
  49. cout << "-W- MpdEmcDigitizer::Init: " << "No EmcPoint array!" << endl;
  50. return kERROR;
  51. }
  52. fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  53. if (!fMcTrackArray) {
  54. cout << "-W- MpdEmcDigitizer::Init: " << "No MCTrack array!" << endl;
  55. return kERROR;
  56. }
  57. // Create and register output array
  58. fDigiArray = new TClonesArray("MpdEmcDigit", 100);
  59. ioman->Register("EmcDigit", "EMC", fDigiArray, kTRUE);
  60. cout << "-I- MpdEmcDigitizer: Intialization successfull" << endl;
  61. return kSUCCESS;
  62. }
  63. //__________________________________________________________________________
  64. void MpdEmcDigitizer::Finish() {
  65. cout << "-I- MpdEmcDigitizer: Finish" << endl;
  66. }
  67. //__________________________________________________________________________
  68. void MpdEmcDigitizer::Exec(Option_t* opt)
  69. {
  70. // Main processing engine
  71. cout << "MpdEmcDigitizer::Exec started" << endl;
  72. static Int_t nSecRows = -1;
  73. if (nSecRows < 0)
  74. nSecRows = fGeoPar->GetNrows() / (fGeoPar->GetNsec() / 2 - 1); // 6 wide, 2 narrow sectors
  75. // Reset output Array
  76. if (!fDigiArray) Fatal("MpdEmcDigitizer::Exec", "No array of digits");
  77. fDigiArray->Delete();
  78. fHitMap.clear();
  79. Int_t nPoints = fPointArray->GetEntriesFast();
  80. for (Int_t iPnt = 0; iPnt < nPoints; ++iPnt) {
  81. FairMCPoint* pnt = (FairMCPoint*) fPointArray->UncheckedAt(iPnt);
  82. Int_t trId = pnt->GetTrackID();
  83. if (trId < 0) continue; // strange case
  84. MpdMCTrack* tr = (MpdMCTrack*) fMcTrackArray->UncheckedAt(trId);
  85. // if (tr->GetMotherId() != -1) continue;
  86. Int_t pdg = tr->GetPdgCode();
  87. Double_t x = pnt->GetX();
  88. Double_t y = pnt->GetY();
  89. Double_t z = pnt->GetZ();
  90. Float_t e = pnt->GetEnergyLoss();
  91. gGeoManager->FindNode(x,y,z);
  92. TString path(gGeoManager->GetPath());
  93. //cout << path << endl;
  94. //if (!path.Contains("sc")) exit(0);
  95. //if (!path.Contains("sc")) continue; // !!! this condition is very bad for Geant4 (not for Geant3)
  96. if (!path.Contains("box")) continue; //
  97. if (path.Contains("_cl_")) gGeoManager->CdUp(); // to tower (channel)
  98. Int_t ip = path.Index("ChH"); // half-EMC
  99. UInt_t ih = TString(path(ip+4,1)).Atoi();
  100. ip = path.Index("Sector"); // sector
  101. Int_t ip1 = path.Index("/",ip);
  102. UInt_t isec = TString(path(ip+8,ip1)).Atoi();
  103. if (path(ip+6) == '2') {
  104. // Narrow sectors
  105. if (isec == 0) isec = 2;
  106. else isec = 6;
  107. } else {
  108. if (isec >= 5) isec += 2;
  109. else if (isec >= 2) ++isec;
  110. }
  111. ip = path.Index("Crate"); // row
  112. ip1 = path.Index("/",ip);
  113. UInt_t irow = TString(path(ip+6,ip1)).Atoi();
  114. ip = path.Index("box");
  115. if (path(ip+4) == '_') --ip; // 1-digit number
  116. ip1 = path.Index("/",ip);
  117. UInt_t ibox = TString(path(ip+6,ip1)).Atoi();
  118. //cout << ih << " " << isec << " " << irow << " " << ibox << endl;
  119. fprintf(lun,"%d %d %d %d %f %f %f %f\n",ih,isec,irow,ibox,TMath::ATan2(y,x),x,y,z);
  120. // Code channel ID
  121. TString tower = "c"; // channel
  122. tower += ibox;
  123. tower += "r"; // row
  124. tower += irow;
  125. tower += "s"; // sector
  126. tower += isec;
  127. tower += "h"; // EMC half (Z <> 0)
  128. tower += ih;
  129. MpdEmcDigit* hit = SearchHit(tower);
  130. if (hit == NULL) {
  131. hit = new((*fDigiArray)[fDigiArray->GetEntriesFast()]) MpdEmcDigit(ih, isec, irow, ibox);
  132. fHitMap[tower] = hit;
  133. Double_t phi, the;
  134. //FindChanPhiZ(phi,z);
  135. FindChanPhiThe(phi,the);
  136. hit->SetPhiCenter(phi);
  137. hit->SetZCenter(the);
  138. Int_t iz = ibox;
  139. if (ih == 1) ++iz; // negative Z
  140. hit->SetChanZId(TMath::Abs(iz));
  141. Int_t iphi = irow + isec * nSecRows;
  142. if (isec >= 7) iphi -= nSecRows; // take into account narrow sector
  143. else if (isec >= 3) iphi -= nSecRows / 2; // take into account narrow sector
  144. hit->SetChanPhiId(iphi);
  145. hit->SetTimeStamp(pnt->GetTime());
  146. }
  147. hit->IncreaseEnergy(e, trId);
  148. hit->SetTimeStamp(TMath::Min(pnt->GetTime(),hit->GetTimeStamp()));
  149. hit->SetNumTracks(hit->GetNumTracks() + 1);
  150. //hit->SetTrackId(trId);
  151. hit->SetPdg(pdg);
  152. //AZ hit->SetZCenter(CalcZCenter(sec, row, mod));
  153. //AZ hit->SetPhiCenter(CalcPhiCenter(sec, supMod, mod));
  154. // Redo timing - take time from the highest energy deposit
  155. if (TMath::Abs(e-hit->GetDy()) < 1.e-6) {
  156. hit->SetTimeStamp(pnt->GetTime()); // Dy - holder of Emax
  157. hit->SetTrackId(iPnt); // point index
  158. }
  159. }
  160. Int_t nHits = fDigiArray->GetEntriesFast();
  161. for (Int_t iHit = 0; iHit < nHits; ++iHit) {
  162. MpdEmcDigit* hit = (MpdEmcDigit*) fDigiArray->UncheckedAt(iHit);
  163. if (hit->GetNumTracks() > 1) {
  164. hit->SetPdg(0);
  165. //hit->SetTrackId(-1);
  166. }
  167. //hit->Print();
  168. }
  169. // Redo track ID numbering
  170. //RedoId();
  171. }
  172. //__________________________________________________________________________
  173. void MpdEmcDigitizer::RedoId(TClonesArray *digis, TClonesArray *mctrs)
  174. {
  175. // Redo track ID numbering
  176. // Take IDs of particles produced outside EMC
  177. static const Double_t rmin = fGeoPar->GetRmin(), rmax = fGeoPar->GetRmax();
  178. TVector3 vert;
  179. Int_t nDigis = digis->GetEntriesFast();
  180. for (Int_t iDig = 0; iDig < nDigis; ++iDig) {
  181. MpdEmcDigit* digi = (MpdEmcDigit*) digis->UncheckedAt(iDig);
  182. map<Int_t,Float_t> contrib = digi->GetContrib();
  183. map<Int_t,Float_t> copy(contrib);
  184. contrib.clear();
  185. for (map<Int_t,Float_t>::iterator it = copy.begin(); it != copy.end(); ++it) {
  186. Int_t id = it->first, idm = -1;
  187. MpdMCTrack *mctr = (MpdMCTrack*) mctrs->UncheckedAt(id);
  188. mctr->GetStartVertex(vert);
  189. Double_t rvert = vert.Pt();
  190. while (rvert > rmin && rvert < rmax) {
  191. // Born inside EMC - find ancestor
  192. idm = mctr->GetMotherId();
  193. if (idm < 0) break;
  194. mctr = (MpdMCTrack*) mctrs->UncheckedAt(idm);
  195. mctr->GetStartVertex(vert);
  196. rvert = vert.Pt();
  197. }
  198. if (idm >= 0) id = idm;
  199. if (contrib.find(id) == contrib.end()) contrib[id] = it->second;
  200. else contrib[id] += it->second;
  201. }
  202. }
  203. }
  204. //__________________________________________________________________________
  205. void MpdEmcDigitizer::FindChanPhiZ(Double_t &phi, Double_t &z)
  206. {
  207. // Compute Z-position and Phi-angle of the tower (channel) center (at inner radius)
  208. TString path(gGeoManager->GetPath());
  209. Int_t ip = path.Last('/');
  210. Int_t ip1 = path.Last('_');
  211. TString volName = path(ip+1,ip1-ip-1);
  212. TGeoVolume *tower = gGeoManager->GetVolume(volName);
  213. TGeoBBox *box = (TGeoBBox*) tower->GetShape();
  214. Double_t xyzL[3] = {0,0,-box->GetDZ()}, xyzM[3];
  215. gGeoManager->LocalToMaster(xyzL,xyzM);
  216. z = xyzM[2];
  217. phi = TMath::ATan2 (xyzM[1],xyzM[0]);
  218. if (phi < 0) phi += TMath::TwoPi();
  219. }
  220. //__________________________________________________________________________
  221. void MpdEmcDigitizer::FindChanPhiThe(Double_t &phi, Double_t &the)
  222. {
  223. // Compute Phi and Theta angles of the tower (channel) center (at inner radius)
  224. TString path(gGeoManager->GetPath());
  225. Int_t ip = path.Last('/');
  226. Int_t ip1 = path.Last('_');
  227. TString volName = path(ip+1,ip1-ip-1);
  228. TGeoVolume *tower = gGeoManager->GetVolume(volName);
  229. TGeoBBox *box = (TGeoBBox*) tower->GetShape();
  230. Double_t xyzL[3] = {0,0,-box->GetDZ()}, xyzM[3];
  231. gGeoManager->LocalToMaster(xyzL,xyzM);
  232. phi = TMath::ATan2 (xyzM[1],xyzM[0]);
  233. if (phi < 0) phi += TMath::TwoPi();
  234. the = TMath::ATan2 (TMath::Sqrt(xyzM[0]*xyzM[0]+xyzM[1]*xyzM[1]), xyzM[2]);
  235. }
  236. //__________________________________________________________________________
  237. MpdEmcDigit* MpdEmcDigitizer::SearchHit(TString tower)
  238. {
  239. // Check if the hit (tower) exists in the map
  240. if (fHitMap.find(tower) == fHitMap.end()) return NULL;
  241. return fHitMap[tower];
  242. }
  243. // -------------------------------------------------------------------------
  244. ClassImp(MpdEmcDigitizer)