MpdTofHitProducer.cxx 8.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdTofHitProducer
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <assert.h>
  8. #include <TMath.h>
  9. #include <TH1D.h>
  10. #include <TH2D.h>
  11. #include <TEfficiency.h>
  12. #include <TRandom2.h>
  13. #include <TClonesArray.h>
  14. #include "MpdMCTrack.h"
  15. #include "FairLogger.h"
  16. #include "MpdTofUtils.h"
  17. #include "MpdTofPoint.h"
  18. #include "MpdTofHit.h"
  19. #include "MpdTofGeoUtils.h"
  20. #include "MpdTofHitProducerQA.h"
  21. #include "MpdTofHitProducer.h"
  22. using namespace std;
  23. ClassImp(MpdTofHitProducer)
  24. //------------------------------------------------------------------------------------------------------------------------
  25. MpdTofHitProducer::MpdTofHitProducer(const char *name, Bool_t useMCdata, Int_t verbose, Bool_t test, const char *flnm, double stripLength)
  26. : MpdTofHitProducerIdeal(name, useMCdata, verbose, test, true, flnm, false), fStripLength(stripLength)
  27. {
  28. pRandom = new TRandom2;
  29. }
  30. //------------------------------------------------------------------------------------------------------------------------
  31. MpdTofHitProducer::~MpdTofHitProducer()
  32. {
  33. delete pRandom;
  34. }
  35. //------------------------------------------------------------------------------------------------------------------------
  36. void MpdTofHitProducer::AddParameters(TString& buf)const
  37. {
  38. MpdTofHitProducerIdeal::AddParameters(buf);
  39. buf += TString::Format(", fTimeSigma=%.5g ns", fTimeSigma);
  40. buf += TString::Format(", fErrX=%.4g cm", fErrX);
  41. buf += TString::Format(", fErrZ=%.4g cm", fErrZ);
  42. }
  43. //------------------------------------------------------------------------------------------------------------------------
  44. InitStatus MpdTofHitProducer::Init()
  45. {
  46. InitStatus status = MpdTofHitProducerIdeal::Init();
  47. MpdTofGeoUtils::Instance()->ParseTGeoManager(pQA, true);//, "geometryTree.txt"); // forced
  48. MpdTofGeoUtils::Instance()->FindNeighborStrips(0.8, pQA, // 0.8 [cm] <--- thresh. distance between neighbor strips
  49. true); // forced
  50. LOG(INFO)<<"[MpdTofHitProducer::Init] Initialization finished succesfully.";
  51. return status;
  52. }
  53. //------------------------------------------------------------------------------------------------------------------------
  54. Bool_t MpdTofHitProducer::IsHitCreated(Double_t value, Int_t gap) // value - distance to the strip edge [cm]
  55. {
  56. constexpr Double_t slope = (0.98 - 0.95)/0.2;
  57. Double_t efficiency = (value > 0.2) ? 0.98 : ( 0.95 + slope*value);
  58. //-------------------------------------
  59. // 99% ---------
  60. // \.
  61. // \.
  62. // \.
  63. // 95% \.
  64. // <-----------|--|
  65. // 0.2 0. value
  66. //-------------------------------------
  67. if(gap == 1 || gap == 3) efficiency /= 2.; // reduce efficiency on outer strip gap
  68. if(pRandom->Rndm() < efficiency) return true;
  69. return false;
  70. }
  71. //------------------------------------------------------------------------------------------------------------------------
  72. Bool_t MpdTofHitProducer::IsCrossHitCreated(Double_t value, Int_t gap) // value - distance to the strip edge [cm]
  73. {
  74. constexpr Double_t slope = (0.3 - 0.0)/0.5;
  75. Double_t efficiency = (value > 0.5) ? 0. : (0.3 - slope*value);
  76. //-------------------------------------
  77. // 30% /
  78. // /
  79. // /
  80. // /
  81. // 0% /
  82. // <-----------|----|
  83. // 0.5 0. value
  84. //-------------------------------------
  85. if(efficiency == 0.) return false;
  86. if(gap == 1 || gap == 3) efficiency /= 2.; // reduce efficiency on outer strip gap
  87. if(pRandom->Rndm() < efficiency) return true;
  88. return false;
  89. }
  90. //------------------------------------------------------------------------------------------------------------------------
  91. void MpdTofHitProducer::Exec(Option_t *option)
  92. {
  93. aTofHits->Clear();
  94. TVector3 mcPosition, hitPosition, hitPosError;
  95. int nSingleHits = 0, nCrossHits = 0;
  96. auto SmearingAlongStrip = [this](const TGeoCombiTrans& gap, const TGeoCombiTrans& centralGap, const TVector3& point, TVector3& hit, TVector3& error)
  97. {
  98. // rotate stip to origin LRS (dX = fStripLength, dY = 0, dZ = strip step)
  99. Double_t local[3], master[3] = {point.X(), point.Y(), point.Z()};
  100. gap.MasterToLocal(master, local);
  101. // smearing the hit position along x axis
  102. double Xsmeared, x = local[0]; size_t nTries = 0;
  103. do
  104. {
  105. Xsmeared = pRandom->Gaus(x, fErrX);
  106. if(++nTries > 100) return false;
  107. }
  108. while(fabs(Xsmeared) > fStripLength/2.); // check strip boundary exit
  109. // shift and rotate back to MRS(shift by Y from gap to central gap)
  110. local[0] = Xsmeared, local[2] = 0.; // smearing along X axis, shift to the strip center along Z axis(at LRS)
  111. centralGap.LocalToMaster(local, master);
  112. hit.SetXYZ(master[0], master[1], master[2]);
  113. // rotate error vector
  114. local[0] = fErrX; local[1] = 0.; local[2] = fErrZ;
  115. centralGap.LocalToMaster(local, master);
  116. error.SetXYZ(master[0], master[1], master[2]);
  117. return true;
  118. };
  119. if(fUseMCData)
  120. {
  121. for(Int_t pointIndex = 0, nTofPoint = aMcPoints->GetEntriesFast(); pointIndex < nTofPoint; pointIndex++ ) // cycle by TOF points
  122. {
  123. auto pPoint = (MpdTofPoint*) aMcPoints->UncheckedAt(pointIndex);
  124. if(fVerbose > 2) pPoint->Print("");
  125. Int_t tid = pPoint->GetTrackID();
  126. Int_t suid = pPoint->GetDetectorID();
  127. Int_t gap = pPoint->GetGap();
  128. Double_t time = pRandom->Gaus(pPoint->GetTime(), fTimeSigma); // default 100 ps
  129. pPoint->Position(mcPosition);
  130. auto strip = MpdTofGeoUtils::Instance()->FindStrip(suid);
  131. auto centralStrip = MpdTofGeoUtils::Instance()->FindStrip(MpdTofPoint::SetCentralGap(suid));
  132. if(! SmearingAlongStrip(strip->fMatrix, centralStrip->fMatrix, mcPosition, hitPosition, hitPosError))
  133. {
  134. strip->Dump(" [MpdTofHitProducer::Exec] -E- Invalid Rotation matrix.");
  135. continue;
  136. }
  137. LStrip::Side_t side;
  138. Double_t distance = strip->MinDistanceToEdge(&mcPosition, side); // [cm]
  139. bool stripFired;
  140. if( (stripFired = IsHitCreated(distance, gap)) ) // simulate hit efficiency (add flag MpdTofUtils::IsSingle)
  141. {
  142. AddHit(MpdTofPoint::ClearGap(suid), hitPosition, hitPosError, pointIndex, tid, time, MpdTofUtils::IsSingle);
  143. nSingleHits++;
  144. }
  145. if(pQA)
  146. {
  147. if(stripFired) pQA->Point2HitSmearingTest(mcPosition, hitPosition);
  148. pQA->HitGapEfficiencyTest(stripFired, distance, gap);
  149. pQA->PositionInsideStripTest(strip->center, mcPosition, hitPosition);
  150. pQA->RotationToOriginTest(strip->fMatrix, mcPosition, hitPosition);
  151. pQA->CentralDetectorTest(suid, hitPosition);
  152. }
  153. if( (stripFired = IsCrossHitCreated(distance, gap)) ) // simulate cross hit efficiency (add flag MpdTofUtils::IsDouble)
  154. {
  155. Int_t crossSuid = (side == LStrip::kRight) ? strip->neighboring[LStrip::kRight] : strip->neighboring[LStrip::kLeft];
  156. if(LStrip::kInvalid == crossSuid) continue; // edge strip on detector
  157. strip = MpdTofGeoUtils::Instance()->FindStrip(crossSuid);
  158. centralStrip = MpdTofGeoUtils::Instance()->FindStrip(MpdTofPoint::SetCentralGap(crossSuid));
  159. if(! SmearingAlongStrip(strip->fMatrix, centralStrip->fMatrix, mcPosition, hitPosition, hitPosError))
  160. {
  161. strip->Dump(" [MpdTofHitProducer::Exec] -E- Invalid Rotation matrix.");
  162. continue;
  163. }
  164. AddHit(MpdTofPoint::ClearGap(crossSuid), hitPosition, hitPosError, pointIndex, tid, time, MpdTofUtils::IsDouble);
  165. nCrossHits++;
  166. }
  167. if(pQA) pQA->FillCrossHitEfficiency(stripFired, distance, gap, mcPosition, hitPosition);
  168. } // cycle by the TOF points
  169. }
  170. else // exp. data used
  171. {
  172. // FIXME: now not realized
  173. //AddHit(Int_t detUID, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
  174. assert(false);
  175. }
  176. MergeHitsOnStrip(); // save only the fastest hit in the strip
  177. int Nhits = CompressHits(); // remove blank slotes
  178. if(fUseMCData && pQA)
  179. {
  180. pQA->FillNPointsHits(aMcPoints->GetEntriesFast(), Nhits);
  181. pQA->PointDistanceTest(aMcPoints);
  182. }
  183. LOG(DEBUG1)<<"[MpdTofHitProducer::Exec] single hits= "<<nSingleHits<<", cross hits= "<<nCrossHits<<", final hits= "<<Nhits;
  184. }
  185. //------------------------------------------------------------------------------------------------------------------------
  186. void MpdTofHitProducer::Finish()
  187. {
  188. if(pQA) pQA->Finish();
  189. }
  190. //--------------------------------------------------------------------------------------------------------------------------------------
  191. void MpdTofHitProducer::SetSeed(UInt_t seed)
  192. {
  193. pRandom->SetSeed(seed);
  194. }
  195. //------------------------------------------------------------------------------------------------------------------------