MpdEtofHitProducer.cxx 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdEtofHitProducer
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <iostream>
  8. #include <fstream>
  9. #include <valarray>
  10. #include <assert.h>
  11. #include "TClonesArray.h"
  12. #include "TRandom.h"
  13. #include "TVector3.h"
  14. #include "TStopwatch.h"
  15. #include "TMath.h"
  16. #include "FairRunAna.h"
  17. #include "FairBaseParSet.h"
  18. #include "FairRuntimeDb.h"
  19. #include "FairMCApplication.h"
  20. #include "FairDetector.h"
  21. #include "FairRootManager.h"
  22. #include "MpdMCTrack.h"
  23. #include "FairGeoVolume.h"
  24. #include "MpdEtofGeoPar.h"
  25. #include "MpdTofPoint.h"
  26. #include "MpdTofHit.h"
  27. #include "MpdTofHitProducerQA.h"
  28. #include "MpdEtof.h"
  29. #include "MpdEtofGeoUtils.h"
  30. #include "MpdEtofHitProducer.h"
  31. using namespace std;
  32. ClassImp(MpdEtofHitProducer)
  33. //------------------------------------------------------------------------------------------------------------------------
  34. MpdEtofHitProducer::MpdEtofHitProducer(const char *name, Bool_t useMCdata, Int_t verbose, Bool_t test, const char *flnm)
  35. : MpdEtofHitProducerIdeal(name, useMCdata, verbose, test, true, flnm), fTimeSigma(0.100), fErrPhi(0.5), fErrR(1./sqrt(12.)), pRandom(new TRandom2)
  36. {
  37. }
  38. //------------------------------------------------------------------------------------------------------------------------
  39. MpdEtofHitProducer::~MpdEtofHitProducer()
  40. {
  41. delete pRandom;
  42. }
  43. //------------------------------------------------------------------------------------------------------------------------
  44. InitStatus MpdEtofHitProducer::Init()
  45. {
  46. // if(fOnlyPrimary) cout<<" Only primary particles are processed!!! \n"; // FIXME NOT used now ADDD
  47. if(fUseMCData)
  48. {
  49. aMcPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("ETOFPoint");
  50. aMcTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
  51. assert(aMcPoints);
  52. assert(aMcTracks);
  53. }
  54. else
  55. {
  56. aExpDigits = (TClonesArray*) FairRootManager::Instance()->GetObject("??????");// FIXME: NOW unknown name
  57. assert(aExpDigits);
  58. }
  59. // Create and register output array
  60. aTofHits = new TClonesArray("MpdTofHit");
  61. FairRootManager::Instance()->Register("ETOFHit", "ETof", aTofHits, kTRUE);
  62. MpdEtofGeoUtils::Instance()->ParseTGeoManager(true); // forced
  63. MpdEtofGeoUtils::Instance()->FindNeighborStrips(0.8, // 0.8 [cm] <--- thresh. distance between neighbor strips
  64. true); // forced
  65. LOG(INFO)<<"MpdTofHitProducer initialization finished succesfully.";
  66. return kSUCCESS;
  67. }
  68. //------------------------------------------------------------------------------------------------------------------------
  69. Bool_t MpdEtofHitProducer::HitExist(Double_t val) // val - rasstojanie do kraja pad
  70. {
  71. constexpr Double_t slope = (0.98 - 0.95)/0.2;
  72. Double_t efficiency = (val > 0.2) ? 0.98 : ( 0.95 + slope*val);
  73. //-------------------------------------
  74. // 99% ---------
  75. // \.
  76. // \.
  77. // \.
  78. // 95% \.
  79. // <-----------|--|
  80. // 0.2 0.
  81. //-------------------------------------
  82. if(pRandom->Rndm() < efficiency) return true;
  83. return false;
  84. }
  85. //------------------------------------------------------------------------------------------------------------------------
  86. Bool_t MpdEtofHitProducer::DoubleHitExist(Double_t val) // val - rasstojanie do kraja pad
  87. {
  88. constexpr Double_t slope = (0.3 - 0.0)/0.5;
  89. Double_t efficiency = (val > 0.5) ? 0. : (0.3 - slope*val);
  90. //-------------------------------------
  91. // 30% /
  92. // /
  93. // /
  94. // /
  95. // 0% /
  96. // <-----------|----|
  97. // 0.5 0.
  98. //-------------------------------------
  99. if(efficiency == 0.) return false;
  100. if(pRandom->Rndm() < efficiency) return HitExist(val);
  101. return false;
  102. }
  103. //------------------------------------------------------------------------------------------------------------------------
  104. void MpdEtofHitProducer::Exec(Option_t *option)
  105. {
  106. /* static const TVector3 XYZ_err(0.1, 0.1, 0.1); // FIXME: dummy now, MUST BE ROTATED!!!!
  107. aTofHits->Clear();
  108. Int_t UID, trackID;
  109. TVector3 pos, XYZ_smeared;
  110. int nSingleHits = 0, nDoubleHits = 0;
  111. if(fUseMCData)
  112. {
  113. for(Int_t pointIndex = 0, nTofPoint = aMcPoints->GetEntriesFast(); pointIndex < nTofPoint; pointIndex++ ) // cycle by TOF points
  114. {
  115. MpdTofPoint *pPoint = (MpdTofPoint*) aMcPoints->UncheckedAt(pointIndex);
  116. if(fVerbose > 2) pPoint->Print("");
  117. trackID = pPoint->GetTrackID();
  118. UID = pPoint->GetDetectorID();
  119. Double_t time = pRandom->Gaus(pPoint->GetTime(), fTimeSigma); // 100 ps
  120. pPoint->Position(pos);
  121. // Calc Rerr, PhiErr onto MRF
  122. const LStrip *pStrip = MpdEtofGeoUtils::Instance()->FindStrip(UID);
  123. double centerPhi = pStrip->center.Phi(); // [rad]
  124. TVector3 posRotated(pos);
  125. posRotated.RotateZ(-centerPhi);
  126. XYZ_smeared.SetXYZ(posRotated.X(), pRandom->Gaus(posRotated.Y(), fErrPhi), posRotated.Z());
  127. XYZ_smeared.RotateZ(centerPhi);
  128. LStrip::Side_t side;
  129. Double_t distance = pStrip->MinDistanceToEdge(&pos, side); // [cm]
  130. bool passed;
  131. if(passed = HitExist(distance)) // check efficiency
  132. {
  133. AddHit(UID, XYZ_smeared, XYZ_err, pointIndex, trackID, time, MpdTofUtils::IsSingle);
  134. nSingleHits++;
  135. if(pHitProducerQA) pHitProducerQA->FillSingleHitPosition(pos, XYZ_smeared);
  136. }
  137. if(pQA) pQA->GetSingleHitEfficiency()->Fill(passed, distance);
  138. if(passed = DoubleHitExist(distance)) // check cross hit
  139. {
  140. Int_t CrossUID = (side == LStrip::kRight) ? pStrip->neighboring[LStrip::kRight] : pStrip->neighboring[LStrip::kLeft];
  141. if(LStrip::kInvalid == CrossUID) continue; // last strip on module
  142. // Calc Rerr, PhiErr onto MRF
  143. pStrip = MpdEtofGeoUtils::Instance()->FindStrip(CrossUID);
  144. // XYZ_smeared.SetXYZ(pRandom->Gaus(posRotated.X(), fErrX), posRotated.Y(), pStrip->center.Z());
  145. // XYZ_smeared.RotateZ(perpPhi);
  146. AddHit(CrossUID, XYZ_smeared, XYZ_err, pointIndex, trackID, time, MpdTofUtils::IsDouble);
  147. nDoubleHits++;
  148. if(pHitProducerQA) pHitProducerQA->FillDoubleHitPosition(pos, XYZ_smeared);
  149. }
  150. if(pQA) pQA->GetDoubleHitEfficiency()->Fill(passed, distance);
  151. } // cycle by the TOF points
  152. }
  153. else // exp. data used
  154. {
  155. // FIXME: now not realized
  156. //AddHit(Int_t detUID, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
  157. assert(false);
  158. }
  159. MergeHitsOnStrip(); // save only the fastest hit in the strip
  160. int nFinally = CompressHits(); // remove blank slotes
  161. cout<<" -I- [MpdEtofHitProducer::Exec] single hits= "<<nSingleHits<<", double hits= "<<nDoubleHits<<", final hits= "<<nFinally<<endl;
  162. */
  163. }
  164. //------------------------------------------------------------------------------------------------------------------------
  165. /*void MpdEtofHitProducer::Finish()
  166. {
  167. if(pQA) pQA->Finish();
  168. }*/
  169. //------------------------------------------------------------------------------------------------------------------------
  170. void MpdEtofHitProducer::SetSeed(UInt_t seed)
  171. {
  172. pRandom->SetSeed(seed);
  173. }
  174. //------------------------------------------------------------------------------------------------------------------------