MpdEtofMatching.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdEtofMatching
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <iostream>
  8. #include <fstream>
  9. #include <map>
  10. #include <TMath.h>
  11. #include <TFile.h>
  12. #include <TRandom2.h>
  13. #include <TStorage.h>
  14. #include <TEfficiency.h>
  15. #include "FairRootManager.h"
  16. #include "MpdMCTrack.h"
  17. #include "FairRunAna.h"
  18. #include "FairLogger.h"
  19. #include "MpdEctKalmanTrack.h"
  20. #include "MpdKalmanFilter.h"
  21. #include "MpdTofPoint.h"
  22. #include "MpdTofHit.h"
  23. #include "MpdTofMatching.h"
  24. #include "MpdTofGeoUtils.h"
  25. #include "MpdEtofGeoUtils.h"
  26. #include "MpdEtof.h"
  27. //#include "LMatchingFilter.h"
  28. #include "MpdTofMatchingQA.h"
  29. #include "MpdEtofMatching.h"
  30. using namespace std;
  31. /*
  32. struct less_by_pointer
  33. {
  34. inline bool operator() (const intervalType& struct1, const intervalType& struct2)
  35. {
  36. return (struct1.value < struct2.value);
  37. }
  38. };
  39. */
  40. ClassImp(MpdEtofMatching)
  41. //------------------------------------------------------------------------------------------------------------------------
  42. MpdEtofMatching::MpdEtofMatching(const char *name, Int_t verbose, Bool_t test, const char *flnm)
  43. : FairTask(name, verbose),
  44. aMcPoints(nullptr), aMcTracks(nullptr), aTofHits(nullptr), aKFectTracks(nullptr), aTofMatchings(nullptr),
  45. fMode(kIntervalTree), fDoTest(test), fUseMCData(false), pRandom(new TRandom2),
  46. fNSmeared(20), fTofZpos(295.2), fTofRmax(140.), fThreshR(15.), fThreshTheta(1.5), //Z- 250.0
  47. pKF(nullptr)
  48. {
  49. pMatchingQA = fDoTest ? new MpdTofMatchingQA(flnm, true) : nullptr;
  50. //pMF = new LMatchingFilter(pMatchingQA, fVerbose);
  51. }
  52. //------------------------------------------------------------------------------------------------------------------------
  53. MpdEtofMatching::~MpdEtofMatching()
  54. {
  55. // delete pMF;
  56. delete pRandom;
  57. delete pMatchingQA;
  58. if (aTofMatchings != nullptr)
  59. {
  60. aTofMatchings->Delete();
  61. delete aTofMatchings;
  62. }
  63. }
  64. //------------------------------------------------------------------------------------------------------------------------
  65. InitStatus MpdEtofMatching::Init()
  66. {
  67. aMcPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("ETOFPoint");
  68. aMcTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
  69. aTofHits = (TClonesArray*) FairRootManager::Instance()->GetObject("ETOFHit");
  70. aKFectTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("EctTrack");
  71. if(aMcPoints && aMcTracks) fUseMCData = true;
  72. if(!aTofHits || !aKFectTracks){ LOG(FATAL)<<"[MpdEtofMatching::Init] Branch not found!"; return kERROR; }
  73. pKF = MpdKalmanFilter::Instance("KF","KF");
  74. // Create and register output array
  75. aTofMatchings = new TClonesArray("MpdTofMatchingData");
  76. FairRootManager::Instance()->Register("ETOFMatching", "ETof", aTofMatchings, kTRUE);
  77. /// pMF->SetContainer(aTofMatchings);
  78. // MpdEtofGeoUtils::Instance()->ParseTGeoManager(fUseMCData, nullptr, false);
  79. // MpdEtofGeoUtils::Instance()->FindNeighborStrips(0.8, nullptr, nullptr, false);// 0.8 [cm] <--- thresh. distance between neighbor strips, (see h1TestDistance histo)
  80. LOG(INFO)<<"MpdEtofMatching initialization finished succesfully.";
  81. return kSUCCESS;
  82. }
  83. //------------------------------------------------------------------------------------------------------------------------
  84. void MpdEtofMatching::Exec(Option_t *option)
  85. {
  86. fDoMCTest = fUseMCData && fDoTest;
  87. /*
  88. // Reset event
  89. aTofMatchings->Clear();
  90. pMF->Reset();
  91. Int_t nTofPoints = -1, nMCTracks = -1;
  92. if(fUseMCData){ nTofPoints = aMcPoints->GetEntriesFast(); nMCTracks = aMcTracks->GetEntriesFast();}
  93. Int_t nTofHits = aTofHits->GetEntriesFast();
  94. Int_t nKFTracks = aKFectTracks->GetEntriesFast();
  95. if(fVerbose) cout<<"\n -I- [MpdEtofMatching::Exec] points= "<<nTofPoints<<", hits= "<<nTofHits<<", mc tracks= "<<nMCTracks<<", kf tracks= "<<nKFTracks<<endl;
  96. // ---------------------------------------------------------------------------------------->>> Sorting & Mapping points to MC tracks
  97. TmmP2T mmMCpoints; // pair< MCtrackID, MpdTofPoint*>
  98. TmmP2T::const_iterator mmMCpointCiter;
  99. if(fUseMCData) MpdTofMatching::MappingMcPoints2McTracks(aMcPoints, mmMCpoints);
  100. // ---------------------------------------------------------------------------------------->>> Mapping hits to stripUIDs
  101. typedef std::map<Int_t, std::pair<MpdTofHit*, Int_t> > TmS2H; // pair< stripUID, pair<MpdTofHit*, hitIndex> >
  102. TmS2H mHits;
  103. for(Int_t hitIndex = 0; hitIndex < nTofHits; hitIndex++) // cycle by tof hits
  104. {
  105. MpdTofHit *pTofHit = (MpdTofHit*) aTofHits->At(hitIndex);
  106. mHits.insert(make_pair( pTofHit->GetDetectorID(), make_pair(pTofHit, hitIndex)));
  107. }
  108. // ---------------------------------------------------------------------------------------->>> Sorting by Pt the KfEctTracks
  109. TsPt tids;
  110. for(Int_t KfIndex = 0; KfIndex < nKFTracks; KfIndex++) // cycle by ECT KF tracks
  111. {
  112. MpdEctKalmanTrack *pKfTrack = (MpdEctKalmanTrack*) aKFectTracks->UncheckedAt(KfIndex);
  113. tids.insert(make_pair(pKfTrack, KfIndex));
  114. }
  115. // ----------------------------------------------------------------------------------------
  116. const intervalTreeType* mDetectorsR = MpdEtofGeoUtils::Instance()->GetDetR();
  117. const intervalTreeType* mDetectorsPhi = MpdEtofGeoUtils::Instance()->GetDetPhi();
  118. vector<intervalType> segmentR, segmentPhi, intersect;
  119. TVector3 hitPosition, Momentum, estPoint;
  120. bool IsInside;
  121. Double_t trackLength, thR, thPhi;
  122. Int_t charge, selectedTracks = 0;
  123. mcInfo MCdata; // The MC run variables. VALID ONLY if fUseMCData = true
  124. //cout<<"\n ------------------------------------------------------------------------------------------------------------->> EVENT";
  125. //mDetectorsR->dump("\n\n ----->>> mDetectorsR INTERVALS");
  126. //mDetectorsPhi->dump("\n\n ----->>> mDetectorsPhi INTERVALS");
  127. TmS2H::const_iterator hitCiter;
  128. for(TsPt::const_iterator it = tids.begin(), itEnd = tids.end(); it != itEnd; it++)// cycle by sorted ECT KF tracks
  129. {
  130. const MpdEctKalmanTrack *pKfTrack = (MpdEctKalmanTrack*) it->first;
  131. Int_t KfIndex = it->second;
  132. if(fUseMCData) MpdTofMatching::GetMcInfo(MCdata, mmMCpoints, pKfTrack, aMcTracks);
  133. size_t TofSide = (pKfTrack->Momentum3().Pz() > 0.) ? 0 : 1;
  134. double TofZpos = (TofSide == 0) ? fTofZpos : -fTofZpos;
  135. if(EstTrackOnPlane(pKfTrack, TofZpos, estPoint, trackLength, Momentum, charge))
  136. {
  137. IsInside = (estPoint.Perp() < fTofRmax) ? true : false;
  138. if(fDoMCTest) pMatchingQA->FilleRest(estPoint.Perp());
  139. }
  140. if(!IsInside) continue; // KF track out ETof Rmax
  141. if(fDoMCTest)
  142. {
  143. if(MCdata.Npoints == 1) // only SINGLE tof point per MCtrack
  144. pMatchingQA->FillDevPoint2EstCyl(MCdata.Position, estPoint, pKfTrack->Momentum());
  145. if(MCdata.TofTouch) // KF track have TOFpoint
  146. {
  147. selectedTracks++;
  148. pMatchingQA->FillSelectedKfTracks(pKfTrack->Momentum3().Eta(), pKfTrack->Momentum());
  149. }
  150. }
  151. // ---------------------------------------------------------------------------------------->>> Looking for overlaping of estPoint & detectors
  152. double estR = estPoint.Perp(), estPhi = estPoint.Phi();
  153. const static double Rerror = 10.; // [cm] // FIXME: should be dependent on the parameters of the KFtrack
  154. const static double PhiError = 0.1; // [rads] // FIXME: should be dependent on the parameters of the KFtrack
  155. segmentR.clear();
  156. segmentPhi.clear();
  157. intersect.clear();
  158. mDetectorsR[TofSide].findOverlapping(estR - Rerror, estR + Rerror, segmentR);
  159. mDetectorsPhi[TofSide].findOverlapping(estPhi - PhiError, estPhi + PhiError, segmentPhi);
  160. //cout<<"\n <<<--R-->>>>> ("<<estR-Rerror<<", "<<estR+Rerror<<")";
  161. //MpdTof::intervalTreeType::dumpIntervals(&segmentR, "\n\n ----->>> mDetectorsR findOverlapping");
  162. //cout<<"\n <<<--Phi-->>>>> ("<<estPhi-PhiError<<", "<<estPhi+PhiError<<")";
  163. //MpdTof::intervalTreeType::dumpIntervals(&segmentPhi, "\n\n ----->>> mDetectorsPhi findOverlapping");
  164. if(!segmentR.empty() && !segmentPhi.empty()) // have overlaped segments both R and Phi
  165. {
  166. // calc. intersection
  167. sort(segmentR.begin(), segmentR.end(), less_by_pointer());
  168. sort(segmentPhi.begin(), segmentPhi.end(), less_by_pointer());
  169. set_intersection(segmentR.begin(), segmentR.end(), segmentPhi.begin(), segmentPhi.end(), std::back_inserter(intersect), less_by_pointer());
  170. for(vector<intervalType>::const_iterator cit = intersect.begin(), citEnd = intersect.end(); cit != citEnd; cit++) // cycle by the overlaped strips
  171. {
  172. Int_t detUID = (*cit).value->volumeUID;
  173. if(fDoMCTest && MCdata.Npoints == 1) // only one tof point per MCtrack
  174. pMatchingQA->FillDevPoint2EstP(MCdata.Position, estPoint, pKfTrack->Momentum());
  175. hitCiter = mHits.find(detUID);
  176. if(hitCiter != mHits.end()) // the estimated strip have hit
  177. {
  178. MpdTofHit *TofHit = hitCiter->second.first;
  179. Int_t hitIndex = hitCiter->second.second;
  180. TofHit->Position(hitPosition);
  181. double devHit, devHitZ, devHitR, devHitPhi;
  182. MpdTofMatching::GetDelta(hitPosition, estPoint, devHit, devHitZ, devHitR, devHitPhi);
  183. bool mcIsSameVolumeUID = (TofHit->GetDetectorID() == MCdata.stripUID);
  184. if(fDoMCTest) pMatchingQA->FillHitDev2EstP(mcIsSameVolumeUID, devHitR, devHitPhi);
  185. if(devHitR < fThreshR && devHitPhi < fThreshTheta) // inside matching window
  186. {
  187. pMF->AddCandidate(MpdTofMatchingData(KfIndex, hitIndex, pKfTrack->GetNofTrHits(), TofHit, MCdata.pid, TofHit->GetFlag(), trackLength, estPoint, Momentum, charge, devHitR, devHitPhi));
  188. //cout<<"\n AddCandidate KfIndex="<<KfIndex<<" hitIndex="<<hitIndex<<" delta="<<delta<<" NmbTrHits="<<pKfTrack->GetNofTrHits()<<" mcHasCand="<<mcHasCand;
  189. }
  190. }
  191. } // cycle by the overlaped strips
  192. } // have overlaped segments both R and Phi
  193. if(fDoMCTest && MCdata.TofTouch) pMatchingQA->FillCandidates(MCdata.HaveTrueCand, MCdata.HaveCand, pKfTrack->Momentum3().Eta(), pKfTrack->Momentum());
  194. } // cycle by KF tracks
  195. if(fDoTest) pMatchingQA->FillTrackPerEvent(nKFTracks, selectedTracks);
  196. double chi2;
  197. Int_t MatchingOK = pMF->Processing(nKFTracks, chi2); // accept candidates
  198. Int_t nEntries = pMF->UpdateContainer(); // remove unmatched candidates
  199. assert(nEntries == MatchingOK);
  200. if(fDoMCTest) pMatchingQA->FillMatchingEfficiency(aTofMatchings, aTofHits, aKFectTracks, aMcTracks);
  201. if(fVerbose)
  202. {
  203. cout<<"\n -I- [MpdEtofMatching::Exec] MatchingOK = "<<MatchingOK;
  204. if(MatchingOK != 0) cout<<" ("<<chi2/MatchingOK<<")";
  205. }
  206. */
  207. }
  208. //------------------------------------------------------------------------------------------------------------------------
  209. void MpdEtofMatching::Finish()
  210. {
  211. if(pMatchingQA) pMatchingQA->Finish();
  212. }
  213. //------------------------------------------------------------------------------------------------------------------------
  214. bool MpdEtofMatching::EstTrackOnPlane(const MpdEctKalmanTrack *tr, Double_t Zetof, TVector3& pos, Double_t& length, TVector3& Mom, Int_t& charge) const
  215. {
  216. MpdEctKalmanTrack tr1(*tr);
  217. TObjArray *hits = tr1.GetHits();
  218. if (hits->GetEntriesFast() == 0)
  219. {
  220. Int_t nh = tr1.GetNofTrHits();
  221. for (Int_t j = 0; j < nh; ++j) hits->Add(tr1.GetTrHits()->UncheckedAt(j));
  222. tr1.SetNofHits(nh);
  223. }
  224. MpdKalmanHit *h = (MpdKalmanHit*) tr1.GetTrHits()->First();
  225. if (h->GetType() == MpdKalmanHit::kFixedZ) tr1.SetType(MpdKalmanTrack::kEndcap);
  226. tr1.SetPos(tr1.GetPosAtHit());
  227. tr1.SetPosNew(tr1.GetPos());
  228. tr1.SetParamNew(*tr1.GetParamAtHit());
  229. tr1.SetParam(*tr1.GetParamAtHit());
  230. tr1.SetWeight(*tr1.GetWeightAtHit());
  231. tr1.SetLength(tr1.GetLengAtHit());
  232. tr1.SetDirection(MpdKalmanTrack::kOutward);
  233. MpdKalmanHit hEnd;
  234. hEnd.SetType(MpdKalmanHit::kFixedZ);
  235. hEnd.SetPos(Zetof); // eTOF Z position, [cm]
  236. if (!pKF->PropagateParamZ(&tr1, &hEnd, kTRUE)) return false;
  237. Double_t Z = tr1.GetPosNew();
  238. if(TMath::IsNaN(Z)) return false;
  239. Double_t Pt_inv = tr1.GetParamNew(4); // 1/Pt
  240. if(Pt_inv == 0.) return false;
  241. length = tr1.GetLength();
  242. Mom = tr1.Momentum3();
  243. charge = tr1.Charge();
  244. pos.SetXYZ(tr1.GetParamNew(0), tr1.GetParamNew(1), Z);
  245. return true;
  246. }
  247. //------------------------------------------------------------------------------------------------------------------------