MpdTpcTrackFollow2Sft.cxx 44 KB

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