MpdZdcDigiProducer.cxx 8.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302
  1. /*************************************************************************************
  2. *
  3. * MpdZdcDigiProducer
  4. * Class to create digital data taken from MpdZdc detector
  5. *
  6. * Author: Elena Litvinenko
  7. * e-mail: litvin@nf.jinr.ru
  8. * Version: 18-Apr-2008
  9. * Modified March 2021 by A.Strijak
  10. *
  11. ************************************************************************************/
  12. #include <iostream>
  13. #include "TClonesArray.h"
  14. #include "FairRootManager.h"
  15. #include "FairRun.h"
  16. #include "FairRunAna.h"
  17. #include "FairRuntimeDb.h"
  18. #include "MpdZdcDigiProducer.h"
  19. #include "MpdZdcDigi.h"
  20. #include "MpdZdcPoint.h"
  21. #include "TROOT.h"
  22. #include "TVectorT.h"
  23. // ----- Default constructor -------------------------------------------
  24. MpdZdcDigiProducer::MpdZdcDigiProducer(const char* name) :
  25. FairTask(name) {
  26. fPointArray=0;
  27. fDigiArray=0;
  28. fGeoPar=0;
  29. // fELossZdc1Value = NULL, fELossZdc2Value = NULLL;
  30. fPix2Mip = 15; // 15 MPPC pixels per MIP
  31. fMIPEnergy = 0.005; // 5 MeV
  32. fMIPNoise = 0.3; // 0.2 MIP noise level
  33. fMIP2GeV = 0.005;
  34. }
  35. // -------------------------------------------------------------------------
  36. // ----- Destructor ----------------------------------------------------
  37. MpdZdcDigiProducer::~MpdZdcDigiProducer() { }
  38. // -------------------------------------------------------------------------
  39. // -------------------------------------------------------------------------
  40. void MpdZdcDigiProducer::SetParContainers()
  41. {
  42. cout << "-I- MpdZdcDigiProducer: SetParContainers started..." << endl;
  43. // Get run and runtime database
  44. FairRunAna* run = FairRunAna::Instance();
  45. if ( ! run ) Fatal("FairMuchDigitize::SetParContainers", "No analysis run");
  46. FairRuntimeDb* rtdb = run->GetRuntimeDb();
  47. if ( ! rtdb ) Fatal("FairMuchDigitize::SetParContainers", "No runtime database");
  48. cout << "-I- MpdZdcDigiProducer: SetParContainers continued..." << endl;
  49. rtdb->activateParIo(rtdb->getFirstInput());
  50. // fGeoPar=( MpdZdcGeoPar*) rtdb->getContainer("MpdZdcGeoPar");
  51. fGeoPar=( MpdZdcGeoPar*) gROOT->FindObject("MpdZdcGeoPar");
  52. fGeoPar->print();
  53. cout << "-I- MpdZdcDigiProducer: SetParContainers finished." << endl;
  54. }
  55. // -------------------------------------------------------------------------
  56. // ----- Public method Init --------------------------------------------
  57. InitStatus MpdZdcDigiProducer::Init() {
  58. cout << "-I- MpdZdcDigiProducer: Init started..." << endl;
  59. fRandom3 = new TRandom3();
  60. // Get RootManager
  61. FairRootManager* ioman = FairRootManager::Instance();
  62. if ( ! ioman ) {
  63. cout << "-E- MpdZdcDigiProducer::Init: "
  64. << "RootManager not instantiated!" << endl;
  65. return kFATAL;
  66. }
  67. // Get input array
  68. fPointArray = (TClonesArray*) ioman->GetObject("ZdcPoint");
  69. if ( ! fPointArray ) {
  70. cout << "-W- MpdZdcDigiProducer::Init: "
  71. << "No ZdcPoint array!" << endl;
  72. return kERROR;
  73. }
  74. // Create and register output array
  75. fDigiArray = new TClonesArray("MpdZdcDigi");
  76. ioman->Register("ZdcDigi","Zdc",fDigiArray,kTRUE);
  77. /*
  78. fELossZdc1Value = new TClonesArray("TParameter<double>");
  79. ioman->Register("ELossZdc1Value","Zdc",fELossZdc1Value,kTRUE);
  80. fELossZdc2Value = new TClonesArray("TParameter<double>");
  81. ioman->Register("ELossZdc2Value","Zdc",fELossZdc2Value,kTRUE);
  82. */
  83. MpdZdcDigiScheme *fDigiScheme = MpdZdcDigiScheme::Instance();
  84. fDigiScheme->Init(fGeoPar,0,kTRUE,2);
  85. cout << "-I- MpdZdcDigiProducer: Intialization successfull" << endl;
  86. return kSUCCESS;
  87. }
  88. // ----- Public method Exec --------------------------------------------
  89. void MpdZdcDigiProducer::Exec(Option_t* opt) {
  90. //#define EDEBUG
  91. #ifdef EDEBUG
  92. static Int_t lEDEBUGcounter=0;
  93. cout << "EDEBUG-- MpdZdcDigiProducer::Exec() started... " << endl;;
  94. #endif
  95. if ( ! fDigiArray ) Fatal("Exec", "No DigiArray");
  96. fDigiArray->Clear();
  97. MpdZdcDigiScheme *pDigiScheme = MpdZdcDigiScheme::Instance();
  98. if (!pDigiScheme)
  99. Fatal("MpdZdcDigiProducer::Exec", "No DigiScheme");
  100. Int_t detID, modID, chanID;
  101. MpdZdcDigiId_t digiID;
  102. //marina
  103. Double_t dEdepSectEv[90][10];
  104. for(Int_t i=0; i<90; i++) {// mod
  105. for(Int_t ii=0;ii<10;ii++) { // section
  106. dEdepSectEv[i][ii] = 0.;
  107. }
  108. }
  109. //end marina
  110. MpdZdcPoint* point = NULL;
  111. map<MpdZdcDigiId_t, Float_t> fDigiIdEnergy;
  112. fDigiIdEnergy.clear();
  113. map<MpdZdcDigiId_t, Float_t>::const_iterator p;
  114. Int_t nPoints = fPointArray->GetEntriesFast();
  115. Double_t e1=0, e2=0;
  116. //cout <<"marina " <<nPoints <<endl;
  117. for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
  118. //marina
  119. point = (MpdZdcPoint*) fPointArray->At(iPoint);
  120. //cout <<"marina 1 " <<point->GetCopyZdc() <<" " <<point->GetCopyMother() <<" " <<point->GetCopy() <<endl;
  121. detID = point->GetCopyZdc();//==1 (z>0), ==2 (z<0)
  122. //modID = point->GetCopyMother(); // modules 1-45
  123. if(detID==1) modID = point->GetCopyMother(); // modules 1-45
  124. else modID = point->GetCopyMother() + 45;//modules 46-90
  125. chanID = (Int_t)((point->GetCopy()-1)/6); //sections 0-9
  126. // cout <<"marina 1 " <<detID <<" " <<modID <<" " <<point->GetCopy() <<' ' <<chanID <<endl;
  127. dEdepSectEv[modID-1][chanID]+=point->GetEnergyLoss();
  128. //end marina
  129. /*
  130. if (detID == 1) {
  131. e1 += point->GetEnergyLoss();
  132. }
  133. else e2 += point->GetEnergyLoss();
  134. */
  135. Int_t pMMcopy=1*(point->GetZ()>0)+2*(point->GetZ()<0);
  136. digiID = pDigiScheme->GetDigiIdFromVolumeData (point->GetDetectorID(), point->GetCopy(), point->GetCopyMother(),pMMcopy);
  137. if ((digiID[0]!=-1)&&(digiID[1]!=-1)) {
  138. if (fDigiIdEnergy.find(digiID)==fDigiIdEnergy.end())
  139. fDigiIdEnergy[digiID] = point->GetEnergyLoss();
  140. else
  141. fDigiIdEnergy[digiID] += point->GetEnergyLoss();
  142. if (pMMcopy==1) {
  143. e1 += point->GetEnergyLoss();
  144. }
  145. else {
  146. e2 += point->GetEnergyLoss();
  147. }
  148. }//if ((digiID[0]!=-1)&&(digiID[1]!=-1))
  149. #ifdef EDEBUG
  150. else {
  151. if (lEDEBUGcounter<100) {
  152. cout << "EDEBUG-- MpdZdcDigiProducer::Exec: Boundary point? : "; point->Print("");
  153. lEDEBUGcounter++;
  154. }
  155. }
  156. #endif
  157. }//for (Int_t iPoint=0; iPoint<nPoints; iPoint++)
  158. //cout <<"marina 2 " <<endl;
  159. /*
  160. TClonesArray& clref1 = *fELossZdc1Value;
  161. new(clref1[0]) TParameter<double>("ELossZdc1",e1);
  162. TClonesArray& clref2 = *fELossZdc2Value;
  163. new(clref2[0]) TParameter<double>("ELossZdc2",e2);
  164. */
  165. e1 = 0;
  166. e2 = 0;
  167. //cout <<"marina 3 " <<endl;
  168. for(Int_t i=0; i<90; i++) {// mod
  169. for(Int_t ii=0;ii<10;ii++) { // section
  170. //cout <<"dEdepSectEv " <<i <<" " <<ii <<" " <<dEdepSectEv[i][ii] <<endl;
  171. if(dEdepSectEv[i][ii]>0) {
  172. if(i<=44) detID = 1;
  173. else detID = 2;
  174. Double_t recEnergy = RecoEnergy(dEdepSectEv[i][ii]);
  175. if(detID==1) {
  176. MpdZdcDigi* digi = AddHit(detID, i+1, ii+1, dEdepSectEv[i][ii]);
  177. digi->ConvertSim();
  178. digi->SetELossReco(recEnergy);
  179. e1 += recEnergy;
  180. }
  181. else {
  182. MpdZdcDigi* digi = AddHit(detID, i-45+1, ii+1, dEdepSectEv[i][ii]);
  183. digi->ConvertSim();
  184. digi->SetELossReco(recEnergy);
  185. e2 += recEnergy;
  186. }
  187. }
  188. }
  189. }//for(Int_t i=0; i<90; i++)
  190. // TClonesArray& clref1 = *fELossZdc1Value;
  191. // new(clref1[0]) TParameter<double>("ELossZdc1",e1);
  192. // TClonesArray& clref2 = *fELossZdc2Value;
  193. // new(clref2[0]) TParameter<double>("ELossZdc2",e2);
  194. /*
  195. for(p=fDigiIdEnergy.begin(); p!=fDigiIdEnergy.end(); ++p) {
  196. pDigiScheme->SplitDigiID((*p).first, detID, modID, chanID);
  197. if ((detID!=-1)&&(chanID!=-1)) {
  198. MpdZdcDigi* digi = AddHit(detID, modID, chanID, (*p).second);
  199. digi->ConvertSim();
  200. #ifdef EDEBUG
  201. if (lEDEBUGcounter<20) {
  202. cout << "EDEBUG-- MpdZdcDigiProducer::Exec: "<< detID<< " " << chanID << " " <<
  203. (*p).second << " " << lEDEBUGcounter << endl;
  204. lEDEBUGcounter++;
  205. }
  206. #endif
  207. }
  208. }
  209. */
  210. #undef EDEBUG
  211. }
  212. // -------------------------------------------------------------------------
  213. // ----- Private method AddDigi --------------------------------------------
  214. MpdZdcDigi* MpdZdcDigiProducer::AddHit(Int_t detID, Int_t modID, Int_t chanID,Float_t energy)
  215. {
  216. TClonesArray& clref = *fDigiArray;
  217. Int_t size = clref.GetEntriesFast();
  218. //cout <<"size " <<size <<endl;
  219. MpdZdcDigi* result = new(clref[size]) MpdZdcDigi(detID,modID,chanID,energy);
  220. //cout <<"result " <<result <<endl;
  221. return result;
  222. }
  223. // ----
  224. Double_t MpdZdcDigiProducer::RecoEnergy (Double_t pfELoss)
  225. {
  226. Double_t energyMIP = pfELoss / fMIPEnergy;
  227. Double_t energyPix = fRandom3->Poisson(energyMIP * fPix2Mip);
  228. Double_t energyMIPSmeared = energyPix / fPix2Mip;
  229. Double_t noise = fRandom3->Gaus(0, fMIPNoise);
  230. energyMIPSmeared += noise;
  231. return energyMIPSmeared * fMIP2GeV;
  232. }
  233. ClassImp(MpdZdcDigiProducer)