MpdEctTrackFollow2Tpc.cxx 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553
  1. // -------------------------------------------------------------------------
  2. // ----- MpdEctTrackFollow2Tpc source file -----
  3. // ----- Created 25/08/08 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdEctTrackFollow2Tpc.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track propagator from MPD End-Cap Tracker (ECT) into TPC
  9. **/
  10. #include "MpdEctTrackFollow2Tpc.h"
  11. #include "MpdKalmanFilter.h"
  12. #include "MpdKalmanHitR.h"
  13. #include "MpdKalmanHitZ.h"
  14. #include "MpdEctKalmanTrack.h"
  15. #include "MpdStrawendcapPoint.h"
  16. #include "MpdEtofPoint.h"
  17. #include "TpcPoint.h"
  18. #include "TpcLheHit.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 MpdEctTrackFollow2Tpc::fgkChi2Cut = 25; //10; //20; //100;
  33. const Double_t MpdEctTrackFollow2Tpc::fgkLayCut = 10;
  34. //__________________________________________________________________________
  35. MpdEctTrackFollow2Tpc::MpdEctTrackFollow2Tpc(const char *name, Int_t iVerbose )
  36. :FairTask(name, iVerbose),
  37. fKHits(new TClonesArray("MpdKalmanHitR")),
  38. fTracks(new TClonesArray("MpdEctKalmanTrack", 100)),
  39. fHistoDir(0x0),
  40. fhLays(new TH1F("hLays","TPC layers",100,0,100)),
  41. fLayPointers(0x0)
  42. {
  43. }
  44. //__________________________________________________________________________
  45. MpdEctTrackFollow2Tpc::~MpdEctTrackFollow2Tpc()
  46. {
  47. delete fKHits;
  48. //delete fTracks;
  49. //delete fTrackCand;
  50. delete [] fLayPointers;
  51. delete fhLays;
  52. }
  53. //__________________________________________________________________________
  54. InitStatus MpdEctTrackFollow2Tpc::Init()
  55. {
  56. return ReInit();
  57. }
  58. //__________________________________________________________________________
  59. InitStatus MpdEctTrackFollow2Tpc::ReInit()
  60. {
  61. fEctPoints =(TClonesArray *) FairRootManager::Instance()->GetObject("STRAWPoint");
  62. fTpcPoints =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcPoint");
  63. fTpcHits = (TClonesArray *) FairRootManager::Instance()->GetObject("LheHit");
  64. fEctTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("EctTrack");
  65. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  66. fTpcTracks = 0x0; //(TClonesArray *) FairRootManager::Instance()->GetObject("TpcPoint"); // just for testing
  67. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  68. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  69. FairRootManager::Instance()->Register("EctTpcTrack", "Ect", fTracks, kTRUE);
  70. fNPass = 1;
  71. return kSUCCESS;
  72. }
  73. //__________________________________________________________________________
  74. void MpdEctTrackFollow2Tpc::Reset()
  75. {
  76. ///
  77. cout << " MpdEctTrackFollow2Tpc::Reset " << endl;
  78. fKHits->Clear();
  79. fTracks->Clear("C");
  80. delete [] fLayPointers;
  81. }
  82. //__________________________________________________________________________
  83. void MpdEctTrackFollow2Tpc::SetParContainers()
  84. {
  85. }
  86. //__________________________________________________________________________
  87. void MpdEctTrackFollow2Tpc::Finish()
  88. {
  89. //Write();
  90. }
  91. //__________________________________________________________________________
  92. void MpdEctTrackFollow2Tpc::Exec(Option_t * option)
  93. {
  94. static int eventCounter = 0;
  95. cout << " EctRec event " << ++eventCounter << endl;
  96. Reset();
  97. // Get TPC Kalman hits
  98. GetTpcHits();
  99. for (Int_t i = 0; i < fNPass; ++i) {
  100. fTracks->Clear();
  101. cout << " Total number of ECT tracks: " << fEctTracks->GetEntriesFast() << endl;
  102. FollowInTPC(); // propagate tracks into TPC
  103. //StoreTracks();
  104. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  105. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  106. }
  107. RemoveDoubles(); // remove double tracks
  108. AddHits(); // add hit objects to tracks
  109. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  110. }
  111. //__________________________________________________________________________
  112. void MpdEctTrackFollow2Tpc::GetTpcHits()
  113. {
  114. /// Get TPC hits near zMax
  115. fhLays->Reset();
  116. Int_t ids[99999] = {0}, idMax = 0;
  117. // !!!!!!!!! Geometry dependent values
  118. //Double_t zMax = 134.5, rMin = 20.5, rMax = 109.5, dR = (rMax-rMin)/50; // dR = 1.78; // old version
  119. Double_t zMax = 124.5, rMin = 29.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.4123; // new version (with dead material)
  120. Int_t nTpc = fTpcHits->GetEntriesFast(), lay0 = -1, layMax = -1;
  121. cout << " TPC hits: " << nTpc << endl;
  122. Double_t tanMax = (rMin + fgkLayCut * dR) / zMax;
  123. for (Int_t i = 0; i < nTpc; ++i) {
  124. TpcLheHit *hit = (TpcLheHit*) fTpcHits->UncheckedAt(i);
  125. //TpcPoint *p = 0x0;
  126. //if (fTpcTracks) p = (TpcPoint*) fTpcTracks->UncheckedAt(hit->GetRefIndex());
  127. //if (p && p->GetTrackID() != 64) continue; // !!!
  128. if (hit->GetZ() < 0 || hit->GetLayer() > fgkLayCut) {
  129. fTpcHits->RemoveAt(i);
  130. continue;
  131. }
  132. if (hit->GetR() / hit->GetZ() > tanMax) {
  133. fTpcHits->RemoveAt(i);
  134. continue;
  135. }
  136. lay0 = hit->GetLayer();
  137. //Int_t id = ((TpcPoint*)fTpcTracks->UncheckedAt(hit->GetRefIndex()))->GetTrackID();
  138. Int_t id = hit->GetTrackID();
  139. cout << i << " " << hit->GetLayer() << " " << hit->GetZ() << " " << id << endl;
  140. if (id > 99998) { fTpcHits->RemoveAt(i); continue; }
  141. idMax = TMath::Max (id, idMax);
  142. ids[id]++;
  143. layMax = TMath::Max (lay0, layMax);
  144. fhLays->Fill(lay0+0.1);
  145. //hit->SetFlag(1);
  146. }
  147. Int_t nid = 0;
  148. for (Int_t i = 0; i <= idMax; ++i) if (ids[i]) ++nid;
  149. fTpcHits->Compress();
  150. cout << " TPC points near Zmax: " << fTpcHits->GetEntriesFast() << " " << nid << " " << fhLays->GetSum() << endl;
  151. fLayPointers = new Int_t [layMax+1];
  152. Int_t ipos = 0;
  153. for (Int_t i = layMax; i > -1; --i) {
  154. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  155. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  156. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  157. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  158. if (cont == 0) {
  159. fLayPointers[i] = -1;
  160. continue;
  161. }
  162. fLayPointers[i] = ipos;
  163. ipos += cont;
  164. }
  165. // Create TPC Kalman hits
  166. Int_t nHits = fTpcHits->GetEntriesFast();
  167. for (Int_t i = 0; i < nHits; ++i) {
  168. TpcLheHit *hit = (TpcLheHit*) fTpcHits->UncheckedAt(i);
  169. MpdKalmanHitR *hitK = new ((*fKHits)[i]) MpdKalmanHitR(hit->GetR(), hit->GetRphi(), hit->GetZ(),
  170. //hit->GetRphiErr(), hit->GetZerr(), lay, j);
  171. hit->GetRphiErr(), hit->GetZerr(), hit->GetLayer(), hit->GetRefIndex());
  172. if (hit->GetTrackID() == 1967) cout << hit->GetTrackID() << " " << hit->GetLayer() << endl;
  173. }
  174. }
  175. //__________________________________________________________________________
  176. void MpdEctTrackFollow2Tpc::FollowInTPC()
  177. {
  178. /// Propagate tracks into TPC
  179. Int_t nCand = fEctTracks->GetEntriesFast(), iok = 0;
  180. Double_t zMax = 134.5, rMin = 20.5, rMax = 109.5, dR = (rMax-rMin)/50; // dR = 1.78; // old version
  181. //Double_t zMax = 124.5, rMin = 29.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.4123; // new version (with dead material)
  182. MpdEctKalmanTrack track;
  183. for (Int_t i = 0; i < nCand; ++i) {
  184. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fEctTracks->UncheckedAt(i);
  185. cout << " Track seed No. " << i << " " << tr->GetNofTrHits() << " " << tr->GetTrackID() << endl;
  186. //if (track->GetTrackID() != 117) continue;
  187. //if (track->GetTrackID() != 10) continue;
  188. //if (track->GetTrackID() != 1105) continue;
  189. //if (track->GetTrackID() != 77) continue;
  190. track = *tr;
  191. track.SetParamNew(*track.GetParam());
  192. track.ReSetWeight();
  193. Int_t nHits = tr->GetNofTrHits();
  194. for (Int_t ih = 0; ih < nHits; ++ih) {
  195. TObject *hit = tr->GetTrHits()->UncheckedAt(ih);
  196. track.GetHits()->Add(hit);
  197. }
  198. MpdKalmanHitZ hit;
  199. hit.SetType(MpdKalmanHit::kFixedZ);
  200. hit.SetZ(zMax);
  201. MpdKalmanFilter::Instance()->PropagateToHit(&track, &hit);
  202. Int_t lay0 = (Int_t) ((track.GetRadNew() - rMin) / dR);
  203. lay0 = TMath::Max (lay0, 0);
  204. cout << lay0 << endl;
  205. if (lay0 <= fgkLayCut) iok = RunKalmanFilter(&track, lay0);
  206. if (iok == -1) cout<<"RunKalmanFilter with error!!!"<<endl;
  207. hit = *((MpdKalmanHitZ*) track.GetHits()->Last());
  208. // Keep only tracks with TPC points
  209. if (hit.GetType() == MpdKalmanHit::kFixedR) new ((*fTracks)[i]) MpdEctKalmanTrack(track);
  210. //new ((*fTracks)[i]) MpdEctKalmanTrack(track);
  211. cout << track.GetNofHits() << endl;
  212. }
  213. }
  214. //__________________________________________________________________________
  215. Int_t MpdEctTrackFollow2Tpc::RunKalmanFilter(MpdEctKalmanTrack *track, Int_t layBeg)
  216. {
  217. /// Run Kalman filter
  218. MpdKalmanHitR *hitOK = 0x0;
  219. Int_t layEnd = -1, dLay = -1, layOK = -1;
  220. Int_t miss = 0;
  221. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  222. TMatrixD param(5,1), paramTmp(5,1);
  223. Double_t saveR = 0;
  224. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  225. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  226. Int_t nLay = GetNofHitsInLayer(lay);
  227. Int_t indx0 = GetHitsInLayer(lay);
  228. Double_t dChi2Min = 1.e+6;
  229. MpdKalmanHitR *hitMin = 0x0;
  230. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  231. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  232. /*
  233. if (trackDir == MpdKalmanTrack::kOutward) {
  234. // !!! This part is due to the simplified hit merging (digitization) -
  235. // different hit position inside one layer - should be removed later
  236. indxBeg = nLay - 1;
  237. indxEnd = -1;
  238. dIndx = -1;
  239. }
  240. */
  241. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  242. MpdKalmanHitR *hit = (MpdKalmanHitR*) fKHits->UncheckedAt(indx0+indx);
  243. // !!! Exact ID match
  244. //if (((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  245. // Exclude used hits
  246. if (hit->GetFlag() != 1) continue;
  247. if (track->GetType() == MpdKalmanTrack::kFixedR) {
  248. // Check if the hit within some window
  249. if (TMath::Abs(hit->GetRphi()-track->GetParamNew(0)) > 9) continue;
  250. if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  251. }
  252. //*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;
  253. //track->Print("");
  254. //hit->Print("");
  255. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) return -1;
  256. //*if (hit->GetTrackID() == -607)*/ cout << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << " " << hit->GetZ() << " " << track->GetPosNew() << endl;
  257. Double_t dChi2 = MpdKalmanFilter::Instance()->CheckHitR(track,hit,pointWeight,param);
  258. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  259. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  260. if (dChi2 < dChi2Min) {
  261. dChi2Min = dChi2;
  262. hitMin = hit;
  263. saveWeight = *track->GetWeight();
  264. saveR = track->GetPosNew();
  265. // temporary storage for the current track
  266. paramTmp = param;
  267. pointWeightTmp = pointWeight;
  268. //cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetR() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() << endl;
  269. //cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  270. //cout << hit->GetRphi() << " " << hit->GetZ() << endl;
  271. //cout << param(0,0) << " " << param(1,0) << endl;
  272. //paramTmp.Print();
  273. }
  274. }
  275. if (dChi2Min < fgkChi2Cut) {
  276. //layOK = lay;
  277. hitOK = hitMin;
  278. track->GetHits()->Add(hitOK);
  279. miss = 0;
  280. // Restore track params at the best hit
  281. track->SetChi2(track->GetChi2()+dChi2Min);
  282. saveWeight += pointWeightTmp;
  283. track->SetWeight(saveWeight);
  284. track->SetPosNew(saveR);
  285. track->SetParamNew(paramTmp);
  286. //cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  287. //paramTmp.Print();
  288. // Check if the accepted hit is the same as the seeded hit
  289. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  290. if (track->GetNofHits() == 1) track->SetParam1();
  291. } else {
  292. ++miss;
  293. //if (miss > 1) return -1;
  294. }
  295. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  296. } // for (Int_t lay = layOK-1; lay >= 0;
  297. return 0;
  298. }
  299. //__________________________________________________________________________
  300. Double_t MpdEctTrackFollow2Tpc::Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit)
  301. {
  302. /// Adjust hit coord. R-Phi to be "around" track R-Phi - to avoid
  303. /// discontinuity around +- Pi
  304. Double_t hitPhi = hit->GetRphi() / hit->GetR();
  305. Double_t phi0 = track->GetParamNew(0) / track->GetPosNew();
  306. return hit->GetR() * MpdKalmanFilter::Instance()->Proxim(phi0, hitPhi);
  307. }
  308. //__________________________________________________________________________
  309. void MpdEctTrackFollow2Tpc::Write()
  310. {
  311. /// Write
  312. TFile histoFile("EctRec.root","RECREATE");
  313. Writedir2current(fHistoDir);
  314. histoFile.Close();
  315. }
  316. //__________________________________________________________________________
  317. void MpdEctTrackFollow2Tpc::Writedir2current( TObject *obj )
  318. {
  319. /// Write
  320. if( !obj->IsFolder() ) obj->Write();
  321. else{
  322. TDirectory *cur = gDirectory;
  323. TDirectory *sub = cur->mkdir(obj->GetName());
  324. sub->cd();
  325. TList *listSub = ((TDirectory*)obj)->GetList();
  326. TIter it(listSub);
  327. while( TObject *obj1=it() ) Writedir2current(obj1);
  328. cur->cd();
  329. }
  330. }
  331. //__________________________________________________________________________
  332. void MpdEctTrackFollow2Tpc::RemoveDoubles()
  333. {
  334. /// Remove double tracks (if number of common hits greater than 50% of hits on track)
  335. Int_t nFound = fTracks->GetEntriesFast(), nOK = 0;
  336. for (Int_t i = 0; i < nFound; ++i) {
  337. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  338. if (track == 0x0) continue;
  339. Int_t nh = track->GetNofHits();
  340. cout << i << " " << nh << " " << ++nOK << endl;
  341. for (Int_t j = i + 1; j < nFound; ++j) {
  342. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  343. if (track1 == 0x0) continue;
  344. Int_t nh1 = track1->GetNofHits();
  345. Int_t nc = NofCommonHits(track, track1);
  346. if (1.*nc/TMath::Min(nh,nh1) < 0.5) continue;
  347. if (nh > nh1) fTracks->RemoveAt(j);
  348. else if (nh < nh1) {
  349. fTracks->RemoveAt(i);
  350. --nOK;
  351. break;
  352. } else {
  353. if (track->GetChi2() > track1->GetChi2()) {
  354. fTracks->RemoveAt(i);
  355. --nOK;
  356. break;
  357. }
  358. fTracks->RemoveAt(j);
  359. }
  360. }
  361. }
  362. fTracks->Compress();
  363. // Remove double tracks (originated from the same ETOF hit)
  364. //*
  365. nFound = fTracks->GetEntriesFast();
  366. for (Int_t i = 0; i < nFound; ++i) {
  367. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  368. if (track == 0x0) continue;
  369. Int_t iTof = track->GetTofIndex();
  370. Int_t nh = track->GetNofHits();
  371. for (Int_t j = i + 1; j < nFound; ++j) {
  372. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  373. if (track1 == 0x0) continue;
  374. Int_t iTof1 = track1->GetTofIndex();
  375. if (iTof1 != iTof) continue;
  376. Int_t nh1 = track1->GetNofHits();
  377. if (nh > nh1) fTracks->RemoveAt(j);
  378. else if (nh < nh1) {
  379. fTracks->RemoveAt(i);
  380. break;
  381. } else {
  382. if (track->GetChi2() > track1->GetChi2()) {
  383. fTracks->RemoveAt(i);
  384. break;
  385. }
  386. fTracks->RemoveAt(j);
  387. }
  388. }
  389. }
  390. fTracks->Compress();
  391. //*/
  392. }
  393. //__________________________________________________________________________
  394. Int_t MpdEctTrackFollow2Tpc::NofCommonHits(MpdEctKalmanTrack* track, MpdEctKalmanTrack* track1)
  395. {
  396. /// Compute number of common hits on 2 tracks
  397. TObjArray *hits = track->GetHits(), *hits1 = track1->GetHits();
  398. MpdKalmanHit *hit = (MpdKalmanHit*) hits->First();
  399. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->First();
  400. Int_t dir = 1;
  401. if (hit->GetLayer() < ((MpdKalmanHit*) hits->Last())->GetLayer()) dir = -1;
  402. Int_t nCom = 0;
  403. while (hit && hit1) {
  404. if (dir * hit->GetLayer() < dir * hit1->GetLayer()) {
  405. hit1 = (MpdKalmanHit*) hits1->After(hit1);
  406. continue;
  407. }
  408. if (hit == hit1) ++nCom;
  409. hit = (MpdKalmanHit*) hits->After(hit);
  410. }
  411. return nCom;
  412. }
  413. //__________________________________________________________________________
  414. void MpdEctTrackFollow2Tpc::AddHits()
  415. {
  416. /// Add hit objects to tracks and compute number of wrongly assigned hits
  417. /// (hits with ID different from ID of starting track)
  418. Int_t nFound = fTracks->GetEntriesFast();
  419. for (Int_t i = 0; i < nFound; ++i) {
  420. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  421. Int_t nHits = track->GetNofHits();
  422. if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  423. track->SetNofHits(nHits);
  424. TClonesArray &trHits = *track->GetTrHits();
  425. SetTrackID(track); // set track ID as ID of majority of its hits
  426. TObjArray *hits = track->GetHits();
  427. Int_t nWrong = 0, motherID = track->GetTrackID(), motherID1 = 0;
  428. cout << i << " " << nHits << " " << motherID << endl;
  429. // Get track mother ID
  430. FairMCTrack *mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(motherID);
  431. while (mctrack->GetMotherId() > 0) {
  432. motherID = mctrack->GetMotherId();
  433. mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  434. }
  435. for (Int_t j = 0; j < nHits; ++j) {
  436. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) hits->UncheckedAt(j);
  437. new (trHits[j]) MpdKalmanHitZ(*hit);
  438. if (hit->GetType() == MpdKalmanHit::kFixedZ) {
  439. // ECT point
  440. Int_t iproj = (hit->GetLayer() - 1) % 3;
  441. if (iproj == 0) cout << " U";
  442. else if (iproj == 1) cout << " R";
  443. else cout << " V";
  444. cout << hit->GetLayer();
  445. MpdStrawendcapPoint *h = (MpdStrawendcapPoint*) fEctPoints->UncheckedAt(hit->GetIndex());
  446. motherID1 = h->GetTrackID();
  447. cout << "-" << motherID1;
  448. } else {
  449. // TPC point
  450. cout << " " << hit->GetLayer();
  451. TpcPoint *h = (TpcPoint*) fTpcPoints->UncheckedAt(hit->GetIndex());
  452. motherID1 = h->GetTrackID();
  453. cout << "-" << motherID1;
  454. ((MpdKalmanHit*)trHits[j])->SetType(MpdKalmanHit::kFixedR);
  455. }
  456. // Get point mother ID
  457. FairMCTrack *mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(motherID1);
  458. while (mctrack->GetMotherId() > 0) {
  459. motherID1 = mctrack->GetMotherId();
  460. mctrack = (FairMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  461. }
  462. if (motherID1 != motherID) nWrong++;
  463. }
  464. cout << "\n" << nWrong << " " << track->GetTrackID() << " " << motherID << endl;
  465. track->SetNofWrong(nWrong);
  466. track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
  467. }
  468. fTracks->Compress();
  469. }
  470. //__________________________________________________________________________
  471. void MpdEctTrackFollow2Tpc::SetTrackID(MpdEctKalmanTrack* track)
  472. {
  473. /// Set track ID as ID of majority of its hits
  474. const Int_t idMax = 99999;
  475. Int_t ids[idMax+1] = {0}, nHits = track->GetNofHits(), locMax = 0;
  476. TObjArray *hits = track->GetHits();
  477. FairMCPoint *p = 0x0;
  478. for (Int_t i = 0; i < nHits; ++i) {
  479. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  480. if (hit->GetType() == MpdKalmanHit::kFixedZ) p = (FairMCPoint*) fEctPoints->UncheckedAt(hit->GetIndex());
  481. else p = (FairMCPoint*) fTpcPoints->UncheckedAt(hit->GetIndex());
  482. Int_t id = p->GetTrackID();
  483. if (id > idMax) exit(0);
  484. ++ids[id];
  485. if (ids[id] > ids[locMax]) locMax = id;
  486. }
  487. if (ids[locMax] > 1) track->SetTrackID(locMax);
  488. }
  489. ClassImp(MpdEctTrackFollow2Tpc);