MpdSftTrackFinderTpc.cxx 44 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144
  1. // -------------------------------------------------------------------------
  2. // ----- MpdSftTrackFinderTpc source file -----
  3. // ----- Created 03/12/08 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdSftTrackFinderTpc.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD "Silicon" Forward Tracker (SFT) and TPC
  9. **/
  10. #include "MpdSftTrackFinderTpc.h"
  11. #include "MpdKalmanFilter.h"
  12. #include "MpdKalmanHit.h"
  13. #include "MpdEctKalmanTrack.h"
  14. #include "MpdTpcKalmanFilter.h"
  15. #include "FairHit.h"
  16. #include "TpcLheHit.h"
  17. #include "FairMCPoint.h"
  18. #include "MpdMCTrack.h"
  19. #include "FairRootManager.h"
  20. #include "FairRun.h"
  21. #include "TMath.h"
  22. #include "TFile.h"
  23. //#include "TLorentzVector.h"
  24. #include "TVector2.h"
  25. //#include "TClonesArray.h"
  26. #include <TRandom.h>
  27. #include <iostream>
  28. //#include <vector>
  29. using std::cout;
  30. using std::endl;
  31. //using std::vector;
  32. // !!! 100 should be used only for track fitting !!!
  33. const Double_t MpdSftTrackFinderTpc::fgkChi2Cut = 100; //10; //100;
  34. FILE *lunSft1 = 0x0; //fopen("error1.dat","w");
  35. //__________________________________________________________________________
  36. MpdSftTrackFinderTpc::MpdSftTrackFinderTpc(const char *name, Int_t iVerbose )
  37. :FairTask(name, iVerbose)
  38. {
  39. fKHits = new TClonesArray("MpdKalmanHit", 100);
  40. fTracks = new TClonesArray("MpdEctKalmanTrack", 100);
  41. fHistoDir = 0x0;
  42. fhLays = new TH1F("hLaysSft","ECT layers",100,0,100);
  43. fLayMinMax[0][0] = fLayMinMax[1][0] = 999;
  44. fLayMinMax[0][1] = fLayMinMax[1][1] = 0;
  45. }
  46. //__________________________________________________________________________
  47. MpdSftTrackFinderTpc::~MpdSftTrackFinderTpc()
  48. {
  49. //delete fKHits;
  50. //delete fTracks;
  51. //delete fTrackCand;
  52. delete [] fLayPointers;
  53. delete fhLays;
  54. }
  55. //__________________________________________________________________________
  56. InitStatus MpdSftTrackFinderTpc::Init()
  57. {
  58. return ReInit();
  59. }
  60. //__________________________________________________________________________
  61. InitStatus MpdSftTrackFinderTpc::ReInit()
  62. {
  63. fSftPoints =(TClonesArray *) FairRootManager::Instance()->GetObject("SftPoint");
  64. fTpcPoints =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcPoint");
  65. fTpcHits =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanHit");
  66. //fTpcTracks = 0x0; //(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  67. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  68. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  69. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  70. FairRootManager::Instance()->Register("EctTrack", "Ect", fTracks, kTRUE);
  71. fNPass = 1;
  72. return kSUCCESS;
  73. }
  74. //__________________________________________________________________________
  75. void MpdSftTrackFinderTpc::Reset()
  76. {
  77. ///
  78. cout << " MpdSftTrackFinderTpc::Reset " << endl;
  79. //fKHits->Clear();
  80. //fTracks->Clear("C");
  81. fKHits->Delete();
  82. fTracks->Delete();
  83. delete [] fLayPointers;
  84. }
  85. //__________________________________________________________________________
  86. void MpdSftTrackFinderTpc::SetParContainers()
  87. {
  88. }
  89. //__________________________________________________________________________
  90. void MpdSftTrackFinderTpc::Finish()
  91. {
  92. //Write();
  93. }
  94. //__________________________________________________________________________
  95. void MpdSftTrackFinderTpc::Exec(Option_t * option)
  96. {
  97. static int eventCounter = 0;
  98. cout << " EctRec event " << ++eventCounter << endl;
  99. Reset();
  100. // Create Kalman hits
  101. MakeKalmanHitsSft();
  102. for (Int_t i = 0; i < fNPass; ++i) {
  103. fTracks->Clear();
  104. GetTrackSeeds(0); // Z > 0
  105. GetTrackSeeds(1); // Z < 0
  106. cout << " Total number of hits for tracking: " << fKHits->GetEntriesFast() << endl;
  107. cout << " Total number of track seeds: " << fTracks->GetEntriesFast() << endl;
  108. DoTracking(i);
  109. //StoreTracks();
  110. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  111. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  112. }
  113. RemoveDoubles(); // remove double tracks
  114. AddHits(); // add hit objects to tracks
  115. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  116. // Remove TPC tracks which are parts of the combined ones
  117. RemoveTpcTracks();
  118. // Propagate to the vertex
  119. //GoToVertex();
  120. }
  121. //__________________________________________________________________________
  122. void MpdSftTrackFinderTpc::MakeKalmanHitsSft()
  123. {
  124. /// Create Kalman hits from SFT points.
  125. Int_t lay, nKH = 0, layMax = 0;
  126. //Double_t phi, r, coord, errR = 0.008, errRphi = 0.008; // 80um in R, 80um in R-Phi
  127. Double_t phi, r, errR = 0.05, errRphi = 0.0170; // coord, 500um in R, 170um in R-Phi
  128. fhLays->Reset();
  129. Int_t nSft = fSftPoints->GetEntriesFast(), iz = 0;
  130. cout << " Number of SFT points: " << nSft << endl;
  131. TVector3 pos;
  132. for (Int_t i = 0; i < nSft; ++i) {
  133. FairMCPoint *p = (FairMCPoint*) fSftPoints->UncheckedAt(i);
  134. p->Position(pos);
  135. Int_t mod = p->GetDetectorID() >> 5;
  136. lay = (p->GetDetectorID() & 15) - 1;
  137. lay += (mod-1) * 50;
  138. //if (lay > 2 && lay < 6) continue;
  139. layMax = TMath::Max (lay, layMax);
  140. iz = (p->GetZ() > 0) ? 0 : 1;
  141. fLayMinMax[iz][0] = TMath::Min (lay, fLayMinMax[iz][0]);
  142. fLayMinMax[iz][1] = TMath::Max (lay, fLayMinMax[iz][1]);
  143. fhLays->Fill(lay+0.1);
  144. // Add error
  145. Double_t dRphi = 0, dR = 0;
  146. gRandom->Rannor(dRphi,dR); // add errors
  147. r = pos.Pt();
  148. phi = pos.Phi();
  149. //MpdKalmanHitZ *hit = new ((*fKHits)[nKH++]) MpdKalmanHitZ(r+dR*errR, phi, 0+dRphi*errRphi,
  150. // p->GetZ(), errR, errRphi, lay, i);
  151. Double_t meas[2] = {0+dRphi*errRphi, r+dR*errR};
  152. Double_t err[2] = {errRphi, errR};
  153. Double_t rphi = meas[1] * phi + meas[0];
  154. phi = rphi / meas[1];
  155. meas[0] = 0.0; //
  156. Double_t cosSin[2] = {TMath::Cos(phi), TMath::Sin(phi)};
  157. //MpdKalmanHit *hit =
  158. new ( (*fKHits)[nKH++]) MpdKalmanHit(lay*1000000, 2,
  159. MpdKalmanHit::kFixedZ, meas, err, cosSin, 9.0, p->GetZ(), i);
  160. }
  161. fKHits->Sort(); // in descending order in layer No. (in abs(Z))
  162. cout << " First and last layers: " << ((MpdKalmanHit*)fKHits->First())->GetLayer() << " "
  163. << ((MpdKalmanHit*)fKHits->Last())->GetLayer() << endl;
  164. fLayPointers = new Int_t [layMax+1];
  165. Int_t ipos = 0;
  166. //for (Int_t i = 0; i <= layMax; ++i) {
  167. for (Int_t i = layMax; i >= 0; --i) {
  168. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  169. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  170. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  171. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  172. if (cont == 0) {
  173. fLayPointers[i] = -1;
  174. continue;
  175. }
  176. fLayPointers[i] = ipos;
  177. ipos += cont;
  178. }
  179. }
  180. //__________________________________________________________________________
  181. void MpdSftTrackFinderTpc::GetTrackSeeds(Int_t iPass)
  182. {
  183. /// Build track seeds from SFT hits and vertex
  184. Int_t nCand = fTracks->GetEntriesFast();//, layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  185. TVector3 vert(0.,0.,0.), pmom;
  186. // Build track candidates from 2 SFT detectors (closest to the vertex)
  187. Int_t nLay[2], indx0[2];
  188. for (Int_t i = 0; i < 2; ++i) {
  189. Int_t lay = (iPass == 0) ? 1 - i : fLayMinMax[iPass][0] + 1 - i;
  190. nLay[i] = GetNofHitsInLayer(lay);
  191. indx0[i] = GetHitsInLayer(lay);
  192. //cout << lay << " " << nLay[i] << endl;
  193. }
  194. /*
  195. // Build track candidates from 2 SFT detectors (farthest from the vertex)
  196. Int_t nLay[2], indx0[2];
  197. for (Int_t lay = 0; lay < 2; ++lay) {
  198. nLay[lay] = GetNofHitsInLayer(layMax-lay);
  199. indx0[lay] = GetHitsInLayer(layMax-lay);
  200. }
  201. */
  202. FairHit hitTmp;
  203. TpcLheHit hit1Tmp;
  204. for (Int_t i = 0; i < nLay[0]; ++i) {
  205. // Loop over hits in the downstream SFT
  206. MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(indx0[0]+i);
  207. Int_t id1 = ((FairMCPoint*)fSftPoints->UncheckedAt(hit1->GetIndex()))->GetTrackID();
  208. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(id1);
  209. if (mctrack->GetMotherId() >= 0) continue; // skip secondary tracks - for test
  210. hitTmp.SetZ(hit1->GetDist());
  211. for (Int_t j = 0; j < nLay[1]; ++j) {
  212. // Loop over hits in the upstream SFT
  213. MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(indx0[1]+j);
  214. Int_t id2 = ((FairMCPoint*)fSftPoints->UncheckedAt(hit2->GetIndex()))->GetTrackID();
  215. if (id2 != id1) continue; // !!! exact match of 2 SFTs
  216. MpdEctKalmanTrack *track =
  217. new ((*fTracks)[nCand++]) MpdEctKalmanTrack(-1, -1, &hitTmp, &hit1Tmp, vert);
  218. track->SetDirection(MpdKalmanTrack::kOutward);
  219. Double_t pt = EvalPt(track, hit1, hit2);
  220. ((FairMCPoint*) fSftPoints->UncheckedAt(hit2->GetIndex()))->Momentum(pmom);
  221. track->SetParam (4, 1./pt); // q/Pt
  222. //track->SetParam (4, 1./pmom.Pt()*TMath::Sign(1.,pt)); // q/Pt - exact value
  223. track->SetTrackID(id1);
  224. EvalCovar(hit1, hit2, track);
  225. cout << nCand << " " << id1 << " " << pmom.Pt() << " " << 1./track->GetParam(4) << " " << track->GetParam(0) << " " << track->GetParam(1) << " " << track->GetParam(2) << " " << track->GetParam(3) << " " << hit2->GetPhi() << endl;
  226. cout << pmom.Phi() << " " << TMath::PiOver2()-pmom.Theta() << endl;
  227. //track->SetParam(2,pmom.Phi()); // exact value
  228. //track->SetParam(3,TMath::PiOver2()-pmom.Theta()); // exact value
  229. track->SetParamNew (*track->GetParam());
  230. track->SetPos(hit2->GetDist());
  231. track->SetPosNew(hit2->GetDist());
  232. if (lunSft1) fprintf(lunSft1,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e \n",track->GetParam(2),MpdKalmanFilter::Instance()->Proxim(track->GetParam(2),pmom.Phi()),track->GetParam(3),MpdKalmanFilter::Instance()->Proxim(track->GetParam(3),TMath::PiOver2()-pmom.Theta()),TMath::Sqrt((*track->Weight2Cov())(2,2)),TMath::Sqrt((*track->Weight2Cov())(3,3)),1./track->GetParam(4),pmom.Pt());
  233. //track->GetHits()->RemoveAt(0);
  234. track->GetHits()->Add(hit2);
  235. track->GetHits()->Compress();
  236. }
  237. }
  238. cout << " Number of SFT track candidates: " << nCand << " " << endl;
  239. //if (lunCpc) fclose(lunCpc);
  240. }
  241. //__________________________________________________________________________
  242. Double_t MpdSftTrackFinderTpc::EvalPt(MpdKalmanTrack *track, const MpdKalmanHit *hit1, const MpdKalmanHit *hit2)
  243. {
  244. /// Evaluate signed track Pt (curvature) assuming the track coming from the
  245. /// primary vertex and fill track parameters
  246. Double_t ph1 = hit1->GetMeas(0) / hit1->GetMeas(1) + hit1->GetPhi();
  247. //Double_t ph2 = hit2->GetRphi() / hit2->GetR() + hit2->GetPhi();
  248. Double_t ph2 = hit2->GetMeas(0) / hit2->GetMeas(1) + hit2->GetPhi();
  249. //Double_t x1 = hit1->GetR() * TMath::Cos(ph1);
  250. Double_t x1 = hit1->GetMeas(1) * TMath::Cos(ph1);
  251. Double_t y1 = hit1->GetMeas(1) * TMath::Sin(ph1);
  252. Double_t x2 = hit2->GetMeas(1) * TMath::Cos(ph2);
  253. Double_t y2 = hit2->GetMeas(1) * TMath::Sin(ph2);
  254. Double_t z1 = hit1->GetDist();
  255. Double_t z2 = hit2->GetDist();
  256. TVector2 vec1(x1-0.,y1-0.);
  257. TVector2 vec2(x2-0.,y2-0.);
  258. TVector2 vec21 = vec1 - vec2;
  259. Double_t cosAlpha = vec2 * vec21 / vec2.Mod() / vec21.Mod();
  260. Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  261. Double_t factor = 0.0015; // 0.3 * 0.5T * 0.01
  262. //factor *= 3; // 1.5T
  263. Double_t charge = ph1 - MpdKalmanFilter::Instance()->Proxim(ph1,ph2);
  264. if (hit1->GetLayer() < hit2->GetLayer()) charge = -charge;
  265. Double_t pt = factor * TMath::Abs(rad);
  266. track->SetParam(0,x2);
  267. track->SetParam(1,y2);
  268. track->SetParam(2,vec21.Phi());
  269. track->SetParam(3,TMath::ATan2(z1-z2,vec21.Mod()));
  270. //return TMath::Min(pt, 25.0) * TMath::Sign(1., -charge);
  271. return TMath::Min(pt, 25.0) * TMath::Sign(1., charge);
  272. }
  273. //__________________________________________________________________________
  274. void MpdSftTrackFinderTpc::EvalCovar(const MpdKalmanHit *hitOut, const MpdKalmanHit *hitIn,
  275. MpdEctKalmanTrack *track)
  276. {
  277. /// Evaluate covariance matrix for track seed
  278. TMatrixD ww(5,5);
  279. //ww(0,0) = hitOut->GetRphiErr() * hitOut->GetRphiErr(); // <RphiRphi>
  280. ww(0,0) = hitOut->GetErr(0) * hitOut->GetErr(0); // <RphiRphi>
  281. //ww(0,0) *= 4.; // extra factor of 2
  282. ww(1,1) = hitOut->GetErr(1) * hitOut->GetErr(1); // <RR>
  283. //Double_t phiOut = hitOut->GetPhi();
  284. //Double_t phiIn = hitIn->GetPhi();
  285. Double_t phiOut = hitOut->GetMeas(0) / hitOut->GetMeas(1) + hitOut->GetPhi();
  286. Double_t phiIn = hitIn->GetMeas(0) / hitIn->GetMeas(1) + hitIn->GetPhi();
  287. Double_t xIn = hitIn->GetMeas(1) * TMath::Cos(phiIn);
  288. Double_t yIn = hitIn->GetMeas(1) * TMath::Sin(phiIn);
  289. Double_t xOut = hitOut->GetMeas(1) * TMath::Cos(phiOut);
  290. Double_t yOut = hitOut->GetMeas(1) * TMath::Sin(phiOut);
  291. Double_t dx = xOut - xIn, dy = yOut - yIn;
  292. Double_t dist2 = dx * dx + dy * dy;
  293. Double_t sinPhi = TMath::Sin (track->GetParam(2));
  294. Double_t cosPhi = TMath::Cos (track->GetParam(2));
  295. Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
  296. Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
  297. ww(2,2) = (pOut * pOut + pIn * pIn) / dist2 * ww(0,0); // <PhiPhi>
  298. //ww(2,2) /= 50.; // extra factor of 2
  299. Double_t tanThe = TMath::Tan(track->GetParam(3));
  300. Double_t dRad = hitOut->GetMeas(1) - hitIn->GetMeas(1);
  301. Double_t denom = dRad * (1.+tanThe*tanThe);
  302. ww(3,3) = ww(1,1) * 2. / denom / denom; // <TheThe>
  303. ww(3,3) *= 4;
  304. //ww(4,4) = (track->GetParam(4)*0.5) * (track->GetParam(4)*0.5); // error 50%
  305. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
  306. ww(4,4) = (track->GetParam(4)*1.) * (track->GetParam(4)*1.); // error 100%
  307. // Error transformation from R-Rphi to X-Y
  308. /*
  309. Double_t cosp = TMath::Cos(phiIn);
  310. Double_t sinp = TMath::Sin(phiIn);
  311. Double_t xx = cosp * cosp * ww(1,1) + sinp * sinp * ww(0,0);
  312. Double_t yy = sinp * sinp * ww(1,1) + cosp * cosp * ww(0,0);
  313. Double_t xy = 2 * sinp * cosp * ww(1,1) * ww(0,0);
  314. ww(0,0) = xx * 1;
  315. ww(1,1) = yy;
  316. ww(0,1) = ww(1,0) = xy;
  317. */
  318. //fWeight->Print();
  319. //fWeight->Invert(); // weight matrix
  320. Int_t iok = 0;
  321. MpdKalmanFilter::Instance()->MnvertLocal(ww.GetMatrixArray(), 5, 5, 5, iok);
  322. track->SetWeight(ww);
  323. //fWeight->Print();
  324. }
  325. //__________________________________________________________________________
  326. void MpdSftTrackFinderTpc::DoTracking(Int_t iPass)
  327. {
  328. /// Run Kalman tracking
  329. Int_t nCand = fTracks->GetEntriesFast(), iok = 0;
  330. for (Int_t i = 0; i < nCand; ++i) {
  331. //for (Int_t i = 20; i < 21; ++i) {
  332. cout << " Track seed No. " << i << endl;
  333. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  334. //if (track->GetTrackID() != 117) continue;
  335. //if (track->GetTrackID() != 10) continue;
  336. //if (track->GetTrackID() != 1105) continue;
  337. //if (track->GetTrackID() != 77) continue;
  338. //(*(track->GetParamNew()))(4,0) = -0.5; // just for check
  339. /*
  340. for (Int_t k = 0; k < 5; ++k) {
  341. for (Int_t j = i; j < 5; ++j) {
  342. if (j == i) (*track->GetWeight())(i,j) /= 100.;
  343. else (*track->GetWeight())(i,j) = (*track->GetWeight())(j,i) = 0.;
  344. }
  345. }
  346. */
  347. // Start ECT tracking from different layers to account for missing hits
  348. const Int_t layMax = 1; //10;
  349. MpdEctKalmanTrack tracks[layMax];
  350. Int_t lay = ((MpdKalmanHit*)track->GetHits()->UncheckedAt(0))->GetLayer() + 1;
  351. for (Int_t j = 0; j < layMax; ++j) {
  352. tracks[j] = *track;
  353. //cout << track->GetParamNew(0) << endl;
  354. //cout << i << " " << lay << " " << tracks[lay].GetNofHits() << " " << tracks[lay].GetChi2() << " " << tracks[lay].GetParam(0) << endl;
  355. if (j > 0) {
  356. if (tracks[j-1].GetNofHits() > 2) // 2 hits from CPC there
  357. lay = ((MpdKalmanHit*)tracks[j-1].GetHits()->UncheckedAt(2))->GetLayer() - 1;
  358. else break; // no more hits found for track
  359. }
  360. iok = RunKalmanFilter(&tracks[j], lay);
  361. //iok = RunKalmanFilter(&tracks[lay], 0);
  362. cout << j << " " << lay << " " << tracks[j].GetNofHits() << " " << tracks[j].GetChi2() << " " << iok << endl;
  363. }
  364. // Select the best track (with max number of hits)
  365. Int_t nHitsMax = tracks[0].GetNofHits(), iMax = 0;
  366. for (Int_t j = 1; j < layMax; ++j) {
  367. Int_t nhits = tracks[j].GetNofHits();
  368. if (nhits == 0) break;
  369. if (nhits > nHitsMax) {
  370. nHitsMax = nhits;
  371. iMax = j;
  372. } else if (nhits == nHitsMax) {
  373. if (tracks[j].GetChi2() < tracks[iMax].GetChi2()) {
  374. iMax = j;
  375. nHitsMax = nhits;
  376. }
  377. }
  378. }
  379. fTracks->RemoveAt(i);
  380. new ((*fTracks)[i]) MpdEctKalmanTrack(tracks[iMax]);
  381. // Refit
  382. if (1) {
  383. //track->Weight2Cov()->Print();
  384. //track->GetParamNew()->Print();
  385. MpdKalmanHit hitTmp;
  386. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  387. hitTmp.SetDist(track->GetPosNew()+TMath::Sign(0.2,track->GetPosNew()));
  388. MpdKalmanFilter::Instance()->PropagateToHit(track, &hitTmp);
  389. track->SetDirection(MpdKalmanTrack::kInward);
  390. //cout << " Chi2 = " << track->GetChi2() << endl;
  391. MpdKalmanFilter::Instance()->Refit(track, 1);
  392. //track->GetParamNew()->Print();
  393. //track->SetDirection(MpdKalmanTrack::kOutward);
  394. //MpdKalmanFilter::Instance()->Refit(track, 1);
  395. }
  396. cout << i << " " << track->GetNofHits() << " " << 1 / track->GetParamNew(4) << endl;
  397. } //for (Int_t i = 0; i < nCand;
  398. FollowInTPC();
  399. }
  400. //__________________________________________________________________________
  401. void MpdSftTrackFinderTpc::FollowInTPC()
  402. {
  403. /// Follow tracks in TPC
  404. // Manage TPC hits
  405. if (fTpcHits->GetEntriesFast() == 0) { cout << " !!! No TPC hits" << endl; return; }
  406. fTpcHits->Sort();
  407. fhLays->Reset();
  408. Int_t nHits = fTpcHits->GetEntriesFast(), layMax = 0;
  409. for (Int_t i = 0; i < nHits; ++i) {
  410. MpdKalmanHit *hit = (MpdKalmanHit*) fTpcHits->UncheckedAt(i);
  411. Int_t lay = hit->GetLayer();
  412. fhLays->Fill(lay+0.1);
  413. layMax = TMath::Max (layMax, lay);
  414. }
  415. delete [] fLayPointers;
  416. fLayPointers = new Int_t [layMax+1];
  417. Int_t ipos = 0;
  418. for (Int_t i = layMax; i > -1; --i) {
  419. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  420. if (cont == 0) {
  421. fLayPointers[i] = -1;
  422. continue;
  423. }
  424. fLayPointers[i] = ipos;
  425. ipos += cont;
  426. }
  427. Int_t nCand = fTracks->GetEntriesFast(), iok = 0;
  428. for (Int_t i = 0; i < nCand; ++i) {
  429. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  430. // Propagate track to TPC first hit (-0.5 cm)
  431. /*
  432. MpdKalmanHit *hit = (MpdKalmanHit*) fTpcHits->Last();
  433. MpdKalmanHitZ hitZ;
  434. hitZ.SetZ (hit->GetZ() - 0.5);
  435. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitZ,kTRUE);
  436. */
  437. /*
  438. // Propagate track to TPC inner shell
  439. MpdKalmanHit hitR;
  440. hitR.SetType(MpdKalmanHit::kFixedR);
  441. hitR.SetDist (20.);
  442. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitR,kTRUE);
  443. // Propagate thru TPC inner shell
  444. PassInnerShell(track);
  445. iok = RunKalmanFilterTPC(track, 0);
  446. */
  447. Int_t nh1 = track->GetNofHits();
  448. iok = RunKalmanFilterTPC(track, layMax);
  449. if (iok < 0 || track->GetNofHits() == nh1) {
  450. fTracks->RemoveAt(i);
  451. continue; // failed propagation or no TPC hits added
  452. }
  453. track->Weight2Cov();
  454. MpdKalmanHit hitR;
  455. hitR.SetType(MpdKalmanHit::kFixedR);
  456. hitR.SetDist(10.);
  457. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitR,kTRUE);
  458. track->SetParam(*track->GetParamNew());
  459. }
  460. fTracks->Compress();
  461. }
  462. //__________________________________________________________________________
  463. Int_t MpdSftTrackFinderTpc::RunKalmanFilter(MpdEctKalmanTrack *track, Int_t layBeg)
  464. {
  465. /// Run Kalman filter
  466. //const Double_t rMin = -28., rMax = 123.; // min and max radial size of ECT - to be taken elsewhere
  467. const Double_t rMin = -28., rMax = 170.; // min and max radial size of ECT - to be taken elsewhere
  468. //cout << fHits->GetEntriesFast() << endl;
  469. //Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  470. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  471. MpdKalmanHit *hitOK = 0x0;
  472. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  473. //Int_t layBeg = layMax, layEnd = -1, dLay = -1, layOK = -1;
  474. Int_t layEnd = -1, dLay = -1;//, layOK = -1;
  475. if (trackDir == MpdKalmanTrack::kOutward) {
  476. layEnd = layMax + 1;
  477. dLay = 1;
  478. }
  479. layEnd = (track->GetPos() > 0) ? fLayMinMax[0][1] + 1 : fLayMinMax[1][1] + 1;
  480. //Int_t indxOK = hits->IndexOf(hitOK);
  481. //Int_t nHits = hits->GetEntriesFast();
  482. Int_t miss = 0, istop = 0;
  483. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  484. TMatrixD param(5,1), paramTmp(5,1);
  485. Double_t saveZ = 0;
  486. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  487. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  488. Int_t nLay = GetNofHitsInLayer(lay);
  489. Int_t indx0 = GetHitsInLayer(lay);
  490. Double_t dChi2Min = 1.e+6;
  491. MpdKalmanHit *hitMin = 0x0;
  492. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  493. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  494. /*
  495. if (trackDir == MpdKalmanTrack::kOutward) {
  496. // !!! This part is due to the simplified hit merging (digitization) -
  497. // different hit position inside one layer - should be removed later
  498. indxBeg = nLay - 1;
  499. indxEnd = -1;
  500. dIndx = -1;
  501. }
  502. */
  503. //for (Int_t indx = 0; indx < nLay; ++indx) {
  504. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  505. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  506. if (track->GetPosNew() * hit->GetDist() < 0) { istop = 1; break; } // hit with opposite Z - finish
  507. // !!! Exact ID match
  508. if (((FairMCPoint*)fSftPoints->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  509. // Exclude used hits
  510. if (hit->GetFlag() != 1) continue;
  511. //cout << " window:" << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << hit->GetZ() << " " << track->GetParamNew(1) << endl;
  512. // Check if the hit within some window (15x15cm for the moment - check!!!)
  513. //if (TMath::Abs(hit->GetRphi()-track->GetParamNew(0)) > 9) continue;
  514. //if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  515. TVector2 dist = GetDistance(track, hit);
  516. if (dist.X() > 15.) continue; // distance in transverse to the tube direction
  517. if (hit->GetNofDim() > 1 && dist.Y() > 25.) continue; // distance in R
  518. // For tests
  519. //FairMCPoint *mc = (FairMCPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  520. //*if (hit->GetTrackID() == 186)*/ cout << " Hit: " << mc->GetTrackID() << " " << hit->GetLayer() << " " << hit->GetR()*TMath::Cos(hit->GetRphi()/hit->GetR()) << " " << hit->GetR()*TMath::Sin(hit->GetRphi()/hit->GetR()) << " " << hit->GetR() << " " << hit->GetZ() << "; Track: " << track->GetParamNew(1)*TMath::Cos(track->GetParamNew(0)/track->GetParamNew(1)) << " " << track->GetParamNew(1)*TMath::Sin(track->GetParamNew(0)/track->GetParamNew(1)) << " " << track->GetParamNew(1) << " " << track->GetPosNew() << "; Dist: " << dist.X() << " " << dist.Y() << endl;
  521. //track->Print("");
  522. //hit->Print("");
  523. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) return -1;
  524. // For tests
  525. //FairMCPoint *mc = (FairMCPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  526. //FairMCPoint *mc = (FairMCPoint*) fSftPoints->UncheckedAt(hit->GetIndex());
  527. //*if (hit->GetTrackID() == 186)*/ cout << " Hit: " << mc->GetTrackID() << " " << hit->GetLayer() << " " << mc->GetX() << " " << mc->GetY() << " " << hit->GetDist() << " " << hit->GetPhi() << "; Track: " << track->GetParamNew(0) << " " << track->GetParamNew(1) << " " << track->GetPosNew() << "; Dist: " << dist.X() << " " << dist.Y() << endl;
  528. //
  529. // For debugging
  530. /*
  531. TVector2 dist0 = GetDistance(track, hit);
  532. cout << dist0.X() << " " << dist0.Y() << endl;
  533. *.
  534. MpdKalmanHitZ hitDbg = *hit;
  535. Double_t xDbg = hit->GetXY(0) * TMath::Cos(hit->GetPhi()) + hit->GetXY(1) * TMath::Sin(hit->GetPhi());
  536. Double_t yDbg = -hit->GetXY(0) * TMath::Sin(hit->GetPhi()) + hit->GetXY(1) * TMath::Cos(hit->GetPhi());
  537. hitDbg.SetRphi(yDbg);
  538. hitDbg.SetR(xDbg);
  539. dist = GetDistance(track, &hitDbg);
  540. cout << dist.X() << endl;
  541. */
  542. Double_t radNew = track->GetRadNew();
  543. if (radNew < rMin || radNew > rMax) return -1;
  544. //Double_t err = hit->GetRphiErr();
  545. //if (track->GetNofHits() == 0) hit->SetRphiErr(0.04);
  546. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(track,hit,pointWeight,param);
  547. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  548. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  549. if (dChi2 < dChi2Min) {
  550. //if (dChi2 < dChi2Min && dist0.X() <= dist.X()) {
  551. //if (dChi2 < dChi2Min && dist.X() <= 0.2) {
  552. dChi2Min = dChi2;
  553. hitMin = hit;
  554. saveWeight = *track->GetWeight();
  555. saveZ = track->GetPosNew();
  556. // temporary storage for the current track
  557. paramTmp = param;
  558. pointWeightTmp = pointWeight;
  559. /*cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetMeas(0) << " " << hitMin->GetMeas(1) << " " << ((FairMCPoint*)fSftPoints->UncheckedAt(hit->GetIndex()))->GetTrackID() << endl;
  560. cout << 1/param(4,0) << " " << endl;*/
  561. }
  562. } // for (Int_t indx = indxBeg; indx != indxEnd;
  563. if (istop) break;
  564. Double_t cut = fgkChi2Cut;
  565. //if (track->GetNofHits() == 0) cut /= 2.;
  566. //if (dChi2Min < fgkChi2Cut) {
  567. if (dChi2Min < cut) {
  568. //layOK = lay;
  569. hitOK = hitMin;
  570. track->GetHits()->Add(hitOK);
  571. miss = 0;
  572. // Restore track params at the best hit
  573. track->SetChi2(track->GetChi2()+dChi2Min);
  574. saveWeight += pointWeightTmp;
  575. track->SetWeight(saveWeight);
  576. track->SetPosNew(saveZ);
  577. track->SetParamNew(paramTmp);
  578. //cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  579. //paramTmp.Print();
  580. // Check if the accepted hit is the same as the seeded hit
  581. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  582. if (track->GetNofHits() == 1) track->SetParam1();
  583. } else {
  584. ++miss;
  585. //if (miss > 1) return -1;
  586. }
  587. //track->GetParamNew()->Print();
  588. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  589. } // for (Int_t lay = layBeg; lay != layEnd;
  590. // Update track params to those at the last accepted hit
  591. if (miss == 0) {
  592. track->SetWeight(saveWeight);
  593. track->SetPos(saveZ);
  594. track->SetParam(paramTmp);
  595. } else {
  596. track->SetPos(track->GetPosNew());
  597. track->SetParam(*track->GetParamNew());
  598. }
  599. return 0;
  600. }
  601. //__________________________________________________________________________
  602. Int_t MpdSftTrackFinderTpc::RunKalmanFilterTPC(MpdEctKalmanTrack *track, Int_t layBeg)
  603. {
  604. /// Run Kalman filter in TPC
  605. // Get TpcKalmanFilter task
  606. MpdTpcKalmanFilter *tpcKF = (MpdTpcKalmanFilter*) FairRun::Instance()->GetTask("TPC Kalman filter");
  607. if (tpcKF == nullptr) Fatal("RunKalmanFilterTPC", "!!! No TPC kalman filter task is found !!! ");
  608. MpdKalmanHit *hitOK = 0x0;
  609. Int_t layEnd = -1, dLay = -1;//, layOK = -1;
  610. Int_t layMax = ((MpdKalmanHit*)fTpcHits->UncheckedAt(0))->GetLayer();
  611. if (track->GetDirection() == MpdKalmanTrack::kOutward) {
  612. layEnd = layMax + 1;
  613. dLay = 1;
  614. }
  615. Int_t miss = 0;
  616. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  617. TMatrixD param(5,1), paramTmp(5,1);
  618. Double_t saveR = 0;
  619. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  620. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  621. Int_t nLay = GetNofHitsInLayer(lay);
  622. if (nLay == 0) continue;
  623. Int_t indx0 = GetHitsInLayer(lay);
  624. Double_t dChi2Min = 1.e+6;
  625. MpdKalmanHit *hitMin = 0x0;
  626. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  627. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  628. /*
  629. if (trackDir == MpdKalmanTrack::kOutward) {
  630. // !!! This part is due to the simplified hit merging (digitization) -
  631. // different hit position inside one layer - should be removed later
  632. indxBeg = nLay - 1;
  633. indxEnd = -1;
  634. dIndx = -1;
  635. }
  636. */
  637. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  638. MpdKalmanHit *hit = (MpdKalmanHit*) fTpcHits->UncheckedAt(indx0+indx);
  639. // !!! Exact ID match
  640. //if (((FairMCPoint*)fTpcPoints->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  641. if (tpcKF->GetHitID(hit) != track->GetTrackID()) continue;
  642. // Exclude used hits
  643. if (hit->GetFlag() != 1) continue;
  644. if (track->GetType() == MpdKalmanTrack::kBarrel) {
  645. // Check if the hit within some window
  646. //if (TMath::Abs(hit->GetMeas(0)-track->GetParamNew(0)) > 9) continue;
  647. //if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  648. if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) / hit->GetPos() > TMath::PiOver4()) continue;
  649. }
  650. //*if (hit->GetTrackID() == 186)*/ cout << " Hit: " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() << " " << hit->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetDetectorID() << " " << hit->GetRphi() << " " << hit->GetR() << " " << hit->GetZ() << " " << dist.X() << " " << dist.Y() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  651. //track->Print("");
  652. //hit->Print("");
  653. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) return -1;
  654. /*Debug
  655. cout << " Hit: " << hit->GetLayer() << " " << hit->GetMeas(0) << " " << hit->GetMeas(1) << " "
  656. << hit->GetDist() << ", Track: " << track->GetParamNew(0) << " " << track->GetParamNew(1)
  657. << " " << track->GetPosNew() << " " << track->GetPos() << endl;
  658. */
  659. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitR(track,hit,pointWeight,param);
  660. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  661. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  662. if (dChi2 < dChi2Min) {
  663. dChi2Min = dChi2;
  664. hitMin = hit;
  665. saveWeight = *track->GetWeight();
  666. saveR = track->GetPosNew();
  667. // temporary storage for the current track
  668. paramTmp = param;
  669. pointWeightTmp = pointWeight;
  670. /*cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetMeas(0) << " " << hitMin->GetMeas(1)
  671. << " " << tpcKF->GetHitID(hit) << endl;
  672. cout << 1/param(4,0) << " " << endl;*/
  673. }
  674. }
  675. if (dChi2Min < fgkChi2Cut) {
  676. //layOK = lay;
  677. hitOK = hitMin;
  678. track->GetHits()->Add(hitOK);
  679. track->SetTpcIndex(track->GetTrackID()); //
  680. miss = 0;
  681. // Restore track params at the best hit
  682. track->SetChi2(track->GetChi2()+dChi2Min);
  683. saveWeight += pointWeightTmp;
  684. track->SetWeight(saveWeight);
  685. track->SetPosNew(saveR);
  686. track->SetParamNew(paramTmp);
  687. //cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  688. //paramTmp.Print();
  689. // Check if the accepted hit is the same as the seeded hit
  690. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  691. if (track->GetNofHits() == 1) track->SetParam1();
  692. } else {
  693. ++miss;
  694. //if (miss > 1) return -1;
  695. }
  696. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  697. } // for (Int_t lay = layBeg; lay != layEnd;
  698. // Update track params to those at the last accepted hit
  699. if (miss == 0) {
  700. track->SetWeight(saveWeight);
  701. track->SetPos(saveR);
  702. track->SetParam(paramTmp);
  703. } else {
  704. track->SetPos(track->GetPosNew());
  705. track->SetParam(*track->GetParamNew());
  706. }
  707. return 0;
  708. }
  709. //__________________________________________________________________________
  710. TVector2 MpdSftTrackFinderTpc::GetDistance(MpdEctKalmanTrack *track, MpdKalmanHit *hit)
  711. {
  712. /// Compute distance between track and hit
  713. Double_t xTr, yTr;
  714. if (track->GetType() == MpdKalmanTrack::kEndcap) {
  715. xTr = track->GetParamNew(0);
  716. yTr = track->GetParamNew(1);
  717. } else {
  718. Double_t rPhi = track->GetParamNew(0);
  719. Double_t r = track->GetPosNew();
  720. Double_t ph = rPhi / r;
  721. xTr = r * TMath::Cos(ph);
  722. yTr = r * TMath::Sin(ph);
  723. }
  724. Double_t phi = hit->GetPhi();
  725. Double_t cosPhi = TMath::Cos(phi);
  726. Double_t sinPhi = TMath::Sin(phi);
  727. // Transform track coordinates to local tube coordinates
  728. Double_t xLoc = xTr * cosPhi + yTr * sinPhi;
  729. Double_t yLoc = -xTr * sinPhi + yTr * cosPhi;
  730. //Double_t xLoc = (xTr - hit->GetXY(0)) * cosPhi + (yTr - hit->GetXY(1)) * sinPhi;
  731. //Double_t yLoc = -(xTr - hit->GetXY(0)) * sinPhi + (yTr - hit->GetXY(1)) * cosPhi;
  732. //TVector2 dist(yLoc-hit->GetRphi(), xLoc-hit->GetR());
  733. TVector2 dist(TMath::Abs(yLoc-hit->GetMeas(0)), TMath::Abs(xLoc-hit->GetMeas(1)));
  734. return dist;
  735. }
  736. //__________________________________________________________________________
  737. Double_t MpdSftTrackFinderTpc::Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit)
  738. {
  739. /// Adjust hit coord. R-Phi to be "around" track R-Phi - to avoid
  740. /// discontinuity around +- Pi
  741. Double_t hitPhi = hit->GetMeas(0) / hit->GetDist();
  742. Double_t phi0 = track->GetParamNew(0) / track->GetPosNew();
  743. return hit->GetDist() * MpdKalmanFilter::Instance()->Proxim(phi0, hitPhi);
  744. }
  745. //__________________________________________________________________________
  746. void MpdSftTrackFinderTpc::Write()
  747. {
  748. /// Write
  749. TFile histoFile("EctRec.root","RECREATE");
  750. Writedir2current(fHistoDir);
  751. histoFile.Close();
  752. }
  753. //__________________________________________________________________________
  754. void MpdSftTrackFinderTpc::Writedir2current( TObject *obj )
  755. {
  756. /// Write
  757. if( !obj->IsFolder() ) obj->Write();
  758. else{
  759. TDirectory *cur = gDirectory;
  760. TDirectory *sub = cur->mkdir(obj->GetName());
  761. sub->cd();
  762. TList *listSub = ((TDirectory*)obj)->GetList();
  763. TIter it(listSub);
  764. while( TObject *obj1=it() ) Writedir2current(obj1);
  765. cur->cd();
  766. }
  767. }
  768. //__________________________________________________________________________
  769. void MpdSftTrackFinderTpc::RemoveDoubles()
  770. {
  771. /// Remove double tracks (if number of common hits greater than 50% of hits on track)
  772. Int_t nFound = fTracks->GetEntriesFast(), nOK = 0;
  773. for (Int_t i = 0; i < nFound; ++i) {
  774. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  775. if (track == 0x0) continue;
  776. Int_t nh = track->GetNofHits();
  777. cout << i << " " << nh << " " << ++nOK << endl;
  778. for (Int_t j = i + 1; j < nFound; ++j) {
  779. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  780. if (track1 == 0x0) continue;
  781. Int_t nh1 = track1->GetNofHits();
  782. Int_t nc = NofCommonHits(track, track1);
  783. if (1.*nc/TMath::Min(nh,nh1) < 0.5) continue;
  784. if (nh > nh1) fTracks->RemoveAt(j);
  785. else if (nh < nh1) {
  786. fTracks->RemoveAt(i);
  787. --nOK;
  788. break;
  789. } else {
  790. if (track->GetChi2() > track1->GetChi2()) {
  791. fTracks->RemoveAt(i);
  792. --nOK;
  793. break;
  794. }
  795. fTracks->RemoveAt(j);
  796. }
  797. }
  798. }
  799. fTracks->Compress();
  800. // Remove double tracks (originated from the same ETOF hit)
  801. /*
  802. nFound = fTracks->GetEntriesFast();
  803. for (Int_t i = 0; i < nFound; ++i) {
  804. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  805. if (track == 0x0) continue;
  806. Int_t iTof = track->GetTofIndex();
  807. Int_t nh = track->GetNofHits();
  808. for (Int_t j = i + 1; j < nFound; ++j) {
  809. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  810. if (track1 == 0x0) continue;
  811. Int_t iTof1 = track1->GetTofIndex();
  812. if (iTof1 != iTof) continue;
  813. Int_t nh1 = track1->GetNofHits();
  814. if (nh > nh1) fTracks->RemoveAt(j);
  815. else if (nh < nh1) {
  816. fTracks->RemoveAt(i);
  817. break;
  818. } else {
  819. if (track->GetChi2() > track1->GetChi2()) {
  820. fTracks->RemoveAt(i);
  821. break;
  822. }
  823. fTracks->RemoveAt(j);
  824. }
  825. }
  826. }
  827. fTracks->Compress();
  828. */
  829. }
  830. //__________________________________________________________________________
  831. Int_t MpdSftTrackFinderTpc::NofCommonHits(MpdEctKalmanTrack* track, MpdEctKalmanTrack* track1)
  832. {
  833. /// Compute number of common hits on 2 tracks
  834. TObjArray *hits = track->GetHits(), *hits1 = track1->GetHits();
  835. MpdKalmanHit *hit = (MpdKalmanHit*) hits->First();
  836. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->First();
  837. Int_t dir = 1;
  838. if (hit->GetLayer() < ((MpdKalmanHit*) hits->Last())->GetLayer()) dir = -1;
  839. Int_t nCom = 0;
  840. while (hit && hit1) {
  841. if (dir * hit->GetLayer() < dir * hit1->GetLayer()) {
  842. hit1 = (MpdKalmanHit*) hits1->After(hit1);
  843. continue;
  844. }
  845. if (hit == hit1) ++nCom;
  846. hit = (MpdKalmanHit*) hits->After(hit);
  847. }
  848. return nCom;
  849. }
  850. //__________________________________________________________________________
  851. void MpdSftTrackFinderTpc::AddHits()
  852. {
  853. /// Add hit objects to tracks and compute number of wrongly assigned hits
  854. /// (hits with ID different from ID of starting track)
  855. MpdTpcKalmanFilter *tpcKF = (MpdTpcKalmanFilter*) FairRun::Instance()->GetTask("TPC Kalman filter");
  856. Int_t nFound = fTracks->GetEntriesFast();
  857. //FairMCPoint *p;
  858. for (Int_t i = 0; i < nFound; ++i) {
  859. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  860. Int_t nHits = track->GetNofHits();
  861. if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  862. track->SetNofHits(nHits);
  863. TClonesArray &trHits = *track->GetTrHits();
  864. //SetTrackID(track); // set track ID as ID of majority of its hits
  865. TObjArray *hits = track->GetHits();
  866. Int_t nWrong = 0, motherID = track->GetTrackID();
  867. cout << i << " " << nHits << " " << motherID << endl;
  868. // Get track mother ID
  869. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  870. while (mctrack->GetMotherId() > 0) {
  871. motherID = mctrack->GetMotherId();
  872. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  873. }
  874. for (Int_t j = 0; j < nHits; ++j) {
  875. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  876. new (trHits[j]) MpdKalmanHit(*hit);
  877. //MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  878. //new (trHits[j]) MpdKalmanHit(*hit);
  879. cout << " " << hit->GetLayer();
  880. //Int_t iproj = (hit->GetLayer() - 1) % 3;
  881. Int_t iproj = -1, motherID1 = -1;
  882. //if (hit->GetType() == MpdKalmanHit::kFixedZ) p = (FairMCPoint*) fSftPoints->UncheckedAt(hit->GetIndex()); // SFT
  883. if (hit->GetType() == MpdKalmanHit::kFixedZ)
  884. motherID1 = ((FairMCPoint*) fSftPoints->UncheckedAt(hit->GetIndex()))->GetTrackID(); // SFT
  885. else {
  886. //p = (FairMCPoint*) fTpcPoints->UncheckedAt(hit->GetIndex()); // TPC
  887. motherID1 = tpcKF->GetHitID(hit);
  888. iproj = -1;
  889. }
  890. if (iproj == 0) cout << "U";
  891. else if (iproj == 1) cout << "R";
  892. else if (iproj > 0) cout << "V";
  893. //Int_t motherID1 = p->GetTrackID();
  894. cout << "-" << motherID1;
  895. // Get point mother ID
  896. MpdMCTrack *mctrack1 = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID1);
  897. while (mctrack1->GetMotherId() > 0) {
  898. motherID1 = mctrack1->GetMotherId();
  899. mctrack1 = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack1->GetMotherId());
  900. }
  901. //if (motherID1 != motherID && hit->GetMult() == 1) nWrong++;
  902. if (motherID1 != motherID) nWrong++;
  903. }
  904. cout << "\nWrongs: " << nWrong << " " << track->GetTrackID() << " " << motherID << endl;
  905. track->SetNofWrong(nWrong);
  906. track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
  907. }
  908. fTracks->Compress();
  909. }
  910. //__________________________________________________________________________
  911. void MpdSftTrackFinderTpc::SetTrackID(MpdEctKalmanTrack* track)
  912. {
  913. /// Set track ID as ID of majority of its hits
  914. const Int_t idMax = 99999;
  915. Int_t ids[idMax+1] = {0}, nHits = track->GetNofHits(), locMax = 0;
  916. TObjArray *hits = track->GetHits();
  917. for (Int_t i = 0; i < nHits; ++i) {
  918. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  919. FairMCPoint *p = (FairMCPoint*) fSftPoints->UncheckedAt(hit->GetIndex());
  920. Int_t id = p->GetTrackID();
  921. if (id > idMax) exit(0);
  922. ++ids[id];
  923. if (ids[id] > ids[locMax]) locMax = id;
  924. }
  925. if (ids[locMax] > 1) track->SetTrackID(locMax);
  926. }
  927. //__________________________________________________________________________
  928. void MpdSftTrackFinderTpc::RemoveTpcTracks()
  929. {
  930. /// Remove TPC tracks which are parts of the combined ones
  931. MpdTpcKalmanFilter *tpcKF = (MpdTpcKalmanFilter*) FairRun::Instance()->GetTask("TPC Kalman filter");
  932. TClonesArray *tpcTracks = tpcKF->GetTracks();
  933. Int_t nTr = fTracks->GetEntriesFast(), nTpc = tpcTracks->GetEntriesFast();
  934. for (Int_t i = 0; i < nTr; ++i) {
  935. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  936. if (tr == nullptr) continue;
  937. //Int_t id = tr->GetTrackID();
  938. if (tr->GetTpcIndex() >= 0) {
  939. // There are TPC hits
  940. for (Int_t j = 0; j < nTpc; ++j) {
  941. MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) tpcTracks->UncheckedAt(j);
  942. if (tpc == nullptr) continue;
  943. if (tpc->GetTrackID() != tr->GetTrackID()) continue;
  944. if (tpc->GetNofTrHits() < tr->GetNofTrHits()) {
  945. tpcTracks->RemoveAt(j);
  946. } else {
  947. cout << tpc->GetNofTrHits() << " " << tr->GetNofTrHits() << " " << tr->GetTrackID()
  948. << " " << tpc->Momentum3().Eta() << endl;
  949. //Fatal("RemoveTpcTrack", "Strange ???");
  950. Warning("RemoveTpcTrack", "Strange ???");
  951. fTracks->RemoveAt(i);
  952. }
  953. }
  954. }
  955. }
  956. tpcTracks->Compress();
  957. fTracks->Compress();
  958. }
  959. //__________________________________________________________________________
  960. void MpdSftTrackFinderTpc::GoToVertex()
  961. {
  962. /// Propagate tracks to the vertex
  963. //return;
  964. Int_t nTr = fTracks->GetEntriesFast();
  965. // Save last parameters
  966. /*
  967. for (Int_t i = 0; i < nTr; ++i) {
  968. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  969. //cout << track->GetPos() << " " << track->GetPosNew() << endl;
  970. track->SetPos(track->GetPosNew());
  971. track->SetParam(*track->GetParamNew());
  972. }
  973. return;
  974. */
  975. TMatrixD param(5,1);
  976. TMatrixDSym weight(5), pointWeight(5);
  977. MpdKalmanHit hitv;
  978. hitv.SetType(MpdKalmanHit::kFixedZ);
  979. hitv.SetErr(1,0.0005); // 50um
  980. hitv.SetErr(0,0.0005); // 50um
  981. //hitv.SetRerr(0.05); // 1mm
  982. //hitv.SetRphiErr(0.05); // 1mm
  983. //hitv.SetPhi(TMath::PiOver2());
  984. hitv.SetCosSin(0,0.0);
  985. hitv.SetCosSin(1,1.0);
  986. hitv.SetNofDim(2);
  987. for (Int_t i = 0; i < nTr; ++i) {
  988. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  989. cout << " Refit0: " << i << " " << track->GetNofTrHits() << " " << track->GetChi2() << " "
  990. << 1./track->GetParam(4) << " " << 1./track->GetParamNew(4) << endl;
  991. //MpdKalmanFilter::Instance()->Refit(track, -1);
  992. cout << " Refit1: " << track->GetNofTrHits() << " " << track->GetChi2() << " "
  993. << 1./track->GetParam(4) << " " << 1./track->GetParamNew(4) << endl;
  994. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitv);
  995. //TVector2 dist = GetDistance(track,&hitv);
  996. //cout << dist.X() << " " << dist.Y() << endl;
  997. cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  998. /*
  999. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(track,&hitv,pointWeight,param);
  1000. track->SetChi2(track->GetChi2()+dChi2);
  1001. weight = *track->GetWeight();
  1002. weight += pointWeight;
  1003. track->SetWeight(weight);
  1004. track->SetParam(param);
  1005. cout << " Vertex: " << 1./track->GetParam(4) << endl;
  1006. */
  1007. }
  1008. }
  1009. //__________________________________________________________________________
  1010. void MpdSftTrackFinderTpc::PassInnerShell(MpdKalmanTrack* track)
  1011. {
  1012. /// Propagate track thru TPC inner shell (include MS)
  1013. const Int_t nR = 7;
  1014. const Double_t rad[nR] = {20.035, 21.07, 22.105, 24.64, 27.15, 28.165, 29.180};
  1015. const Double_t dx[nR] = {0.0031, 0.0017, 0.0031, 0.0003, 0.001, 0.0017, 0.0013};
  1016. MpdKalmanHit hit;
  1017. hit.SetType(MpdKalmanHit::kFixedR);
  1018. for (Int_t i = 0; i < nR; ++i) {
  1019. hit.SetDist(rad[i]);
  1020. MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kTRUE);
  1021. // Add multiple scattering
  1022. TMatrixDSym *cov = track->Weight2Cov();
  1023. Double_t th = track->GetParamNew(3);
  1024. Double_t cosTh = TMath::Cos(th);
  1025. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[i]);
  1026. //cov->Print();
  1027. (*cov)(2,2) += (angle2 / cosTh / cosTh );
  1028. (*cov)(3,3) += angle2;
  1029. //cov->Print();
  1030. Int_t iok = 0;
  1031. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1032. track->SetWeight(*cov);
  1033. }
  1034. }
  1035. ClassImp(MpdSftTrackFinderTpc);