MpdEctTrackFinderTofTpc.cxx 34 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880
  1. // -------------------------------------------------------------------------
  2. // ----- MpdEctTrackFinderTofTpc source file -----
  3. // ----- Created 17/06/08 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdEctTrackFinderTof.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD End-Cap Tracker (ECT) using seeds built from TOF, TPC
  9. ** and vertex points
  10. **/
  11. #include "MpdEctTrackFinderTofTpc.h"
  12. #include "MpdKalmanFilter.h"
  13. #include "MpdKalmanHitZ.h"
  14. #include "MpdEctKalmanTrack.h"
  15. #include "MpdStrawendcapPoint.h"
  16. #include "MpdEtofPoint.h"
  17. #include "TpcLheHit.h"
  18. #include "TpcLheKalmanTrack.h"
  19. #include "FairMCTrack.h"
  20. #include "FairRootManager.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. const Double_t MpdEctTrackFinderTofTpc::fgkChi2Cut = 20; //20; //100;
  33. FILE *lun = 0x0; //fopen("error.dat","w");
  34. //__________________________________________________________________________
  35. MpdEctTrackFinderTofTpc::MpdEctTrackFinderTofTpc(const char *name, Int_t iVerbose )
  36. :FairTask(name, iVerbose)
  37. {
  38. fKHits = new TClonesArray("MpdKalmanHitZ", 100);
  39. fTracks = new TClonesArray("MpdEctKalmanTrack", 100);
  40. fHistoDir = 0x0;
  41. fhLays = new TH1F("hLays","ECT layers",100,0,100);
  42. fLayPointers = 0x0;
  43. fTPC = new TObjArray(300);
  44. }
  45. //__________________________________________________________________________
  46. MpdEctTrackFinderTofTpc::~MpdEctTrackFinderTofTpc()
  47. {
  48. //delete fKHits;
  49. //delete fTracks;
  50. //delete fTrackCand;
  51. delete [] fLayPointers;
  52. delete fhLays;
  53. delete fTPC;
  54. }
  55. //__________________________________________________________________________
  56. InitStatus MpdEctTrackFinderTofTpc::Init()
  57. {
  58. return ReInit();
  59. }
  60. //__________________________________________________________________________
  61. InitStatus MpdEctTrackFinderTofTpc::ReInit()
  62. {
  63. fTpcHits = (TClonesArray *) FairRootManager::Instance()->GetObject("LheHit");
  64. fEctHits =(TClonesArray *) FairRootManager::Instance()->GetObject("STRAWPoint");
  65. fTofHits =(TClonesArray *) FairRootManager::Instance()->GetObject("ETOFPoint");
  66. //fTpcTracks = 0x0; //(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  67. fTpcTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("TpcPoint"); // just for testing
  68. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  69. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  70. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  71. FairRootManager::Instance()->Register("EctTrack", "Ect", fTracks, kTRUE);
  72. fNPass = 1;
  73. return kSUCCESS;
  74. }
  75. //__________________________________________________________________________
  76. void MpdEctTrackFinderTofTpc::Reset()
  77. {
  78. ///
  79. cout << " MpdEctTrackFinderTof::Reset " << endl;
  80. fKHits->Clear();
  81. fTracks->Clear();
  82. fTPC->Clear();
  83. delete [] fLayPointers;
  84. }
  85. //__________________________________________________________________________
  86. void MpdEctTrackFinderTofTpc::SetParContainers()
  87. {
  88. }
  89. //__________________________________________________________________________
  90. void MpdEctTrackFinderTofTpc::Finish()
  91. {
  92. //Write();
  93. }
  94. //__________________________________________________________________________
  95. void MpdEctTrackFinderTofTpc::Exec(Option_t * option)
  96. {
  97. static int eventCounter = 0;
  98. cout << " EctRec event " << ++eventCounter << endl;
  99. Reset();
  100. // Create Kalman hits
  101. MakeKalmanHits();
  102. for (Int_t i = 0; i < fNPass; ++i) {
  103. fTracks->Clear();
  104. GetTrackSeeds(i);
  105. cout << " Total number of hits for tracking: " << fKHits->GetEntriesFast() << endl;
  106. cout << " Total number of track seeds: " << fTracks->GetEntriesFast() << endl;
  107. DoTracking(i);
  108. //StoreTracks();
  109. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  110. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  111. }
  112. RemoveDoubles(); // remove double tracks
  113. AddHits(); // add hit objects to tracks
  114. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  115. }
  116. //__________________________________________________________________________
  117. void MpdEctTrackFinderTofTpc::MakeKalmanHits()
  118. {
  119. /// Create Kalman hits from ECT points.
  120. fhLays->Reset();
  121. Int_t nHits = fEctHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0;
  122. Double_t phi, r, coord, errR = 0.2, errRphi = 0.02; // 2mm in R, 200um in R-Phi
  123. //Double_t phi, r, coord, errR = 0.2, errRphi = 0.005; // 2mm in R, 50um in R-Phi
  124. for (Int_t ih = 0; ih < nHits; ++ih) {
  125. MpdStrawendcapPoint *h = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(ih);
  126. //if (h->GetTrackID() != 1265) continue; // !!!
  127. //if (h->GetTrackID() != 64) continue; // !!!
  128. //if (h->GetTrackID() != 1574) continue; // !!!
  129. //if (h->GetTrackID() != 1427) continue; // !!!
  130. //if (h->GetTrackID() != 64 && h->GetTrackID() != 92) continue; // !!!
  131. phi = TMath::ATan2 (h->GetY(), h->GetX()); // tube Phi - radial case
  132. phi = h->GetPhi(); // tube Phi
  133. lay = h->GetDetectorID() / 1000;
  134. // Extrapolate track to Z = Ztube
  135. Double_t dZ = h->GetZ() - h->GetTrackZ();
  136. Double_t dt = 0.; // dZ / h->GetPz();
  137. //if (TMath::Abs(h->GetPz()) > 1.e-6 && h->GetPz() * h->GetZ() > 0) dt = dZ / h->GetPz();
  138. Double_t xNew = h->GetTrackX() + dt * h->GetPx();
  139. Double_t yNew = h->GetTrackY() + dt * h->GetPy();
  140. //cout << dZ << " " << h->GetTrackX() << " " << xNew << " " << h->GetTrackY() << " " << yNew << " " << lay << endl;
  141. //Double_t zNew = h->GetTrackZ() + dt * h->GetPz(); // just for cross-check
  142. // Transform to the tube local coordinate system
  143. Double_t cosPhi = TMath::Cos(phi);
  144. Double_t sinPhi = TMath::Sin(phi);
  145. //Double_t xLoc = h->GetX() * cosPhi + h->GetY() * sinPhi; // cross-check
  146. //Double_t yLoc = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  147. Double_t xLoc = xNew * cosPhi + yNew * sinPhi;
  148. Double_t yLoc = -xNew * sinPhi + yNew * cosPhi;
  149. //Double_t xLoc = (xNew - h->GetX()) * cosPhi + (yNew - h->GetY()) * sinPhi;
  150. //Double_t yLoc = -(xNew - h->GetX()) * sinPhi + (yNew - h->GetY()) * cosPhi;
  151. //r = xNew * xNew + yNew * yNew;
  152. //r = TMath::Sqrt (r);
  153. //r = TMath::Abs(xLoc);
  154. r = xLoc;
  155. //cout << xLoc << " " << yLoc << " " << r << " " << h->GetPz() << endl;
  156. coord = yLoc;
  157. // Add error
  158. Double_t dRphi = 0, dR = 0;
  159. gRandom->Rannor(dRphi,dR); // add errors
  160. MpdKalmanHitZ *hit = new ((*fKHits)[nKH++]) MpdKalmanHitZ(r+dR*errR, phi, coord+dRphi*errRphi,
  161. h->GetTrackZ(), errR, errRphi, lay, ih);
  162. hit->SetXY(h->GetX(), h->GetY());
  163. hit->SetSin(sinPhi);
  164. hit->SetCos(cosPhi);
  165. // Add second measurement - just for test at the moment
  166. //!!!
  167. //hit->SetNofDim(2);
  168. //!!!
  169. //lay = hit.GetLayer();
  170. layMax = TMath::Max (lay, layMax);
  171. fhLays->Fill(lay+0.1);
  172. }
  173. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << endl;
  174. fKHits->Sort(); // in ascending order in abs(Z)
  175. fLayPointers = new Int_t [layMax+1];
  176. Int_t ipos = 0;
  177. for (Int_t i = 0; i <= layMax; ++i) {
  178. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  179. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  180. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  181. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  182. if (cont == 0) {
  183. fLayPointers[i] = -1;
  184. continue;
  185. }
  186. fLayPointers[i] = ipos;
  187. ipos += cont;
  188. }
  189. }
  190. //__________________________________________________________________________
  191. void MpdEctTrackFinderTofTpc::GetTrackSeeds(Int_t iPass)
  192. {
  193. /// Build track seeds from TOF hits, TPC hits at max Z and vertex
  194. Int_t tmp[99999] = {0}, idMax = 0, nID = 0;
  195. GetTpcPoints(); // select TPC points near zMax
  196. //GetTofPoints(); //
  197. Int_t nTpc = fTPC->GetEntriesFast();
  198. Int_t nTof = fTofHits->GetEntriesFast();
  199. Int_t nCand = 0;
  200. //const Double_t dphiMax = 0.2, dslMax = 0.1; // cuts dPhi and dSlope
  201. //const Double_t dphiMax = 0.25, dslMax = 0.1; // cuts dPhi and dSlope
  202. const Double_t dphiMax = 0.5, dslMax = 0.1; // cuts dPhi and dSlope
  203. TVector3 vert(0.,0.,0.);
  204. for (Int_t itof = 0; itof < nTof; ++itof) {
  205. MpdEtofPoint *tof = (MpdEtofPoint*) fTofHits->UncheckedAt(itof);
  206. //if (tof->GetTrackID() != 64) continue; // !!!
  207. //if (tof->GetTrackID() != 1265) continue; // !!!
  208. Double_t rTof = TMath::Sqrt (tof->GetX() * tof->GetX() + tof->GetY() * tof->GetY());
  209. Double_t slope0 = rTof / tof->GetZ();
  210. Double_t phTof = TMath::ATan2 (tof->GetY(), tof->GetX());
  211. for (Int_t itpc = 0; itpc < nTpc; ++itpc) {
  212. TpcLheHit *tpc = (TpcLheHit*) fTPC->UncheckedAt(itpc);
  213. // For testing
  214. TpcPoint *p = 0x0;
  215. Int_t idTpc = -1;
  216. if (fTpcTracks) p = (TpcPoint*) fTpcTracks->UncheckedAt(tpc->GetRefIndex());
  217. if (p) idTpc = p->GetTrackID();
  218. //if (idTpc != tof->GetTrackID()) continue; // !!! for test - keep only correct combinations
  219. Double_t rTpc = tpc->GetR();
  220. if (rTpc > rTof) continue;
  221. Double_t phTpc = tpc->GetRphi() / rTpc;
  222. Double_t dphi = MpdKalmanFilter::Instance()->Proxim(phTof, phTpc) - phTof;
  223. if (TMath::Abs(dphi) > dphiMax) continue;
  224. Double_t slope = tpc->GetR() / tpc->GetZ();
  225. if (TMath::Abs(slope-slope0) > dslMax) continue;
  226. if (lun && idTpc == tof->GetTrackID()) fprintf(lun,"%10.3e %10.3e %10.3e %10.3e\n",rTpc,rTof,dphi,slope-slope0);
  227. if (tof->GetTrackID() < 99999) { tmp[tof->GetTrackID()]++; idMax = TMath::Max(idMax,tof->GetTrackID()); }
  228. MpdEctKalmanTrack *track = new ((*fTracks)[nCand++]) MpdEctKalmanTrack(itof, itpc, tof, tpc, vert);
  229. track->SetDirection(MpdKalmanTrack::kOutward);
  230. track->SetPos(tpc->GetZ());
  231. track->SetPosNew(tpc->GetZ());
  232. EvalParams(tof, tpc, track);
  233. // Eval. track Phi - linear extrapolation of 2 estimates
  234. //track->SetParam(2, phTof); // phi
  235. Double_t dx = tof->GetX() - rTpc*TMath::Cos(phTpc);
  236. Double_t dy = tof->GetY() - rTpc*TMath::Sin(phTpc);
  237. Double_t dr = TMath::Sqrt (dx * dx + dy * dy);
  238. Double_t phTofTpc = TMath::ATan2 (dy, dx);
  239. Double_t ph = phTpc - 1.0 * (MpdKalmanFilter::Instance()->Proxim(phTofTpc,phTof)-phTofTpc)/(rTof-dr) * rTpc;
  240. track->SetParam(2, ph); // phi
  241. //track->SetParam(2, TMath::ATan2(p->GetPy(),p->GetPx())); // exact value - for test
  242. Double_t th = TMath::ATan2(tof->GetZ()-tpc->GetZ(),dr);
  243. track->SetParam(3, th); // dip angle
  244. //track->SetParam(3, TMath::ATan2(p->GetPz(),TMath::Sqrt(p->GetPx()*p->GetPx()+p->GetPy()*p->GetPy()))); // exact value - for test
  245. track->SetParamNew(*track->GetParam());
  246. track->GetParam()->Print();
  247. EvalCovar(tof, tpc, track);
  248. //if (lun && idTpc == tof->GetTrackID()) fprintf(lun,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",track->GetParam(2),ph,track->GetParam(3),th,dphi,TMath::Sqrt((*track->GetWeight())(2,2)),TMath::Sqrt((*track->GetWeight())(3,3)));
  249. cout << nCand-1 << " " << track->GetTrackID() << " " << idTpc << endl;
  250. }
  251. }
  252. for (Int_t j = 0; j <= idMax; ++j) if (tmp[j] > 0) ++nID;
  253. cout << " Number of ECT track candidates: " << nCand << " " << nID << endl;
  254. if (lun) fclose(lun);
  255. //exit(0);
  256. }
  257. //__________________________________________________________________________
  258. void MpdEctTrackFinderTofTpc::GetTpcPoints()
  259. {
  260. /// Get TPC hits near zMax
  261. Int_t ids[99999] = {0}, idMax = 0;
  262. Double_t zMax = 134.5, rMin = 20.5, rMax = 109.5, dR = (rMax-rMin)/50; // dR = 1.78; // old version
  263. Int_t nTpc = fTpcHits->GetEntriesFast(), lay0 = -1;
  264. cout << " TPC hits: " << nTpc << endl;
  265. Double_t dz = 0.;
  266. //fTpcHits->Sort();
  267. for (Int_t i = 0; i < nTpc; ++i) {
  268. TpcLheHit *hit = (TpcLheHit*) fTpcHits->UncheckedAt(i);
  269. TpcPoint *p = 0x0;
  270. if (fTpcTracks) p = (TpcPoint*) fTpcTracks->UncheckedAt(hit->GetRefIndex());
  271. //if (p && p->GetTrackID() != 64) continue; // !!!
  272. if (hit->GetZ() < 0) continue;
  273. if (hit->GetLayer() != lay0) {
  274. lay0 = hit->GetLayer();
  275. dz = zMax / hit->GetR() * dR / 2.;
  276. }
  277. if (hit->GetZ() < zMax - dz) continue;
  278. //Int_t id = ((TpcPoint*)fTpcTracks->UncheckedAt(hit->GetRefIndex()))->GetTrackID();
  279. Int_t id = hit->GetTrackID();
  280. cout << i << " " << hit->GetLayer() << " " << hit->GetZ() << " " << dz << " " << id << endl;
  281. fTPC->Add(hit);
  282. if (id > 99998) continue;
  283. idMax = TMath::Max (id, idMax);
  284. ids[id]++;
  285. }
  286. Int_t nid = 0;
  287. for (Int_t i = 0; i <= idMax; ++i) if (ids[i]) ++nid;
  288. cout << " TPC points for seeding: " << fTPC->GetEntriesFast() << " " << nid << endl;
  289. }
  290. //__________________________________________________________________________
  291. void MpdEctTrackFinderTofTpc::EvalParams(const MpdEtofPoint *tof, const TpcLheHit *tpc,
  292. MpdEctKalmanTrack *track)
  293. {
  294. /// Evaluate track parameters
  295. // Evaluate signed track Pt (curvature) assuming the track coming from the
  296. // primary vertex
  297. Double_t rTof = TMath::Sqrt (tof->GetX() * tof->GetX() + tof->GetY() * tof->GetY());
  298. Double_t phTof = TMath::ATan2 (tof->GetY(), tof->GetX());
  299. Double_t phTpc = tpc->GetRphi() / tpc->GetR();
  300. TVector2 vecTof(rTof*TMath::Cos(phTof)-0.,rTof*TMath::Sin(phTof)-0.);
  301. TVector2 vecTpc(tpc->GetR()*TMath::Cos(phTpc)-0.,tpc->GetR()*TMath::Sin(phTpc)-0.);
  302. TVector2 dvec = vecTof - vecTpc;
  303. Double_t cosAlpha = vecTpc * dvec / vecTpc.Mod() / dvec.Mod();
  304. Double_t rad = vecTof.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  305. Double_t factor = 0.0015; // 0.3 * 0.5T * 0.01
  306. //Double_t charge = phTof - MpdKalmanFilter::Instance()->Proxim(phTof,phTpc);
  307. Double_t charge = MpdKalmanFilter::Instance()->Proxim(phTof,phTpc) - phTof;
  308. if (rTof < tpc->GetR()) charge = -charge;
  309. Double_t pt = factor * TMath::Abs(rad) * TMath::Sign(1., -charge);
  310. track->SetParam(0, vecTpc.X()); // X
  311. track->SetParam(1, vecTpc.Y()); // Y
  312. track->SetParam(4, 1./pt); // 1/pt
  313. }
  314. //__________________________________________________________________________
  315. void MpdEctTrackFinderTofTpc::EvalCovar(const MpdEtofPoint *tof, const TpcLheHit *tpc,
  316. MpdEctKalmanTrack *track)
  317. {
  318. /// Evaluate covariance matrix for the track seed
  319. static const Double_t xErrTof = 2. / TMath::Sqrt(12.), yErrTof = xErrTof;
  320. //static const Double_t xErrTof = 0.02, yErrTof = xErrTof;
  321. TMatrixDSym weight(5);
  322. //weight(0,0) = xErrTof * xErrTof; // <xx>
  323. weight(0,0) = tpc->GetRphiErr() * tpc->GetRphiErr(); // <xx>
  324. weight(0,0) *= 4.; // extra factor of 2
  325. //weight(1,1) = yErrTof * yErrTof; // <yy>
  326. weight(1,1) = tpc->GetRphiErr() * tpc->GetRphiErr(); // <yy>
  327. weight(1,1) *= 4.; // extra factor of 2
  328. Double_t phTpc = tpc->GetRphi() / tpc->GetR();
  329. Double_t xTpc = tpc->GetR() * TMath::Cos(phTpc);
  330. Double_t yTpc = tpc->GetR() * TMath::Sin(phTpc);
  331. Double_t xTof = tof->GetX(), yTof = tof->GetY();
  332. Double_t dx = xTof - xTpc, dy = yTof - yTpc;
  333. Double_t dist2 = dx * dx + dy * dy;
  334. Double_t sinPhi = TMath::Sin (track->GetParam(2));
  335. Double_t cosPhi = TMath::Cos (track->GetParam(2));
  336. Double_t pTpc = TMath::Cos(phTpc) * cosPhi + TMath::Sin(phTpc) * sinPhi;
  337. weight(2,2) = pTpc * pTpc * tpc->GetRphiErr() * tpc->GetRphiErr();
  338. weight(2,2) += xErrTof * xErrTof;
  339. weight(2,2) /= dist2; // <PhiPhi>
  340. //weight(2,2) *= 4.; // extra factor of 2
  341. weight(2,2) *= 2.; // extra factor of sqrt(2)
  342. Double_t tanThe = TMath::Tan (track->GetParam(3));
  343. Double_t rTof = TMath::Sqrt (tof->GetX() * tof->GetX() + tof->GetY() * tof->GetY());
  344. Double_t dR = rTof - tpc->GetR();
  345. Double_t denom = dR * (1. + tanThe * tanThe);
  346. weight(3,3) = tpc->GetZerr() * tpc->GetZerr() / denom / denom;
  347. Double_t dz = tof->GetZ() - tpc->GetZ();
  348. Double_t dThdR = dz / (dR * dR + dz * dz);
  349. weight(3,3) += dThdR * dThdR * xErrTof * xErrTof; // <TheThe>
  350. //weight(3,3) *= 4.; // extra factor of 2
  351. weight(3,3) *= 8.; // extra factor of 2
  352. //weight(4,4) = (track->GetParam(4)*0.5) * (track->GetParam(4)*0.5); // error 50%
  353. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
  354. weight(4,4) = (track->GetParam(4)*1.) * (track->GetParam(4)*1.); // error 100%
  355. //weight(4,4) = (track->GetParam(4)*2.) * (track->GetParam(4)*2.); // error 200%
  356. //weight(4,4) = (track->GetParam(4)*4.) * (track->GetParam(4)*4.); // error 400%
  357. weight.Print();
  358. //fWeight->Invert(); // weight matrix
  359. Int_t iok = 0;
  360. MpdKalmanFilter::Instance()->MnvertLocal(weight.GetMatrixArray(), 5, 5, 5, iok);
  361. //fWeight->Print();
  362. track->SetWeight(weight);
  363. }
  364. //__________________________________________________________________________
  365. void MpdEctTrackFinderTofTpc::DoTracking(Int_t iPass)
  366. {
  367. /// Run Kalman tracking
  368. Int_t nCand = fTracks->GetEntriesFast(), iok = 0;
  369. for (Int_t i = 0; i < nCand; ++i) {
  370. cout << " Track seed No. " << i << endl;
  371. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  372. //if (track->GetTrackID() != 117) continue;
  373. //if (track->GetTrackID() != 10) continue;
  374. //if (track->GetTrackID() != 1105) continue;
  375. //if (track->GetTrackID() != 77) continue;
  376. //(*(track->GetParamNew()))(4,0) = -0.5; // just for check
  377. /*
  378. for (Int_t k = 0; k < 5; ++k) {
  379. for (Int_t j = i; j < 5; ++j) {
  380. if (j == i) (*track->GetWeight())(i,j) /= 100.;
  381. else (*track->GetWeight())(i,j) = (*track->GetWeight())(j,i) = 0.;
  382. }
  383. }
  384. */
  385. // Start ECT tracking from different layers to account for missing hits
  386. const Int_t laySkip = 10; //10;
  387. Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  388. MpdEctKalmanTrack tracks[laySkip];
  389. for (Int_t j = 0; j < laySkip; ++j) {
  390. tracks[j] = *track;
  391. //cout << track->GetParamNew(0) << endl;
  392. //cout << i << " " << lay << " " << tracks[lay].GetNofHits() << " " << tracks[lay].GetChi2() << " " << tracks[lay].GetParam(0) << endl;
  393. Int_t lay = 1;
  394. if (j > 0 && tracks[j-1].GetNofHits() > 0)
  395. lay = ((MpdKalmanHit*)tracks[j-1].GetHits()->First())->GetLayer() + 1;
  396. iok = RunKalmanFilter(&tracks[j], lay);
  397. //iok = RunKalmanFilter(&tracks[lay], 0);
  398. //cout << i << " " << lay << " " << tracks[lay].GetNofHits() << " " << tracks[lay].GetChi2() << endl;
  399. }
  400. // Select the best track (with max number of hits)
  401. Int_t nHitsMax = tracks[0].GetNofHits(), iMax = 0;
  402. for (Int_t j = 1; j < laySkip; ++j) {
  403. Int_t nhits = tracks[j].GetNofHits();
  404. if (nhits > nHitsMax) {
  405. nHitsMax = nhits;
  406. iMax = j;
  407. } else if (nhits == nHitsMax) {
  408. if (tracks[j].GetChi2() < tracks[iMax].GetChi2()) {
  409. iMax = j;
  410. nHitsMax = nhits;
  411. }
  412. }
  413. }
  414. *track = tracks[iMax];
  415. track->GetParamNew()->Print();
  416. if (0) {
  417. track->SetParam(*track->GetParamNew());
  418. track->SetParam(4,track->GetParamNew(4)*1.5);
  419. track->SetParamNew(*track->GetParam());
  420. MpdKalmanFilter::Instance()->Refit(track, -1);
  421. MpdKalmanFilter::Instance()->Refit(track, 1);
  422. }
  423. //iok = RunKalmanFilter(track);
  424. cout << i << " " << track->GetNofHits() << endl;
  425. if (0 && iok == 0 && track->GetNofHits() < 0) {
  426. //*
  427. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) track->GetHits()->UncheckedAt(0);
  428. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  429. track->SetPos(hit->GetPos()-dir);
  430. track->SetPosNew(track->GetPos());
  431. //for (Int_t j = 0; j < 5; ++j) track->SetParam(j,track->GetParam1()[j]);
  432. //track->SetParamNew(*track->GetParam());
  433. //MpdKalmanHitZ hitTmp(*hit);
  434. //hitTmp.SetZ(track->GetPosNew()-dir);
  435. //MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  436. MpdKalmanFilter::Instance()->Refit(track,1);
  437. //*/
  438. /*
  439. track->Print("");
  440. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  441. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) track->GetHits()->UncheckedAt(0);
  442. MpdKalmanHitZ hitTmp(*hit);
  443. hitTmp.SetZ(track->GetPosNew()+dir);
  444. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  445. track->Print("");
  446. MpdKalmanFilter::Instance()->Refit(track,-1);
  447. track->Print("");
  448. hitTmp.SetZ(track->GetPosNew()-dir);
  449. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  450. track->Print("");
  451. MpdKalmanFilter::Instance()->Refit(track,1);
  452. track->Print("");
  453. */
  454. }
  455. }
  456. }
  457. //__________________________________________________________________________
  458. Int_t MpdEctTrackFinderTofTpc::RunKalmanFilter(MpdEctKalmanTrack *track, Int_t layBeg)
  459. {
  460. /// Run Kalman filter
  461. const Double_t rMin = 28., rMax = 123.; // min and max radial size of ECT - to be taken elsewhere
  462. //cout << fHits->GetEntriesFast() << endl;
  463. Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  464. MpdKalmanHitZ *hitOK = 0x0;
  465. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  466. //Int_t layBeg = layMax, layEnd = -1, dLay = -1, layOK = -1;
  467. Int_t layEnd = -1, dLay = -1, layOK = -1;
  468. if (trackDir == MpdKalmanTrack::kOutward) {
  469. layEnd = layMax + 1;
  470. dLay = 1;
  471. }
  472. //Int_t indxOK = hits->IndexOf(hitOK);
  473. //Int_t nHits = hits->GetEntriesFast();
  474. Int_t miss = 0;
  475. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  476. TMatrixD param(5,1), paramTmp(5,1);
  477. Double_t saveZ = 0;
  478. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  479. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  480. Int_t nLay = GetNofHitsInLayer(lay);
  481. Int_t indx0 = GetHitsInLayer(lay);
  482. Double_t dChi2Min = 1.e+6;
  483. MpdKalmanHitZ *hitMin = 0x0;
  484. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  485. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  486. /*
  487. if (trackDir == MpdKalmanTrack::kOutward) {
  488. // !!! This part is due to the simplified hit merging (digitization) -
  489. // different hit position inside one layer - should be removed later
  490. indxBeg = nLay - 1;
  491. indxEnd = -1;
  492. dIndx = -1;
  493. }
  494. */
  495. //for (Int_t indx = 0; indx < nLay; ++indx) {
  496. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  497. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) fKHits->UncheckedAt(indx0+indx);
  498. // !!! Exact ID match
  499. //if (((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  500. // Exclude used hits
  501. if (hit->GetFlag() != 1) continue;
  502. //cout << " window:" << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << hit->GetZ() << " " << track->GetParamNew(1) << endl;
  503. // Check if the hit within some window (15x15cm for the moment - check!!!)
  504. //if (TMath::Abs(hit->GetRphi()-track->GetParamNew(0)) > 9) continue;
  505. //if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  506. TVector2 dist = GetDistance(track, hit);
  507. if (dist.X() > 15.) continue; // distance in transverse to the tube direction
  508. if (hit->GetNofDim() > 1 && dist.Y() > 15.) continue; // distance in R
  509. //*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;
  510. //track->Print("");
  511. //hit->Print("");
  512. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) return -1;
  513. //*if (hit->GetTrackID() == -607)*/ cout << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << " " << hit->GetZ() << " " << track->GetPosNew() << endl;
  514. //
  515. // For debugging
  516. /*
  517. TVector2 dist0 = GetDistance(track, hit);
  518. cout << dist0.X() << " ";
  519. MpdKalmanHitZ hitDbg = *hit;
  520. Double_t xDbg = hit->GetXY(0) * TMath::Cos(hit->GetPhi()) + hit->GetXY(1) * TMath::Sin(hit->GetPhi());
  521. Double_t yDbg = -hit->GetXY(0) * TMath::Sin(hit->GetPhi()) + hit->GetXY(1) * TMath::Cos(hit->GetPhi());
  522. hitDbg.SetRphi(yDbg);
  523. hitDbg.SetR(xDbg);
  524. dist = GetDistance(track, &hitDbg);
  525. cout << dist.X() << endl;
  526. */
  527. Double_t radNew = track->GetRadNew();
  528. if (radNew < rMin || radNew > rMax) return -1;
  529. //Double_t err = hit->GetRphiErr();
  530. //if (track->GetNofHits() == 0) hit->SetRphiErr(0.04);
  531. Double_t dChi2 = MpdKalmanFilter::Instance()->CheckHitZ(track,hit,pointWeight,param);
  532. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  533. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  534. if (dChi2 < dChi2Min) {
  535. //if (dChi2 < dChi2Min && dist0.X() <= dist.X()) {
  536. //if (dChi2 < dChi2Min && dist.X() <= 0.2) {
  537. dChi2Min = dChi2;
  538. hitMin = hit;
  539. saveWeight = *track->GetWeight();
  540. saveZ = track->GetPosNew();
  541. // temporary storage for the current track
  542. paramTmp = param;
  543. pointWeightTmp = pointWeight;
  544. //cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetR() << endl;
  545. //cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  546. //cout << hit->GetRphi() << " " << hit->GetZ() << endl;
  547. //cout << param(0,0) << " " << param(1,0) << endl;
  548. //paramTmp.Print();
  549. }
  550. }
  551. Double_t cut = fgkChi2Cut;
  552. //if (track->GetNofHits() == 0) cut /= 2.;
  553. //if (dChi2Min < fgkChi2Cut) {
  554. if (dChi2Min < cut) {
  555. //layOK = lay;
  556. hitOK = hitMin;
  557. track->GetHits()->Add(hitOK);
  558. miss = 0;
  559. // Restore track params at the best hit
  560. track->SetChi2(track->GetChi2()+dChi2Min);
  561. saveWeight += pointWeightTmp;
  562. track->SetWeight(saveWeight);
  563. track->SetPosNew(saveZ);
  564. track->SetParamNew(paramTmp);
  565. //cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  566. //paramTmp.Print();
  567. // Check if the accepted hit is the same as the seeded hit
  568. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  569. if (track->GetNofHits() == 1) track->SetParam1();
  570. } else {
  571. ++miss;
  572. //if (miss > 1) return -1;
  573. }
  574. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  575. } // for (Int_t lay = layOK-1; lay >= 0;
  576. return 0;
  577. }
  578. //__________________________________________________________________________
  579. TVector2 MpdEctTrackFinderTofTpc::GetDistance(MpdEctKalmanTrack *track, MpdKalmanHitZ *hit)
  580. {
  581. /// Compute distance between track and hit
  582. Double_t xTr, yTr;
  583. if (track->GetType() == MpdKalmanTrack::kFixedZ) {
  584. xTr = track->GetParamNew(0);
  585. yTr = track->GetParamNew(1);
  586. } else {
  587. Double_t rPhi = track->GetParamNew(0);
  588. Double_t r = track->GetPosNew();
  589. Double_t ph = rPhi / r;
  590. xTr = r * TMath::Cos(ph);
  591. yTr = r * TMath::Sin(ph);
  592. }
  593. Double_t phi = hit->GetPhi();
  594. Double_t cosPhi = TMath::Cos(phi);
  595. Double_t sinPhi = TMath::Sin(phi);
  596. // Transform track coordinates to local tube coordinates
  597. Double_t xLoc = xTr * cosPhi + yTr * sinPhi;
  598. Double_t yLoc = -xTr * sinPhi + yTr * cosPhi;
  599. //Double_t xLoc = (xTr - hit->GetXY(0)) * cosPhi + (yTr - hit->GetXY(1)) * sinPhi;
  600. //Double_t yLoc = -(xTr - hit->GetXY(0)) * sinPhi + (yTr - hit->GetXY(1)) * cosPhi;
  601. TVector2 dist(TMath::Abs(yLoc-hit->GetRphi()), TMath::Abs(xLoc-hit->GetR()));
  602. return dist;
  603. }
  604. //__________________________________________________________________________
  605. Double_t MpdEctTrackFinderTofTpc::Proxim(Double_t phi0, Double_t phi)
  606. {
  607. /// Adjust angle phi to be "around" phi0 - to avoid discontinuity around +- Pi
  608. Double_t dPhi = phi0 - phi;
  609. if (TMath::Abs(dPhi) > TMath::Pi()) phi += TMath::Pi() * 2 * TMath::Sign(1.,dPhi);
  610. return phi;
  611. }
  612. //__________________________________________________________________________
  613. Double_t MpdEctTrackFinderTofTpc::Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit)
  614. {
  615. /// Adjust hit coord. R-Phi to be "around" track R-Phi - to avoid
  616. /// discontinuity around +- Pi
  617. Double_t hitPhi = hit->GetRphi() / hit->GetR();
  618. Double_t phi0 = track->GetParamNew(0) / track->GetPosNew();
  619. return hit->GetR() * Proxim(phi0, hitPhi);
  620. }
  621. //__________________________________________________________________________
  622. void MpdEctTrackFinderTofTpc::Write()
  623. {
  624. /// Write
  625. TFile histoFile("EctRec.root","RECREATE");
  626. Writedir2current(fHistoDir);
  627. histoFile.Close();
  628. }
  629. //__________________________________________________________________________
  630. void MpdEctTrackFinderTofTpc::Writedir2current( TObject *obj )
  631. {
  632. /// Write
  633. if( !obj->IsFolder() ) obj->Write();
  634. else{
  635. TDirectory *cur = gDirectory;
  636. TDirectory *sub = cur->mkdir(obj->GetName());
  637. sub->cd();
  638. TList *listSub = ((TDirectory*)obj)->GetList();
  639. TIter it(listSub);
  640. while( TObject *obj1=it() ) Writedir2current(obj1);
  641. cur->cd();
  642. }
  643. }
  644. //__________________________________________________________________________
  645. void MpdEctTrackFinderTofTpc::RemoveDoubles()
  646. {
  647. /// Remove double tracks (if number of common hits greater than 50% of hits on track)
  648. Int_t nFound = fTracks->GetEntriesFast(), nOK = 0;
  649. for (Int_t i = 0; i < nFound; ++i) {
  650. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  651. if (track == 0x0) continue;
  652. Int_t nh = track->GetNofHits();
  653. cout << i << " " << nh << " " << ++nOK << endl;
  654. for (Int_t j = i + 1; j < nFound; ++j) {
  655. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  656. if (track1 == 0x0) continue;
  657. Int_t nh1 = track1->GetNofHits();
  658. Int_t nc = NofCommonHits(track, track1);
  659. if (1.*nc/TMath::Min(nh,nh1) < 0.5) continue;
  660. if (nh > nh1) fTracks->RemoveAt(j);
  661. else if (nh < nh1) {
  662. fTracks->RemoveAt(i);
  663. --nOK;
  664. break;
  665. } else {
  666. if (track->GetChi2() > track1->GetChi2()) {
  667. fTracks->RemoveAt(i);
  668. --nOK;
  669. break;
  670. }
  671. fTracks->RemoveAt(j);
  672. }
  673. }
  674. }
  675. fTracks->Compress();
  676. // Remove double tracks (originated from the same TPC hit)
  677. //*
  678. nFound = fTracks->GetEntriesFast();
  679. for (Int_t i = 0; i < nFound; ++i) {
  680. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  681. if (track == 0x0) continue;
  682. Int_t iTof = track->GetTpcIndex();
  683. Int_t nh = track->GetNofHits();
  684. for (Int_t j = i + 1; j < nFound; ++j) {
  685. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  686. if (track1 == 0x0) continue;
  687. Int_t iTof1 = track1->GetTpcIndex();
  688. if (iTof1 != iTof) continue;
  689. Int_t nh1 = track1->GetNofHits();
  690. if (nh > nh1) fTracks->RemoveAt(j);
  691. else if (nh < nh1) {
  692. fTracks->RemoveAt(i);
  693. break;
  694. } else {
  695. if (track->GetChi2() > track1->GetChi2()) {
  696. fTracks->RemoveAt(i);
  697. break;
  698. }
  699. fTracks->RemoveAt(j);
  700. }
  701. }
  702. }
  703. fTracks->Compress();
  704. //*/
  705. }
  706. //__________________________________________________________________________
  707. Int_t MpdEctTrackFinderTofTpc::NofCommonHits(MpdEctKalmanTrack* track, MpdEctKalmanTrack* track1)
  708. {
  709. /// Compute number of common hits on 2 tracks
  710. TObjArray *hits = track->GetHits(), *hits1 = track1->GetHits();
  711. MpdKalmanHit *hit = (MpdKalmanHit*) hits->First();
  712. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->First();
  713. Int_t nCom = 0;
  714. while (hit && hit1) {
  715. if (hit->GetLayer() > hit1->GetLayer()) {
  716. hit1 = (MpdKalmanHit*) hits1->After(hit1);
  717. continue;
  718. }
  719. if (hit == hit1) ++nCom;
  720. hit = (MpdKalmanHit*) hits->After(hit);
  721. }
  722. return nCom;
  723. }
  724. //__________________________________________________________________________
  725. void MpdEctTrackFinderTofTpc::AddHits()
  726. {
  727. /// Add hit objects to tracks and compute number of wrongly assigned hits
  728. /// (hits with ID different from ID of starting track)
  729. Int_t nFound = fTracks->GetEntriesFast();
  730. for (Int_t i = 0; i < nFound; ++i) {
  731. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  732. Int_t nHits = track->GetNofHits();
  733. if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  734. track->SetNofHits(nHits);
  735. TClonesArray &trHits = *track->GetTrHits();
  736. SetTrackID(track); // set track ID as ID of majority of its hits
  737. TObjArray *hits = track->GetHits();
  738. Int_t nWrong = 0, motherID = track->GetTrackID();
  739. cout << i << " " << nHits << " " << motherID << endl;
  740. // Get track mother ID
  741. FairMCTrack *mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(motherID);
  742. while (mctrack->GetMotherId() > 0) {
  743. motherID = mctrack->GetMotherId();
  744. mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  745. }
  746. for (Int_t j = 0; j < nHits; ++j) {
  747. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) hits->UncheckedAt(j);
  748. new (trHits[j]) MpdKalmanHitZ(*hit);
  749. Int_t iproj = (hit->GetLayer() - 1) % 3;
  750. if (iproj == 0) cout << " U";
  751. else if (iproj == 1) cout << " R";
  752. else cout << " V";
  753. cout << hit->GetLayer();
  754. MpdStrawendcapPoint *h = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  755. Int_t motherID1 = h->GetTrackID();
  756. cout << "-" << motherID1;
  757. // Get point mother ID
  758. FairMCTrack *mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(motherID1);
  759. while (mctrack->GetMotherId() > 0) {
  760. motherID1 = mctrack->GetMotherId();
  761. mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  762. }
  763. if (motherID1 != motherID) nWrong++;
  764. }
  765. cout << "\n" << nWrong << " " << track->GetTrackID() << " " << motherID << endl;
  766. track->SetNofWrong(nWrong);
  767. track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
  768. }
  769. fTracks->Compress();
  770. }
  771. //__________________________________________________________________________
  772. void MpdEctTrackFinderTofTpc::SetTrackID(MpdEctKalmanTrack* track)
  773. {
  774. /// Set track ID as ID of majority of its hits
  775. const Int_t idMax = 99999;
  776. Int_t ids[idMax+1] = {0}, nHits = track->GetNofHits(), locMax = 0;
  777. TObjArray *hits = track->GetHits();
  778. for (Int_t i = 0; i < nHits; ++i) {
  779. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) hits->UncheckedAt(i);
  780. MpdStrawendcapPoint *p = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  781. Int_t id = p->GetTrackID();
  782. if (id > idMax) exit(0);
  783. ++ids[id];
  784. if (ids[id] > ids[locMax]) locMax = id;
  785. }
  786. if (ids[locMax] > 1) track->SetTrackID(locMax);
  787. }
  788. ClassImp(MpdEctTrackFinderTofTpc);