MpdFfdHitProducer.cxx 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdFfdPoint
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <iostream>
  8. #include <TClonesArray.h>
  9. #include <TGraph.h>
  10. #include <TFile.h>
  11. #include <TRandom3.h>
  12. #include "FairRootManager.h"
  13. #include "FairDetector.h"
  14. #include "FairBaseParSet.h"
  15. #include "FairGeoNode.h"
  16. #include "FairRunAna.h"
  17. #include "FairRun.h"
  18. #include "FairRuntimeDb.h"
  19. #include "MpdMCTrack.h"
  20. #include "MpdFfdGeoPar.h"
  21. #include "MpdFfdHit.h"
  22. #include "MpdFfdPoint.h"
  23. #include "MpdFfdHitProducer.h"
  24. using namespace std;
  25. ClassImp(MpdFfdHitProducer)
  26. TGraph* MpdFfdHitProducer::gPMTeff_pdf = nullptr;
  27. //------------------------------------------------------------------------------------------------------------------------
  28. MpdFfdHitProducer::MpdFfdHitProducer(const char* name, bool useMCdata, bool timeWindow, Int_t verbose, const char* flnm)
  29. : FairTask(name, verbose), fUseMCData(useMCdata), fUseTimeWindow(timeWindow), fFlnm(flnm)
  30. {
  31. fDoTest = (flnm != nullptr);
  32. if(fDoTest)
  33. {
  34. Add(hSuids = new TH1D("FFD_suids", "suid;entries", 200, -0.5, 199.5));
  35. Add(hOpYield = new TH2D("FFD_opYield", "Yield of op(per quartz 1 cm) vs #beta;#beta;N op", 1000, 0.1, 1.01, 1000, 0.5, 10000.5));
  36. Add(hOpYieldPion = new TH2D("FFD_opYieldPion", "Yield of op(per quartz 1 cm) vs #beta;#beta;N op", 1000, 0.1, 1.01, 1000, 0.5, 10000.5));
  37. Add(hOpYieldProton = new TH2D("FFD_opYieldProton", "Yield of op(per quartz 1 cm) vs #beta;#beta;N op", 1000, 0.1, 1.01, 1000, 0.5, 10000.5));
  38. Add(hOpYieldElectron = new TH2D("FFD_opYieldElectron", "Yield of op(per quartz 1 cm) vs #beta;#beta;N op", 1000, 0.1, 1.01, 1000, 0.5, 10000.5));
  39. Add(hPeYield = new TH2D("FFD_peYield", "Yield of op(per quartz 1 cm) vs #beta;#beta;N pe", 1000, 0.1, 1.01, 1000, 0.5, 1000.5));
  40. Add(hPeYieldPion = new TH2D("FFD_peYieldPion", "Yield of op(per quartz 1 cm) vs #beta;#beta;N pe", 1000, 0.1, 1.01, 1000, 0.5, 1000.5));
  41. Add(hPeYieldProton = new TH2D("FFD_peYieldProton", "Yield of op(per quartz 1 cm) vs #beta;#beta;N pe", 1000, 0.1, 1.01, 1000, 0.5, 1000.5));
  42. Add(hPeYieldElectron = new TH2D("FFD_peYieldElectron", "Yield of op(per quartz 1 cm) vs #beta;#beta;N pe", 1000, 0.1, 1.01, 1000, 0.5, 1000.5));
  43. Add(hOpTrash = new TH2D("FFD_opTrashRatio", "trash op ratio;N op; trash op, %", 1000, 0.5, 1000000.5, 1000, 0., 50.));
  44. Add(hOccup = new TH2D("FFD_Occup", ";occupancy;suid", 100, -0.5, 99.5, 200, -0.5, 199.5));
  45. Add(hChannelsEW = new TH2D("FFD_ChannelsEW", ";channels, east;channels, west", 100, -0.5, 99.5, 100, -0.5, 99.5));
  46. Add(hTimeWindows = new TH2D("FFD_TimeWindows", ";t_{mean}, ns;#Deltat, ns", 1000, 0., 10., 1000, 0., 1.));
  47. Add(hTimeSize = new TH2D("FFD_TimeSize", ";t_{mean}, ns;n pe", 1000, 0., 100., 1000, 0.5, 1000.5));
  48. Add(hXY = new TH2D("FFD_XY", "; X, cm; Y, cm", 1000, -20., 20., 1000, -20., 20.));
  49. Add(hXmap = new TH2D("FFD_Xmap", "; X, cm; suid", 1000, -20., 20., 200, -0.5, 199.5));
  50. Add(hYmap = new TH2D("FFD_Ymap", "; Y, cm; suid", 1000, -20., 20., 200, -0.5, 199.5));
  51. Add(hZmap = new TH2D("FFD_Zmap", "; Z, cm; suid", 1000, -150., 150., 200, -0.5, 199.5));
  52. Add(hXcenter = new TH2D("FFD_Xcenter", "; X, cm; suid", 1000, -20., 20., 200, -0.5, 199.5));
  53. Add(hYcenter = new TH2D("FFD_Ycenter", "; Y, cm; suid", 1000, -20., 20., 200, -0.5, 199.5));
  54. Add(hCenter = new TH2D("FFD_XYcenterShift", "; X, cm; Y, cm", 1000, -50., 50., 1000, -50., 50.));
  55. Add(hPMTeff = new TEfficiency("FFD_PMTeff", ";op energy, eV;PMT efficiency", 100, 1., 9.));
  56. }
  57. }
  58. //------------------------------------------------------------------------------------------------------------------------
  59. MpdFfdHitProducer::~MpdFfdHitProducer()
  60. {
  61. fList.Delete();
  62. }
  63. //------------------------------------------------------------------------------------------------------------------------
  64. InitStatus MpdFfdHitProducer::Init()
  65. {
  66. if(fUseMCData)
  67. {
  68. aMcPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("FFDPoint");
  69. aMcTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
  70. assert(aMcPoints);
  71. assert(aMcTracks);
  72. }
  73. else
  74. {
  75. aExpDigits = (TClonesArray*) FairRootManager::Instance()->GetObject("??????");// FIXME: NOW unknown name
  76. assert(aExpDigits);
  77. }
  78. // Create and register output array
  79. aFfdHits = new TClonesArray("MpdFfdHit");
  80. FairRootManager::Instance()->Register("FfdHit", "Ffd", aFfdHits, kTRUE);
  81. LOG(INFO)<<"[MpdFfdHitProducer::Init] Initialization finished succesfully.";
  82. return kSUCCESS;
  83. }
  84. //------------------------------------------------------------------------------------------------------------------------
  85. void MpdFfdHitProducer::Exec(Option_t* opt)
  86. {
  87. // Cleanup containers
  88. aFfdHits->Clear();
  89. mmOccup.clear();
  90. if(fUseTimeWindow)
  91. {
  92. for(size_t i = 0; i < 160; i++) fData[i].clear();
  93. }
  94. auto CenterTest = [this](MpdFfdPoint* point, Int_t suid, MpdFfdHit* hit, const TVector3& padCenter)
  95. {
  96. mmOccup.insert(make_pair(suid, hit));
  97. TVector3 pos; point->Position(pos);
  98. double padCenterX = padCenter.X(), padCenterY = padCenter.Y();
  99. hCenter->Fill(padCenterX - pos.X(), padCenterY - pos.Y());
  100. hXcenter->Fill(padCenterX, suid);
  101. hYcenter->Fill(padCenterY, suid);
  102. };
  103. if(fUseMCData) // input data from MC MpdFfdPoint
  104. {
  105. size_t nOpTotal = 0, nOpTrash = 0;
  106. for(Int_t pointIndex = 0, nFfdPoint = aMcPoints->GetEntriesFast(); pointIndex < nFfdPoint; pointIndex++ ) // cycle by FFD points
  107. {
  108. auto pPoint = (MpdFfdPoint*) aMcPoints->UncheckedAt(pointIndex);
  109. Int_t suid = pPoint->GetDetectorID(); // [1,160]
  110. Int_t tid = pPoint->GetTrackID();
  111. double beta = pPoint->GetBeta(), nOp = pPoint->GetEntries();
  112. nOpTotal += nOp;
  113. if(! pPoint->IsClosed()) // invalid parent track parameters, skip MpdFfdPoint data.
  114. {
  115. nOpTrash += nOp;
  116. continue;
  117. }
  118. if(fDoTest)
  119. {
  120. hOpYield->Fill(beta, nOp);
  121. hSuids->Fill(suid);
  122. double x = pPoint->GetX(), y = pPoint->GetY();
  123. hXY->Fill(x, y);
  124. hXmap->Fill(x, suid);
  125. hYmap->Fill(y, suid);
  126. hZmap->Fill(pPoint->GetZ(), suid);
  127. }
  128. if(fVerbose > 2) pPoint->Print("");
  129. size_t nOp2_8 = 0, nPe = 0;
  130. if(pPoint->GetMode() == MpdFfdPoint::kCherenkovPhoton)
  131. {
  132. for(const auto& entry : pPoint->GetData()) // cycle by op
  133. {
  134. double energy = entry.first;
  135. if(2. < energy && energy < 8.) nOp2_8++;; // <<<----------- [2,8] eV cut
  136. bool pass = IsPeCreated(energy); // simulate PMT efficiency
  137. if(pass)
  138. {
  139. nPe++;
  140. if(fUseTimeWindow) AddEntry(entry.second, suid, pPoint, pointIndex);
  141. }
  142. hPMTeff->Fill(pass, energy);
  143. } // cycle by op
  144. }
  145. else // kPhotoElectron mode, saved only pe
  146. {
  147. nPe = pPoint->GetEntries();
  148. if(fUseTimeWindow)
  149. {
  150. for(const auto& entry : pPoint->GetData()) AddEntry(entry.second, suid, pPoint, pointIndex); // cycle by pe
  151. }
  152. }
  153. if(! fUseTimeWindow && nPe >= fNpeThresh) // create FFD hit at fast mode
  154. {
  155. TVector3 padCenter = GetPadCenter(suid);
  156. auto hit = AddHit(pPoint, padCenter, TVector3(fErrXY, fErrXY, fErrZ), pointIndex, nPe, 0);
  157. if(fDoTest) CenterTest(pPoint, suid, hit, padCenter);
  158. }
  159. if(fDoTest)
  160. {
  161. hPeYield->Fill(beta, nPe);
  162. auto mcTrack = (const MpdMCTrack*) aMcTracks->UncheckedAt(tid);
  163. Int_t pid = mcTrack->GetPdgCode();
  164. if(2212 == pid || -2212 == pid) { if(nOp2_8)hOpYieldProton->Fill(beta, nOp2_8); hPeYieldProton->Fill(beta, nPe); }
  165. else if(211 == pid || -211 == pid){ if(nOp2_8)hOpYieldPion->Fill(beta, nOp2_8); hPeYieldPion->Fill(beta, nPe); }
  166. else if(11 == pid || -11 == pid) { if(nOp2_8)hOpYieldElectron->Fill(beta, nOp2_8); hPeYieldElectron->Fill(beta, nPe); }
  167. assert(! pPoint->GetData().empty());
  168. double Tmin = 1.E10, Tmax = -1., Tmean = 0.;
  169. for(const auto& entry : pPoint->GetData()) // cycle by pe
  170. {
  171. double time = entry.second;
  172. if(Tmin > time) Tmin = time;
  173. if(time > Tmax) Tmax = time;
  174. Tmean += time;
  175. }
  176. Tmean /= pPoint->GetData().size();
  177. hTimeWindows->Fill(Tmean, Tmax - Tmin);
  178. hTimeSize->Fill(Tmean, nPe);
  179. }
  180. } // cycle by FFD points
  181. if(fUseTimeWindow)
  182. {
  183. Int_t nWindowSize = fTimeWindow * fRoIBins/fRoISize; nWindowSize = ((UInt_t)nWindowSize < fRoIBins) ? nWindowSize : fRoIBins - 1;
  184. for(size_t suid = 0; suid < 160; suid++)
  185. {
  186. auto data = fData[suid];
  187. for(auto iter = data.begin(), itEnd = data.end(); iter != itEnd; iter++) // cycle by non-zero deposit time bins
  188. {
  189. Tintegrals mIntegrals; // key = MpdFfdPoint index, v
  190. double timeWindowEnd = GetTime(iter->first) + fTimeWindow;
  191. size_t bin;
  192. auto itWindowEnd = GetTimeBin(timeWindowEnd, bin) ? data.upper_bound(bin) : itEnd;
  193. for(auto it = iter; it != itWindowEnd; it++) // cycle by window bins
  194. {
  195. auto entry = mIntegrals.find(it->second.pointIndex);
  196. if(entry != mIntegrals.end()) entry->second.nPe += it->second.nPe;
  197. else mIntegrals.insert(make_pair(it->second.pointIndex, it->second));
  198. }
  199. size_t Integral = 0, maxDeposit = 0, maxIndex = -1;
  200. for(const auto& entry : mIntegrals) // cycle by MpdFfdPoint deposit into this time window
  201. {
  202. size_t nPe = entry.second.nPe;
  203. Integral += nPe;
  204. if(nPe > maxDeposit)
  205. {
  206. maxDeposit = nPe;
  207. maxIndex = entry.first;
  208. }
  209. }
  210. if(Integral >= fNpeThresh) // create FFD hit at TimeWindow mode
  211. {
  212. TVector3 padCenter = GetPadCenter(suid+1); // [0,159] -> [1,160]
  213. auto pPoint = (MpdFfdPoint*) aMcPoints->UncheckedAt(maxIndex);
  214. auto hit = AddHit(pPoint, padCenter, TVector3(fErrXY, fErrXY, fErrZ), maxIndex, maxDeposit, mIntegrals.size(), mIntegrals);
  215. if(fDoTest) CenterTest(pPoint, suid, hit, padCenter);
  216. iter = itWindowEnd; // integrate pe by only ONCE time
  217. if(iter == itEnd) break;
  218. }
  219. } // cycle by time bins
  220. } // cycle by suid
  221. } // fUseTimeWindow
  222. if(fDoTest)
  223. {
  224. if(nOpTotal) hOpTrash->Fill(nOpTotal, (100.* nOpTrash) / nOpTotal); // [%]
  225. size_t nE = 0, nW = 0;
  226. for(auto iter = mmOccup.begin(); iter != mmOccup.end(); iter = mmOccup.upper_bound(iter->first)) // cycle by unique channel hits
  227. {
  228. size_t suid = iter->first;
  229. hOccup->Fill(mmOccup.count(suid), suid);
  230. if(suid <= 80) nE++; else nW++;
  231. }
  232. hChannelsEW->Fill(nE, nW);
  233. }
  234. }
  235. else // input data from experimental MpdFfdDigit
  236. {
  237. // FIXME: now not realized
  238. //AddHit(Int_t detUID, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
  239. assert(false);
  240. }
  241. LOG(DEBUG1)<<"[MpdFfdHitProducer::Exec] FFD hits = "<<aFfdHits->GetEntriesFast();
  242. }
  243. //------------------------------------------------------------------------------------------------------------------------
  244. void MpdFfdHitProducer::AddEntry(double time, size_t suid, const MpdFfdPoint* ptr, size_t index)
  245. {
  246. assert(suid <= 160);
  247. suid--; // [1,160] -> [0,159]
  248. size_t timeBin;
  249. if(GetTimeBin(time, timeBin))
  250. {
  251. bool found = false;
  252. auto range = fData[suid].equal_range(timeBin);
  253. for(auto it = range.first; it != range.second; it++)
  254. {
  255. if(it->second.pointIndex == index)
  256. {
  257. it->second.nPe++;
  258. found = true;
  259. break;
  260. }
  261. }
  262. if(!found) fData[suid].insert(std::make_pair(timeBin, PointData(index, ptr, 1))); // insert first time deposit for pair(timebin, MpdFfdPoint)
  263. }
  264. }
  265. //------------------------------------------------------------------------------------------------------------------------
  266. TVector3 MpdFfdHitProducer::GetPadCenter(Int_t suid) // [1,161]
  267. {
  268. assert(1 <= suid && suid <= 161);
  269. static double *X = nullptr, *Y = nullptr;
  270. auto calcCenters = [this](double x, double y, size_t index, const char* comment = nullptr, std::ostream& os = std::cout)
  271. {
  272. TVector3 vect(x, y, 0.);
  273. TRotation rot; rot.RotateZ(45. * TMath::DegToRad());
  274. vect = rot * vect;
  275. X[index] = vect.X();
  276. Y[index] = vect.Y();
  277. if(fDoTest && fVerbose > 1)
  278. {
  279. if(comment) os<<comment;
  280. os<<"\n suid="<<(index + 1)<<" ("<<X[index]<<","<<Y[index]<<")";
  281. }
  282. };
  283. if(X == nullptr) // init by first time.
  284. {
  285. const double posX[] = { 7.9, 14.5, -7.9, -14.5, 0., 0., 0., 0., 6.6, -6.6, 6.6, -6.6, 6.6, -6.6, 6.6, -6.6, 13.2, -13.2, 13.2, -13.2 }; //[cm]
  286. const double posY[] = { 0., 0., 0., 0., 7.9, 14.5, -7.9, -14.5, 6.6, 6.6, -6.6, -6.6, 13.2, 13.2, -13.2, -13.2, 6.6, 6.6, -6.6, -6.6 }; //[cm]
  287. X = new double[160];
  288. Y = new double[160];
  289. size_t index = -1; // [0,159]
  290. for(size_t side=0;side<2;side++) // cycle by sides
  291. {
  292. for(size_t k=0;k<20;k++) // cycle by detectors
  293. {
  294. calcCenters(posX[k] +1.4, posY[k] +1.4, ++index, "\n-------------------------");
  295. calcCenters(posX[k] -1.4, posY[k] +1.4, ++index);
  296. calcCenters(posX[k] -1.4, posY[k] -1.4, ++index);
  297. calcCenters(posX[k] +1.4, posY[k] -1.4, ++index);
  298. }
  299. }
  300. } // init by first time
  301. size_t index = suid - 1; // [1,161] -> [0,159]
  302. double Z = (index < 80) ? 140.7+2.55 : -140.7-2.55;
  303. return TVector3(X[index], Y[index], Z);
  304. }
  305. //------------------------------------------------------------------------------------------------------------------------
  306. void MpdFfdHitProducer::Finish()
  307. {
  308. if(fFlnm.Length() > 0)
  309. {
  310. LOG(DEBUG2)<<"[MpdFfdHitProducer::Finish] Update "<<fFlnm.Data()<<" file. ";
  311. auto ptr = gFile;
  312. TFile file(fFlnm.Data(), "RECREATE");
  313. fList.Write();
  314. file.Close();
  315. gFile = ptr;
  316. }
  317. }
  318. //------------------------------------------------------------------------------------------------------------------------
  319. void MpdFfdHitProducer::SetParContainers()
  320. {
  321. // Get run and runtime database
  322. FairRunAna* run = FairRunAna::Instance();
  323. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  324. FairRuntimeDb* db = run->GetRuntimeDb();
  325. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  326. // Get √ geometry parameter container
  327. db->getContainer("MpdFfdGeoPar");
  328. }
  329. //------------------------------------------------------------------------------------------------------------------------
  330. bool MpdFfdHitProducer::IsPeCreated(double energy) // [eV]
  331. {
  332. if(gPMTeff_pdf == nullptr) // init by first fime
  333. {
  334. const size_t n = 11;
  335. const double length[n] =
  336. { 155., 170., 200., 250., 300., 350., 400., 450., 500., 550., 600.}; // op wavelength [nm]
  337. const double eff[n] =
  338. { 0.20, 0.20, 0.18, 0.16, 0.18, 0.21, 0.22, 0.16, 0.11, 0.035, 0.015}; // pe creating efficiency
  339. gPMTeff_pdf = new TGraph(n, length, eff);
  340. }
  341. if(energy < 2. || energy > 8.) return false; // [2,8] eV fast cut
  342. double wavelength = 1239.8/energy;
  343. double efficiency = gPMTeff_pdf->Eval(wavelength);
  344. return (gRandom->Rndm() < efficiency);
  345. }
  346. //------------------------------------------------------------------------------------------------------------------------
  347. MpdFfdHit* MpdFfdHitProducer::AddHit(const MpdFfdPoint *point, const TVector3& pos, const TVector3& dpos, Int_t refIndex, size_t npe, Int_t flag)
  348. {
  349. MpdFfdHit *pHit = new ((*aFfdHits)[aFfdHits->GetEntriesFast()]) MpdFfdHit(point->GetDetectorID(), pos, dpos, refIndex, point->GetTime(), npe, flag);
  350. pHit->AddLink(FairLink(1, point->GetTrackID())); // LS: key value = 1 for mc track index
  351. return pHit;
  352. }
  353. //------------------------------------------------------------------------------------------------------------------------
  354. MpdFfdHit* MpdFfdHitProducer::AddHit(const MpdFfdPoint *point, const TVector3& pos, const TVector3& dpos, Int_t refIndex, size_t npe, Int_t flag, const Tintegrals& map)
  355. {
  356. MpdFfdHit *pHit = new ((*aFfdHits)[aFfdHits->GetEntriesFast()]) MpdFfdHit(point->GetDetectorID(), pos, dpos, refIndex, point->GetTime(), npe, flag);
  357. for(const auto& entry : map)
  358. {
  359. pHit->AddLink(FairLink(1, entry.second.pPoint->GetTrackID())); // LS: key value = 1 for mc track index
  360. pHit->AddLink(FairLink(2, entry.second.pointIndex)); // LS: key value = 2 for mc point index
  361. }
  362. return pHit;
  363. }
  364. //------------------------------------------------------------------------------------------------------------------------