MpdItsToTpcMatching.cxx 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635
  1. #include "MpdItsToTpcMatching.h"
  2. #include "MpdKalmanFilter.h"
  3. #include "MpdKalmanTrack.h"
  4. #include "MpdMCTrack.h"
  5. #include "MpdItsHit5spd.h"
  6. #include "MpdStsPoint.h"
  7. #include "MpdItsKalmanTrack.h"
  8. #include "MpdTpcKalmanTrack.h"
  9. #include "MpdVector.h"
  10. #include "MpdVectorFinder.h"
  11. #include "FairRun.h"
  12. #include "FairRunAna.h"
  13. #include "FairRuntimeDb.h"
  14. #include "FairRootManager.h"
  15. #include <iostream>
  16. #include <map>
  17. #include <set>
  18. #include <string>
  19. using std::cout;
  20. using std::endl;
  21. std::map<Int_t,MpdKalmanTrack*> fmap;
  22. std::ofstream fout("match.dat");
  23. std::ofstream fout_v_match("log_match.txt");
  24. //__________________________________________________________________________
  25. MpdItsToTpcMatching::MpdItsToTpcMatching(const char *name, Int_t iVerbose )
  26. :FairTask(name, iVerbose),
  27. fExact(0)
  28. //fNPass(2),
  29. {
  30. // AZ fItsTracks = new TClonesArray("MpdVector", 100); /// not needed?
  31. fTracks = new TClonesArray("MpdItsKalmanTrack", 100);
  32. fTracks1 = new TClonesArray("MpdTpcKalmanTrack", 100);
  33. fTracksRefit = new TClonesArray("MpdItsKalmanTrack", 100);
  34. fTpcTracksRefit = new TClonesArray("MpdItsKalmanTrack", 100); /// was tpc
  35. }
  36. //__________________________________________________________________________
  37. MpdItsToTpcMatching::~MpdItsToTpcMatching()
  38. {
  39. }
  40. //__________________________________________________________________________
  41. InitStatus MpdItsToTpcMatching::Init()
  42. {
  43. //return ReInit();
  44. if (ReInit() != kSUCCESS) return kERROR;
  45. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  46. LOG(INFO) << __func__ << (rtdb);
  47. //FillGeoScheme(); /// TODO
  48. return kSUCCESS;
  49. }
  50. //__________________________________________________________________________
  51. /*
  52. void MpdItsToTpcMatching::FillGeoScheme()
  53. {
  54. }
  55. */
  56. //__________________________________________________________________________
  57. InitStatus MpdItsToTpcMatching::ReInit()
  58. {
  59. fItsPoints = (TClonesArray *) FairRootManager::Instance()->GetObject("StsPoint");
  60. fItsTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("ItsTrack");
  61. //if (fItsPoints == 0x0 || fItsHits == 0x0) return kERROR;
  62. fTpcTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  63. fMCTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  64. fKHits = (TClonesArray *) FairRootManager::Instance()->GetObject("ItsKHits");
  65. ///FairRootManager::Instance()->Register("ItsTrackTo27", "ItsCopy", fTracks1, kTRUE); /// originally fTracks
  66. FairRootManager::Instance()->Register("ItsTrackRefit", "ItsRefit", fTracksRefit, kTRUE); /// originally fTracks
  67. FairRootManager::Instance()->Register("TpcTrackRefit", "TpcRefit", fTpcTracksRefit, kTRUE); /// originally fTracks
  68. return kSUCCESS;
  69. }
  70. //__________________________________________________________________________
  71. void MpdItsToTpcMatching::SetParContainers()
  72. {
  73. // Get run and runtime database
  74. FairRunAna* run = FairRunAna::Instance();
  75. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  76. FairRuntimeDb* db = run->GetRuntimeDb();
  77. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  78. // Get STS geometry parameter container
  79. db->getContainer("MpdStsGeoPar");
  80. }
  81. //__________________________________________________________________________
  82. void MpdItsToTpcMatching::Finish()
  83. {
  84. //Write();
  85. //if (lunErr) fclose(lunErr);
  86. }
  87. //__________________________________________________________________________
  88. void MpdItsToTpcMatching::Reset()
  89. {
  90. cout << " MpdItsToTpcMatching::Reset " << endl;
  91. fout_v_match << " MpdItsToTpcMatching::Reset " << endl;
  92. fTracks1->Delete(); // AZ
  93. fTpcTracksRefit->Delete(); // AZ
  94. fTracks->Delete();
  95. fTracksRefit->Delete();
  96. }
  97. //__________________________________________________________________________
  98. void MpdItsToTpcMatching::GetItsTracks(multimap <Float_t, MpdItsKalmanTrack*> &multimapIts,
  99. multimap <Float_t, MpdItsKalmanTrack*> &multimapItsPhi)
  100. {
  101. /// Get ITS tracks linked to their outer hit positions (rphi and z)
  102. Int_t n = fItsTracks->GetEntriesFast();
  103. //Float_t phi;
  104. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  105. for (Int_t i = 0; i < n; i++) {
  106. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fItsTracks->UncheckedAt(i);
  107. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetTrHits()->First();
  108. TVector3 pos = geo->GlobalPos(hit);
  109. cout << pos.Pt() << " -aaaaaaaaaa- " << pos.Z() << endl;
  110. multimapIts.insert(pair<Float_t,MpdItsKalmanTrack*>(pos.Z(),track));
  111. }
  112. }
  113. //__________________________________________________________________________
  114. void MpdItsToTpcMatching::RefitItsTo27(multimap <Float_t, MpdItsKalmanTrack*> &multimapIts,
  115. multimap <Float_t, MpdItsKalmanTrack*> &multimapItsPhi)
  116. {
  117. /// Refit uses MpdKalmanTrack while fItsTracks has MpdItsKalmanTrack
  118. Int_t n = fItsTracks->GetEntriesFast();
  119. Float_t phi;
  120. cout << "fItsTracks " << n << endl;
  121. fout_v_match << "fItsTracks " << n << endl;
  122. cout << "refit ITS" << endl;
  123. fout_v_match << "refit ITS" << endl;
  124. for (Int_t i = 0; i < n; i++) {
  125. MpdItsKalmanTrack *temp = (MpdItsKalmanTrack*)fItsTracks->UncheckedAt(i);
  126. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) new ((*fTracks)[i]) MpdItsKalmanTrack(*temp);
  127. ///cout << "Track id's " << temp->GetTrackID() << " " << track->GetTrackID() << endl;
  128. ///cout << "refitting nofhits noftrhits getentries nofits" << track->GetNofHits() << " " << track->GetNofTrHits() << " " << track->GetTrHits()->GetEntriesFast() << " " << track->GetNofIts() << endl;
  129. track->SetDirection(MpdKalmanTrack::kOutward);
  130. MpdKalmanFilter::Instance()->Refit(track, -1, 0); // TODO check params /// -1
  131. /// Propagate ITS to 27 Outward
  132. MpdKalmanHit hitTmp;
  133. hitTmp.SetType(MpdKalmanHit::kFixedR);
  134. hitTmp.SetDist(27.0);
  135. Bool_t ok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hitTmp, kTRUE, kFALSE); // doesn't create GetMeas() coordinates
  136. /// writing to ItsTracksRefit
  137. new ((*fTracksRefit)[i]) MpdItsKalmanTrack(*track);
  138. TMatrixD *parNew = track->GetParamNew();
  139. /// (*parNew)(1,0) -- theta, (*parNew)(0,0) / track->GetPosNew() -- phi
  140. phi = (*parNew)(0,0) / track->GetPosNew();
  141. //AZ - save ITS track
  142. fmap[track->GetTrackID()] = track;
  143. //AZ
  144. multimapIts.insert(pair<Float_t, MpdItsKalmanTrack*>((*parNew)(1,0) , track));
  145. ///multimapItsPhi.insert(pair<Float_t, MpdItsKalmanTrack*>(phi, track));
  146. multimapItsPhi.insert(pair<Float_t, MpdItsKalmanTrack*>((*parNew)(0,0), track));
  147. if (TMath::Abs(TMath::Pi() - TMath::Abs(phi)) < TMath::Pi() / 9.0) {
  148. multimapItsPhi.insert(pair<Float_t, MpdItsKalmanTrack*>((*parNew)(0,0) - TMath::Sign(1, phi) * 2 * TMath::Pi() * track->GetPosNew(), track));
  149. ///multimapItsPhi.insert(pair<Float_t, MpdItsKalmanTrack*>(phi - TMath::Sign(1, phi) * 2 * TMath::Pi(), track)); /// transverse, phi, duplicate hit
  150. }
  151. }
  152. }
  153. //__________________________________________________________________________
  154. void MpdItsToTpcMatching::RefitTpcTo27(multimap <Float_t, MpdTpcKalmanTrack*> &multimapTpc,
  155. multimap <Float_t, MpdTpcKalmanTrack*> &multimapTpcPhi)
  156. {
  157. // Propagate TPC track ro R = 27 cm
  158. Int_t n = fTpcTracks->GetEntriesFast();
  159. fTpcIndSet.clear();
  160. Float_t phi;
  161. cout << "fTpcTracks " << n << endl;
  162. cout << "refit TPC" << endl;
  163. fout_v_match << "fTpcTracks " << n << endl;
  164. fout_v_match << "refit TPC" << endl;
  165. fout_v_match << "id, dphi, drphi, dz, pt, tpc-its" << endl;
  166. for (Int_t i = 0; i < n; i++) {
  167. MpdTpcKalmanTrack *temp = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(i);
  168. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) new ((*fTracks1)[i]) MpdTpcKalmanTrack(*temp); /// TODO explain why this is needed
  169. /// TObject.UniqueID is used to identify track after it is refitted
  170. /// @temp is original track, @track is to be refitted
  171. track->SetUniqueID(i+1);
  172. //AZ track->SetDirection(MpdKalmanTrack::kInward);
  173. track->SetDirection(MpdKalmanTrack::kOutward);
  174. /// Propagate TPC track to 27 cm inward (fPos == 27)
  175. MpdKalmanHit hitTmp;
  176. hitTmp.SetType(MpdKalmanHit::kFixedR);
  177. //AZ hitTmp.SetPos(27.0);
  178. hitTmp.SetPos(track->GetPos()); // AZ: fPos ~= 27 cm
  179. /// see MpdTrackFinderIts5spd::GetTrackSeeds
  180. track->SetParamNew(*track->GetParam());
  181. track->SetPos(track->GetPosNew());
  182. track->ReSetWeight();
  183. TMatrixDSym w = *track->GetWeight(); // save current weight matrix
  184. Bool_t ok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hitTmp, kFALSE, kFALSE);
  185. track->SetWeight(w); // restore original weight matrix (near TPC inner shell)
  186. track->GetHits()->Clear();
  187. track->SetDirection(MpdKalmanTrack::kInward); //AZ
  188. //track->SetChi2Its(track->GetChi2()); // temporary storage
  189. //track->SetChi2(0.); /// commented 27.6.2019
  190. //Int_t j = track->GetTrHits()->GetEntriesFast() - 1;
  191. TMatrixD *parNew = track->GetParamNew();
  192. //cout << (*parNew)(0,0) << " " << (*parNew)(1,0) << " " << (*parNew)(2,0) << " " << (*parNew)(3,0) << " " << (*parNew)(4,0) << " p " << track->GetPosNew() << endl;
  193. /// (*parNew)(1,0) -- theta, (*parNew)(0,0) / track->GetPosNew() -- phi
  194. phi = (*parNew)(0,0) / track->GetPosNew();
  195. // AZ - write parameters
  196. if (fmap.find(track->GetTrackID()) != fmap.end()) {
  197. Int_t id = track->GetTrackID();
  198. fout << id << " " << phi-fmap[id]->GetParamNew(0)/track->GetPosNew() << " " << (*parNew)(0,0) - fmap[id]->GetParamNew(0) << " "
  199. << (*parNew)(1,0)-fmap[id]->GetParamNew(1) << " " << track->Pt() << endl;
  200. fout_v_match << id << " " << phi-fmap[id]->GetParamNew(0)/track->GetPosNew() << " " << (*parNew)(0,0) - fmap[id]->GetParamNew(0) << " "
  201. << (*parNew)(1,0)-fmap[id]->GetParamNew(1) << " " << track->Pt() << endl;
  202. }
  203. //AZ
  204. multimapTpc.insert(pair<Float_t, MpdTpcKalmanTrack*>((*parNew)(1,0) , track));
  205. multimapTpcPhi.insert(pair<Float_t, MpdTpcKalmanTrack*>((*parNew)(0,0), track)); /// was (phi, track)
  206. ///multimapTpcPhi.insert(pair<Float_t, MpdTpcKalmanTrack*>(phi, track));
  207. if (TMath::Abs(TMath::Pi() - phi) < TMath::Pi() / 9.0)
  208. ///multimapTpcPhi.insert(pair<Float_t, MpdTpcKalmanTrack*>(phi - TMath::Sign(1, phi) * 2 * TMath::Pi(), track));
  209. multimapTpcPhi.insert(pair<Float_t, MpdTpcKalmanTrack*>((*parNew)(0,0) - TMath::Sign(1, phi) * 2 * TMath::Pi() * track->GetPosNew(), track));
  210. fTpcIndSet.insert(i); // store TPC track indices
  211. }
  212. }
  213. //__________________________________________________________________________
  214. void MpdItsToTpcMatching::Exec(Option_t * option)
  215. {
  216. fout_v_match << "- track matching -" << endl;
  217. fmap.clear();
  218. //if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  219. Reset();
  220. // MpdKalmanFilter function Refit
  221. //Bool_t MpdKalmanFilter::Refit(MpdKalmanTrack *track, Int_t iDir, Bool_t updLeng)
  222. multimap <Float_t, MpdItsKalmanTrack*> multimapIts, multimapItsPhi;
  223. multimap <Float_t, MpdTpcKalmanTrack*> multimapTpc, multimapTpcPhi;
  224. RefitItsTo27(multimapIts, multimapItsPhi);
  225. //GetItsTracks(multimapIts, multimapItsPhi);
  226. RefitTpcTo27(multimapTpc, multimapTpcPhi);
  227. //Int_t n = fItsTracks->GetEntriesFast();
  228. multimap<Float_t, MpdItsKalmanTrack*>::iterator itits;
  229. multimap<Float_t, MpdTpcKalmanTrack*>::iterator ittpc;
  230. Float_t epsz = 2.5; //0.5;/// TODO should it depend on pt as well?
  231. Float_t epsphi = 2.5; //0.5;//0.2 for phi, 0.5 for rphi;
  232. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07};
  233. Int_t k = 0, tpck = 0;
  234. multimap <Float_t, std::tuple<MpdItsKalmanTrack*, MpdItsKalmanTrack*, MpdTpcKalmanTrack*>> multimapMatch; /// 1 was tpc
  235. /// added 9.10.2019
  236. set<MpdItsKalmanTrack*> tracksWithoutMatch;
  237. //________________________________________________MAIN MATCHING LOOP_____________________________________________________________
  238. for (itits = multimapIts.begin(); itits != multimapIts.end(); ++itits) {
  239. TMatrixD *parNew = (itits->second)->GetParamNew();
  240. ///phi = (*parNew)(0,0) / (itits->second)->GetPosNew();
  241. Float_t phi = (*parNew)(0,0);/// phi became rphi here
  242. //cout << "par " << itits->first << " " << phi << endl;
  243. multimap<Float_t, MpdTpcKalmanTrack*>::iterator itlow = multimapTpc.lower_bound(itits->first - epsz);
  244. multimap<Float_t, MpdTpcKalmanTrack*>::iterator ittop = multimapTpc.upper_bound(itits->first + epsz);
  245. multimap<Float_t, MpdTpcKalmanTrack*>::iterator itlowphi = multimapTpcPhi.lower_bound(phi - epsphi);
  246. multimap<Float_t, MpdTpcKalmanTrack*>::iterator ittopphi = multimapTpcPhi.upper_bound(phi + epsphi);
  247. /// getting window for possible tpc track matches
  248. set<MpdTpcKalmanTrack*> setz, setphi, intersect;
  249. for (multimap<Float_t, MpdTpcKalmanTrack*>::iterator itr = itlow; itr != ittop; ++itr) {
  250. setz.insert((*itr).second);
  251. }
  252. for (multimap<Float_t, MpdTpcKalmanTrack*>::iterator itr = itlowphi; itr != ittopphi; ++itr) {
  253. setphi.insert((*itr).second);
  254. }
  255. set_intersection(setz.begin(), setz.end(), setphi.begin(), setphi.end(), std::inserter(intersect, intersect.begin()));
  256. for (set<MpdTpcKalmanTrack*>::iterator it = intersect.begin(); it != intersect.end(); ++it) {
  257. MpdItsKalmanTrack* temp = new MpdItsKalmanTrack(*(*it)); /// TpcKalmanTrack is converted to ItsKalmanTrack
  258. //TMatrixD *parNew1 = temp->GetParamNew();
  259. // ITS hits are sorted "inward"
  260. TMatrixD param(5,1);
  261. TMatrixDSym weight(5), pointWeight(5);
  262. /// from MpvVectorFinder::AddHits()
  263. TClonesArray &trHits = *temp->GetTrHits();
  264. Int_t lastIndx = trHits.GetEntriesFast();
  265. //temp->SetChi2Its((itits->second)->GetChi2Its());
  266. temp->SetChi2Its(temp->GetChi2()); // tpc chi2 for now
  267. temp->SetChi2(0);
  268. tpck++;
  269. // AZ
  270. TVector3 mom3 = temp->Momentum3(), norm;
  271. mom3.SetMag(1.0);
  272. TString mass2 = "0.0194797849"; // pion mass squared
  273. /*
  274. if (fMCTracks) {
  275. // Get particle mass - ideal PID
  276. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(temp->GetTrackID());
  277. TParticlePDG *pdgP = TDatabasePDG::Instance()->GetParticle(mctrack->GetPdgCode());
  278. if (pdgP) {
  279. Double_t mass = pdgP->Mass();
  280. if (mass < 0.1 || mass > 0.25) {
  281. // Electrons or heavier than pions
  282. mass2 = "";
  283. mass2 += mass*mass;
  284. }
  285. }
  286. }
  287. */
  288. //AZ
  289. Int_t hitindex = 0;
  290. /// loop over each hit from its track
  291. for (Int_t i = 0; i < (itits->second)->GetHits()->GetEntriesFast(); i++) {
  292. MpdKalmanHit* hitTmp = ((MpdKalmanHit*) (itits->second)->GetHits()->UncheckedAt(i));
  293. ///cout << "temp " << temp->GetNofHits() << " " << temp->GetNofIts() << " " << temp->GetNofTrHits()<< endl;
  294. Bool_t ok = MpdKalmanFilter::Instance()->PropagateToHit(temp, hitTmp, kFALSE, kTRUE); // propagate tpc track to its hits
  295. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(temp, hitTmp, pointWeight, param);
  296. /// this is correct. you still propagate track to a given hit, even if hit is not added to the track afterwards because of Chi2
  297. if (dChi2 <= 20) {
  298. /// from MpdItsKalmanTrack::Refit() ???
  299. /// if hit is not included in track then track parameters and weight matrix are not updated
  300. temp->SetChi2(temp->GetChi2()+dChi2); // (*it)
  301. weight = *temp->GetWeight(); // (*it)
  302. weight += pointWeight;
  303. temp->SetWeight(weight); // (*it)
  304. temp->SetParamNew(param); // (*it)
  305. /// from AddHits()
  306. temp->GetHits()->Add(hitTmp); // this changes GetNofHits() result
  307. new (trHits[lastIndx+hitindex]) MpdKalmanHit(*hitTmp);
  308. hitindex++;
  309. } else {
  310. ///temp->SetChi2(temp->GetChi2() - dChi2); /// if hit is not added then Chi2 does not increase
  311. }
  312. //cout << ok /*<< " chi2 " << (*it)->GetChi2() */<< " ";
  313. //cout << " " << param(1,0) << endl;
  314. // AZ - Add multiple scattering in the sensor
  315. norm = MpdKalmanFilter::Instance()->GetGeo()->Normal(hitTmp);
  316. //Double_t step = 0.005 / TMath::Abs(norm * mom3) * 4.0; // extra factor 4. - possible overlaps
  317. Double_t step = thick[hitTmp->GetLayer()] / TMath::Abs(norm * mom3) * 4.0; // extra factor 4. - possible overlaps
  318. Double_t x0 = 9.36; // rad. length
  319. TMatrixDSym *cov = temp->Weight2Cov();
  320. Double_t th = temp->GetParamNew(3);
  321. Double_t cosTh = TMath::Cos(th);
  322. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(temp, x0, step, mass2);
  323. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  324. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  325. (*cov)(3,3) += angle2;
  326. Int_t iok = 0;
  327. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  328. temp->SetWeight(*cov);
  329. }
  330. /// swap GetChi2 - tpc, GetChi2Its - its
  331. Float_t c = temp->GetChi2Its();
  332. temp->SetChi2Its(temp->GetChi2());
  333. temp->SetChi2(c);
  334. ///quality
  335. Float_t qual = - (temp->GetNofHits() + (100.0 - TMath::Min(temp->GetChi2Its(), 100.0)) / 101.0); /// was temp->GetChi2() + temp->GetChi2Its()
  336. /// !!!! IMPORTANT - adding tuple to multimap
  337. /// tuple: 0 - propagated track, 1 - its track matched, 2 - tpc track matched
  338. multimapMatch.insert(pair<Float_t, std::tuple<MpdItsKalmanTrack*, MpdItsKalmanTrack*, MpdTpcKalmanTrack*> >(qual , std::make_tuple(temp, itits->second, *it)));
  339. } // for (set<MpdTpcKalmanTrack*>::iterator it
  340. if (intersect.size() == 0) {
  341. /// add its track standalone to final structure anyway through additional set
  342. k++;
  343. //(itits->second)->SetNofIts(5);
  344. tracksWithoutMatch.insert(itits->second);
  345. }
  346. } // for (itits = multimapIts.begin();
  347. ///sets for its and tpc tracks
  348. std::set<MpdTpcKalmanTrack*> tpcUnique;
  349. std::set<MpdItsKalmanTrack*> itsUnique;
  350. /// tuple: 0 - propagated track, 1 - its track matched, 2 - tpc track matched
  351. Int_t matchedCount = 0, qualTr = 0;
  352. Int_t usedItsHits = 0; // overall amount of ITS hits used when matching
  353. Int_t wrongMatch = 0;
  354. cout << "multimapMatch.size " << multimapMatch.size() << endl;
  355. multimap<Float_t, std::tuple<MpdItsKalmanTrack*, MpdItsKalmanTrack*, MpdTpcKalmanTrack*> >::iterator itr = multimapMatch.begin();
  356. MpdVectorFinder *vectorFinder = (MpdVectorFinder*) FairRun::Instance()->GetTask("MpdVectorFinder");
  357. // Select ITS+TPC or ITS-only tracks
  358. for ( ; itr != multimapMatch.end(); ++itr) {
  359. std::tuple<MpdItsKalmanTrack*, MpdItsKalmanTrack*, MpdTpcKalmanTrack*> &tupl = itr->second;
  360. ///cout << itr->first << " " << std::get<0>(tupl) << " " << std::get<1>(tupl) << " "<< std::get<2>(tupl) << endl;
  361. if (!tpcUnique.count(std::get<2>(tupl)) && !itsUnique.count(std::get<1>(tupl))) {
  362. tpcUnique.insert(std::get<2>(tupl));
  363. itsUnique.insert(std::get<1>(tupl));
  364. usedItsHits += std::get<0>(tupl)->GetNofHits(); //AZ
  365. if (std::get<0>(tupl)->GetNofHits() > 2) { ///good quality
  366. /// tuple: 0 - propagated track, 1 - its track matched, 2 - tpc track matched
  367. qualTr++;
  368. if (std::get<2>(tupl)->GetTrackID() != std::get<1>(tupl)->GetTrackID())
  369. wrongMatch++; //AZ 17.04.2020- track match for tracks with different id's
  370. //AZ }
  371. MpdItsKalmanTrack* track = std::get<0>(tupl); /// matched tpc+its track
  372. vectorFinder->GoToBeamLine(track);
  373. /*
  374. track->SetParam(*track->GetParamNew());
  375. track->SetPos(track->GetPosNew());
  376. Double_t pos = track->GetPos();
  377. TMatrixD par = *track->GetParam();
  378. TMatrixDSym cov = *track->Weight2Cov();
  379. Double_t leng = track->GetLength();
  380. TString nodeNew = track->GetNodeNew();
  381. //cout << " 1: " << nodeNew << ", " << track->GetNode() << endl;
  382. // Go to beam pipe
  383. MpdKalmanHit hit;
  384. hit.SetType(MpdKalmanHit::kFixedR);
  385. hit.SetPos(2.9); // fPipeR is calculated with geometry, which is not present in matching class
  386. Bool_t iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  387. if (iok != 1) {
  388. // Restore track
  389. track->SetParam(par);
  390. track->SetParamNew(par);
  391. track->SetCovariance(cov);
  392. track->ReSetWeight();
  393. track->SetPos(pos);
  394. track->SetPosNew(pos);
  395. track->SetLength(leng);
  396. //track->SetNode(node);
  397. //cout << " 2: " << nodeNew << ", " << track->GetNode() << endl;
  398. track->SetNodeNew(nodeNew);
  399. } else {
  400. // Add multiple scattering
  401. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  402. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  403. TMatrixDSym* pcov = track->Weight2Cov();
  404. Double_t th = track->GetParamNew(3);
  405. Double_t cosTh = TMath::Cos(th);
  406. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  407. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  408. (*pcov)(3,3) += angle2;
  409. Int_t ok = 0;
  410. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  411. track->SetWeight(*pcov);
  412. }
  413. cov = *track->Weight2Cov();
  414. hit.SetPos(0.);
  415. hit.SetMeas(0,track->GetParam(2)); // track Phi
  416. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  417. //Double_t pos = ((MpdKalmanHit*)track->GetHits()->Last())->GetPos();
  418. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  419. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  420. if (iok != 1) MpdKalmanFilter::Instance()->FindPca(track, vert);
  421. //track->SetPos(pos); // restore position
  422. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  423. */
  424. track->SetNofIts(track->GetNofHits());
  425. //AZ usedItsHits += track->GetNofHits();
  426. //AZ if (track->GetNofHits() > 2) {
  427. /// writing to TpcTracksRefit
  428. fTpcIndSet.erase(track->GetUniqueID()-1); // exclude matched TPC track
  429. track->SetUniqueID(std::get<1>(tupl)->GetTrackID()); /// ITS track ID
  430. //cout << "its and tpc track id equal " << track->GetTrackID() << " " << track->GetUniqueID() << endl;
  431. new ((*fTpcTracksRefit)[matchedCount]) MpdItsKalmanTrack(*track); /// was Tpc
  432. ///cout << "TPC + ITS track hits " << track->GetNofTrHits() << endl;
  433. //delete track; //AZ
  434. } else {
  435. /// refit its track anyway without matching to tpc track
  436. /// this track *std::get<1>(tupl) is refitted outward to 27.0, while i need here not refitted, original its track
  437. MpdItsKalmanTrack *temp = (MpdItsKalmanTrack*) fItsTracks->UncheckedAt(std::get<1>(tupl)->GetUniqueID() - 1); /// index is uniqueid - 1
  438. if (temp->GetNofTrHits() < 4) {
  439. // Exclude short ITS tracks
  440. tracksWithoutMatch.insert(std::get<1>(tupl));
  441. continue;
  442. }
  443. MpdItsKalmanTrack *track = new MpdItsKalmanTrack(*temp); //AZ
  444. vectorFinder->GoToBeamLine(track);
  445. /// writing to TpcTracksRefit
  446. MpdItsKalmanTrack *trtr = // AZ - 17.04.2020
  447. //AZ new ((*fTpcTracksRefit)[matchedCount]) MpdItsKalmanTrack(*temp); /// *std::get<1>(tupl)
  448. new ((*fTpcTracksRefit)[matchedCount]) MpdItsKalmanTrack(*track); //AZ
  449. trtr->SetChi2(temp->GetChi2Its()); //AZ
  450. //AZ-010121 trtr->SetUniqueID(std::get<1>(tupl)->GetTrackID());
  451. trtr->SetUniqueID(std::get<1>(tupl)->GetTrackID()+1); //AZ - to avoid 0
  452. delete track;
  453. }
  454. matchedCount++;
  455. ///cout << "ITS only track hits " << temp->GetNofTrHits() << " uniqueid " << std::get<1>(tupl)->GetUniqueID() - 1 << " " << temp->GetChi2Its() << endl;
  456. //AZ delete track;
  457. } /// if (!tpcUnique.count(std::get<2>(tupl)) && !itsUnique.count(std::get<1>(tupl))) {
  458. } // for ( ; itr != multimapMatch.end();
  459. //for (multimap<Float_t, std::tuple<MpdItsKalmanTrack*, MpdItsKalmanTrack*, MpdTpcKalmanTrack*>>::iterator itr = multimapMatch.begin(); itr != multimapMatch.end(); ++itr) {
  460. for (itr = multimapMatch.begin(); itr != multimapMatch.end(); ++itr) {
  461. if (!itsUnique.count(std::get<1>(itr->second))) {
  462. tracksWithoutMatch.insert(std::get<1>(itr->second));
  463. itsUnique.insert(std::get<1>(itr->second));
  464. }
  465. delete std::get<0>(itr->second);
  466. }
  467. Int_t i = matchedCount - 1;
  468. for (set<MpdItsKalmanTrack*>::iterator ittemp = tracksWithoutMatch.begin(); ittemp != tracksWithoutMatch.end(); ittemp++) {
  469. //AZ new ((*fTpcTracksRefit)[i]) MpdItsKalmanTrack(*(*ittemp));
  470. MpdItsKalmanTrack *temp = (MpdItsKalmanTrack*) fItsTracks->UncheckedAt((*ittemp)->GetUniqueID() - 1); //AZ
  471. Int_t nhits = temp->GetNofTrHits();
  472. /*
  473. if (nhits < 4) {
  474. // Exclude unmatched tracks with less than 4 ITS hits - reset "Use" flag
  475. TClonesArray *hits = temp->GetTrHits();
  476. for (Int_t ih = 0; ih < nhits; ++ih) {
  477. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(ih);
  478. MpdKalmanHit *hitK = (MpdKalmanHit*) fKHits->UncheckedAt(hit->GetUniqueID()-1);
  479. hitK->SetFlag(0);
  480. }
  481. continue;
  482. }
  483. */
  484. if (nhits < 6) {
  485. // Unmatched tracks - reset "Use" flag for reusing in KF-tracking
  486. TClonesArray *hits = temp->GetTrHits();
  487. for (Int_t ih = 0; ih < nhits; ++ih) {
  488. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(ih);
  489. MpdKalmanHit *hitK = (MpdKalmanHit*) fKHits->UncheckedAt(hit->GetUniqueID()-1);
  490. hitK->SetFlag(0);
  491. }
  492. if (nhits < 4) continue; // do not store short tracks (< 4 hits)
  493. }
  494. MpdItsKalmanTrack *track = new MpdItsKalmanTrack(*temp); //AZ
  495. vectorFinder->GoToBeamLine(track);
  496. i++;
  497. //AZ new ((*fTpcTracksRefit)[i]) MpdItsKalmanTrack(*temp);
  498. new ((*fTpcTracksRefit)[i]) MpdItsKalmanTrack(*track); //AZ
  499. delete track; //AZ
  500. }
  501. // Add TPC-only tracks to the same container
  502. Int_t nout = fTpcTracksRefit->GetEntriesFast();
  503. for (set<Int_t>::iterator sit = fTpcIndSet.begin(); sit != fTpcIndSet.end(); ++sit) {
  504. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(*sit);
  505. new ((*fTpcTracksRefit)[nout++]) MpdItsKalmanTrack(*tr);
  506. }
  507. cout << "fItsTracks size " << fItsTracks->GetEntriesFast() << endl;
  508. cout << "Its tracks without match: " << k << endl;
  509. cout << "Overall possible tpc+its tracks: " << tpck << endl;
  510. cout << "Overall matched tracks: " << matchedCount << endl;
  511. cout << "Overall good matched tracks: " << qualTr << endl;
  512. cout << "Overall used ITS hits when matching: " << usedItsHits << endl;
  513. cout << "wrongMatch: " << wrongMatch << endl;
  514. cout << "leftover unmatched TPC tracks: " << fTpcIndSet.size() << endl;
  515. fout_v_match << "fItsTracks size " << fItsTracks->GetEntriesFast() << endl;
  516. fout_v_match << "Its tracks without match: " << k << endl;
  517. fout_v_match << "Overall possible tpc+its tracks: " << tpck << endl;
  518. fout_v_match << "Overall matched tracks: " << matchedCount << endl;
  519. fout_v_match << "Overall good matched tracks: " << qualTr << endl;
  520. fout_v_match << "wrongMatch: " << wrongMatch << endl;
  521. fout_v_match << "- track matching done -" << endl;
  522. fout_v_match.close();
  523. }
  524. //__________________________________________________________________________
  525. ClassImp(MpdItsToTpcMatching)