MpdEctTrackFinderTpc.cxx 50 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353
  1. // -------------------------------------------------------------------------
  2. // ----- MpdEctTrackFinderTpc source file -----
  3. // ----- Created 28/03/08 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdEctTrackFinderTpc.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD End-Cap Tracker (ECT) using seeds from TPC
  9. **/
  10. #include "MpdEctTrackFinderTpc.h"
  11. #include "MpdKalmanFilter.h"
  12. #include "MpdKalmanHit.h"
  13. #include "MpdEctKalmanTrack.h"
  14. #include "MpdStrawendcapGeoPar.h"
  15. #include "MpdStrawendcapPoint.h"
  16. #include "MpdDchPoint.h"
  17. #include "MpdTgemPoint.h"
  18. #include "MpdTpcKalmanTrack.h"
  19. #include "FairGeoNode.h"
  20. #include "MpdMCTrack.h"
  21. #include "FairRootManager.h"
  22. #include "FairRunAna.h"
  23. #include "FairRuntimeDb.h"
  24. #include <TGeoManager.h>
  25. #include <TGeoTube.h>
  26. #include "TMath.h"
  27. #include "TVector2.h"
  28. #include <TRandom.h>
  29. #include <TFile.h>
  30. #include <iostream>
  31. //#include <vector>
  32. using std::cout;
  33. using std::endl;
  34. //using std::set;
  35. const Double_t MpdEctTrackFinderTpc::fgkChi2Cut = 10; //10 //20; //100;
  36. //__________________________________________________________________________
  37. MpdEctTrackFinderTpc::MpdEctTrackFinderTpc(const char *name, Int_t iVerbose )
  38. :FairTask(name, iVerbose)
  39. {
  40. fKHits = new TClonesArray("MpdKalmanHit", 100);
  41. fTracks = new TClonesArray("MpdEctKalmanTrack", 100);
  42. fTrackCand = new TClonesArray("MpdEctKalmanTrack", 100);
  43. fHistoDir = 0x0;
  44. fhLays = new TH1F("hLays0","ECT layers",fgkNlays+1,0,fgkNlays+1);
  45. fLayPointers = 0x0;
  46. fDetType = 1;
  47. fMirror = kTRUE; //kFALSE;
  48. fExact = 0;
  49. if (fExact) fMirror = kFALSE;
  50. }
  51. //__________________________________________________________________________
  52. MpdEctTrackFinderTpc::~MpdEctTrackFinderTpc()
  53. {
  54. //delete fKHits;
  55. //delete fTracks;
  56. //delete fTrackCand;
  57. delete [] fLayPointers;
  58. delete fhLays;
  59. }
  60. //__________________________________________________________________________
  61. InitStatus MpdEctTrackFinderTpc::Init()
  62. {
  63. ReInit();
  64. if (fDetType == 1) InitGeo();
  65. TGeoVolume *inW = gGeoManager->GetVolume("tpc01InWall");
  66. TGeoTube *tube = (TGeoTube*) inW->GetShape();
  67. fZtpc = tube->GetDZ();
  68. return kSUCCESS;
  69. }
  70. //__________________________________________________________________________
  71. InitStatus MpdEctTrackFinderTpc::ReInit()
  72. {
  73. fEctPoints = 0x0; //(TClonesArray *) FairRootManager::Instance()->GetObject("STRAWPoint");
  74. fEctHits =(TClonesArray *) FairRootManager::Instance()->GetObject("STRAWPoint");
  75. if (fEctHits == 0x0) {
  76. fEctHits =(TClonesArray *) FairRootManager::Instance()->GetObject("TgemPoint");
  77. fDetType = 0;
  78. }
  79. if (fEctHits == 0x0) {
  80. fEctHits =(TClonesArray *) FairRootManager::Instance()->GetObject("DchPoint");
  81. fDetType = 2;
  82. }
  83. if (fEctHits == 0x0) Fatal("ReInit", "No ECT points found !!!");
  84. fTpcTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  85. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  86. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  87. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  88. FairRootManager::Instance()->Register("EctTrack", "Ect", fTracks, kTRUE);
  89. FairRootManager::Instance()->Register("EctKalmanHit", "EctHit", fKHits, kFALSE);
  90. fNPass = 2; //2;
  91. return kSUCCESS;
  92. }
  93. //__________________________________________________________________________
  94. void MpdEctTrackFinderTpc::Reset()
  95. {
  96. ///
  97. cout << " MpdEctTrackFinderTpc::Reset " << endl;
  98. //fKHits->Clear("C");
  99. fKHits->Delete();
  100. fTracks->Delete();
  101. fTrackCand->Delete();
  102. if (fLayPointers) delete [] fLayPointers;
  103. fLayPointers = 0x0;
  104. }
  105. //__________________________________________________________________________
  106. void MpdEctTrackFinderTpc::SetParContainers()
  107. {
  108. // Get run and runtime database
  109. FairRunAna* run = FairRunAna::Instance();
  110. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  111. FairRuntimeDb* db = run->GetRuntimeDb();
  112. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  113. // Get ECT geometry parameter container
  114. db->getContainer("MpdStrawendcapGeoPar");
  115. }
  116. //__________________________________________________________________________
  117. void MpdEctTrackFinderTpc::Finish()
  118. {
  119. //Write();
  120. }
  121. //__________________________________________________________________________
  122. void MpdEctTrackFinderTpc::Exec(Option_t * option)
  123. {
  124. static int eventCounter = 0;
  125. cout << " EctRec event " << ++eventCounter << endl;
  126. Reset();
  127. // Create Kalman hits
  128. MakeKalmanHits();
  129. if (fKHits->GetEntriesFast() == 0) return;
  130. for (Int_t i = 0; i < fNPass; ++i) {
  131. fTrackCand->Delete();
  132. GetTrackSeeds(i);
  133. Int_t nHitsOk = 0, nHits = fKHits->GetEntriesFast();
  134. for (Int_t j = 0; j < nHits; ++j) {
  135. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(j);
  136. if (hit->GetFlag() >= 0) ++nHitsOk;
  137. }
  138. cout << " Total number of hits for tracking: " << nHits << ", good: " << nHitsOk << endl;
  139. cout << " Total number of track seeds: " << fTrackCand->GetEntriesFast() << endl;
  140. DoTracking(i);
  141. RemoveDoubles(); // remove double (clone) tracks
  142. StoreTracks(i);
  143. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  144. //if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  145. ExcludeHits(); // exclude used hits
  146. }
  147. SelectTracks(0); // do track selection and compute shared hit multiplicities
  148. AddHits(); // add hit objects to tracks
  149. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  150. }
  151. //__________________________________________________________________________
  152. void MpdEctTrackFinderTpc::InitGeo()
  153. {
  154. // Initialize detector geometry
  155. // Create configuration
  156. cout << " Creating configuration ... " << endl;
  157. map<TString, FairGeoNode*> layers;
  158. map<TString, FairGeoNode*> tubes;
  159. map<TString, FairGeoNode*> modules;
  160. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  161. // Create mother volume
  162. FairGeoNode *cave = new FairGeoNode();
  163. cave->setVolumeType(kFairGeoTopNode);
  164. MpdStrawendcapGeoPar *geoPar = (MpdStrawendcapGeoPar*) rtdb->getContainer("MpdStrawendcapGeoPar");
  165. TObjArray* passNodes = geoPar->GetGeoPassiveNodes();
  166. Int_t nVol = passNodes->GetEntriesFast();
  167. // Get volumes
  168. Int_t nTripl = 0, nMod = 0, nTube = 0;
  169. frMinMax[0] = frMinMax[2] = -1.;
  170. for (Int_t i = 0; i < nVol; ++i) {
  171. FairGeoNode* passVol = (FairGeoNode*) (passNodes->UncheckedAt(i));
  172. TString name = passVol->GetName();
  173. TArrayD* params = passVol->getParameters();
  174. //cout << i << " " << name << " " << (*params)[0] << " " << (*params)[1] << " "
  175. // << (*params)[2] << " " << passVol->GetListOfDaughters()->GetEntriesFast() << endl;
  176. if (name.Contains("layer")) {
  177. layers.insert(pair<TString, FairGeoNode*>(name,passVol));
  178. if (frMinMax[0] < 0) {
  179. frMinMax[0] = (*params)[0];
  180. frMinMax[1] = (*params)[1];
  181. }
  182. Ssiz_t pos = name.First("#") + 1;
  183. nTripl = TMath::Max (nTripl, (TString(name(pos,name.Length()))).Atoi());
  184. } else if (name.Contains("tube")) {
  185. tubes.insert(pair<TString, FairGeoNode*>(name,passVol));
  186. if (frMinMax[2] < 0) frMinMax[2] = (*params)[1] * 10.; // mm
  187. Ssiz_t pos = name.First("#") + 1;
  188. nTube = TMath::Max (nTube, (TString(name(pos,name.Length()))).Atoi());
  189. } else if (name.Contains("mod")) {
  190. modules.insert(pair<TString, FairGeoNode*>(name,passVol));
  191. Ssiz_t pos = name.First("#") + 1;
  192. nMod = TMath::Max (nMod, (TString(name(pos,name.Length()))).Atoi());
  193. }
  194. /*
  195. CbmGeoTransform *transf = passVol->getLabTransform();
  196. CbmGeoVector trans = transf->getTranslation();
  197. //cout << trans << endl;
  198. CbmGeoRotation rot = transf->getRotMatrix();
  199. //rot.print();
  200. passVol->setMother(cave);
  201. passVol->calcLabTransform();
  202. */
  203. }
  204. //fNofLays = nTripl * 3;
  205. //fNofLays *= 2;
  206. fNofTubes = nTube;
  207. // Set mother volumes and find out views ordering
  208. map<TString,FairGeoNode*>::iterator it;
  209. Int_t ishft[3] = {0};
  210. Double_t zzz[3] = {0.};
  211. for (it = layers.begin(); it != layers.end(); ++it) {
  212. FairGeoNode *lay = it->second;
  213. lay->setMother(cave);
  214. TString name = lay->GetName();
  215. if (name == "stt01layerradial#1") zzz[0] = lay->getLabTransform()->getTranslation().getZ();
  216. else if (name == "stt01layerleft#1") zzz[1] = lay->getLabTransform()->getTranslation().getZ();
  217. else if (name == "stt01layerright#1") zzz[2] = lay->getLabTransform()->getTranslation().getZ();
  218. //cout << lay->getLabTransform()->getTranslation() << endl;
  219. //lay->getLabTransform()->getRotMatrix().print();
  220. }
  221. TMath::Sort(3,zzz,ishft,kFALSE);
  222. // Get transforms
  223. for (it = tubes.begin(); it != tubes.end(); ++it) {
  224. FairGeoNode *tube = it->second;
  225. TString layName = tube->GetName();
  226. //if (layName.Contains("stt03") || layName.Contains("stt04")) continue; // skip z<0 for the moment !!!
  227. layName.ReplaceAll("tube","layer");
  228. Ssiz_t pos = layName.First("#") + 1;
  229. Int_t group, itube, ind = 0;
  230. sscanf(&layName[4],"%d",&group);
  231. --group;
  232. itube = (TString(layName(pos,layName.Length()))).Atoi() - 1;
  233. //cout << lay << " " << module << " " << itube << endl;
  234. if (itube >= fgkNtube) {
  235. cout << " Too many tubes !!! " << itube << endl;
  236. exit(0);
  237. }
  238. if (layName.Contains("left")) ind = 1;
  239. else if (layName.Contains("right")) ind = 2;
  240. // Loop over triplets
  241. for (Int_t i = 0; i < nTripl; ++i) {
  242. layName.Remove(pos);
  243. layName += (i+1);
  244. tube->setMother(layers[layName]);
  245. Int_t ilay = i * 3 + group * nTripl * 3 + ishft[ind];
  246. if (ilay >= fgkNlays) {
  247. cout << " Too many layers !!! " << ilay << endl;
  248. exit(0);
  249. }
  250. //cout << tube->GetName() << " " << lay << endl;
  251. FairGeoTransform *tran = tube->calcLabTransform();
  252. //cout << tran->getTranslation() << endl;
  253. //tran->getRotMatrix().print();
  254. fTransf[ilay][itube] = tran;
  255. //transf[ilay][itube] = layers[lay]->getLabTransform();
  256. if (itube == 0) fZplanes[ilay] = layers[layName]->getLabTransform()->getTranslation().getZ();
  257. }
  258. } // for (it = tubes.begin();
  259. }
  260. //__________________________________________________________________________
  261. void MpdEctTrackFinderTpc::MakeKalmanHits()
  262. {
  263. /// Create Kalman hits from ECT points.
  264. if (fDetType == 0) {
  265. // TGEM
  266. MakeKalmanHitsTgem(); // TGEM
  267. return;
  268. } else if (fDetType == 2) {
  269. // DCH
  270. MakeKalmanHitsDch(); // DCH
  271. return;
  272. }
  273. Int_t nHits = fEctHits->GetEntriesFast(), lay = 0, nKH = 0, itube = 0;
  274. Double_t phi, coord, errR = 0.2, errRphi = 0.02; //r, 2mm in R, 200um in R-Phi
  275. /*Double_t firstHit[100][1000] = {{0},{0}};
  276. Int_t indxTube[100][1000];
  277. for (Int_t i = 0; i < 100; ++i) {
  278. for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
  279. }*/
  280. Double_t firstHit[fgkNlays+1][fgkNtube] = {{0},{0}};
  281. Int_t indxTube[fgkNlays+1][fgkNtube];
  282. Int_t iend = fgkNlays + 1;
  283. for (Int_t i = 0; i < iend; ++i) {
  284. for (Int_t j = 0; j < fgkNtube; ++j) indxTube[i][j] = -1;
  285. }
  286. for (Int_t ih = 0; ih < nHits; ++ih) {
  287. MpdStrawendcapPoint *h = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(ih);
  288. //if (h->GetZ() < 0) continue; // !!! exclude one side for now
  289. //if (h->GetTrackID() != 1063) continue; // !!!
  290. phi = h->GetPhi(); // tube Phi
  291. lay = h->GetDetectorID() / 1000; // + 1;
  292. if (h->GetZ() < 0) lay += fgkNlays2;;
  293. itube = h->GetDetectorID() % 1000; // tube number
  294. // Extrapolate track to Z = Ztube
  295. Double_t dZ = h->GetZ() - h->GetTrackZ();
  296. Double_t dt = 0.; // dZ / h->GetPz();
  297. if (TMath::Abs(h->GetPz()) > 1.e-6 && h->GetPz() * h->GetZ() > 0) dt = dZ / h->GetPz();
  298. Double_t xNew = h->GetTrackX() + dt * h->GetPx();
  299. Double_t yNew = h->GetTrackY() + dt * h->GetPy();
  300. //cout << dZ << " " << h->GetTrackX() << " " << xNew << " " << h->GetTrackY() << " " << yNew << " " << lay << endl;
  301. //Double_t zNew = h->GetTrackZ() + dt * h->GetPz(); // just for cross-check
  302. // Transform to the rotated coordinate system
  303. Double_t cosSin[2] = {TMath::Cos(phi), TMath::Sin(phi)}; // phi - tube direction
  304. //Double_t xLoc = h->GetX() * cosPhi + h->GetY() * sinPhi; // cross-check
  305. //Double_t yLoc = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  306. //Double_t xRot = xNew * cosSin[0] + yNew * cosSin[1];
  307. Double_t yRot = -xNew * cosSin[1] + yNew * cosSin[0];
  308. //Double_t xLoc = (xNew - h->GetX()) * cosPhi + (yNew - h->GetY()) * sinPhi;
  309. //Double_t yLoc = -(xNew - h->GetX()) * sinPhi + (yNew - h->GetY()) * cosPhi;
  310. //r = xNew * xNew + yNew * yNew;
  311. //r = TMath::Sqrt (r);
  312. //r = TMath::Abs(xLoc);
  313. //r = xRot;
  314. //cout << xRot << " " << yRot << " " << r << " " << h->GetPz() << endl;
  315. coord = yRot;
  316. Double_t yWire = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  317. Double_t dY = TMath::Abs (yWire - coord);
  318. // Add error
  319. Double_t dRphi = 0, dR = 0;
  320. gRandom->Rannor(dRphi,dR);
  321. Double_t meas[2] = {coord+dRphi*errRphi, 0};
  322. Double_t err[2] = {errRphi, errR};
  323. MpdKalmanHit *hit = new ((*fKHits)[nKH]) MpdKalmanHit(lay*1000000+itube, 1, MpdKalmanHit::kFixedZ,
  324. meas, err, cosSin, 9., h->GetZ(), ih);
  325. //hit->SetXY(h->GetX(), h->GetY());
  326. // Add second measurement - just for test at the moment
  327. //!!!
  328. //hit->SetNofDim(2);
  329. //!!!
  330. // Check if multiple hits per tube
  331. if (1 && indxTube[lay][itube] >= 0) {
  332. // Multiple hits
  333. MpdKalmanHit *hitOld = (MpdKalmanHit*) fKHits->UncheckedAt(indxTube[lay][itube]); // previous hit
  334. if (dY < firstHit[lay][itube]) {
  335. // Closer to anode wire - remove previous hit
  336. firstHit[lay][itube] = dY;
  337. const TArrayI *inds = hitOld->Index();
  338. Int_t nOld = inds->GetSize();
  339. for (Int_t io = 0; io < nOld; ++io) hit->SetIndex((*inds)[io]);
  340. fKHits->RemoveAt(indxTube[lay][itube]);
  341. if (fMirror) fKHits->RemoveAt(indxTube[lay][itube]+1);
  342. indxTube[lay][itube] = nKH;
  343. } else {
  344. hitOld->SetIndex(hit->GetIndex());
  345. fKHits->RemoveAt(nKH); // remove last hit
  346. }
  347. } else {
  348. firstHit[lay][itube] = dY;
  349. indxTube[lay][itube] = nKH;
  350. }
  351. if (fMirror && fKHits->UncheckedAt(nKH)) {
  352. // Add mirror hit
  353. // wire point transverse coordinate in rotated system:
  354. //yRot = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  355. //Double_t dy = hit->GetRphi() - yRot;
  356. Double_t dy = hit->GetMeas(0) - yWire;
  357. MpdKalmanHit *hitM = new ((*fKHits)[++nKH]) MpdKalmanHit(*hit);
  358. //hitM->SetRphi(yRot-dy);
  359. hitM->SetMeas(0, yWire-dy);
  360. hitM->SetMirror();
  361. }
  362. ++nKH;
  363. }
  364. if (nKH == 0) return;
  365. fKHits->Compress();
  366. fKHits->Sort(); // in descending order in abs(Z)
  367. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  368. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << endl;
  369. fLayPointers = new Int_t [layMax+1];
  370. fhLays->Reset();
  371. /*
  372. Int_t ipos = 0;
  373. for (Int_t i = 0; i <= layMax; ++i) {
  374. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  375. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  376. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  377. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  378. if (cont == 0) {
  379. fLayPointers[i] = -1;
  380. continue;
  381. }
  382. fLayPointers[i] = ipos;
  383. ipos += cont;
  384. }
  385. */
  386. for (Int_t i = 0; i <= layMax; ++i) fLayPointers[i] = -1;
  387. nKH = fKHits->GetEntriesFast();
  388. for (Int_t i = 0; i < nKH; ++i) {
  389. lay = ((MpdKalmanHit*) fKHits->UncheckedAt(i))->GetLayer();
  390. fhLays->Fill(lay+0.1);
  391. if (fLayPointers[lay] < 0) fLayPointers[lay] = i;
  392. }
  393. }
  394. //__________________________________________________________________________
  395. void MpdEctTrackFinderTpc::MakeKalmanHitsTgem()
  396. {
  397. /// Create Kalman hits from TGEM points.
  398. fhLays->Reset();
  399. Int_t nHits = fEctHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0, itube = 0;
  400. //Double_t phi, r, coord, errR = 0.2, errRphi = 0.02; // 2mm in R, 200um in R-Phi
  401. //Double_t phi, r, coord, errR = 0.01, errRphi = 0.01; // 100um in R, 100um in R-Phi
  402. static const Double_t phi[3] = {0, 60*TMath::DegToRad(), -60*TMath::DegToRad()}, errR = 0.01, errRphi = 0.15 / TMath::Sqrt(12); // 100um in R, 1.5mm pitch
  403. Double_t firstHit[100][1000] = {{0},{0}};
  404. Int_t indxTube[100][1000];
  405. for (Int_t i = 0; i < 100; ++i) {
  406. for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
  407. }
  408. for (Int_t ih = 0; ih < nHits; ++ih) {
  409. MpdTgemPoint *h = (MpdTgemPoint*) fEctHits->UncheckedAt(ih);
  410. if (h->GetZ() < 0) continue; // !!! exclude one side for now
  411. lay = h->GetDetectorID() / 10000;
  412. Int_t lay3 = (lay-1) % 3;
  413. Double_t cosSin[2] = {TMath::Cos(phi[lay3]), TMath::Sin(phi[lay3])}; // phi - strip direction
  414. Double_t u = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  415. // Add error
  416. Double_t dRphi = 0, dR = 0;
  417. gRandom->Rannor(dRphi,dR);
  418. Double_t meas[2] = {u+dRphi*errRphi, 0};
  419. Double_t err[2] = {errRphi, errR};
  420. //MpdKalmanHit *hit =
  421. new ((*fKHits)[nKH]) MpdKalmanHit(h->GetDetectorID()*100, 1,
  422. MpdKalmanHit::kFixedZ, meas, err, cosSin, 9., h->GetZ(), ih);
  423. // Add second measurement - just for test at the moment
  424. //!!!
  425. //hit->SetNofDim(2);
  426. //!!!
  427. //lay = hit.GetLayer();
  428. layMax = TMath::Max (lay, layMax);
  429. //fhLays->Fill(lay+0.1);
  430. // Check if multiple hits per tube
  431. if (0 && indxTube[lay][itube] >= 0) {
  432. // Multiple hits
  433. if (TMath::Abs(u) < firstHit[lay][itube]) {
  434. // Closer to anode wire - remove previous hit
  435. firstHit[lay][itube] = TMath::Abs(u);
  436. fKHits->RemoveAt(indxTube[lay][itube]);
  437. indxTube[lay][itube] = nKH;
  438. } else fKHits->RemoveAt(nKH); // remove last hit
  439. } else {
  440. firstHit[lay][itube] = TMath::Abs(u);
  441. indxTube[lay][itube] = nKH;
  442. fhLays->Fill(lay+0.1);
  443. }
  444. ++nKH;
  445. }
  446. fKHits->Compress();
  447. fKHits->Sort(); // in descending order in abs(Z)
  448. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << " "
  449. << ((MpdKalmanHit*)fKHits->First())->GetDist() << endl;
  450. fLayPointers = new Int_t [layMax+1];
  451. /*
  452. Int_t ipos = 0;
  453. for (Int_t i = 0; i <= layMax; ++i) {
  454. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  455. if (cont == 0) {
  456. fLayPointers[i] = -1;
  457. continue;
  458. }
  459. fLayPointers[i] = ipos;
  460. ipos += cont;
  461. }
  462. */
  463. for (Int_t i = 0; i <= layMax; ++i) fLayPointers[i] = -1;
  464. nKH = fKHits->GetEntriesFast();
  465. for (Int_t i = 0; i < nKH; ++i) {
  466. lay = ((MpdKalmanHit*) fKHits->UncheckedAt(i))->GetLayer();
  467. if (fLayPointers[lay] < 0) fLayPointers[lay] = i;
  468. }
  469. }
  470. //__________________________________________________________________________
  471. void MpdEctTrackFinderTpc::MakeKalmanHitsDch()
  472. {
  473. /// Create Kalman hits from DCH points.
  474. fhLays->Reset();
  475. Int_t nHits = fEctHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0, itube = 0;
  476. if (nHits == 0) return;
  477. static const Double_t errR = 0.01, errRphi = 0.01; // 100um in R, 100um in R-Phi
  478. //static const Double_t errR = 0.01, errRphi = 0.02; // 100um in R, 200um in R-Phi
  479. Double_t firstHit[100][1000] = {{0},{0}};
  480. Int_t indxTube[100][1000];
  481. for (Int_t i = 0; i < 100; ++i) {
  482. for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
  483. }
  484. for (Int_t ih = 0; ih < nHits; ++ih) {
  485. MpdDchPoint *h = (MpdDchPoint*) fEctHits->UncheckedAt(ih);
  486. if (h->GetZ() < 0) continue; // !!! exclude one side for now
  487. lay = h->GetDetectorID();
  488. Double_t cosSin[2] = {TMath::Cos(h->GetPhi()), TMath::Sin(h->GetPhi())}; // phi - wire direction
  489. Double_t u = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  490. // Add error
  491. Double_t dRphi = 0, dR = 0;
  492. gRandom->Rannor(dRphi,dR);
  493. Double_t meas[2] = {u+dRphi*errRphi, 0};
  494. Double_t err[2] = {errRphi, errR};
  495. //MpdKalmanHit *hit =
  496. new ((*fKHits)[nKH]) MpdKalmanHit(h->GetDetectorID()*1000000, 1,
  497. MpdKalmanHit::kFixedZ, meas, err, cosSin, 9., h->GetZ(), ih);
  498. // Add second measurement - just for test at the moment
  499. //!!!
  500. //hit->SetNofDim(2);
  501. //!!!
  502. //lay = hit.GetLayer();
  503. layMax = TMath::Max (lay, layMax);
  504. //fhLays->Fill(lay+0.1);
  505. // Check if multiple hits per tube
  506. if (0 && indxTube[lay][itube] >= 0) {
  507. // Multiple hits
  508. if (TMath::Abs(u) < firstHit[lay][itube]) {
  509. // Closer to anode wire - remove previous hit
  510. firstHit[lay][itube] = TMath::Abs(u);
  511. fKHits->RemoveAt(indxTube[lay][itube]);
  512. indxTube[lay][itube] = nKH;
  513. } else fKHits->RemoveAt(nKH); // remove last hit
  514. } else {
  515. firstHit[lay][itube] = TMath::Abs(u);
  516. indxTube[lay][itube] = nKH;
  517. fhLays->Fill(lay+0.1);
  518. }
  519. ++nKH;
  520. }
  521. fKHits->Compress();
  522. fKHits->Sort(); // in descending order in abs(Z)
  523. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << " "
  524. << ((MpdKalmanHit*)fKHits->First())->GetDist() << endl;
  525. fLayPointers = new Int_t [layMax+1];
  526. for (Int_t i = 0; i <= layMax; ++i) fLayPointers[i] = -1;
  527. nKH = fKHits->GetEntriesFast();
  528. for (Int_t i = 0; i < nKH; ++i) {
  529. lay = ((MpdKalmanHit*) fKHits->UncheckedAt(i))->GetLayer();
  530. if (fLayPointers[lay] < 0) fLayPointers[lay] = i;
  531. }
  532. }
  533. //__________________________________________________________________________
  534. void MpdEctTrackFinderTpc::GetTrackSeeds(Int_t iPass)
  535. {
  536. /// Build ECT track seeds from TPC tracks (apply some angular selection cuts?)
  537. Int_t nTpcTracks = fTpcTracks->GetEntriesFast();
  538. cout << " Seed tracks: " << nTpcTracks << endl;
  539. Int_t nCand = 0;
  540. for (Int_t itr = 0; itr < nTpcTracks; ++itr) {
  541. MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(itr);
  542. if (tpc == 0x0 || tpc->GetUniqueID() == 1) continue;
  543. if (TMath::Abs((*tpc->GetParamAtHit())(3,0)) < 0.3) continue; // track in central region
  544. MpdKalmanTrack tmp(*tpc);
  545. tmp.SetParam(*tpc->GetParamAtHit());
  546. if (tmp.Momentum() < 0.07) continue; // !!! track with low momentum
  547. //if (tpc->GetParam(3) < 0) continue; // !!! going in backward hemisphere
  548. MpdEctKalmanTrack *track = new ((*fTrackCand)[nCand++]) MpdEctKalmanTrack(itr, *tpc);
  549. //cout << nCand-1 << " " << track->GetTrackID() << endl;
  550. //cout << track->GetHits()->GetEntriesFast() << endl;
  551. //track->GetParam()->Print();
  552. //track->GetCovariance()->Print();
  553. //track->GetWeight()->Print();
  554. /*AZ 27-12-2011
  555. Double_t distFirst = ((MpdKalmanHit*) track->GetHits()->First())->GetDist();
  556. Double_t distLast = ((MpdKalmanHit*) track->GetHits()->Last())->GetDist();
  557. track->SetPos(TMath::Min(distFirst,distLast));
  558. track->SetPosNew(track->GetPos());
  559. Bool_t ok = MpdKalmanFilter::Instance()->Refit(track, -1); // from last to first hit
  560. //if (!ok || track->GetParamNew(3) < 0) {
  561. if (!ok) {
  562. // propagation failure or going in backward hemisphere
  563. fTrackCand->RemoveLast();
  564. --nCand;
  565. continue;
  566. }
  567. */
  568. //AZ 27-12-2011
  569. /*track->SetPos(((MpdKalmanHit*) track->GetHits()->First())->GetDist());
  570. track->SetPosNew(track->GetPos());
  571. track->SetLength(track->GetLengAtHit());*/
  572. //
  573. track->GetHits()->Clear();
  574. track->SetChi2(0.);
  575. }
  576. cout << " Number of ECT track candidates: " << nCand << endl;
  577. }
  578. //__________________________________________________________________________
  579. void MpdEctTrackFinderTpc::DoTracking(Int_t iPass)
  580. {
  581. /// Run Kalman tracking
  582. Int_t nCand = fTrackCand->GetEntriesFast(), iok = 0;
  583. for (Int_t i = 0; i < nCand; ++i) {
  584. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  585. //Int_t id = track->GetTrackID();
  586. //cout << " Track seed No. " << i << ", ID: " << track->GetTrackID() << endl;
  587. //if (track->GetTrackID() != 77) continue;
  588. //(*(track->GetParamNew()))(4,0) = -0.5; // just for check
  589. /* "Standalone" tracking
  590. for (Int_t k = 0; k < 5; ++k) {
  591. for (Int_t j = k; j < 5; ++j) {
  592. if (j == k) (*track->GetWeight())(k,j) /= 100.;
  593. else (*track->GetWeight())(k,j) = (*track->GetWeight())(j,k) = 0.;
  594. }
  595. }
  596. */
  597. // Propagate to TPC end-plate
  598. MpdKalmanHit hit;
  599. hit.SetType(MpdKalmanHit::kFixedZ);
  600. //Double_t zEnd = 150.; // get it from geometry !!!
  601. //hit.SetPos(TMath::Sign(zEnd,track->GetParam(3)));
  602. hit.SetPos(TMath::Sign(fZtpc,track->GetParam(3)));
  603. //Double_t thick = 0.06; // material thickness in rad. lengths
  604. Double_t thick = 0.20; // material thickness in rad. lengths - v7
  605. MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  606. PassWall(track,thick);
  607. // Start ECT tracking from different layers to account for missing hits
  608. MpdEctKalmanTrack tracks[10];
  609. Int_t layMax = 1; // TGEM or DCH
  610. //Int_t layMax = 2; // TGEM or DCH
  611. if (fDetType == 1) layMax = 10; // ECT
  612. for (Int_t j = 0; j < layMax; ++j) {
  613. tracks[j] = *track;
  614. //cout << track->GetParamNew(0) << endl;
  615. //cout << i << " " << lay << " " << tracks[lay].GetNofHits() << " " << tracks[lay].GetChi2() << " " << tracks[lay].GetParam(0) << endl;
  616. //Int_t lay = 1;
  617. Int_t lay = j + 1;
  618. if (hit.GetPos() < 0) lay += fgkNlays2;
  619. if (j > 0 && tracks[j-1].GetNofHits() > 0)
  620. lay = ((MpdKalmanHit*)tracks[j-1].GetHits()->First())->GetLayer() + 1;
  621. iok = RunKalmanFilter(&tracks[j], lay);
  622. //iok = RunKalmanFilter(&tracks[lay], 0);
  623. //cout << i << " " << lay << " " << tracks[j].GetNofHits() << " " << tracks[j].GetChi2() << endl;
  624. }
  625. // Select the best track (with max number of hits)
  626. Int_t nHitsMax = tracks[0].GetNofHits(), iMax = 0;
  627. for (Int_t j = 1; j < layMax; ++j) {
  628. Int_t nhits = tracks[j].GetNofHits();
  629. if (nhits > nHitsMax) {
  630. nHitsMax = nhits;
  631. iMax = j;
  632. } else if (nhits == nHitsMax) {
  633. if (tracks[j].GetChi2() < tracks[iMax].GetChi2()) {
  634. iMax = j;
  635. nHitsMax = nhits;
  636. }
  637. }
  638. }
  639. //*track = tracks[iMax];
  640. fTrackCand->RemoveAt(i);
  641. new ((*fTrackCand)[i]) MpdEctKalmanTrack(tracks[iMax]);
  642. if (0 && iok == 0 && track->GetNofHits() < 0) {
  643. //*
  644. MpdKalmanHit *pHit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(0);
  645. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  646. track->SetPos(pHit->GetPos()-dir);
  647. track->SetPosNew(track->GetPos());
  648. //for (Int_t j = 0; j < 5; ++j) track->SetParam(j,track->GetParam1()[j]);
  649. //track->SetParamNew(*track->GetParam());
  650. //MpdKalmanHitZ hitTmp(*pHit);
  651. //hitTmp.SetZ(track->GetPosNew()-dir);
  652. //MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  653. MpdKalmanFilter::Instance()->Refit(track,1);
  654. //*/
  655. /*
  656. track->Print("");
  657. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  658. MpdKalmanHitZ *hitZ = (MpdKalmanHitZ*) track->GetHits()->UncheckedAt(0);
  659. MpdKalmanHitZ hitTmp(*hitZ);
  660. hitTmp.SetZ(track->GetPosNew()+dir);
  661. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  662. track->Print("");
  663. MpdKalmanFilter::Instance()->Refit(track,-1);
  664. track->Print("");
  665. hitTmp.SetZ(track->GetPosNew()-dir);
  666. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  667. track->Print("");
  668. MpdKalmanFilter::Instance()->Refit(track,1);
  669. track->Print("");
  670. */
  671. }
  672. } // for (Int_t i = 0; i < nCand;
  673. }
  674. //__________________________________________________________________________
  675. void MpdEctTrackFinderTpc::PassWall(MpdEctKalmanTrack *track, Double_t thick)
  676. {
  677. /// Propagate track thru TPC end plate (include MS)
  678. // Add multiple scattering
  679. TMatrixDSym *cov = track->Weight2Cov();
  680. Double_t th = track->GetParamNew(3);
  681. Double_t cosTh = TMath::Cos(th);
  682. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick);
  683. (*cov)(2,2) += (angle2 / cosTh / cosTh );
  684. (*cov)(3,3) += angle2;
  685. Int_t iok = 0;
  686. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  687. track->SetWeight(*cov);
  688. }
  689. //__________________________________________________________________________
  690. Int_t MpdEctTrackFinderTpc::RunKalmanFilter(MpdEctKalmanTrack *track, Int_t layBeg)
  691. {
  692. /// Run Kalman filter
  693. Double_t rMin = 49.2, rMax = 110.5; // min and max radial size of TGEM - to be taken elsewhere
  694. Double_t dX = 0.009; // 0.9% of X0 (TGEM layer thickness)
  695. if (fDetType == 1) {
  696. // ECT
  697. rMin = 49.2;
  698. rMax = 110.5;
  699. dX = 0.001;
  700. } else if (fDetType == 2) {
  701. // DCH
  702. rMin = 9.2;
  703. rMax = 135.4;
  704. dX = 0.001;
  705. }
  706. if (track->GetRadNew() > rMax) return -1;
  707. //cout << fHits->GetEntriesFast() << endl;
  708. //Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  709. Int_t layMax = fgkNlays2;
  710. MpdKalmanHit *hitOK = 0x0;
  711. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  712. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1;
  713. Int_t layEnd = -1, dLay = -1;//, layOK = -1;
  714. //if (track->GetPosNew() < 0) layEnd = fgkNlays2;
  715. if (track->GetPosNew() < 0) layMax = 2*fgkNlays2;
  716. if (trackDir == MpdKalmanTrack::kOutward) {
  717. layEnd = layMax + 1;
  718. dLay = 1;
  719. //if (track->GetPosNew() > 0 && layMax > fgkNlays2) layEnd -= fgkNlays2;
  720. }
  721. //Int_t indxOK = hits->IndexOf(hitOK);
  722. //Int_t nHits = hits->GetEntriesFast();
  723. Int_t miss = 0;
  724. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  725. TMatrixD param(5,1), paramTmp(5,1);
  726. Double_t saveZ = 0., saveLeng = 0.;
  727. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  728. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  729. //if (lay > 30) rMin = 38.;
  730. //else rMin = 33.;
  731. Int_t nLay = GetNofHitsInLayer(lay);
  732. Int_t indx0 = GetHitsInLayer(lay);
  733. Double_t dChi2Min = 1.e+6;
  734. MpdKalmanHit *hitMin = 0x0;
  735. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  736. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  737. /*
  738. if (trackDir == MpdKalmanTrack::kOutward) {
  739. // !!! This part is due to the simplified hit merging (digitization) -
  740. // different hit position inside one layer - should be removed later
  741. indxBeg = nLay - 1;
  742. indxEnd = -1;
  743. dIndx = -1;
  744. }
  745. */
  746. //for (Int_t indx = 0; indx < nLay; ++indx) {
  747. Double_t radNew = -1;
  748. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  749. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  750. // Exclude used hits
  751. if (hit->GetFlag() < 0) continue;
  752. // !!! Exact ID match
  753. Int_t itmp = 0;
  754. //if (fExact && ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  755. if (fExact && HitMotherId(hit,track->GetTrackID(),itmp) != track->GetTrackID()) continue;
  756. if (TMath::Abs(hit->GetPos()-track->GetPosNew()) > 1.e-3) {
  757. //cout << " window:" << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << hit->GetZ() << " " << track->GetParamNew(1) << endl;
  758. // Check if the hit within some window (15x25cm for the moment - check!!!)
  759. //if (TMath::Abs(hit->GetRphi()-track->GetParamNew(0)) > 9) continue;
  760. //if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  761. TVector2 dist = GetDistance(track, hit);
  762. //if (track->GetTrackID() == 1335) cout << ((TpcLheKalmanTrack*)fTpcTracks->UncheckedAt(track->GetTpcIndex()))->GetNofTrHits() << " " << dist.X() << " " << dist.Y() << " " << hit->GetLayer() << endl;
  763. if (dist.X() > 15.) continue; // distance in transverse to the tube direction
  764. if (hit->GetNofDim() > 1 && dist.Y() > 25.) continue; // distance in R
  765. //*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;
  766. //track->Print("");
  767. //hit->Print("");
  768. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) return -1;
  769. dist = GetDistance(track, hit);
  770. //if (track->GetTrackID() == 1335) cout << track->GetTrackID() << " " << dist.X() << " " << dist.Y() << " " << endl;
  771. //*if (hit->GetTrackID() == -607)*/ cout << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << " " << hit->GetZ() << " " << track->GetPosNew() << endl;
  772. //
  773. // For debugging
  774. /*
  775. TVector2 dist0 = GetDistance(track, hit);
  776. cout << dist0.X() << " ";
  777. MpdKalmanHitZ hitDbg = *hit;
  778. Double_t xDbg = hit->GetXY(0) * TMath::Cos(hit->GetPhi()) + hit->GetXY(1) * TMath::Sin(hit->GetPhi());
  779. Double_t yDbg = -hit->GetXY(0) * TMath::Sin(hit->GetPhi()) + hit->GetXY(1) * TMath::Cos(hit->GetPhi());
  780. hitDbg.SetRphi(yDbg);
  781. hitDbg.SetR(xDbg);
  782. dist = GetDistance(track, &hitDbg);
  783. cout << dist.X() << endl;
  784. */
  785. //
  786. //Double_t radNew = track->GetRadNew();
  787. if (radNew < 0) radNew = track->GetRadNew();
  788. //if (radNew < rMin || radNew > rMax) return -1;
  789. if (radNew < rMin) break; // next layer
  790. else if (radNew > rMax) return -1;
  791. //Double_t err = hit->GetRphiErr();
  792. //if (track->GetNofHits() == 0) hit->SetRphiErr(0.04);
  793. // Add multiple scattering
  794. if (TMath::Abs(hit->GetPos()-track->GetPos()) > 1.e-1) PassWall(track, dX);
  795. }
  796. if (fDetType == 1) {
  797. // ECT
  798. Double_t phi = TMath::ATan2 (track->GetParamNew(1),track->GetParamNew(0));
  799. Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi,hit->GetPhi()) - phi;
  800. if (TMath::Abs(dPhi) > 1) continue;
  801. }
  802. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(track,hit,pointWeight,param);
  803. //if (track->GetTrackID() == 1335) cout << dChi2 << endl;
  804. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  805. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  806. if (dChi2 < dChi2Min) {
  807. //if (dChi2 < dChi2Min && dist0.X() <= dist.X()) {
  808. //if (dChi2 < dChi2Min && dist.X() <= 0.2) {
  809. dChi2Min = dChi2;
  810. hitMin = hit;
  811. saveWeight = *track->GetWeight();
  812. saveZ = track->GetPosNew();
  813. saveLeng = track->GetLength();
  814. // temporary storage for the current track
  815. paramTmp = param;
  816. pointWeightTmp = pointWeight;
  817. //cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetR() << endl;
  818. //cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  819. //cout << hit->GetRphi() << " " << hit->GetZ() << endl;
  820. //cout << param(0,0) << " " << param(1,0) << endl;
  821. //paramTmp.Print();
  822. }
  823. } //for (Int_t indx = indxBeg; indx != indxEnd;
  824. Double_t cut = fgkChi2Cut;
  825. //if (track->GetNofHits() == 0) cut /= 2.;
  826. //if (dChi2Min < fgkChi2Cut) {
  827. if (dChi2Min < cut) {
  828. //layOK = lay;
  829. hitOK = hitMin;
  830. track->GetHits()->Add(hitOK);
  831. //miss = 0;
  832. // Restore track params at the best hit
  833. track->SetChi2(track->GetChi2()+dChi2Min);
  834. saveWeight += pointWeightTmp;
  835. track->SetWeight(saveWeight);
  836. track->SetPosNew(saveZ);
  837. track->SetLength(saveLeng);
  838. track->SetParamNew(paramTmp);
  839. //if (track->GetTrackID() == 1335) cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  840. //paramTmp.Print();
  841. // Check if the accepted hit is the same as the seeded hit
  842. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  843. if (track->GetNofHits() == 1) track->SetParam1();
  844. //if (((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() == track->GetTrackID()) track->fLengAtHit = saveLeng;
  845. } else if (fDetType == 1 && radNew > rMin) {
  846. // ECT
  847. //++miss;
  848. //if (miss > 1) return -1;
  849. // Check whether tracks goes through a tube
  850. Int_t layM1 = lay - 1;
  851. if (nLay == 0) {
  852. // No hits in the layer - extrapolate track
  853. MpdKalmanHit hitTmp;
  854. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  855. hitTmp.SetDist(fZplanes[layM1]);
  856. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp)) return -1;
  857. }
  858. FairGeoVector vec(track->GetParamNew(0)*10,track->GetParamNew(1)*10,track->GetPosNew()*10); // mm
  859. //cout << vec.getX() << " " << vec.getY() << " " << vec.getZ() << endl;
  860. Int_t iTube1 = 0, iter = 0, idTubes = 0;//, iTubeOk = 0;
  861. Int_t iTube = Int_t (TMath::ATan2 (vec.getY(),vec.getX()) / TMath::TwoPi() * fNofTubes);
  862. if (iTube < 0) iTube += fNofTubes;
  863. iTube = iTube % fNofTubes;
  864. FairGeoVector vecLoc = fTransf[layM1][iTube]->transTo(vec);
  865. //cout << vecLoc.getX() << " " << vecLoc.getY() << " " << vecLoc.getZ() << " " << iTube << endl;
  866. //cout << transf[lay][0]->getTranslation() << endl;
  867. //transf[lay][0]->getRotMatrix().print();
  868. /*if (TMath::Abs(vecLoc.getZ()) > z2) cout << vecLoc << " " << ph << endl;
  869. if (TMath::Abs(vecLoc.getZ()) > z2) {
  870. iTube = nTubes / 2;
  871. vecLoc = transf[lay][iTube]->transTo(vec);
  872. cout << vecLoc;
  873. }*/
  874. Double_t dyMin = 0;
  875. while (1) {
  876. FairGeoVector vecLoc1 = fTransf[layM1][(iTube+1)%fNofTubes]->transTo(vec);
  877. Double_t dy = vecLoc.getY() - vecLoc1.getY();
  878. idTubes = TMath::Nint (vecLoc.getY() / dy);
  879. iTube1 = iTube + idTubes;
  880. if (iTube1 < 0) iTube1 += fNofTubes;
  881. iTube1 = iTube1 % fNofTubes;
  882. vecLoc1 = fTransf[layM1][iTube1]->transTo(vec);
  883. dyMin = TMath::Abs(vecLoc1.getY());
  884. //iTubeOk = iTube1;
  885. //cout << vecLoc1.getX() << " " << vecLoc1.getY() << " " << vecLoc1.getZ() << " " << iTubeOk << " " << dyMin << " " << iter << endl;
  886. if (dyMin <= frMinMax[2] || (iter && TMath::Abs(idTubes) < 2 && TMath::Abs(dyMin/dy) < 0.51)) break;
  887. iTube = iTube1;
  888. vecLoc = vecLoc1;
  889. ++iter;
  890. }
  891. //if (dyMin > 4.) ++miss; // too far from the tube
  892. //if (dyMin < frMinMax[2]+0.5) ++miss; // too far from the tube
  893. if (dyMin < frMinMax[2]-0.2) ++miss; // too far from the tube
  894. //cout << dyMin << " " << dChi2Min << endl;
  895. }
  896. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  897. } // for (Int_t lay = layOK-1; lay >= 0;
  898. track->SetMisses(miss);
  899. return 0;
  900. }
  901. //__________________________________________________________________________
  902. TVector2 MpdEctTrackFinderTpc::GetDistance(MpdEctKalmanTrack *track, MpdKalmanHit *hit)
  903. {
  904. /// Compute distance between track and hit
  905. Double_t xTr, yTr;
  906. if (track->GetType() == MpdKalmanTrack::kEndcap) {
  907. xTr = track->GetParamNew(0);
  908. yTr = track->GetParamNew(1);
  909. } else {
  910. Double_t rPhi = track->GetParamNew(0);
  911. Double_t r = track->GetPosNew();
  912. Double_t ph = rPhi / r;
  913. xTr = r * TMath::Cos(ph);
  914. yTr = r * TMath::Sin(ph);
  915. }
  916. Double_t phi = hit->GetPhi();
  917. Double_t cosPhi = TMath::Cos(phi);
  918. Double_t sinPhi = TMath::Sin(phi);
  919. // Transform track coordinates to local tube coordinates
  920. Double_t xLoc = xTr * cosPhi + yTr * sinPhi;
  921. Double_t yLoc = -xTr * sinPhi + yTr * cosPhi;
  922. //Double_t xLoc = (xTr - hit->GetXY(0)) * cosPhi + (yTr - hit->GetXY(1)) * sinPhi;
  923. //Double_t yLoc = -(xTr - hit->GetXY(0)) * sinPhi + (yTr - hit->GetXY(1)) * cosPhi;
  924. TVector2 dist(TMath::Abs(yLoc-hit->GetMeas(0)), TMath::Abs(xLoc-hit->GetMeas(1)));
  925. //cout << dist.X() << " " << dist.Y() << endl;
  926. return dist;
  927. }
  928. //__________________________________________________________________________
  929. Double_t MpdEctTrackFinderTpc::Proxim(Double_t phi0, Double_t phi)
  930. {
  931. /// Adjust angle phi to be "around" phi0 - to avoid discontinuity around +- Pi
  932. Double_t dPhi = phi0 - phi;
  933. if (TMath::Abs(dPhi) > TMath::Pi()) phi += TMath::Pi() * 2 * TMath::Sign(1.,dPhi);
  934. return phi;
  935. }
  936. //__________________________________________________________________________
  937. Double_t MpdEctTrackFinderTpc::Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit)
  938. {
  939. /// Adjust hit coord. R-Phi to be "around" track R-Phi - to avoid
  940. /// discontinuity around +- Pi
  941. Fatal("Proxim", " !!! Not implemented !!!");
  942. /*
  943. Double_t hitPhi = hit->GetRphi() / hit->GetR();
  944. Double_t phi0 = track->GetParamNew(0) / track->GetPosNew();
  945. return hit->GetR() * Proxim(phi0, hitPhi);
  946. */
  947. return 0.;
  948. }
  949. //__________________________________________________________________________
  950. void MpdEctTrackFinderTpc::Write()
  951. {
  952. /// Write
  953. TFile histoFile("EctRec.root","RECREATE");
  954. Writedir2current(fHistoDir);
  955. histoFile.Close();
  956. }
  957. //__________________________________________________________________________
  958. void MpdEctTrackFinderTpc::Writedir2current( TObject *obj )
  959. {
  960. /// Write
  961. if( !obj->IsFolder() ) obj->Write();
  962. else{
  963. TDirectory *cur = gDirectory;
  964. TDirectory *sub = cur->mkdir(obj->GetName());
  965. sub->cd();
  966. TList *listSub = ((TDirectory*)obj)->GetList();
  967. TIter it(listSub);
  968. while( TObject *obj1=it() ) Writedir2current(obj1);
  969. cur->cd();
  970. }
  971. }
  972. //__________________________________________________________________________
  973. void MpdEctTrackFinderTpc::SelectTracks(Int_t iPass)
  974. {
  975. /// Do track selection and compute shared hit multiplicities
  976. //if (iPass) return;
  977. Int_t nFound = fTracks->GetEntriesFast(), nHitsTot = fKHits->GetEntriesFast();
  978. Int_t *index = new Int_t [nHitsTot];
  979. for (Int_t i = 0; i < nHitsTot; ++i) index[i] = 0;
  980. for (Int_t i = 0; i < nFound; ++i) {
  981. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  982. Int_t nHits = track->GetNofHits();
  983. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  984. if (nHits == 0) continue;
  985. TObjArray *hits = track->GetHits();
  986. for (Int_t j = 0; j < nHits; ++j) {
  987. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  988. Int_t ind = fKHits->IndexOf(hit);
  989. if (index[ind]) hit->SetFlag(TMath::Abs(hit->GetFlag())+10);
  990. ++index[ind];
  991. }
  992. }
  993. //fTracks->Compress();
  994. // Remove ghost tracks (with too many shared hits)
  995. nFound = fTracks->GetEntriesFast();
  996. //Int_t *nh = new Int_t [nFound];
  997. Double_t *etas = new Double_t [nFound];
  998. Int_t *order = new Int_t [nFound];
  999. for (Int_t i = 0; i < nFound; ++i) {
  1000. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1001. //nh[i] = track->GetNofHits();
  1002. etas[i] = TMath::Abs (track->Momentum3().Eta());
  1003. }
  1004. //TMath::Sort(nFound,nh,order,kFALSE); // in ascending order
  1005. TMath::Sort(nFound,etas,order,kTRUE); // in descending order
  1006. for (Int_t i = 0; i < nFound; ++i) {
  1007. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(order[i]);
  1008. Int_t nHits = track->GetNofHits();
  1009. TObjArray *hits = track->GetHits();
  1010. Int_t nShared = 0;
  1011. for (Int_t j = 0; j < nHits; ++j) {
  1012. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1013. if (hit->GetFlag()%1000 >= 10) ++nShared;
  1014. }
  1015. if (1.*nShared / nHits >= 0.5) {
  1016. // Remove track
  1017. // Update hit sharing
  1018. for (Int_t j = 0; j < nHits; ++j) {
  1019. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1020. if (hit->GetFlag()%1000 >= 10) hit->SetFlag(hit->GetFlag()-10);
  1021. else {
  1022. if (!fMirror) {
  1023. hit->SetFlag(TMath::Abs(hit->GetFlag())); // use hit again
  1024. } else {
  1025. // Check mirror hit
  1026. Int_t nLay = GetNofHitsInLayer(hit->GetLayer());
  1027. Int_t indx0 = GetHitsInLayer(hit->GetLayer());
  1028. for (Int_t indx = 0; indx < nLay; ++indx) {
  1029. MpdKalmanHit *h = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  1030. if (h == hit) continue;
  1031. if (h->GetDetectorID() == hit->GetDetectorID()) {
  1032. if (index[indx0+indx]) break; // do not use hit again
  1033. hit->SetFlag(TMath::Abs(hit->GetFlag())); // use hit again
  1034. h->SetFlag(TMath::Abs(h->GetFlag())); // use hit again
  1035. break;
  1036. }
  1037. }
  1038. }
  1039. }
  1040. }
  1041. fTracks->RemoveAt(order[i]);
  1042. cout << " Removing track: " << nHits << " " << nShared << endl;
  1043. }
  1044. }
  1045. fTracks->Compress();
  1046. //delete [] nh;
  1047. delete [] etas;
  1048. delete [] index;
  1049. delete [] order;
  1050. }
  1051. //__________________________________________________________________________
  1052. void MpdEctTrackFinderTpc::RemoveDoubles()
  1053. {
  1054. /// Remove double tracks (if number of common hits greater than 50% of hits on track)
  1055. Int_t nFound = fTrackCand->GetEntriesFast(), nOK = 0;
  1056. for (Int_t i = 0; i < nFound; ++i) {
  1057. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  1058. if (track == 0x0) continue;
  1059. Int_t nh = track->GetNofHits();
  1060. if (nh == 0) continue;
  1061. //TObjArray *hits = track->GetHits();
  1062. //cout << i << " " << nh << " " << ++nOK << endl;
  1063. for (Int_t j = i + 1; j < nFound; ++j) {
  1064. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(j);
  1065. if (track1 == 0x0) continue;
  1066. Int_t nh1 = track1->GetNofHits();
  1067. if (nh1 == 0) continue;
  1068. //TObjArray *hits1 = track1->GetHits();
  1069. Int_t nc = NofCommonHits(track, track1);
  1070. if (1.*nc/TMath::Min(nh,nh1) < 0.5) continue;
  1071. if (nh > nh1+5) fTrackCand->RemoveAt(j);
  1072. else if (nh < nh1-5) {
  1073. fTrackCand->RemoveAt(i);
  1074. --nOK;
  1075. break;
  1076. } else {
  1077. if (track->GetChi2()/nh > track1->GetChi2()/nh1) {
  1078. fTrackCand->RemoveAt(i);
  1079. --nOK;
  1080. break;
  1081. }
  1082. fTrackCand->RemoveAt(j);
  1083. }
  1084. }
  1085. }
  1086. fTrackCand->Compress();
  1087. }
  1088. //__________________________________________________________________________
  1089. Int_t MpdEctTrackFinderTpc::NofCommonHits(MpdEctKalmanTrack* track, MpdEctKalmanTrack* track1, Int_t dir)
  1090. {
  1091. /// Compute number of common hits on 2 tracks
  1092. TObjArray *hits = track->GetHits(), *hits1 = track1->GetHits();
  1093. MpdKalmanHit *hit = (MpdKalmanHit*) hits->First();
  1094. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->First();
  1095. Int_t nCom = 0;
  1096. while (hit && hit1) {
  1097. if (hit->GetLayer() > hit1->GetLayer()) {
  1098. hit1 = (MpdKalmanHit*) hits1->After(hit1);
  1099. continue;
  1100. }
  1101. if (hit == hit1) ++nCom;
  1102. hit = (MpdKalmanHit*) hits->After(hit);
  1103. }
  1104. return nCom;
  1105. }
  1106. //__________________________________________________________________________
  1107. void MpdEctTrackFinderTpc::AddHits()
  1108. {
  1109. /// Add hit objects to tracks and compute number of wrongly assigned hits
  1110. /// (hits with ID different from ID of starting TPC track)
  1111. Int_t nFound = fTracks->GetEntriesFast();
  1112. for (Int_t i = 0; i < nFound; ++i) {
  1113. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1114. Int_t nHits = track->GetNofHits();
  1115. if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1116. cout << i << " " << nHits << " " << track->GetTrackID() << " " << track->GetChi2() << " "
  1117. << track->Phi()*TMath::RadToDeg() << " " << track->Momentum3().Eta() << " "
  1118. << 1./track->GetParam(4) << " " << track->GetMisses() << endl;
  1119. track->SetNofHits(nHits);
  1120. TClonesArray &trHits = *track->GetTrHits();
  1121. TObjArray *hits = track->GetHits();
  1122. Int_t nWrong = 0, nMirr = 0, motherID = track->GetTrackID(), motherID1 = 0;
  1123. // Get track mother ID
  1124. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  1125. while (mctrack->GetMotherId() >= 0) {
  1126. motherID = mctrack->GetMotherId();
  1127. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1128. }
  1129. for (Int_t j = 0; j < nHits; ++j) {
  1130. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1131. new (trHits[j]) MpdKalmanHit(*hit);
  1132. Int_t iproj = (hit->GetLayer() - 1) % 3;
  1133. if (iproj == 0) cout << " R";
  1134. else if (iproj == 1) cout << " U";
  1135. else cout << " V";
  1136. cout << (hit->GetLayer()-1) % fgkNlays2 + 1;
  1137. Int_t hitId = HitMotherId(hit,motherID,motherID1);
  1138. //cout << "-" << motherID1;
  1139. cout << "-" << hitId;
  1140. if (hit->Index()->GetSize() > 1) cout << "-" << hit->Index()->GetSize();
  1141. if (hit->IsMirror()) cout << "(M)";
  1142. if (motherID1 != motherID) ++nWrong;
  1143. else if (hit->IsMirror()) ++nMirr;
  1144. if (hit->GetFlag() > -1) hit->SetFlag(-hit->GetFlag());
  1145. }
  1146. cout << "\n" << nWrong << " " << nMirr << endl;
  1147. track->SetNofWrong(nWrong);
  1148. track->SetMirrors(nMirr);
  1149. track->SetLastLay();
  1150. //track->GetParam()->Print();
  1151. }
  1152. fTracks->Compress();
  1153. }
  1154. //__________________________________________________________________________
  1155. Int_t MpdEctTrackFinderTpc::HitMotherId(MpdKalmanHit* hit, Int_t idM, Int_t &id1)
  1156. {
  1157. /// Check if hit has the same mother ID as idM
  1158. Int_t nOver = hit->Index()->GetSize(), idHit = 0;
  1159. for (Int_t i = 0; i < nOver; ++i) {
  1160. FairMCPoint *h = (FairMCPoint*) fEctHits->UncheckedAt(hit->GetIndex(i));
  1161. id1 = idHit = h->GetTrackID();
  1162. // Get point mother ID
  1163. MpdMCTrack* mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(id1);
  1164. while (mctrack->GetMotherId() >= 0) {
  1165. id1 = mctrack->GetMotherId();
  1166. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1167. }
  1168. if (id1 == idM) return idHit;
  1169. }
  1170. return idHit;
  1171. }
  1172. //__________________________________________________________________________
  1173. void MpdEctTrackFinderTpc::StoreTracks(Int_t iPass)
  1174. {
  1175. /// Transfer tracks from fTrackCand to fTracks
  1176. static const Int_t nHitMin = 10, nHitMinAbs = 5;
  1177. static const Double_t etaMin = 1.3, etaMax = 1.5;
  1178. Int_t nFound = fTracks->GetEntriesFast(), nCand = fTrackCand->GetEntriesFast();
  1179. for (Int_t i = 0; i < nCand; ++i) {
  1180. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  1181. if (fDetType == 1) {
  1182. // Selection criteria for ECT
  1183. if (track->GetNofHits() < nHitMinAbs) continue;
  1184. Double_t eta = TMath::Abs(track->Momentum3().Eta());
  1185. //if (iPass == 0 && track->GetChi2() / (track->GetNofHits()-5) > 3.0) continue;
  1186. //if (iPass == 0 && eta < etaMin) continue;
  1187. if (eta < etaMin) continue;
  1188. if (eta > etaMax && track->GetNofHits() < nHitMin) continue;
  1189. // !!!!
  1190. //if (track->GetMisses() >= track->GetNofHits()) continue;
  1191. // !!!
  1192. } else if (fDetType == 2 && iPass == 0) {
  1193. // Selection criteria for DCH
  1194. if (track->GetNofHits() < 11) continue;
  1195. }
  1196. track->Weight2Cov();
  1197. //SetTrackID(track);
  1198. new ((*fTracks)[nFound++]) MpdEctKalmanTrack(*track);
  1199. //fTpcTracks->RemoveAt(track->GetTpcIndex());
  1200. fTpcTracks->UncheckedAt(track->GetTpcIndex())->SetUniqueID(1); // for one-job running of TPC and ECT tracking
  1201. }
  1202. }
  1203. //__________________________________________________________________________
  1204. void MpdEctTrackFinderTpc::ExcludeHits()
  1205. {
  1206. /// Exclude hits, already used for tracking, from consideration during the next passes
  1207. Int_t nReco = fTracks->GetEntriesFast();
  1208. cout << " nReco: " << nReco << endl;
  1209. for (Int_t i = 0; i < nReco; ++i) {
  1210. MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
  1211. Int_t nhitsKF = track->GetNofHits();
  1212. TObjArray *hits = track->GetHits();
  1213. for (Int_t j = 0; j < nhitsKF; ++j) {
  1214. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1215. hit->SetFlag(-TMath::Abs(hit->GetFlag()));
  1216. if (fMirror) {
  1217. // Exclude mirror hits as well
  1218. Int_t nLay = GetNofHitsInLayer(hit->GetLayer());
  1219. Int_t indx0 = GetHitsInLayer(hit->GetLayer());
  1220. for (Int_t indx = 0; indx < nLay; ++indx) {
  1221. MpdKalmanHit *h = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  1222. if (h == hit) continue;
  1223. if (h->GetDetectorID() == hit->GetDetectorID()) { h->SetFlag(-TMath::Abs(h->GetFlag())); break; }
  1224. }
  1225. }
  1226. }
  1227. }
  1228. }
  1229. ClassImp(MpdEctTrackFinderTpc);