LMatchingFilter.h 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849
  1. //----------------------------------------------------------------------------------------------------------------------------------------
  2. #ifndef __MPD_LMATCHING_FILTER_H
  3. #define __MPD_LMATCHING_FILTER_H 1
  4. #include<assert.h>
  5. #include<vector>
  6. #include<map>
  7. //#include<debug/map>
  8. #include<set>
  9. #include<iostream>
  10. #include<TClonesArray.h>
  11. #include<TH2D.h>
  12. #include "FairLogger.h"
  13. #include "FairMCTrack.h"
  14. #include "MpdTofMatchingQA.h"
  15. #include "MpdTofMatchingData.h"
  16. #include "MpdTofHit.h"
  17. #include "MpdTofPoint.h"
  18. //#define DoTest
  19. struct CandHistory
  20. {
  21. enum kFlags{ kNothing = 0, k1thRank = (1 << 1), kNoCompetitor = (1 << 2), kHaveCompetitor = (1 << 3) };
  22. Int_t flag;
  23. std::vector<MpdTofMatchingData*> vCompetitors;
  24. MpdTofMatchingData *pCandidate;
  25. void AddFlag(Int_t f){ flag = flag | f ;}
  26. void ResetFlag(Int_t f)
  27. {
  28. Clear();
  29. flag = f;
  30. }
  31. void Clear()
  32. {
  33. flag = kNothing;
  34. vCompetitors.clear();
  35. pCandidate = nullptr;
  36. }
  37. void Print()
  38. {
  39. switch(flag)
  40. {
  41. case k1thRank | kNoCompetitor:
  42. std::cout<<"\n k1thRank & kNoCompetitor pair<"<<pCandidate->GetKFTrackIndex()<<", "<<pCandidate->GetTofHitIndex()<<"> deltaZ="<<pCandidate->GetDelta1()<<" deltaPhi="<<pCandidate->GetDelta2();
  43. break;
  44. case k1thRank | kHaveCompetitor:
  45. std::cout<<"\n k1thRank & kHaveCompetitor pair<"<<pCandidate->GetKFTrackIndex()<<", "<<pCandidate->GetTofHitIndex()<<"> deltaZ="<<pCandidate->GetDelta1()<<" deltaPhi="<<pCandidate->GetDelta2();
  46. for(int i=0; i < vCompetitors.size();i++)
  47. {
  48. MpdTofMatchingData* ptr = vCompetitors[i];
  49. std::cout<<"\n\t competitiors pair<"<<ptr->GetKFTrackIndex()<<", "<<ptr->GetTofHitIndex()<<"> deltaZ="<<ptr->GetDelta1()<<" deltaPhi="<<ptr->GetDelta2();
  50. }
  51. break;
  52. };
  53. }
  54. };
  55. //----------------------------------------------------------------------------------------------------------------------------------------
  56. class LMatchingFilter
  57. {
  58. public:
  59. template<typename T> using TmmIntT = std::multimap<Int_t, T>;
  60. using TmmCand = TmmIntT<MpdTofMatchingData*>;
  61. using candIter = TmmCand::iterator;
  62. using linkIter = TmmIntT<candIter>::iterator;
  63. private:
  64. TClonesArray *aMdata; // Matching data container;
  65. TmmIntT<MpdTofMatchingData> mcTrueMdata, mcMaybeMdata; // pair<kfTrackIndex, MpdTofMatchingData>, filled by mcFindTrueMachings method, if used MC run data.
  66. TmmIntT<MpdTofMatchingData*> mmCandidate, mmCandidateSnapshot, mmAccepted, mmAcceptedSnapshot; // pair<KfIndex, MpdTofMatchingData*>
  67. TmmIntT<candIter> mmLinks, mmLinksSnapshot; // pair<hitIndex, candIter>
  68. TmmIntT<Int_t> mmRanks; // pair<number of candidate hits for current KFtrack, KfIndex>
  69. MpdTofMatchingQA *pQA;
  70. int fVerbose;
  71. bool isOwner;
  72. bool isMCdata; // true, if used MC run data. (TofHits and KalmanTracks have links to MC origins.)
  73. double fChi2;
  74. CandHistory candHist;
  75. // for debug
  76. void Commit(void);
  77. void Status(const char* comment = nullptr, std::ostream& os = std::cout) const;
  78. void PrintChanges(const TmmIntT<candIter>& src, const TmmIntT<candIter>& checked, std::ostream& os = std::cout) const;
  79. void PrintChanges(const TmmCand& src, const TmmCand& checked, std::ostream& os = std::cout) const; // cycle by src entries, print missed into the checked container entries.
  80. bool CheckAdequacy(const char* comment = nullptr); // check mmCandidate & mmLinks entries accordance; return true if pass test
  81. Int_t ProcessSingleTracks(void);
  82. void ProcessNRankTrack(Int_t rank);
  83. void UpdateRanks();
  84. void RemoveLink(Int_t hitIndex, Int_t KfIndex); // remove link to pair<KfIndex, hitIndex>
  85. void RemoveCandidate(Int_t hitIndex, Int_t KfIndex);
  86. candIter FindCandidate(TmmCand *map, Int_t hitIndex, Int_t KfIndex); // find candidate iter to map pair<KfIndex, hitIndex>
  87. void AcceptCandidate(candIter itCand); // move candidate to accepted, clean links and cross candidates
  88. candIter FindClosest1hRankCandidate(const std::pair<linkIter, linkIter>& range, Double_t& minDelta);
  89. double FindMinDeltaForTrack(Int_t KfIndex, Int_t disabledHitIndex, candIter& minCand, bool& IsWithoutCompetition); // Search through all the KfIndex candidates except disabledHitIndex
  90. inline bool IsSameRank(Int_t rank, Int_t KfIndex) // return true, if KfIndex cantidate have the "rank" rank
  91. {
  92. auto range = mmRanks.equal_range(rank);
  93. for(auto iter = range.first; iter != range.second; ++iter) if(iter->second == KfIndex) return true;
  94. return false;
  95. }
  96. inline Int_t GetMinRank()const { auto iter = mmRanks.begin(); assert(iter != mmRanks.end()); return iter->first; }
  97. inline Int_t GetMinRankAfter(Int_t rank)const{ auto iter = mmRanks.upper_bound(rank); if(iter != mmRanks.end()) return iter->first; else return -1; }
  98. public:
  99. LMatchingFilter(MpdTofMatchingQA*, int verbose = 0);
  100. ~LMatchingFilter();
  101. void SetContainer(TClonesArray *array);
  102. void AddCandidate(const MpdTofMatchingData& data);
  103. void Reset();
  104. Int_t Processing(Int_t nKFTracks, double& chi2);
  105. Int_t UpdateContainer();
  106. void Dump(const char* comment = nullptr, std::ostream& os = std::cout) const;
  107. void PrintRanks(const char* comment = nullptr, std::ostream& os = std::cout) const;
  108. // method used MC data
  109. template<typename kfTrackT>
  110. void mcFindTrueMachings( TClonesArray *aPoints, TClonesArray *aHits, TClonesArray *aMCTracks, TClonesArray *aKfTracks, bool print = false); // work correctly only if used MC run data.
  111. bool mcIsTrueMatching(const MpdTofMatchingData*);
  112. bool mcIsMaybeMatching(const MpdTofMatchingData*);
  113. };
  114. //----------------------------------------------------------------------------------------------------------------------------------------
  115. LMatchingFilter::LMatchingFilter(MpdTofMatchingQA *ptr, int verbose)
  116. : pQA(ptr), fVerbose(verbose), isOwner(true), isMCdata(false)
  117. {
  118. aMdata = new TClonesArray("MpdTofMatchingData", 1000);
  119. }
  120. //----------------------------------------------------------------------------------------------------------------------------------------
  121. LMatchingFilter::~LMatchingFilter()
  122. {
  123. if(isOwner) delete aMdata;
  124. }
  125. //----------------------------------------------------------------------------------------------------------------------------------------
  126. void LMatchingFilter::AddCandidate(const MpdTofMatchingData& data)
  127. {
  128. MpdTofMatchingData *pData = new ((*aMdata) [aMdata->GetEntriesFast()]) MpdTofMatchingData(data); // save data to the TClonesArray container
  129. candIter it = mmCandidate.insert(std::make_pair(pData->GetKFTrackIndex(), pData)); // insert matching candidate into the mmCandidate
  130. mmLinks.insert(std::make_pair(pData->GetTofHitIndex(), it)); // insert back link to matching candidate
  131. }
  132. //----------------------------------------------------------------------------------------------------------------------------------------
  133. Int_t LMatchingFilter::UpdateContainer()
  134. {
  135. // select accepted candidates
  136. std::set<MpdTofMatchingData*> sAccepted;
  137. for(const auto& cit : mmAccepted) sAccepted.insert(cit.second);
  138. for(int entry = 0, Nentries = aMdata->GetEntriesFast(); entry < Nentries; entry++) // cycle by MpdTofMatchingData at TClonesArray container
  139. {
  140. MpdTofMatchingData *pData = (MpdTofMatchingData*) aMdata->At(entry);
  141. if(sAccepted.find(pData) == sAccepted.end()) // no exist -> not accepted
  142. aMdata->RemoveAt(entry);
  143. }
  144. aMdata->Compress();
  145. return aMdata->GetEntriesFast();
  146. }
  147. //----------------------------------------------------------------------------------------------------------------------------------------
  148. LMatchingFilter::candIter LMatchingFilter::FindCandidate(TmmCand *map, Int_t hitIndex, Int_t KfIndex)
  149. {
  150. auto rangeCands = map->equal_range(KfIndex);
  151. for(candIter iter = rangeCands.first; iter != rangeCands.second; ++iter) // cycle by hits for this KfIndex
  152. {
  153. if(hitIndex == iter->second->GetTofHitIndex()) return iter;
  154. }
  155. assert(false);
  156. return map->end();
  157. }
  158. #include"MpdDetectorList.h"
  159. //----------------------------------------------------------------------------------------------------------------------------------------
  160. template<typename kfTrackT>
  161. void LMatchingFilter::mcFindTrueMachings(TClonesArray *aPoints, TClonesArray *aHits, TClonesArray *aMCTracks, TClonesArray *aKfTracks, bool print)
  162. {
  163. // Cleanup or initialize container
  164. mcTrueMdata.clear();
  165. mcMaybeMdata.clear();
  166. TmmIntT<kfTrackT*> mTracks;
  167. TmmIntT<MpdTofHit*> mHits;
  168. std::multimap<kfTrackT*, int> mTracksInv; // pair <Mpd?KalmanTrack*, kftrackIndex>
  169. std::map<MpdTofHit*, int> mHitsInv; // pair <MpdTofHit*, hitIndex>
  170. // Sorting MpdTpcKalmanTracks by tid
  171. Int_t nKFTracks = aKfTracks->GetEntriesFast();
  172. for(size_t KfIndex = 0; KfIndex < nKFTracks; KfIndex++) // cycle by MpdTpcKalmanTrack*
  173. {
  174. auto pKfTrack = (kfTrackT*) aKfTracks->UncheckedAt(KfIndex);
  175. Int_t tid = pKfTrack->GetTrackID();
  176. mTracks.insert(make_pair(tid, pKfTrack));
  177. mTracksInv.insert(make_pair(pKfTrack, KfIndex));
  178. }
  179. // Sorting MpdTofHits by tid
  180. Int_t nHits = aHits->GetEntriesFast();
  181. for(size_t hitIndex = 0; hitIndex < nHits; hitIndex++) // cycle by MpdTofHit*
  182. {
  183. auto pHit = (MpdTofHit*) aHits->UncheckedAt(hitIndex);
  184. auto pPoint = (MpdTofPoint*) aPoints->UncheckedAt(pHit->GetRefIndex());
  185. Int_t tid = pPoint->GetTrackID();
  186. mHits.insert(make_pair(tid, pHit));
  187. mHitsInv.insert(make_pair(pHit, hitIndex));
  188. }
  189. // Fill "True" & "maybe" matching pairs
  190. Int_t tid = 0; // dummy value, init later
  191. for(auto itTrack = mTracks.begin(), itTrackEnd = mTracks.end(); itTrack != itTrackEnd; itTrack = mTracks.upper_bound(tid)) // cycle by MpdKalmanTrack* with unique tid
  192. {
  193. tid = itTrack->first;
  194. auto mcTrack = (const FairMCTrack*) aMCTracks->UncheckedAt(itTrack->second->GetTrackID());
  195. if(! mcTrack->GetNPoints(kTOF)) continue; // don't have Tof hit
  196. auto itHit = mHits.find(tid);
  197. if(itHit == mHits.end()) continue; // don't have hit for same tid;
  198. size_t counterT = mTracks.count(tid);
  199. size_t counterH = mHits.count(tid);
  200. // container = mcTrueMdata, if relation one to one (both kftracks to mctrack and hits to mctrack)
  201. auto container = ( counterT == 1 && counterH == 1 ) ? &mcTrueMdata : &mcMaybeMdata;
  202. for(size_t t = 0; t < counterT; t++, itTrack++)
  203. {
  204. kfTrackT *track = itTrack->second;
  205. Int_t kfTrackIndex = (mTracksInv.find(track))->second;
  206. for(size_t h = 0; h < counterH; h++, itHit++)
  207. {
  208. MpdTofHit *hit = itHit->second;
  209. Int_t hitIndex = (mHitsInv.find(hit))->second;
  210. container->insert(make_pair(kfTrackIndex, MpdTofMatchingData(kfTrackIndex, hitIndex)));
  211. }
  212. }
  213. }
  214. if(print)
  215. {
  216. std::cout<<"\n\t True Matching data: ----------------------------------------------------------------------";
  217. for(const auto& it : mcTrueMdata) it.second.Print();
  218. std::cout<<"\n\t Maybe Matching data: ----------------------------------------------------------------------";
  219. for(const auto& it : mcMaybeMdata) it.second.Print();
  220. }
  221. }
  222. //----------------------------------------------------------------------------------------------------------------------------------------
  223. bool LMatchingFilter::mcIsTrueMatching(const MpdTofMatchingData *pSample)
  224. {
  225. auto it = mcTrueMdata.find(pSample->GetKFTrackIndex()); // pair<kfindex, MdataT>
  226. if(it == mcTrueMdata.end()) return false; // no one pair for same kfindex exist;
  227. return (it->second == (*pSample)); // check by equality operator
  228. return false;
  229. }
  230. //----------------------------------------------------------------------------------------------------------------------------------------
  231. bool LMatchingFilter::mcIsMaybeMatching(const MpdTofMatchingData *pSample)
  232. {
  233. Int_t kfindex = pSample->GetKFTrackIndex();
  234. auto it = mcMaybeMdata.find(kfindex); // pair<kfindex, MdataT>
  235. if(it == mcMaybeMdata.end()) return false; // no one pair for same kfindex exist;
  236. int counter = mcMaybeMdata.count(kfindex);
  237. for(int i=0; i<counter;i++, it++) // cycle by pairs for same kfindex
  238. {
  239. if(it->second == (*pSample)) return true;
  240. }
  241. return false;
  242. }
  243. //----------------------------------------------------------------------------------------------------------------------------------------
  244. Int_t LMatchingFilter::Processing(Int_t nKFTracks, double& chi2)
  245. {
  246. fChi2 = 0.;
  247. Int_t candNmb = mmCandidate.size();
  248. #ifdef DoTest
  249. std::cout<<"\n -I- LMatchingFilter::Processing --------------------------------------------------------- mmCandidate.size()="<<mmCandidate.size()<<std::flush;
  250. #endif
  251. Int_t nAccepted, iterNmb = 0;
  252. if(mmCandidate.empty()) goto finish;
  253. UpdateRanks();
  254. #ifdef DoTest
  255. PrintRanks(" FIRST ITER ");
  256. Commit();
  257. Dump();
  258. #endif
  259. ProcessSingleTracks(); // process only single hit (1th rank) candidate tracks
  260. iterNmb++;
  261. #ifdef DoTest
  262. Status("AFTER FIRST ProcessSingleTracks");
  263. #endif
  264. if(mmCandidate.empty()) goto finish;
  265. newIteration:
  266. iterNmb++;
  267. UpdateRanks();
  268. #ifdef DoTest
  269. PrintRanks(" new ITER ");
  270. int nn = GetMinRankAfter(10); if (nn>0) std::cout<<"\n MInRank = "<<nn<<" count()= "<<mmRanks.count(nn);
  271. if( mmCandidate.size() <=20) Dump();
  272. #endif
  273. nAccepted = ProcessSingleTracks();
  274. if(nAccepted == 0)
  275. {
  276. ProcessNRankTrack(GetMinRank());
  277. }
  278. if(mmCandidate.empty()) goto finish;
  279. if(iterNmb > 1000)
  280. {
  281. FairLogger::GetLogger()->Warning(MESSAGE_ORIGIN, " <TofMatchingFilter::Processing> Too many tries.");
  282. goto finish;
  283. }
  284. goto newIteration;
  285. finish:
  286. chi2 = fChi2;
  287. if(pQA) pQA->FillCandidateNumber(candNmb, iterNmb);
  288. #ifdef DoTest
  289. std::cout<<"\n -I- [TofMatchingFilter::Processing] FINAL SIZE: mmLinks="<<mmLinks.size()<<" mmCandidate="<<mmCandidate.size()<<" mmRanks="<<mmRanks.size()<<"\n";
  290. #endif
  291. if(fVerbose) std::cout<<" -I- [TofMatchingFilter::Processing] Finished with "<<iterNmb<<" iterations.\n";
  292. return mmAccepted.size();
  293. }
  294. //----------------------------------------------------------------------------------------------------------------------------------------
  295. int LMatchingFilter::ProcessSingleTracks() // only single hit (1th rank) candidate tracks processing
  296. {
  297. auto range = mmRanks.equal_range(1);
  298. if(std::distance(range.first, range.second) == 0) return 0; // no 1th rank candidates
  299. #ifdef DoTest
  300. std::cout<<"\n LMatchingFilter::ProcessSingleTracks_size="<<std::distance(range.first, range.second);
  301. #endif
  302. Int_t MatchingOK = 0;
  303. for(auto iter = range.first; iter != range.second; ++iter) // cycle by ONLY 1th rank candidates
  304. {
  305. Int_t KfIndex = iter->second;
  306. candIter itCand = mmCandidate.find(KfIndex);
  307. #ifdef DoTest
  308. std::cout<<"\nProcessSingleTracks_KfIndex="<<KfIndex;
  309. #endif
  310. if(itCand == mmCandidate.end()) continue; // already removed by competitor track candidate
  311. #ifdef DoTest
  312. std::cout<<" EXIST";
  313. #endif
  314. Int_t hitIndex = itCand->second->GetTofHitIndex();
  315. Int_t counter = mmLinks.count(hitIndex); // number of competitor tracks for this hit
  316. if(counter == 1)// don't have competitor tracks
  317. {
  318. /// candHist.ResetFlag(CandHistory< MdataT>::k1thRank | CandHistory< MdataT>::kNoCompetitor);
  319. /// candHist.pCandidate = itCand->second;
  320. AcceptCandidate(itCand);
  321. MatchingOK++;
  322. }
  323. else // accept best(closest) 1th rank candidate
  324. {
  325. Double_t delta;
  326. candIter it = FindClosest1hRankCandidate(mmLinks.equal_range(hitIndex), delta);
  327. /// candHist.pCandidate = it->second;
  328. AcceptCandidate(it);
  329. MatchingOK++;
  330. }
  331. }
  332. #ifdef DoTest
  333. std::cout<<"\n LMatchingFilter::ProcessSingleTracks 1th rank cands="<<std::distance(range.first, range.second)<<" MatchingOK="<<MatchingOK;
  334. #endif
  335. return MatchingOK;
  336. }
  337. //----------------------------------------------------------------------------------------------------------------------------------------
  338. double LMatchingFilter::FindMinDeltaForTrack(Int_t KfIndex, Int_t disabledHitIndex, candIter& minCand, bool& IsWithoutCompetition)
  339. {
  340. Double_t delta, minDelta = 1.e+10; // big value
  341. auto rangeCands = mmCandidate.equal_range(KfIndex);
  342. minCand = mmCandidate.end();
  343. IsWithoutCompetition = true;
  344. assert(std::distance(rangeCands.first, rangeCands.second) != 0);
  345. for(candIter iter = rangeCands.first; iter != rangeCands.second; ++iter) // cycle by candidates for this KfIndex
  346. {
  347. int hitindex = iter->second->GetTofHitIndex();
  348. int nmbTrHits = iter->second->GetNmbTrHits();
  349. if(disabledHitIndex != hitindex) // skip <KfIndex, disabledHitIndex> candidate
  350. {
  351. delta = iter->second->GetDelta();
  352. if(delta < minDelta) // new local minimum
  353. {
  354. IsWithoutCompetition = (1 == mmLinks.count(hitindex)) ? true : false;
  355. minCand = iter;
  356. minDelta = delta;
  357. }
  358. }
  359. }
  360. assert(minCand != mmCandidate.end()); // don't have local minimum ????
  361. return minDelta;
  362. }
  363. //----------------------------------------------------------------------------------------------------------------------------------------
  364. void LMatchingFilter::ProcessNRankTrack(Int_t rank) // process Nth rank ONE track
  365. {
  366. auto rangeRanks = mmRanks.equal_range(rank);
  367. if(std::distance(rangeRanks.first, rangeRanks.second) == 0) return; // no Nth rank candidates
  368. // select random(first) KfIndex from Nrank list
  369. Int_t KfIndex = rangeRanks.first->second; // track index for processing
  370. #ifdef DoTest
  371. std::cout<<"\n ---------------ProcessNRankTrack---- rank="<<rank<<" mmRanks.count(rank)="<<std::distance(rangeRanks.first, rangeRanks.second)<<" KfIndex="<<KfIndex<<" mmCandidate.count(KfIndex)="<<mmCandidate.count(KfIndex);
  372. #endif
  373. bool IsWithoutCompetition = true;
  374. candIter cand, bestCand;
  375. Double_t minDelta = FindMinDeltaForTrack(KfIndex, -100, bestCand, IsWithoutCompetition);
  376. int maxKfTrackWeight = bestCand->second->GetNmbTrHits();
  377. if(IsWithoutCompetition)
  378. {
  379. assert(cand != mmCandidate.end());
  380. AcceptCandidate(bestCand);
  381. }
  382. else
  383. {
  384. Int_t hitIndex = bestCand->second->GetTofHitIndex();
  385. auto rangeLinks = mmLinks.equal_range(hitIndex);
  386. for(linkIter it = rangeLinks.first; it != rangeLinks.second; ++it) // cycle by tracks for this hitIndex
  387. {
  388. if(it->second == bestCand) continue; // skip himself
  389. Int_t kfindex = it->second->first;
  390. if(IsSameRank(rank, kfindex)) // check only same rank candidates
  391. {
  392. Double_t delta = FindMinDeltaForTrack(kfindex, -100, cand, IsWithoutCompetition);
  393. int kftrackweight = cand->second->GetNmbTrHits();
  394. if(kftrackweight > maxKfTrackWeight) // new track have bigger weight -> set as best candidate
  395. {
  396. bestCand = cand;
  397. minDelta = delta;
  398. }
  399. else if(kftrackweight == maxKfTrackWeight)
  400. {
  401. if(delta < minDelta) // new track have same weight but smaller delta -> set as best candidate
  402. {
  403. bestCand = cand;
  404. minDelta = delta;
  405. }
  406. }
  407. }
  408. } // cycle by tracks for this hitIndex
  409. assert(cand != mmCandidate.end());
  410. AcceptCandidate(bestCand);
  411. }
  412. }
  413. //----------------------------------------------------------------------------------------------------------------------------------------
  414. LMatchingFilter::candIter LMatchingFilter::FindClosest1hRankCandidate(const std::pair< linkIter, linkIter>& range, Double_t& minDelta)
  415. {
  416. /////////// candHist.ResetFlag(CandHistory< MdataT>::k1thRank | CandHistory< MdataT>::kHaveCompetitor);
  417. candIter retvalue = mmCandidate.end();
  418. minDelta = 1.e+10; // big value
  419. for(linkIter it = range.first; it != range.second; it++) // cycle by competitor candidates
  420. {
  421. Double_t delta = it->second->second->GetDelta();
  422. ///////////// candHist.vCompetitors.push_back(it->second->second);
  423. if(delta < minDelta) // best (closest) candidate
  424. {
  425. retvalue = it->second;
  426. minDelta = delta;
  427. }
  428. }
  429. assert(retvalue != mmCandidate.end());
  430. return retvalue;
  431. }
  432. //----------------------------------------------------------------------------------------------------------------------------------------
  433. void LMatchingFilter::AcceptCandidate(candIter itCand)
  434. {
  435. bool isTrueMatching = mcIsTrueMatching(itCand->second);
  436. bool isMaybeTrueMatching = mcIsMaybeMatching(itCand->second);
  437. #ifdef DoTest
  438. if(! (isTrueMatching || isMaybeTrueMatching) ) // mismatching accepted !!!
  439. {
  440. /// candHist.Print();
  441. //itCand->second->Print();
  442. }
  443. Commit();
  444. #endif
  445. // Accept candidate
  446. Int_t acceptedKfIndex = itCand->first, acceptedHitIndex = itCand->second->GetTofHitIndex();
  447. mmAccepted.insert(make_pair(acceptedKfIndex, itCand->second)); // pair< KfIndex, MdataT*>
  448. double delta = itCand->second->GetDelta();
  449. fChi2 += delta*delta;
  450. auto rangeCand = mmCandidate.equal_range(acceptedKfIndex);
  451. for(candIter iter = rangeCand.first, iterEnd = rangeCand.second; iter != iterEnd; iter++) // all candidates for acceptedKfIndex
  452. {
  453. int hitindex = iter->second->GetTofHitIndex();
  454. // remove all links to hits for this KfIndex (track erase)
  455. RemoveLink(hitindex, acceptedKfIndex);
  456. }
  457. auto rangeLinks = mmLinks.equal_range(acceptedHitIndex);
  458. for(linkIter iter = rangeLinks.first, iterEnd = rangeLinks.second; iter != iterEnd; iter++) // all candidates for acceptedHitIndex
  459. {
  460. int kfindex = iter->second->first;
  461. // remove all candidate tracks for this hitIndex (hit erase)
  462. RemoveCandidate(acceptedHitIndex, kfindex);
  463. }
  464. // remove all candidates for this KfIndex
  465. rangeCand = mmCandidate.equal_range(acceptedKfIndex);
  466. mmCandidate.erase(rangeCand.first, rangeCand.second);
  467. // remove all links for this hitIndex hit (hit erase)
  468. rangeLinks = mmLinks.equal_range(acceptedHitIndex);
  469. mmLinks.erase(rangeLinks.first, rangeLinks.second);
  470. #ifdef DoTest
  471. bool pass = CheckAdequacy("AFTER ACCEPT");
  472. if(!pass)
  473. {
  474. Status(" change AFTER ACCEPT"); std::cout<<std::flush;
  475. }
  476. assert(pass);
  477. #else
  478. assert(CheckAdequacy("AFTER ACCEPT"));
  479. #endif
  480. }
  481. //----------------------------------------------------------------------------------------------------------------------------------------
  482. bool LMatchingFilter::CheckAdequacy(const char* comment)
  483. {
  484. TmmIntT<candIter> mmLinksTmp;
  485. // for(const auto& entry : mmCandidate) mmLinksTmp.insert(make_pair(entry.second->GetTofHitIndex(), mmCandidate.find(entry.first)));
  486. for(candIter it = mmCandidate.begin(), itEnd = mmCandidate.end(); it != itEnd; it++) mmLinksTmp.insert(make_pair(it->second->GetTofHitIndex(), it));
  487. bool equalFront = true;// (mmLinksTmp == mmLinks); // iterators point to another mmaps, so operator == return always true
  488. if(mmLinksTmp.size() == mmLinks.size())
  489. {
  490. for(const auto& iter : mmLinksTmp)
  491. {
  492. bool founded = false;
  493. auto rangeLinks = mmLinks.equal_range(iter.first);
  494. for(auto it = rangeLinks.first; it != rangeLinks.second; ++it) // cycle by tracks for this hitIndex
  495. {
  496. if(it->second->second == iter.second->second) // MdataT*(mmLinksTmp) == MdataT* (mmLinks)
  497. {
  498. founded = true; break;
  499. }
  500. }
  501. if(false == founded){ equalFront = false; break; }
  502. }
  503. }
  504. else equalFront = false;
  505. TmmCand mmCandidateTmp;
  506. for(const auto& entry : mmLinks) mmCandidateTmp.insert(make_pair(entry.second->first, entry.second->second));
  507. bool equalBack = true; // (mmCandidateTmp == mmCandidate); // because of a different order filling the entries may be mismatch, so operator == may be return wrong false
  508. if(mmCandidateTmp.size() == mmCandidate.size())
  509. {
  510. for(const auto& iter : mmCandidateTmp)
  511. {
  512. bool founded = false;
  513. auto range = mmCandidate.equal_range(iter.first);
  514. for(auto it = range.first; it != range.second; ++it) // cycle by candidates for same KfIndex
  515. {
  516. if(it->second == iter.second) // MdataT*(mmCandidateTmp) == MdataT* (mmCandidate)
  517. {
  518. founded = true; break;
  519. }
  520. }
  521. if(false == founded){ equalFront = false; break; }
  522. }
  523. }
  524. else equalBack = false;
  525. if( !equalFront || !equalBack) // print report, if it's problems
  526. {
  527. std::cout<<"\n CheckAdequacy: ";
  528. if(comment) std::cout<<comment;
  529. std::cout<<" mmCandidate.size="<<mmCandidate.size()<<" mmLinks.size="<<mmLinks.size()<<" (mmLinks ---> mmCandidate) = "<<equalBack<<", (mmCandidate ---> mmLinks) = "<<equalFront;
  530. }
  531. return (equalFront && equalBack);
  532. }
  533. //----------------------------------------------------------------------------------------------------------------------------------------
  534. void LMatchingFilter::RemoveLink(Int_t hitIndex, Int_t KfIndex)
  535. {
  536. auto range = mmLinks.equal_range(hitIndex); // links for same hitIndex
  537. for(linkIter it = range.first, itEnd = range.second; it != itEnd; ) // iterating multimap && erasing
  538. {
  539. if(it->second->first == KfIndex) mmLinks.erase(it++);
  540. else ++it;
  541. }
  542. }
  543. //----------------------------------------------------------------------------------------------------------------------------------------
  544. void LMatchingFilter::RemoveCandidate(Int_t hitIndex, Int_t KfIndex)
  545. {
  546. auto range = mmCandidate.equal_range(KfIndex);
  547. for(candIter it = range.first, itEnd = range.second; it != itEnd; ) // iterating multimap && erasing
  548. {
  549. if(it->second->GetTofHitIndex() == hitIndex) mmCandidate.erase(it++);
  550. else ++it;
  551. }
  552. }
  553. //----------------------------------------------------------------------------------------------------------------------------------------
  554. void LMatchingFilter::UpdateRanks()
  555. {
  556. mmRanks.clear();
  557. for(auto it = mmCandidate.begin(), end = mmCandidate.end(); it != end; it = mmCandidate.upper_bound(it->first)) // cycle by unique KfIndex
  558. {
  559. mmRanks.insert(make_pair(mmCandidate.count(it->first ), it->first)); // pair <count(KfIndex), KfIndex>
  560. }
  561. }
  562. //----------------------------------------------------------------------------------------------------------------------------------------
  563. void LMatchingFilter::Reset()
  564. {
  565. mmCandidate.clear();
  566. mmAccepted.clear();
  567. mmLinks.clear();
  568. candHist.Clear();
  569. }
  570. //----------------------------------------------------------------------------------------------------------------------------------------
  571. void LMatchingFilter::SetContainer(TClonesArray *array)
  572. {
  573. aMdata = array;
  574. isOwner = false;
  575. }
  576. //----------------------------------------------------------------------------------------------------------------------------------------
  577. void LMatchingFilter::Commit()
  578. {
  579. mmCandidateSnapshot = mmCandidate;
  580. mmAcceptedSnapshot = mmAccepted;
  581. //mmLinksSnapshot = mmLinks; ERROR: DON'T USE THIS: mmlinks used iterators to mmCandidate, so it may be singular iterator!!!!
  582. mmLinksSnapshot.clear();
  583. candIter iter;
  584. for(const auto& it : mmLinks)
  585. {
  586. Int_t hitIndex = it.first;
  587. iter = FindCandidate(&mmCandidateSnapshot, hitIndex, it.second->second->GetKFTrackIndex()); // mmLinksSnapshot used iterators to mmCandidateSnapshot NOW
  588. mmLinksSnapshot.insert(make_pair(hitIndex, iter));
  589. }
  590. }
  591. //----------------------------------------------------------------------------------------------------------------------------------------
  592. void LMatchingFilter::Status(const char *comment, std::ostream& os) const
  593. {
  594. os<<"\n -------------------------------------------------STATUS----- "; if(comment != nullptr) os<<comment;
  595. os<<"\n mmCandidate added:";
  596. PrintChanges(mmCandidate, mmCandidateSnapshot, os);
  597. os<<"\n mmCandidate removed:";
  598. PrintChanges(mmCandidateSnapshot, mmCandidate, os);
  599. os<<"\n mmAccepted added:";
  600. PrintChanges(mmAccepted, mmAcceptedSnapshot, os);
  601. os<<"\n mmAccepted removed:";
  602. PrintChanges(mmAcceptedSnapshot, mmAccepted, os);
  603. os<<"\n mmLinks added:";
  604. PrintChanges(mmLinks, mmLinksSnapshot, os);
  605. os<<"\n mmLinks removed:";
  606. PrintChanges(mmLinksSnapshot, mmLinks, os);
  607. os<<"\n ------------------------------------------------------------"<<std::flush;
  608. }
  609. //----------------------------------------------------------------------------------------------------------------------------------------
  610. void LMatchingFilter::PrintChanges(const TmmCand& src, const TmmCand& checked, std::ostream& os) const
  611. {
  612. if(src == checked)
  613. {
  614. os<<"\t nothing.";
  615. return;
  616. }
  617. Int_t KfIndex, hitIndex; bool exist;
  618. for(const auto& it : src)
  619. {
  620. KfIndex = it.first;
  621. auto range = checked.equal_range(KfIndex); // cands for same KfIndex
  622. exist = false;
  623. hitIndex = it.second->GetTofHitIndex();
  624. for(auto iter = range.first; iter != range.second; ++iter) // cycle by same KfIndex entries from checked container
  625. {
  626. if(hitIndex == iter->second->GetTofHitIndex()) // same hitIndex found
  627. {
  628. exist = true;
  629. break;
  630. }
  631. }
  632. if(!exist) os<<"\n KfIndex="<<KfIndex<<", hitIndex="<<hitIndex;
  633. }
  634. }
  635. //------------------------------------------------------------------------------------------------------------------------
  636. void LMatchingFilter::PrintChanges(const TmmIntT<candIter>& src, const TmmIntT<candIter>& checked, std::ostream& os) const
  637. {
  638. if(src == checked)
  639. {
  640. os<<"\t nothing.";
  641. return;
  642. }
  643. Int_t KfIndex, hitIndex; bool exist;
  644. for(const auto& it : src)
  645. {
  646. hitIndex = it.first;
  647. auto range = checked.equal_range(hitIndex); // links for same hitIndex
  648. exist = false;
  649. KfIndex = it.second->second->GetKFTrackIndex();
  650. for(auto iter = range.first; iter != range.second; ++iter) // cycle by same hitIndex entries from checked container
  651. {
  652. if(KfIndex == iter->second->second->GetKFTrackIndex()) // same KfIndex found
  653. {
  654. exist = true;
  655. break;
  656. }
  657. }
  658. if(!exist) os<<"\n KfIndex="<<KfIndex<<", hitIndex="<<hitIndex;
  659. }
  660. }
  661. //------------------------------------------------------------------------------------------------------------------------
  662. void LMatchingFilter::Dump(const char* comment, std::ostream& os) const
  663. {
  664. os<<"\n ---------------------- DUMP --------------- "; if(comment) os<<comment; os<<" ------------ entries: "<<mmCandidate.size();
  665. Int_t KfIndex, hitIndex; bool exist; int n = 0;
  666. for(auto it = mmCandidate.begin(); it != mmCandidate.end(); )
  667. {
  668. KfIndex = it->first;
  669. os<<"\n "<<++n<<") KfIndex="<<KfIndex<<" --- ";
  670. int counter = mmCandidate.count(KfIndex);
  671. for(int hit=0; hit<counter; hit++)
  672. {
  673. hitIndex = it->second->GetTofHitIndex();
  674. os<<" hitIndex="<< hitIndex;
  675. auto iter = mmLinks.find(hitIndex);
  676. if(iter != mmLinks.end())
  677. {
  678. int counter2 = mmLinks.count(hitIndex);
  679. if(counter2 > 1)os<<" (";
  680. for(int track=0; track<counter2;track++)
  681. {
  682. if(counter2 > 1)os<<" "<<iter->second->second->GetKFTrackIndex()<<" ";
  683. ++iter;
  684. }
  685. if(counter2 > 1)os<<") ";
  686. }
  687. ++it;
  688. }
  689. }
  690. os<<"\n ---------------------- DUMP --------------------------------------"<<std::flush;
  691. }
  692. //------------------------------------------------------------------------------------------------------------------------
  693. void LMatchingFilter::PrintRanks(const char* comment, std::ostream& os) const
  694. {
  695. os<<std::endl;
  696. if(comment) os<<comment;
  697. os<<"\n mmRanks.count(n) 1= "<<mmRanks.count(1)<<", 2= "<<mmRanks.count(2)<<", 3= "<<mmRanks.count(3)<<", 4= "<<mmRanks.count(4)<<", 5= "<<mmRanks.count(5)
  698. <<", 6= "<<mmRanks.count(6)<<", 7= "<<mmRanks.count(7)<<", 8= "<<mmRanks.count(8)<<", 9= "<<mmRanks.count(9)<<", 10= "<<mmRanks.count(10)<<" cand. total: "<< mmCandidate.size();
  699. }
  700. //------------------------------------------------------------------------------------------------------------------------
  701. #endif