MpdTofMatching.cxx 31 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921
  1. //------------------------------------------------------------------------------------------------------------------------
  2. /// \class MpdTofMatching
  3. ///
  4. /// \brief
  5. /// \author Sergei Lobastov (LHE, JINR, Dubna)
  6. //------------------------------------------------------------------------------------------------------------------------
  7. #include <iostream>
  8. #include <fstream>
  9. #include <limits>
  10. #include <TMath.h>
  11. #include <TFile.h>
  12. #include <TRandom2.h>
  13. #include <TRotation.h>
  14. #include <TStorage.h>
  15. #include <TEfficiency.h>
  16. #include "FairRootManager.h"
  17. #include "MpdMCTrack.h"
  18. #include "FairRunAna.h"
  19. #include "FairLogger.h"
  20. #include "MpdEctKalmanTrack.h"
  21. #include "MpdTofPoint.h"
  22. #include "MpdTof.h"
  23. #include "MpdKalmanFilter.h"
  24. #include "MpdTofMatchingQA.h"
  25. #include "MpdTofGeoUtils.h"
  26. #include "MpdTofMatching.h"
  27. #ifdef _OPENMP
  28. #include "omp.h"
  29. #include <sys/time.h>
  30. #endif
  31. using namespace std;
  32. //------------------------------------------------------------------------------------------------------------------------
  33. bool mcInfo::GetPosition(TVector3& pos, Int_t suid)const
  34. {
  35. for(const auto& cit : vPoints) // cycle by vector<MpdTofPoint*> (selected kftrack -> mctrack -> mcpoints)
  36. {
  37. if(suid == cit->GetDetectorID()){ cit->Position(pos); return true; }
  38. }
  39. return false;
  40. }
  41. //------------------------------------------------------------------------------------------------------------------------
  42. pair<Int_t, Int_t> mcInfo::GetClosestPointOnDetector(TVector3& pos, TVector3& mom, const TVector3& from, Int_t duid)const
  43. {
  44. pair<Int_t,Int_t> found = {-1,-1}; // {suid, mctid}
  45. if(vPoints.empty()) return found;
  46. double dev = 1.e+10;
  47. TVector3 point, P;
  48. for(const auto& cit : vPoints) // cycle by vector<MpdTofPoint*> (selected kftrack -> mctrack -> mcpoints)
  49. {
  50. if(MpdTofPoint::IsSameDetector(duid, cit->GetSuid()))
  51. {
  52. cit->Position(point);
  53. double newdev = (point - from).Mag();
  54. if( newdev < dev)
  55. {
  56. dev = newdev;
  57. pos = point;
  58. cit->Momentum(mom);
  59. found.first = cit->GetSuid();
  60. found.second = cit->GetTrackID();
  61. }
  62. }
  63. }
  64. return found;
  65. }
  66. //------------------------------------------------------------------------------------------------------------------------
  67. bool mcInfo::IsSameStrip(Int_t suid)const
  68. {
  69. if(vPoints.empty()) return false;
  70. for(const auto& cit : vPoints) // cycle by vector<MpdTofPoint*> (selected kftrack -> mctrack -> mcpoints)
  71. {
  72. if(MpdTofPoint::IsSameStrip(cit->GetDetectorID(),suid)) return true;
  73. }
  74. return false;
  75. }
  76. //------------------------------------------------------------------------------------------------------------------------
  77. const MpdTofPoint* mcInfo::FindPoint(Int_t tid, Int_t suid)const
  78. {
  79. if(vPoints.empty()) return nullptr;
  80. for(const auto& cit : vPoints) // cycle by vector<MpdTofPoint*> (selected kftrack -> mctrack -> mcpoints)
  81. {
  82. if(cit->GetTrackID() == tid)
  83. {
  84. if(-1 == suid) return cit;
  85. else if(cit->GetDetectorID() == suid) return cit;
  86. }
  87. }
  88. return nullptr;
  89. }
  90. //------------------------------------------------------------------------------------------------------------------------
  91. //------------------------------------------------------------------------------------------------------------------------
  92. //------------------------------------------------------------------------------------------------------------------------
  93. LWeightMatrix::LWeightMatrix(MpdTofMatchingQA *ptr)
  94. : pData(nullptr), pMatchingQA(ptr)
  95. {
  96. }
  97. //------------------------------------------------------------------------------------------------------------------------
  98. LWeightMatrix::~LWeightMatrix()
  99. {
  100. delete []pData;
  101. }
  102. //------------------------------------------------------------------------------------------------------------------------
  103. void LWeightMatrix::SetWeight(size_t trackId, size_t hitId, double value) // [0,...N-1]
  104. {
  105. assert(trackId<fNTracks);
  106. assert(hitId<fNHits);
  107. *(pData + trackId * fNHits + hitId) = value;
  108. }
  109. //------------------------------------------------------------------------------------------------------------------------
  110. void LWeightMatrix::AddWeight(size_t trackId, size_t hitId, double value)
  111. {
  112. assert(trackId<fNTracks);
  113. assert(hitId<fNHits);
  114. *(pData + trackId * fNHits + hitId) += value;
  115. }
  116. //------------------------------------------------------------------------------------------------------------------------
  117. void LWeightMatrix::Normalize(bool copyNormWeight)
  118. {
  119. for(size_t hitId = 0; hitId<fNHits; hitId++)
  120. {
  121. double sum = SumWeightsForHit(hitId);
  122. if(sum > 0.)
  123. {
  124. for(size_t trackId = 0; trackId < fNTracks; trackId++)
  125. {
  126. double norm_weight = GetWeight(trackId, hitId) / sum;
  127. SetWeight(trackId, hitId, norm_weight);
  128. if(copyNormWeight && norm_weight > 0.)
  129. {
  130. auto it = mCandidates.find(hash(trackId, hitId));
  131. assert(it != mCandidates.end());
  132. it->second.SetNormWeight(norm_weight);
  133. }
  134. }
  135. }
  136. }
  137. }
  138. //------------------------------------------------------------------------------------------------------------------------
  139. void LWeightMatrix::Boost1thRankCand(double mult)
  140. {
  141. if(mult == 1.) return; // do nothing
  142. for(size_t trackId = 0; trackId < fNTracks; trackId++)
  143. {
  144. size_t nHits = 0, singleHitIndex; double singleHitWeight;
  145. for(size_t hitId = 0; hitId<fNHits; hitId++)
  146. {
  147. double weight = GetWeight(trackId, hitId);
  148. if(weight > 0.){ nHits++; singleHitIndex = hitId; singleHitWeight = weight; }
  149. if(nHits > 1) break;
  150. }
  151. // boost 1th rank candidate.
  152. if(nHits == 1) SetWeight(trackId, singleHitIndex, singleHitWeight*mult);
  153. }
  154. }
  155. //------------------------------------------------------------------------------------------------------------------------
  156. void LWeightMatrix::Reset(size_t tracks, size_t hits)
  157. {
  158. if(tracks*hits > fNTracks*fNHits) // buffer size increased --> MUST recreate buffer
  159. {
  160. delete []pData;
  161. pData = new double[tracks*hits];
  162. }
  163. fNTracks = tracks;
  164. fNHits = hits;
  165. // cleanup
  166. memset(pData, 0, sizeof(double)*tracks*hits);
  167. mCandidates.clear();
  168. fmT2MC.clear();
  169. }
  170. //------------------------------------------------------------------------------------------------------------------------
  171. double LWeightMatrix::SumWeightsForHit(size_t hitId)const
  172. {
  173. assert(hitId<fNHits);
  174. double sum = 0.;
  175. for(size_t trackId = 0; trackId < fNTracks; trackId++) sum += GetWeight(trackId, hitId);
  176. return sum;
  177. }
  178. //------------------------------------------------------------------------------------------------------------------------
  179. mcInfo* LWeightMatrix::saveMCinfo(Int_t trackId, const mcInfo& mc)
  180. {
  181. auto ret = fmT2MC.insert(make_pair(trackId, mc));
  182. return &(ret.first->second);
  183. }
  184. //------------------------------------------------------------------------------------------------------------------------
  185. const mcInfo* LWeightMatrix::getMCinfo(Int_t trackId)const
  186. {
  187. auto it = fmT2MC.find(trackId);
  188. if(it == fmT2MC.end()) return nullptr;
  189. return &(it->second);
  190. }
  191. //----------------------------------------------------------------------------------------------------------------------------------------
  192. MpdTofMatchingData* LWeightMatrix::AddCandidate(const MpdTofMatchingData& data)
  193. {
  194. const long key = hash(data.GetKFTrackIndex(), data.GetTofHitIndex());
  195. auto ret = mCandidates.insert(make_pair(key, data));
  196. return &(ret.first->second);
  197. }
  198. //----------------------------------------------------------------------------------------------------------------------------------------
  199. const MpdTofMatchingData* LWeightMatrix::FindCandidate(size_t trackId, size_t hitId)const
  200. {
  201. auto it = mCandidates.find(hash(trackId, hitId));
  202. if(it != mCandidates.cend()) return &(it->second);
  203. return nullptr;
  204. }
  205. //------------------------------------------------------------------------------------------------------------------------
  206. void LWeightMatrix::Print(const char* comment, const TClonesArray* aKfTracks, const TClonesArray* aHits, const TmmT2H* mmT2H, std::ostream& os)const
  207. {
  208. auto prec = os.precision(5);
  209. os<<"\n [LWeightMatrix::Print]-------------------------------------------------------------------->>> ";
  210. if(comment != nullptr) os<<comment;
  211. for(size_t trackId = 0; trackId < fNTracks; trackId++) // cycle by kf track index
  212. {
  213. auto pKfTrack = (MpdTpcKalmanTrack*) aKfTracks->UncheckedAt(trackId);
  214. auto mctid = pKfTrack->GetTrackID();
  215. os<<"\nKfTid="<<trackId<<" tid="<<mctid;
  216. // find max weight hit index
  217. size_t iMaxWeight = -1; double maxWeight = 0.;
  218. for(size_t hitId = 0; hitId<fNHits; hitId++) // cycle by hit index
  219. {
  220. double weight = GetWeight(trackId, hitId);
  221. if(weight > maxWeight)
  222. {
  223. maxWeight = weight;
  224. iMaxWeight = hitId;
  225. }
  226. }
  227. auto info = fmT2MC.find(trackId); // pair<kf trackIndex, mcInfo>
  228. int nPoints = 0;
  229. if(info != fmT2MC.end())
  230. {
  231. nPoints = info->second.Npoints();
  232. os<<" ("<<nPoints<<")\t";
  233. }
  234. if(mmT2H && info != fmT2MC.end())
  235. {
  236. auto range = mmT2H->equal_range(info->second.trackIndex); // pair< trackID, pair<MpdTofHit*, hitIndex> >
  237. for (auto i = range.first; i != range.second; ++i) os<<" hid="<<i->second.second;
  238. }
  239. os<<"\t\t";
  240. if(aHits)
  241. for(size_t hitId = 0; hitId<fNHits; hitId++) // cycle by hit index
  242. {
  243. double weight = GetWeight(trackId, hitId);
  244. auto pHit = (MpdTofHit*) aHits->At(hitId);
  245. Int_t suid = pHit->GetDetectorID();
  246. auto cand = FindCandidate(trackId, hitId); // MpdTofMatchingData
  247. Int_t sector, detector, gap, strip;
  248. MpdTofPoint::ParseSuid(suid, sector, detector, gap, strip);
  249. bool IsTrueMatching = pHit->CheckTid(mctid);
  250. if(weight > 0. || IsTrueMatching)
  251. {
  252. os<<" hitId="<<hitId<<",w=["<<weight<<"]";
  253. os<<"{"<<sector<<","<<detector<<","<<gap<<","<<strip<<"}";
  254. if(cand) os<<"<"<<cand->GetDelta()<<">";
  255. }
  256. if(IsTrueMatching) os<<" *T*";
  257. if(iMaxWeight == hitId) // hit with max weight
  258. {
  259. if(IsTrueMatching) os<<"<<< ";
  260. else if(nPoints != 0) os<<" @@MISMATCH@@";
  261. }
  262. }
  263. }
  264. os<<"\n [LWeightMatrix::Print]--------------------------------------------------------------------<<< ";
  265. os.precision(prec);
  266. }
  267. //------------------------------------------------------------------------------------------------------------------------
  268. void LWeightMatrix::CountFoundCandidates()
  269. {
  270. if(pMatchingQA)
  271. {
  272. size_t nFoundTrueMatch = 0, nFoundMisMatch = 0;
  273. for(auto cand = mCandidates.cbegin(), candEnd = mCandidates.cend(); cand != candEnd; cand++)
  274. {
  275. pMatchingQA->FillParameter(cand->second.fIsTrueMatching, cand->second.GetWeight(), cand->second.GetNormWeight());
  276. if(cand->second.fIsTrueMatching) nFoundTrueMatch++;
  277. else nFoundMisMatch++;
  278. }
  279. fNFoundTrueMatch += nFoundTrueMatch;
  280. fNFoundMisMatch += nFoundMisMatch;
  281. }
  282. }
  283. //------------------------------------------------------------------------------------------------------------------------
  284. void LWeightMatrix::CountCandidates(const TmPt& tids)
  285. {
  286. size_t nCandTofTracks = 0, nCandTrueMatch = 0;
  287. for(const auto& iter : tids) // cycle by TPC KF tracks(Pt descending)
  288. {
  289. Int_t KfIndex = iter.second;
  290. bool HadTofSignal = false, IsTrueMatching = false;
  291. for(size_t hitId = 0; hitId < fNHits; hitId++) // cycle by hits of KfIndex track
  292. {
  293. long uid = hash(KfIndex, hitId);
  294. auto it = mCandidates.find(uid);// pair< hash(trackId, hitId), MpdTofMatchingData >
  295. if(it != mCandidates.end())
  296. {
  297. if(it->second.fHadTofSignal) HadTofSignal = true;
  298. if(it->second.fIsTrueMatching) IsTrueMatching = true;
  299. }
  300. }
  301. if(HadTofSignal) nCandTofTracks++;
  302. if(IsTrueMatching) nCandTrueMatch++;
  303. } // cycle by TPC KF tracks
  304. fNCandTrueMatch += nCandTrueMatch;
  305. fNCandTofTracks += nCandTofTracks;
  306. }
  307. //------------------------------------------------------------------------------------------------------------------------
  308. void LWeightMatrix::Process(TmPt tids, const TClonesArray *aKfTracks)
  309. {
  310. CountCandidates(tids);
  311. for(const auto& iter : tids) // cycle by TPC KF tracks(Pt descending)
  312. {
  313. Int_t KfIndex = iter.second;
  314. // auto pKfTrack = (const MpdTpcKalmanTrack*) iter.first;
  315. // auto mcInfo = getMCinfo(KfIndex);
  316. size_t hitIndex[50];
  317. size_t rank = GetRank(KfIndex, hitIndex);
  318. double bestParameter = 1.E+5; // big value
  319. TmCandidate::iterator bestCand;
  320. size_t bestHitId = -1;
  321. for(size_t i = 0; i < rank; i++)
  322. {
  323. auto it = mCandidates.find(hash(KfIndex, hitIndex[i])); // pair<hash(kf trackId, hitId), MpdTofMatchingData>
  324. if(it != mCandidates.end())
  325. {
  326. double delta = it->second.GetDelta();
  327. if(delta < bestParameter)
  328. {
  329. bestCand = it;
  330. bestParameter = delta;
  331. bestHitId = hitIndex[i];
  332. }
  333. }
  334. else assert(false);// == "ERROR: MpdTofMatchingData candidate MUST exist."
  335. }
  336. if( bestHitId != -1) // accepted candidate
  337. {
  338. // set flag of accepted candidate.
  339. bestCand->second.fBestParameter = true;
  340. // erase this hit at all competitive candidates
  341. RemoveHit(bestHitId);
  342. }
  343. } // cycle by TPC KF tracks
  344. // clean rejected candidates
  345. for(auto cand = mCandidates.cbegin(), candEnd = mCandidates.cend(); cand != candEnd; ) // pair< hash(kftrackId, hitId), MpdTofMatchingData >
  346. {
  347. // fill debug histos
  348. if(pMatchingQA) pMatchingQA->FillHitDeviation(cand->second);
  349. if(cand->second.fBestParameter == false)
  350. {
  351. cand = mCandidates.erase(cand);
  352. }
  353. else cand++;
  354. }
  355. CountFoundCandidates();
  356. }
  357. //------------------------------------------------------------------------------------------------------------------------
  358. void LWeightMatrix::MoveEntries(TClonesArray *aTofMatchings)const
  359. {
  360. for(const auto& cit : mCandidates)
  361. {
  362. new ((*aTofMatchings) [aTofMatchings->GetEntriesFast()]) MpdTofMatchingData(cit.second);
  363. }
  364. }
  365. //------------------------------------------------------------------------------------------------------------------------
  366. ClassImp(MpdTofMatching)
  367. //------------------------------------------------------------------------------------------------------------------------
  368. MpdTofMatching::MpdTofMatching(const char *name, Int_t verbose, Bool_t test, const char *flnm)
  369. : FairTask(name, verbose), fDoTest(test)
  370. {
  371. pRandom = new TRandom2;
  372. if(fDoTest)
  373. {
  374. pQA = new MpdTofMatchingQA(flnm, false);
  375. pTimer = new TStopwatch;
  376. }
  377. fWeights = new LWeightMatrix(pQA);
  378. }
  379. //------------------------------------------------------------------------------------------------------------------------
  380. MpdTofMatching::~MpdTofMatching()
  381. {
  382. delete pTimer;
  383. delete fWeights;
  384. delete pRandom;
  385. delete pQA;
  386. }
  387. //------------------------------------------------------------------------------------------------------------------------
  388. InitStatus MpdTofMatching::Init()
  389. {
  390. LOG(DEBUG2)<<"[MpdTofMatching::Init] Begin initialization.";
  391. aMcPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("TOFPoint");
  392. aMcTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
  393. aTofHits = (TClonesArray*) FairRootManager::Instance()->GetObject("TOFHit");
  394. aTPCkfTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  395. aECTkfTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("EctTrack");
  396. if(aMcPoints && aMcTracks) fIsMcRun = true;
  397. if(!aTofHits || !aTPCkfTracks){ LOG(ERROR)<<"Branch not found! (TOFHit or TpcKalmanTrack)"; return kERROR; }
  398. pKF = MpdKalmanFilter::Instance("KF","KF");
  399. // Create and register output array
  400. aTofMatchings = new TClonesArray("MpdTofMatchingData");
  401. FairRootManager::Instance()->Register("TOFMatching", "Tof", aTofMatchings, kTRUE);
  402. MpdTofGeoUtils::Instance()->ParseTGeoManager(nullptr, false);
  403. MpdTofGeoUtils::Instance()->FindNeighborStrips(0.8, nullptr, false); // 0.8 [cm] <--- thresh. distance between neighbor strips, (see h1TestDistance histo)
  404. LOG(DEBUG)<<"[MpdTofMatching::Init] Initialization finished succesfully.";
  405. return kSUCCESS;
  406. }
  407. //------------------------------------------------------------------------------------------------------------------------
  408. void MpdTofMatching::Finish()
  409. {
  410. if(pQA) pQA->Finish();
  411. if(pTimer){ Printf("\nMpdTofMatching TStopwatch report: "); pTimer->Print("m");}
  412. }
  413. //------------------------------------------------------------------------------------------------------------------------
  414. void MpdTofMatching::Exec(Option_t *option)
  415. {
  416. if(pTimer) pTimer->Start(false);
  417. fDoMCtest = fDoTest && fIsMcRun;
  418. // Reset event
  419. aTofMatchings->Clear();
  420. Int_t nTofPoints = -1, nMCTracks = -1;
  421. if(fIsMcRun){ nTofPoints = aMcPoints->GetEntriesFast(); nMCTracks = aMcTracks->GetEntriesFast();}
  422. Int_t nTofHits = aTofHits->GetEntriesFast();
  423. Int_t nKFTracks = aTPCkfTracks->GetEntries();
  424. Int_t nECTtracks = (aECTkfTracks) ? aECTkfTracks->GetEntriesFast() : 0;
  425. LOG(DEBUG2)<<"[MpdTofMatching::Exec] points= "<<nTofPoints<<", hits= "<<nTofHits<<", mc tracks= "<<nMCTracks<<", kf tracks= "<<nKFTracks;
  426. // ---------------------------------------------------------------------------------------->>> Select (! propagated to Etof) & sort by Pt the KfTpcTracks
  427. std::set<Int_t> sEndcapTracks;
  428. if(aECTkfTracks != nullptr)
  429. {
  430. for(Int_t i = 0; i < nECTtracks; i++) // cycle by ECT KF tracks
  431. {
  432. auto track = (MpdEctKalmanTrack*) aECTkfTracks->UncheckedAt(i);
  433. if(track->IsFromTpc()) sEndcapTracks.insert(track->GetTpcIndex());
  434. }
  435. }
  436. // ----------------------------------------------------------------------------------------
  437. TmPt mPt;
  438. for(Int_t index = 0; index < nKFTracks; index++) // cycle by TPC KF tracks
  439. {
  440. if(aECTkfTracks != nullptr && sEndcapTracks.find(index) != sEndcapTracks.end()) continue; // postpone for matching with ETof
  441. auto pKfTrack = (MpdTpcKalmanTrack*) aTPCkfTracks->UncheckedAt(index);
  442. mPt.insert(make_pair(pKfTrack, index));
  443. }
  444. // ---------------------------------------------------------------------------------------->>> Sorting & Mapping points to MC tracks
  445. TmmT2P mmT2P; // pair< MCtrackID, MpdTofPoint*>
  446. if(fDoMCtest) MapMcT2P(aMcPoints, mmT2P);
  447. // ---------------------------------------------------------------------------------------->>> Mapping hits to detectorUIDs
  448. TmmD2H mmD2H; // pair< duid, pair<MpdTofHit*, hitIndex> >
  449. TmS2H mS2H; // pair< suid, pair<MpdTofHit*, hitIndex> >
  450. TmmT2H mmT2H; // pair< tid, pair<MpdTofHit*, hitIndex> >
  451. vector<Int_t> v;
  452. for(Int_t hitIndex = 0; hitIndex < nTofHits; hitIndex++) // cycle by tof hits
  453. {
  454. auto pTofHit = (MpdTofHit*) aTofHits->At(hitIndex);
  455. Int_t suid = pTofHit->GetDetectorID();
  456. mmD2H.insert(make_pair( (suid & 0xFFFF0000), make_pair(pTofHit, hitIndex))); // convert suid to duid (reset stripID to 0 and gap to 0)
  457. mS2H.insert(make_pair(suid, make_pair(pTofHit, hitIndex)));
  458. pTofHit->getLinks(MpdTofUtils::mcTrackIndex, v); // vector<tid>
  459. for(auto it=v.begin(), itEnd = v.end(); it != itEnd; it++) mmT2H.insert(make_pair(*it, make_pair(pTofHit, hitIndex)));
  460. }
  461. // ----------------------------------------------------------------------------------------
  462. if(fDoMCtest) pQA->FillIdealMatching(mPt, mmT2H, aMcTracks);
  463. // resize matrix & clean by zero weight
  464. fWeights->Reset(mPt.size(), mS2H.size());
  465. // fill matching weight matrix.
  466. FillWeights(mPt, mmD2H, mS2H, mmT2P);
  467. if(fDoMCtest) pQA->FillNtracks(aMcTracks, fNTofTracksEvent); // MUST be call after FillWeights(calc. fNTofTracksEvent)
  468. //fWeights->Print("BEFORE", aTPCkfTracks, aTofHits, &mmT2H);
  469. // fWeights->Boost1thRankCand(1000.);
  470. // fWeights->Normalize();
  471. ///fWeights->Print("AfterNormalize", aTPCkfTracks, aTofHits, &mmT2H);
  472. fWeights->Process(mPt, aTPCkfTracks);
  473. //fWeights->Print("AFTER Process", aTPCkfTracks, aTofHits, &mmT2H);
  474. fWeights->MoveEntries(aTofMatchings);
  475. //MpdTof::Dump(aMcPoints, aTofHits, aMcTracks);
  476. if(fDoMCtest) pQA->FillMatchingEfficiency(mPt, aTofMatchings, aTofHits, aMcTracks);
  477. //pQA->tidsTofTouch.Print("\n tidsTofTouch ------------------------>>>");
  478. set<int> matchedTracks; // MCtid
  479. for(Int_t i = 0, N = aTofMatchings->GetEntriesFast(); i < N; i++)
  480. {
  481. auto match = (MpdTofMatchingData*) aTofMatchings->At(i);
  482. int KFtid = match->GetKFTrackIndex();
  483. auto pKfTrack = (MpdTpcKalmanTrack*) aTPCkfTracks->UncheckedAt(KFtid);
  484. matchedTracks.insert(pKfTrack->GetTrackID());
  485. }
  486. set<int> kalmanTracks; // MCtid
  487. for(Int_t i = 0, N = aTPCkfTracks->GetEntriesFast(); i < N; i++)
  488. {
  489. auto pKfTrack = (MpdTpcKalmanTrack*) aTPCkfTracks->UncheckedAt(i);
  490. kalmanTracks.insert(pKfTrack->GetTrackID());
  491. }
  492. TVector3 momentum;
  493. for(Int_t tid = 0, Ntracks = aMcTracks->GetEntriesFast(); tid < Ntracks; tid++) // cycle by mc tracks
  494. {
  495. //if(()->GetNPoints(kTOF)) NmcTracks++;
  496. auto track = (MpdMCTrack*) aMcTracks->UncheckedAt(tid);
  497. if(track->GetMotherId() == -1) // primary mc track
  498. {
  499. bool IsKalman = (kalmanTracks.find(tid) != kalmanTracks.end());
  500. bool IsMatched = (matchedTracks.find(tid) != matchedTracks.end());
  501. track->GetMomentum(momentum);
  502. if(pQA) pQA->FillTrackEfficiency(IsKalman, IsMatched, track->GetPdgCode(), momentum);
  503. }
  504. }
  505. Report(fVerbose);
  506. if(pTimer) pTimer->Stop();
  507. }
  508. //------------------------------------------------------------------------------------------------------------------------
  509. void MpdTofMatching::Report(Int_t verbose, ostream& os)const
  510. {
  511. auto prec = os.precision(3);
  512. if(verbose >0)
  513. {
  514. Int_t nMatchings = aTofMatchings->GetEntries();
  515. os<<"\n -I- [MpdTofMatching::Exec] nTofTracks="<<fNTofTracksEvent<<", Matchings = "<<nMatchings;
  516. if(fNTofTracksEvent != 0) os<<" ("<<((double) nMatchings)/fNTofTracksEvent*100.<<"%)";
  517. }
  518. if(verbose >1)
  519. {
  520. auto cand = fWeights->GetMatchCand();
  521. auto found = fWeights->GetMatchFound();
  522. if(fNTofTracksRun != 0)
  523. {
  524. os<<"\n Cand.: tofTracks = ("<<((double) cand.first)/fNTofTracksRun*100.<<"%), TrueMatch = ("<<((double) cand.second)/fNTofTracksRun*100.<<"%)";
  525. os<<"; Found: TrueMatch = ("<<((double) found.first)/fNTofTracksRun*100.<<"%), MisMatch = ("<<((double) found.second)/fNTofTracksRun*100.<<"%)";
  526. }
  527. }
  528. os.precision(prec);
  529. }
  530. //------------------------------------------------------------------------------------------------------------------------
  531. #include "MpdTofGeoUtils.h"
  532. //------------------------------------------------------------------------------------------------------------------------
  533. void MpdTofMatching::FillWeights(const TmPt& mPt, const TmmD2H& mmD2H, const TmS2H& mS2H, const TmmT2P& mmT2P)
  534. {
  535. const double Zerror = 10.; // 10. [cm] // FIXME: should be dependent on the parameters of the KFtrack
  536. const double PhiError = 0.1;// 0.1 // [rads] // FIXME: should be dependent on the parameters of the KFtrack
  537. //if(pQA) pQA->tidsTofTouch.Reset();
  538. mcInfo *pMCdata = {};
  539. double trLength; //, trackLength;
  540. //Int_t charge;
  541. vector<Tinterval> list; // list of strips overlapped probe point
  542. fNTofTracksEvent = 0;
  543. for(const auto& it : mPt) // cycle by TPC KF tracks
  544. {
  545. auto pKfTrack = (const MpdTpcKalmanTrack*) it.first;
  546. Int_t KfIndex = it.second;
  547. //Int_t tid = pKfTrack->GetTrackID();
  548. if(fDoMCtest)
  549. {
  550. pMCdata = fWeights->saveMCinfo(KfIndex, GetMcInfo(mmT2P, pKfTrack, aMcTracks));
  551. pQA->FillKfTracks(pMCdata->TofTouch, pKfTrack->Momentum3().Eta(), pKfTrack->Momentum());
  552. }
  553. if(pMCdata && pMCdata->TofTouch) fNTofTracksEvent++;
  554. TVector3 estPointOnCyl = EstTrackOnR(pKfTrack); // Estimate point on cylinder
  555. if(MpdTofGeoUtils::Instance()->IsPointInsideDetectors(estPointOnCyl, list, Zerror, PhiError)) // !empty detector list
  556. {
  557. if(pQA) pQA->FillPointInsideDetector(estPointOnCyl);
  558. MpdTpcKalmanTrack ReFittedTrack(RefitTrack(pKfTrack));
  559. for(const auto& detector : list) // cycle by overlapped detector list
  560. {
  561. Int_t duid = detector.value->volumeUID;
  562. auto detCenter = detector.value->center;
  563. auto detPerp = detector.value->perp;
  564. // Propagate the KF track to the Tof detector plane; return true, if propagation successful
  565. TVector3 estPointOnPlate, hitPosition, Momentum;
  566. if(0 == EstTrackOnPlate(ReFittedTrack, detCenter, detPerp, estPointOnPlate, trLength, Momentum)) // Estimate point on detector plate
  567. {
  568. /* // Propagate the KF marker tracks to the detector plane; return true, if all propagations successful
  569. bool retvalue = true;
  570. TVector3 est[4], P[4];
  571. for(size_t i=0;i<4;i++) retvalue = retvalue && EstMarkerTrackOnPlate(i, pKfTrack, detCenter, detPerp, est[i], P[i]);
  572. if(pQA) pQA->FillSmearedPoints(estPointOnPlate, est, P);
  573. */
  574. if(pMCdata && pQA)
  575. {
  576. TVector3 pos, mom;
  577. if(pMCdata->GetClosestPointOnDetector(pos, mom, estPointOnPlate, duid).first != -1)
  578. pQA->FillMcPointDeviation(pos, estPointOnPlate, mom, Momentum, ReFittedTrack.GetNofHits());
  579. }
  580. auto range = mmD2H.equal_range(duid);
  581. for(auto iter = range.first; iter != range.second; iter++) // cycle by estimated detector hits
  582. {
  583. MpdTofHit *pHit = iter->second.first;
  584. Int_t hitIndex = iter->second.second;
  585. pHit->Position(hitPosition);
  586. double devHit, devHitZ, devHitR, devHitPhi;
  587. MpdTof::GetDelta(hitPosition, estPointOnPlate, devHit, devHitZ, devHitR, devHitPhi);
  588. devHitZ = abs(devHitZ);
  589. devHitPhi = abs(devHitPhi);
  590. if(devHitZ < fThreshZ && devHitPhi < fThreshPhi) // inside matching window
  591. {
  592. double weight = 1./(estPointOnPlate - hitPosition).Mag();
  593. fWeights->AddWeight(KfIndex, hitIndex, weight);
  594. auto ptr = fWeights->AddCandidate(MpdTofMatchingData(KfIndex, hitIndex, weight, pHit, trLength, pKfTrack->GetNofTrHits(), Momentum, estPointOnPlate));
  595. if(pMCdata)
  596. {
  597. ptr->fIsTrueMatching = pMCdata->IsSameStrip(pHit->GetDetectorID());
  598. ptr->fHadTofSignal = pMCdata->TofTouch;
  599. }
  600. } // inside matching window
  601. }
  602. }
  603. }
  604. }
  605. } // cycle by TPC KF tracks
  606. fNTofTracksRun += fNTofTracksEvent;
  607. }
  608. //------------------------------------------------------------------------------------------------------------------------
  609. mcInfo MpdTofMatching::GetMcInfo(const TmmT2P& mmT2P, const MpdKalmanTrack* pKfTrack, TClonesArray* aMcTracks)
  610. {
  611. mcInfo mcData;
  612. mcData.trackIndex = pKfTrack->GetTrackID();
  613. auto mcTrack = (const MpdMCTrack*) aMcTracks->UncheckedAt(mcData.trackIndex);
  614. mcData.pid = mcTrack->GetPdgCode();
  615. auto track = const_cast<MpdMCTrack*>(mcTrack);
  616. track->GetMomentum(mcData.vertMomentum);
  617. track->GetStartVertex(mcData.vertPosition);
  618. auto range = mmT2P.equal_range(mcData.trackIndex); // pair< MCtrackID, MpdTofPoint*>
  619. for(auto it = range.first; it != range.second; ++it) mcData.vPoints.push_back(it->second); // cycle by mcPoints for trackIndex mcTrack
  620. mcData.TofTouch = ! mcData.vPoints.empty();
  621. return mcData;
  622. }
  623. //------------------------------------------------------------------------------------------------------------------------
  624. void MpdTofMatching::MapMcT2P(TClonesArray *aPoints, TmmT2P& mmT2P)const // TmmT2P = pair< MCtrackID, MpdTofPoint*>
  625. {
  626. for(Int_t index = 0, nTofPoints = aPoints->GetEntriesFast(); index < nTofPoints; index++) // cycle by MpdTofPoint
  627. {
  628. auto mcPoint = (MpdTofPoint*) aPoints->UncheckedAt(index);
  629. Int_t tid = mcPoint->GetTrackID();
  630. Double_t time = mcPoint->GetTime();
  631. auto iter = mmT2P.find(tid);
  632. if(iter != mmT2P.end()) // same trackID already inserted, insert to position(sorting by time)
  633. {
  634. int count = mmT2P.count(tid);
  635. for(int i = 0; i < count; i++, iter++) // cycle by hits with same trackID
  636. {
  637. if(time < iter->second->GetTime())
  638. {
  639. mmT2P.insert(iter, make_pair(tid, mcPoint));
  640. break;
  641. }
  642. if(i == count-1) mmT2P.insert(++iter, make_pair(tid, mcPoint)); // insert to last
  643. }
  644. }
  645. else mmT2P.insert(make_pair(tid, mcPoint));
  646. } // cycle by MpdTofPoint
  647. }
  648. //------------------------------------------------------------------------------------------------------------------------
  649. TVector3 MpdTofMatching::EstTrackOnR(const MpdTpcKalmanTrack *tr)const
  650. {
  651. MpdKalmanHit hEnd;
  652. hEnd.SetType(MpdKalmanHit::kFixedR);
  653. MpdTpcKalmanTrack tr1(*tr);
  654. TObjArray *hits = tr1.GetHits();
  655. if(hits->GetEntriesFast() == 0)
  656. {
  657. Int_t nh = tr1.GetNofTrHits();
  658. for (Int_t j = 0; j < nh; ++j) hits->Add(tr1.GetTrHits()->UncheckedAt(j));
  659. tr1.SetNofHits(nh);
  660. }
  661. tr1.SetParam(*tr1.GetParamAtHit());
  662. tr1.SetParamNew(*tr1.GetParamAtHit());
  663. tr1.SetWeight(*tr1.GetWeightAtHit());
  664. tr1.SetLength(tr1.GetLengAtHit());
  665. auto h = (MpdKalmanHit*) tr1.GetTrHits()->First();
  666. tr1.SetPos(tr1.GetPosAtHit());
  667. if(h->GetType() == MpdKalmanHit::kFixedZ) tr1.SetType(MpdKalmanTrack::kEndcap);
  668. tr1.SetPosNew(tr1.GetPos());
  669. hEnd.SetPos(fTofBarrelRadius); // barrel TOF radius, [cm]
  670. pKF->PropagateToHit(&tr1, &hEnd, kTRUE);
  671. TVector3 pos(tr1.GetPosNew(), 0.,0.);
  672. pos.SetPhi(tr1.GetParamNew(0)/tr1.GetPosNew());
  673. pos.SetZ(tr1.GetParamNew(1));
  674. return pos;
  675. }
  676. //------------------------------------------------------------------------------------------------------------------------
  677. MpdTpcKalmanTrack MpdTofMatching::RefitTrack(const MpdTpcKalmanTrack *tr)const
  678. {
  679. MpdTpcKalmanTrack tr1(*tr);
  680. tr1.SetDirection(MpdKalmanTrack::kOutward);
  681. TObjArray *hits = tr1.GetHits();
  682. if(hits->GetEntriesFast() == 0)
  683. {
  684. Int_t nh = tr1.GetNofTrHits();
  685. for (Int_t j = 0; j < nh; ++j) hits->Add(tr1.GetTrHits()->UncheckedAt(j));
  686. tr1.SetNofHits(nh);
  687. }
  688. tr1.SetParam(*tr1.GetParamAtHit());
  689. tr1.SetParamNew(*tr1.GetParamAtHit());
  690. tr1.SetWeight(*tr1.GetWeightAtHit());
  691. tr1.SetLength(tr1.GetLengAtHit());
  692. tr1.SetPos(tr1.GetPosAtHit());
  693. tr1.SetPosNew(tr1.GetPos());
  694. return tr1;
  695. }
  696. //------------------------------------------------------------------------------------------------------------------------
  697. int MpdTofMatching::EstTrackOnPlate(const MpdTpcKalmanTrack& tr, const TVector3& point, const TVector3& perp, TVector3& pos, Double_t& length, TVector3& P) const
  698. {
  699. MpdTpcKalmanTrack tr1(tr);
  700. Double_t plane[6] = {point.X(), point.Y(), point.Z(),perp.X(), perp.Y(), perp.Z()};
  701. if (!pKF->PropagateParamP(&tr1,plane,kTRUE)) return -1;
  702. Double_t R = tr1.GetPosNew();
  703. if(TMath::IsNaN(R)) return -2;
  704. length = tr1.GetLength();
  705. P = tr1.Momentum3();
  706. pos.SetXYZ(R, 0., 0.);
  707. pos.SetPhi(tr1.GetParamNew(0)/R);
  708. pos.SetZ(tr1.GetParamNew(1));
  709. return 0;
  710. }
  711. //------------------------------------------------------------------------------------------------------------------------
  712. bool MpdTofMatching::EstMarkerTrackOnPlate(size_t mode, const MpdTpcKalmanTrack* tr, const TVector3& point, const TVector3& perp, TVector3& estPoP, TVector3& P) const
  713. {
  714. const static Double_t fError_dZ = 0.03, fError_dPhi = 0.03; // [rads] msc
  715. auto tr1 = RefitTrack(tr);
  716. TMatrixD& param = *tr1.GetParamNew();
  717. switch(mode)
  718. {
  719. case 0:
  720. param(3,0) += fError_dZ;
  721. break;
  722. case 1:
  723. param(2,0) -= fError_dPhi;
  724. break;
  725. case 2:
  726. param(3,0) -= fError_dZ;
  727. break;
  728. case 3:
  729. param(2,0) += fError_dPhi;
  730. break;
  731. default:;
  732. };
  733. tr1.SetParam(param);
  734. Double_t plane[6] = {point.X(), point.Y(), point.Z(),perp.X(), perp.Y(), perp.Z()};
  735. if (!pKF->PropagateParamP(&tr1,plane,kTRUE)) return false;
  736. Double_t R = tr1.GetPosNew();
  737. if(TMath::IsNaN(R)) return false;
  738. P = tr1.Momentum3();
  739. estPoP.SetXYZ(R, 0., 0.);
  740. estPoP.SetPhi(tr1.GetParamNew(0)/R);
  741. estPoP.SetZ(tr1.GetParamNew(1));
  742. return true;
  743. }
  744. //------------------------------------------------------------------------------------------------------------------------
  745. void MpdTofMatching::Print(const MpdTpcKalmanTrack *track, const char* comment, ostream& os)const
  746. {
  747. if(comment != nullptr) os<<comment;
  748. Double_t rad = track->GetPosNew();
  749. Double_t phi = track->GetParamNew(2);
  750. double X = rad * TMath::Cos(phi); // X
  751. double Y = rad * TMath::Sin(phi); // Y
  752. double Z = track->GetParamNew(1); // Z
  753. TVector3 P = track->Momentum3();
  754. os<<"Mom=("<<P.X()<<","<<P.Y()<<","<<P.Z()<<";"<<P.Perp()<<","<<P.Mag()<<") Pos=("<<X<<","<<Y<<","<<Z<<"; "<<sqrt(X*X+Y*Y)<<","<<sqrt(X*X+Y*Y+Z*Z)<<")";
  755. }
  756. //------------------------------------------------------------------------------------------------------------------------