MpdTrackFinderIts5spd.cxx 45 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216
  1. // -------------------------------------------------------------------------
  2. // ----- MpdTrackFinderIts source file -----
  3. // ----- Created 21/07/09 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. // ----- Modified by VK for ideal hit producer and 5 spd-layers geometry -------
  6. /** MpdTrackFinderIts.h
  7. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  8. **
  9. ** Track finder in MPD Inner Tracking System (ITS) using seeds from TPC
  10. **/
  11. #include "MpdStsGeoPar.h"
  12. #include "MpdTrackFinderIts5spd.h"
  13. #include "MpdKalmanFilter.h"
  14. #include "MpdKalmanGeoScheme.h"
  15. #include "MpdKalmanHit.h"
  16. #include "MpdKalmanTrack.h"
  17. //#include "MpdKalmanStripHit.h"
  18. #include "MpdItsHit5spd.h"
  19. #include "MpdStsPoint.h"
  20. #include "MpdItsKalmanTrack.h"
  21. #include "MpdTpcKalmanFilter.h"
  22. #include "MpdTpcKalmanTrack.h"
  23. #include "MpdMCTrack.h"
  24. #include "FairGeoNode.h"
  25. #include "FairMCPoint.h"
  26. #include "FairRootManager.h"
  27. #include "FairRun.h"
  28. #include "FairRunAna.h"
  29. #include "FairRuntimeDb.h"
  30. #include "TGeoMatrix.h"
  31. #include "TGeoBBox.h"
  32. #include "TGeoNode.h"
  33. #include "TGeoTube.h"
  34. #include "TGeoManager.h"
  35. #include "TMath.h"
  36. #include "TFile.h"
  37. //#include "TLorentzVector.h"
  38. #include "TVector2.h"
  39. #include "TClonesArray.h"
  40. #include <TRandom.h>
  41. #include <algorithm>
  42. #include <iostream>
  43. #include <map>
  44. using std::cout;
  45. using std::endl;
  46. //using std::map;
  47. const Double_t MpdTrackFinderIts5spd::fgkChi2Cut = 20; //20; //100;
  48. //__________________________________________________________________________
  49. MpdTrackFinderIts5spd::MpdTrackFinderIts5spd(Bool_t useVector, const char *name, Int_t iVerbose )
  50. : FairTask(name, iVerbose),
  51. fExact(0), //1),
  52. fGeo(0),
  53. fUseVector(useVector),
  54. fTpcKF(nullptr)
  55. {
  56. if (!fUseVector) fKHits = new TClonesArray("MpdKalmanHit", 100);
  57. fTracks = new TClonesArray("MpdItsKalmanTrack", 100);
  58. fHistoDir = 0x0;
  59. fhLays = new TH1F("hLaysITS","ITS layers",12,0,12);
  60. fLayPointers = 0x0;
  61. }
  62. //__________________________________________________________________________
  63. MpdTrackFinderIts5spd::~MpdTrackFinderIts5spd()
  64. {
  65. //delete fKHits;
  66. //delete fTracks;
  67. //delete fTrackCand;
  68. delete [] fLayPointers;
  69. delete fhLays;
  70. }
  71. //__________________________________________________________________________
  72. InitStatus MpdTrackFinderIts5spd::Init()
  73. {
  74. //----- Get pipe radius
  75. fPipeR = ((TGeoTube*)gGeoManager->GetVolume("pipe1")->GetShape())->GetRmax();
  76. cout << "************** Pipe radius: " << fPipeR << " cm " << endl;
  77. //return ReInit();
  78. if (ReInit() != kSUCCESS) return kERROR;
  79. //----- Read database
  80. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  81. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  82. // cout << geoPar << endl;
  83. TString volName = "sts01 ", path = "";
  84. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  85. // cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  86. //----- Process modular geometry
  87. fGeo = 1;
  88. Int_t nLay = 5;
  89. //Double_t safety = 0.005;
  90. for (Int_t i = 0; i < nLay; ++i) {
  91. fNladders[i] = fNsectors[i] = 0;
  92. volName = "sts01ladder";
  93. volName += (i+1);
  94. path = "/cave_1/sts01_0/" + volName;
  95. path += "_";
  96. TString path0 = path;
  97. //----- Loop over all ladders to find the one with the smallest radius
  98. //fRad[i] = 999.9;
  99. fRad[i] = 0.;
  100. Int_t nladd = 0;
  101. for (Int_t jlad = 1; jlad < 99; ++jlad) {
  102. path = path0;
  103. path += jlad;
  104. if (!gGeoManager->CheckPath(path)) break;
  105. gGeoManager->cd(path);
  106. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  107. //TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  108. Double_t xyzL[3] = {0}, xyzM[3];
  109. gGeoManager->LocalToMaster(xyzL,xyzM);
  110. Double_t rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  111. fRad[i] += rad;
  112. ++nladd;
  113. /*AZ - 3.12.2017
  114. fRad[i] = TMath::Min (fRad[i],rad);
  115. xyzL[0] = box->GetDX();
  116. gGeoManager->LocalToMaster(xyzL,xyzM);
  117. rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  118. fRad[i] = TMath::Min (fRad[i],rad);
  119. xyzL[0] = -box->GetDX();
  120. gGeoManager->LocalToMaster(xyzL,xyzM);
  121. rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  122. fRad[i] = TMath::Min (fRad[i],rad);
  123. */
  124. }
  125. cout << " *** Layer # " << i << " min_R= " << fRad[i] << endl;
  126. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  127. if (ladd == nullptr) { nLay = i; break; }
  128. fRad[i] /= nladd; // mean radius
  129. /*
  130. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  131. //safety = -box->GetDY();
  132. //safety = box->GetDY();
  133. safety = 2 * box->GetDY(); // new
  134. fRad[i] -= safety;
  135. */
  136. cout << " +++ Layer # " << i << " safety min_R= " << fRad[i] << endl;
  137. }
  138. FillGeoScheme();
  139. //----- Get pipe radius
  140. //fPipeR = ((TGeoTube*)gGeoManager->GetVolume("pipe1")->GetShape())->GetRmax();
  141. //cout << "************** Pipe radius: " << fPipeR << " cm " << endl;
  142. //----- Get cables
  143. TObjArray *vols = gGeoManager->GetListOfVolumes();
  144. Int_t nvols = vols->GetEntriesFast(), ncables = 0;
  145. for (Int_t i = 0; i < nvols; ++i) {
  146. TGeoVolume *vol = (TGeoVolume*) vols->UncheckedAt(i);
  147. TString cable = TString(vol->GetName());
  148. if (!(cable.Contains("sts") && cable.Contains("cable"))) continue;
  149. // cout << cable << endl;
  150. ++ncables;
  151. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  152. TString lad = cable;
  153. Int_t ip = lad.Index("cable");
  154. lad.Replace(ip,lad.Length()-ip+1,"");
  155. lad.Replace(lad.Length()-2,1,"");
  156. lad += "_1/";
  157. path = "/cave_1/sts01_0/" + lad + cable;
  158. path += "_1";
  159. gGeoManager->cd(path);
  160. Double_t xyzL[3] = {0}, xyzM[3];
  161. gGeoManager->LocalToMaster(xyzL,xyzM);
  162. // cout << xyzM[2] - box->GetDZ() << " " << xyzM[2] + box->GetDZ() << " " << box->GetDZ() << endl;
  163. if (xyzM[2] - box->GetDZ() > 0) {
  164. Int_t lay = TString(lad(lad.Length()-4,1)).Atoi();
  165. fCables[lay-1].insert(pair<Double_t,Double_t>(xyzM[2] - box->GetDZ(),xyzM[2] + box->GetDZ()));
  166. }
  167. }
  168. cout << "-I- MpdTrackFinderIts5spd: Intialization finished successfully" << endl;
  169. return kSUCCESS;
  170. }
  171. //__________________________________________________________________________
  172. void MpdTrackFinderIts5spd::FillGeoScheme()
  173. {
  174. //----- Fill Kalman filter geometry manager info
  175. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  176. TGeoVolume *vSts = gGeoManager->GetVolume("sts01");
  177. //----- Loop over layers
  178. for (Int_t layer = 1; layer < 999; ++layer) {
  179. TString sladder = "sts01ladder";
  180. sladder += layer;
  181. TGeoVolume *vLay = gGeoManager->GetVolume(sladder);
  182. if (vLay == 0x0 && layer == 1) continue; // no first layer
  183. if (vLay == 0x0) break;
  184. sladder += "_";
  185. //----- Loop over ladders
  186. for (Int_t ladder = 1; ladder < 999; ++ladder) {
  187. TString sladder1 = sladder;
  188. sladder1 += ladder;
  189. TGeoNode *node = vSts->FindNode(sladder1);
  190. if (node == 0x0) break;
  191. ++fNladders[layer-1];
  192. TGeoVolume *vLad = node->GetVolume();
  193. //cout << "Node: " << layer << " " << vLad->GetNodes()->GetEntriesFast() << " " << vLad->GetNdaughters() << endl;
  194. Int_t nDaught = vLad->GetNdaughters(), detID = -1, detIDsts = -1;
  195. TObjArray *daught = vLad->GetNodes();
  196. // for (Int_t j = 0; j < vLad->GetNdaughters(); ++j)
  197. // cout << "Name: " << ((TGeoNode*)(vLad->GetNodes()->UncheckedAt(j)))->GetName() << endl;
  198. //----- Loop over ladder daughters
  199. Int_t iZ = 0;
  200. for (Int_t det = 0; det < nDaught; ++det) {
  201. TString sdet1 = ((TGeoNode*)(daught->UncheckedAt(det)))->GetName();
  202. // if (!sdet1.Contains("sensor") && !sdet1.Contains("sector")) continue;
  203. if (ladder == 1) ++fNsectors[layer-1];
  204. Int_t det1 = TString(sdet1(sdet1.Index("_")+1,2)).Atoi();
  205. // Int_t secType = -1;
  206. // if (sdet1.Contains("sector")) secType = TString(sdet1(sdet1.Index("_")-2,1)).Atoi();
  207. ++iZ;
  208. Int_t side = 0.;
  209. detIDsts = fHitSts.SetDetId(layer, ladder, det1);
  210. detID = fHitSts.SetDetId(layer, ladder, iZ);
  211. //cout << "det1= " << det1 << " " << iZ << " layer= " << layer << " ladder=" << ladder << " detID= " << detID << endl;
  212. fId2Id[layer-1].insert(pair<Int_t,Int_t>(detIDsts,detID)); // detIDsts == detID for pixel detectors (AZ)
  213. detID += 1000000 * (layer-1);
  214. TString detName = sladder1 + "/" + sdet1 + "#";
  215. detName += side;
  216. geo->SetDetId(detName, detID);
  217. TString path = "/cave_1/sts01_0/" + detName(0,detName.Length()-2);
  218. //cout << "DetName: " << detName << " Path: " << path << endl;
  219. gGeoManager->cd(path);
  220. geo->SetPath(detID, path);
  221. node = gGeoManager->GetCurrentNode();
  222. // cout << node->GetName() << " " << detID << endl;
  223. TGeoVolume *vol = node->GetVolume();
  224. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  225. TVector2 size(box->GetDX(), box->GetDZ());
  226. geo->SetSize(detID, size);
  227. // cout << "##### Node: " << node->GetName() << " DX= " << box->GetDX() << " DZ=" << box->GetDZ() << endl;
  228. Double_t xyzL[3] = {0}, xyzM[3], vecM[3];
  229. xyzL[1] = 1;
  230. gGeoManager->LocalToMasterVect(xyzL,vecM);
  231. xyzL[1] = 0;
  232. gGeoManager->LocalToMaster(xyzL,xyzM);
  233. Double_t sign = TMath::Sign (1.,xyzM[0]*vecM[0]+xyzM[1]*vecM[1]);
  234. TVector3 pos(xyzM[0], xyzM[1], xyzM[2]);
  235. geo->SetGlobalPos(detID, pos);
  236. TVector3 norm(sign*vecM[0], sign*vecM[1], sign*vecM[2]);
  237. geo->SetNormal(detID, norm);
  238. }
  239. }
  240. }
  241. }
  242. //__________________________________________________________________________
  243. InitStatus MpdTrackFinderIts5spd::ReInit()
  244. {
  245. fItsPoints = (TClonesArray *) FairRootManager::Instance()->GetObject("StsPoint");
  246. fItsHits = (TClonesArray *) FairRootManager::Instance()->GetObject("StsHit");
  247. if (fItsPoints == 0x0 || fItsHits == 0x0) return kERROR;
  248. fTpcTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  249. fMCTracks = (TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  250. if (fUseVector) fKHits = (TClonesArray *) FairRootManager::Instance()->GetObject("ItsKHits"); // use from VectorFinder
  251. FairRootManager::Instance()->Register("ItsTrackKF", "Its", fTracks, kTRUE);
  252. fNPass = 1;
  253. return kSUCCESS;
  254. }
  255. //__________________________________________________________________________
  256. void MpdTrackFinderIts5spd::Reset()
  257. {
  258. cout << " MpdTrackFinderIts5spd::Reset " << endl;
  259. if (!fUseVector) fKHits->Delete();
  260. fTracks->Delete();
  261. delete [] fLayPointers;
  262. fLayPointers = nullptr;
  263. }
  264. //__________________________________________________________________________
  265. void MpdTrackFinderIts5spd::SetParContainers()
  266. {
  267. //----- Get run and runtime database
  268. FairRunAna* run = FairRunAna::Instance();
  269. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  270. FairRuntimeDb* db = run->GetRuntimeDb();
  271. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  272. //----- Get STS geometry parameter container
  273. db->getContainer("MpdStsGeoPar");
  274. }
  275. //__________________________________________________________________________
  276. void MpdTrackFinderIts5spd::Finish()
  277. {
  278. // Write();
  279. }
  280. //__________________________________________________________________________
  281. void MpdTrackFinderIts5spd::Exec(Option_t * option)
  282. {
  283. static int eventCounter = 0;
  284. cout << " ##### MpdTrackFinderIts5spd \n ItsRec event " << ++eventCounter << endl;
  285. Reset();
  286. //----- Create Kalman hits
  287. if (fItsHits->GetEntriesFast() == 0) return;
  288. MakeKalmanHits();
  289. for (Int_t i = 0; i < fNPass; ++i) {
  290. fTracks->Clear();
  291. GetTrackSeeds(i);
  292. cout << " Total number of hits for tracking: " << fKHits->GetEntriesFast() << endl;
  293. cout << " Total number of track seeds: " << fTracks->GetEntriesFast() << endl;
  294. DoTracking(i);
  295. // StoreTracks();
  296. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  297. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  298. }
  299. AddHits(); // add hit objects to tracks
  300. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  301. }
  302. //__________________________________________________________________________
  303. void MpdTrackFinderIts5spd::MakeKalmanHits()
  304. {
  305. //----- Create Kalman hits from ITS hits.
  306. fhLays->Reset();
  307. Int_t nHits = fItsHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0;
  308. //Double_t r, z, xloc, errZ = 0.00012, errX = 0.00023; // 120um in Z, 23um in R-Phi (local X)
  309. //Double_t r, z, xloc, errZ = 0.012, errX = 0.0023; // 120um in Z, 23um in R-Phi (local X)
  310. //Double_t r, z, xloc, errZ = 0.12, errX = 0.0023; // 1.2mm in Z, 23um in R-Phi (local X)
  311. Double_t r, z, xloc, errZ = 0.001, errX = 0.001; // 10um in Z, 10um in R-Phi (local X)
  312. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  313. if (!fUseVector) {
  314. for (Int_t ih = 0; ih < nHits; ++ih) {
  315. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(ih);
  316. r = TMath::Sqrt (h->GetX() * h->GetX() + h->GetY() * h->GetY());
  317. z = h->GetZ();
  318. xloc = h->GetLocalX();
  319. lay = h->Layer() - 1;
  320. // cout <<" n= " << ih << " r= " << r << " z= " << z << " layer= " << h->Layer() << " xloc= " << xloc << endl;
  321. //cout << h->GetDetectorID() << " " << geo->Path(h->GetDetectorID()+1000000*lay) << endl;
  322. // Get local Z
  323. TString path = geo->Path(h->GetDetectorID()+1000000*lay);
  324. gGeoManager->cd(path);
  325. Double_t posLoc[3] = {0}, posMas[3] = {h->GetX(), h->GetY(), h->GetZ()};
  326. gGeoManager->MasterToLocal(posMas,posLoc);
  327. //----- Add error
  328. Double_t dX = 0, dZ = 0;
  329. gRandom->Rannor(dX,dZ);
  330. // if (errZ > 2) dZ = 0.0; // 1-D case
  331. //dZ = 0.0; // 1-D case
  332. //Double_t meas[2] = {xloc+dX*errX, z+dZ*errZ};
  333. Double_t meas[2] = {xloc+dX*errX, posLoc[2]+dZ*errZ};
  334. Double_t err[2] = {errX, errZ};
  335. Double_t cossin[2] = {TMath::Cos(fStereoA[0]), TMath::Sin(fStereoA[0])};
  336. MpdKalmanHit *hit = 0x0;
  337. // cout << h->GetDetectorID() << " " << fId2Id[lay][h->GetDetectorID()] << endl;
  338. if (fGeo && h->GetUniqueID()) {
  339. hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+fId2Id[lay][h->GetDetectorID()], 1,
  340. MpdKalmanHit::kFixedP, meas, err, cossin, 0., r, ih);
  341. // cout << " DetID = " << h->GetDetectorID() << " " << lay*1000000+fId2Id[lay][h->GetDetectorID()] << endl;
  342. }
  343. hit->SetUniqueID(0);
  344. hit->SetNofDim(2);
  345. }
  346. }
  347. fKHits->Sort(); // in descending order in R
  348. layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  349. nKH = fKHits->GetEntriesFast();
  350. for (Int_t ih = 0; ih < nKH; ++ih) {
  351. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(ih);
  352. //if (hit->GetFlag() < 0) hit->SetFlag(hit->GetFlag() | MpdKalmanHit::kUsed);
  353. fhLays->Fill(hit->GetLayer()+0.1);
  354. }
  355. for (Int_t k = 0; k <= layMax; ++k)
  356. cout << "Layer= " << k+1 << " Number of hits=" << fhLays->GetBinContent(k+1,0) << endl;
  357. cout << " Max layer = " << layMax << " Number of hits= " << fKHits->GetEntriesFast() << endl;
  358. for (lay = 0; lay < fgkNlays; ++lay) {
  359. fHitMapRphi[lay].clear();
  360. fHitMapZ[lay].clear();
  361. }
  362. for (Int_t ihit = 0; ihit < nKH; ++ihit) {
  363. MpdKalmanHit *hk = (MpdKalmanHit*) fKHits->UncheckedAt(ihit);
  364. lay = hk->GetLayer();
  365. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(hk->GetIndex());
  366. TVector3 pos;
  367. h->Position(pos);
  368. z = pos.Z();
  369. fHitMapZ[lay].insert(pair<Double_t,Int_t>(z,ihit));
  370. Double_t phi = pos.Phi();
  371. //if (phi < 0) phi += TMath::TwoPi();
  372. Double_t rphi = fRad[lay] * phi;
  373. fHitMapRphi[lay].insert(pair<Double_t,Int_t>(rphi,ihit));
  374. Double_t rphiMax = fRad[lay] * TMath::Pi();
  375. if (TMath::Abs(rphi-rphiMax) < 3.0) {
  376. // Add overlap 3 cm near phi=180
  377. if (rphi < 0) fHitMapRphi[lay].insert(pair<Double_t,Int_t>(rphi+2*rphiMax,ihit)); // overlap 3 cm near phi=180
  378. else fHitMapRphi[lay].insert(pair<Double_t,Int_t>(rphi-2*rphiMax,ihit)); // overlap 3 cm near phi=180
  379. }
  380. }
  381. // for (Int_t ihit = 0; ihit < fKHits->GetEntriesFast(); ++ihit) {
  382. // cout << "ITS hit number= " << ihit << " Position= "<< ((MpdKalmanHit*)fKHits->UncheckedAt(ihit))->GetPos() <<
  383. // " Distance= " << ((MpdKalmanHit*)fKHits->UncheckedAt(ihit))->GetDist() << endl;
  384. // }
  385. fLayPointers = new Int_t [layMax+1];
  386. Int_t ipos = 0;
  387. for (Int_t i = layMax; i > -1; --i) {
  388. //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
  389. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  390. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  391. Int_t cont = TMath::Nint (fhLays->GetBinContent(i+1,0));
  392. if (cont == 0) {
  393. fLayPointers[i] = -1;
  394. continue;
  395. }
  396. fLayPointers[i] = ipos;
  397. ipos += cont;
  398. }
  399. }
  400. //__________________________________________________________________________
  401. void MpdTrackFinderIts5spd::GetTrackSeeds(Int_t iPass)
  402. {
  403. /// Build ITS track seeds from TPC tracks
  404. Int_t nTpcTracks = fTpcTracks->GetEntriesFast();
  405. cout << " Seed tracks: " << nTpcTracks << endl;
  406. Int_t nCand = 0;
  407. MpdKalmanHit hit;
  408. hit.SetType(MpdKalmanHit::kFixedR);
  409. for (Int_t itr = 0; itr < nTpcTracks; ++itr) {
  410. MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(itr);
  411. //tpc->GetParam()->Print();
  412. //MpdEctKalmanTrack *track = new ((*fTracks)[nCand++]) MpdEctKalmanTrack(itr, *tpc);
  413. MpdItsKalmanTrack *track = new ((*fTracks)[nCand++]) MpdItsKalmanTrack(*tpc);
  414. track->SetUniqueID(itr+1);
  415. /*
  416. track->SetParam(*track->GetParamAtHit());
  417. track->SetParamNew(*track->GetParamAtHit());
  418. track->SetPosNew(track->GetPos());
  419. track->SetWeight(*track->GetWeightAtHit());
  420. track->SetLength(track->GetLengAtHit());
  421. */
  422. hit.SetPos(track->GetPos());
  423. track->SetParamNew(*track->GetParam());
  424. track->SetPos(track->GetPosNew());
  425. track->ReSetWeight();
  426. TMatrixDSym w = *track->GetWeight(); // save current weight matrix
  427. MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  428. track->SetWeight(w); // restore original weight matrix (near TPC inner shell)
  429. //cout << nCand-1 << " " << track->GetTrackID() << endl;
  430. //cout << track->GetHits()->GetEntriesFast() << " " << track->GetTrHits()->GetEntriesFast() << endl;
  431. track->GetHits()->Clear();
  432. track->SetChi2Its(track->GetChi2()); // temporary storage
  433. track->SetChi2(0.);
  434. track->SetDirection(MpdKalmanTrack::kInward); //AZ - 6.12.2017
  435. }
  436. cout << " Number of ITS track candidates: " << nCand << endl;
  437. }
  438. //__________________________________________________________________________
  439. void MpdTrackFinderIts5spd::DoTracking(Int_t iPass)
  440. {
  441. //----- Run Kalman tracking
  442. Double_t vert[3] = {0.0,0.0,0.0};
  443. Int_t nCand = fTracks->GetEntriesFast(), iok = 0;
  444. Int_t lay0 = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  445. for (Int_t i = 0; i < nCand; ++i) {
  446. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);
  447. // cout << " Track seed No. " << i << ", ID: " << track->GetTrackID() << ", Hits: " << track->GetNofTrHits() << endl;
  448. //for (Int_t j = 0; j < track->GetNofTrHits(); ++j) {
  449. //MpdKalmanHit *h = (MpdKalmanHit* )track->GetTrHits()->UncheckedAt(j);
  450. // cout << "Track hit number= " << j << " Distance= " << h->GetDist() << " Layer= " << h->GetLayer() << endl;
  451. //}
  452. iok = RunKalmanFilterMod(track, lay0); // modular geometry
  453. if (iok == -1) {
  454. fTracks->RemoveAt(i);
  455. continue;
  456. }
  457. // if (track->GetNofHits() == 0) continue; // no hits added
  458. //------ Propagate track to the beam line
  459. track->SetParam(*track->GetParamNew());
  460. track->SetPos(track->GetPosNew());
  461. Double_t pos = track->GetPos();
  462. TMatrixD par = *track->GetParam();
  463. TMatrixDSym cov = *track->Weight2Cov();
  464. Double_t leng = track->GetLength();
  465. TString nodeNew = track->GetNodeNew();
  466. //cout << " 1: " << nodeNew << ", " << track->GetNode() << " " << track->GetHits()->GetEntriesFast() << endl;
  467. // Go to beam pipe
  468. MpdKalmanHit hit;
  469. hit.SetType(MpdKalmanHit::kFixedR);
  470. hit.SetPos(fPipeR);
  471. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  472. if (iok != 1) {
  473. // Restore track
  474. track->SetParam(par);
  475. track->SetParamNew(par);
  476. track->SetCovariance(cov);
  477. track->ReSetWeight();
  478. track->SetPos(pos);
  479. track->SetPosNew(pos);
  480. track->SetLength(leng);
  481. //track->SetNode(node);
  482. //cout << " 2: " << nodeNew << ", " << track->GetNode() << " " << pos << endl;
  483. track->SetNodeNew(nodeNew);
  484. } else {
  485. // Add multiple scattering
  486. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  487. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  488. TMatrixDSym* pcov = track->Weight2Cov();
  489. Double_t th = track->GetParamNew(3);
  490. Double_t cosTh = TMath::Cos(th);
  491. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  492. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  493. (*pcov)(3,3) += angle2;
  494. Int_t ok = 0;
  495. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  496. track->SetWeight(*pcov);
  497. }
  498. cov = *track->Weight2Cov();
  499. hit.SetPos(0.);
  500. hit.SetMeas(0,track->GetParam(2)); // track Phi
  501. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  502. //Double_t pos = ((MpdKalmanHit*)track->GetHits()->Last())->GetPos();
  503. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  504. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  505. if (iok != 1) {
  506. track->SetNodeNew(""); // for FindPca
  507. MpdKalmanFilter::Instance()->FindPca(track, vert);
  508. track->SetCovariance(cov); //AZ-040121
  509. }
  510. //track->SetPos(pos); // restore position
  511. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  512. //track->GetCovariance()->Print();
  513. //cout << " 3: " << track->GetNodeNew() << ", " << track->GetNode() << endl;
  514. } // for (Int_t i = 0; i < nCand;
  515. fTracks->Compress();
  516. }
  517. //__________________________________________________________________________
  518. Int_t MpdTrackFinderIts5spd::RunKalmanFilterMod(MpdItsKalmanTrack *track, Int_t layBeg)
  519. {
  520. /// Run Kalman filter in the modular geometry (might not work when propagating outward!!!)
  521. const Int_t maxBr = 24, maxBrLay = 2000; // max number of branches
  522. //const Int_t maxBr = 100, maxBrLay = 1000; // max number of branches
  523. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07}; // layer thicknesses
  524. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  525. //cout << fHits->GetEntriesFast() << endl;
  526. //Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  527. //Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  528. //MpdKalmanHit *hitOK = 0x0;
  529. MpdKalmanHit hitTmp;
  530. //MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  531. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1;
  532. Int_t layEnd = -1, dLay = -1;//, layOK = -1;
  533. // if (trackDir == MpdKalmanTrack::kOutward) {
  534. // layEnd = layMax + 1;
  535. // dLay = 1;
  536. // }
  537. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  538. TMatrixD param(5,1), paramTmp(5,1);
  539. Double_t quality[maxBrLay] = {0.0}; //,saveZ = 0.0, saveLeng = 0.0, dChi2Min = 0.0, posNew = 0.0;
  540. // cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  541. MpdItsKalmanTrack branch[maxBr];
  542. Int_t nBr = 1, nBr0 = 0, indx_lay[maxBrLay] = {0}, maxInLay = 0, ok = 0, nAdd = 0;
  543. branch[0] = *track;
  544. TString mass2 = "0.0194797849"; // pion mass squared
  545. /*
  546. if (fMCTracks) {
  547. // Get particle mass - ideal PID
  548. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(track->GetTrackID());
  549. TParticlePDG *pdgP = TDatabasePDG::Instance()->GetParticle(mctrack->GetPdgCode());
  550. if (pdgP) {
  551. Double_t mass = pdgP->Mass();
  552. if (mass < 0.1 || mass > 0.25) {
  553. // Electrons or heavier than pions
  554. mass2 = "";
  555. mass2 += mass*mass;
  556. }
  557. }
  558. }
  559. */
  560. // cout << "##### Cycle start ##### " << endl;
  561. // cout << " layBeg, LayEnd: " << layBeg << " " << layEnd << endl;
  562. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  563. //Int_t nLay = GetNofHitsInLayer(lay);
  564. //Int_t indx0 = GetHitsInLayer(lay);
  565. //MpdKalmanHit *hitMin = 0x0;
  566. // cout << " lay, nLay, indx: " << lay << " " << nLay << " " << indx0 << endl;
  567. //Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  568. MpdItsKalmanTrack branchLay[maxBrLay];
  569. nBr0 = nBr;
  570. nBr = 0;
  571. // cout << " Point#1: " << endl;
  572. for (Int_t ibr = 0; ibr < nBr0; ++ibr) {
  573. MpdItsKalmanTrack *curTr = &branch[ibr], *branchTr = nullptr;
  574. // Check for missing layer (2 sides)
  575. Int_t lastLay = 4, frozen = 0; //VK
  576. MpdKalmanHit *h = (MpdKalmanHit*) curTr->GetHits()->Last();
  577. if (h) lastLay = h->GetLayer();
  578. if (TMath::Abs(lay - lastLay) > 2) frozen = 1;
  579. ok = nAdd = 0;
  580. //Int_t firstHit = -1, skipHit = -1;
  581. //MpdItsKalmanTrack trackBr[4] = {*curTr, *curTr, *curTr, *curTr};
  582. map<Int_t,MpdItsKalmanTrack> trackBr;
  583. multimap<Int_t,Int_t> hitsInWin;
  584. Bool_t navOK = NavigateToLayer(lay, curTr, trackBr, hitsInWin); // navigate to layer
  585. //AZ if (!navOK || trackBr.size() == 0) { ok = -1; continue; }
  586. if (!navOK) { ok = -1; continue; }
  587. map<Int_t,MpdItsKalmanTrack>::iterator it;
  588. // Debug
  589. /*
  590. cout << " Forks: " << trackBr.size() << " " << hitsInWin.size() << endl;
  591. for (it = trackBr.begin(); it != trackBr.end(); ++it) {
  592. cout << " detID " << it->first << " " << fHitSts.GetLadder(it->first%1000000) << " "
  593. << fHitSts.GetSensor(it->first%1000000) << " "
  594. << it->second.GetParamNew(0) << " " << it->second.GetParamNew(1) << endl;
  595. }
  596. */
  597. for (it = trackBr.begin(); it != trackBr.end(); ++it) {
  598. //if (frozen) break;
  599. //Double_t leng = curTr->GetLength(),
  600. Double_t step = 0;
  601. branchTr = &it->second;
  602. /*
  603. step = branchTr->GetLength() - leng;
  604. step = 0.005 * 2; // 50um
  605. if (step > 1.e-4) {
  606. // Crossing silicon layer - add mult. scat. in the sensor
  607. Double_t x0 = 9.36; // rad. length
  608. TMatrixDSym *cov = branchTr->Weight2Cov();
  609. Double_t th = branchTr->GetParamNew(3);
  610. Double_t cosTh = TMath::Cos(th);
  611. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(branchTr, x0, step, mass2);
  612. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  613. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  614. (*cov)(3,3) += angle2;
  615. Int_t iok = 0;
  616. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  617. branchTr->SetWeight(*cov);
  618. }
  619. */
  620. TVector3 mom3 = branchTr->Momentum3(), norm;
  621. mom3.SetMag(1.0);
  622. typedef multimap<Int_t,Int_t>::iterator mmit;
  623. pair<mmit,mmit> ret = hitsInWin.equal_range(it->first);
  624. map<Int_t,Int_t>::iterator it1;
  625. // Loop over hits in one sensor
  626. for (it1 = ret.first; it1 != ret.second; ++it1) {
  627. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(it1->second);
  628. if (it1 == ret.first) {
  629. norm = geo->Normal(hit);
  630. //step = 0.005 / TMath::Abs(norm * mom3) * 4.0; // extra factor 4. - possible overlaps
  631. step = thick[hit->GetLayer()] / TMath::Abs(norm * mom3) * 4.0; // extra factor 4. - possible overlaps
  632. }
  633. // Exclude used hits
  634. if (hit->GetFlag() & MpdKalmanHit::kUsed) continue;
  635. if (frozen) continue;
  636. // !!! Exact ID match
  637. if (fExact && TrackID(hit) != track->GetTrackID()) continue;
  638. //cout << lay << " " << hit->GetLayer() << " " << branchTr->GetParamNew(0) << " " << hit->GetMeas(0) << " "
  639. // << branchTr->GetParamNew(1) << " " << hit->GetMeas(1) << endl;
  640. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(branchTr,hit,pointWeight,param);
  641. //cout << lay << " " << dChi2 << " " << track->GetTrackID() << " " << TrackID(hit) << endl;
  642. if (TMath::Abs(dChi2) < fgkChi2Cut) {
  643. //if (dChi2 < fgkChi2Cut && lay % 2 == 0) {
  644. if (nBr < maxBrLay) {
  645. branchLay[nBr] = *branchTr;
  646. branchLay[nBr].GetHits()->Add(hit);
  647. branchLay[nBr].SetChi2(branchTr->GetChi2()+dChi2);
  648. TMatrixDSym w = *branchTr->GetWeight();
  649. w += pointWeight;
  650. branchLay[nBr].SetWeight(w);
  651. branchLay[nBr].SetPosNew(branchTr->GetPosNew());
  652. //if (hit->GetType() != MpdKalmanHit::kFixedP) branchLay[nBr].SetPosNew(curTr->GetPosNew());
  653. //else branchLay[nBr].SetPosNew(posNew); // modular geometry
  654. branchLay[nBr].SetLength(branchTr->GetLength());
  655. branchLay[nBr].SetParamNew(param);
  656. // Save track params at last hit
  657. /*AZ 5-jan-2018
  658. branchLay[nBr].SetLengAtHit(branchTr->GetLength());
  659. branchLay[nBr].SetParamAtHit(param);
  660. branchLay[nBr].SetWeightAtHit(*branchLay[nBr].GetWeight());
  661. */
  662. // Add multiple scattering in the sensor
  663. //*
  664. Double_t x0 = 9.36; // rad. length
  665. TMatrixDSym *cov = branchLay[nBr].Weight2Cov();
  666. Double_t th = branchLay[nBr].GetParamNew(3);
  667. Double_t cosTh = TMath::Cos(th);
  668. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(&branchLay[nBr], x0, step, mass2);
  669. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  670. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  671. (*cov)(3,3) += angle2;
  672. Int_t iok = 0;
  673. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  674. branchLay[nBr].SetWeight(*cov);
  675. //*/
  676. ++nBr;
  677. ++nAdd;
  678. }
  679. }
  680. } // for (it1 = ret.first; it1 != ret.second; - loop over hits in one sensor
  681. } // for (it = trackBr.begin(); - loop over forks
  682. //cout << " Point#3: " << endl;
  683. // Add branch with missing layer
  684. if (nAdd == 0 && nBr < maxBrLay && ok != -1) {
  685. hitTmp.SetType(MpdKalmanHit::kFixedR);
  686. hitTmp.SetPos(fRad[lay]);
  687. if (!MpdKalmanFilter::Instance()->PropagateToHit(curTr,&hitTmp,kTRUE,kTRUE)) { }
  688. // Add multiple scattering in the sensor
  689. //Double_t x0 = 9.36, step = 0.005 * 4.0; // rad. length
  690. Double_t x0 = 9.36, step = thick[lay] * 4.0; // rad. length
  691. TMatrixDSym *cov = curTr->Weight2Cov();
  692. Double_t th = curTr->GetParamNew(3);
  693. Double_t cosTh = TMath::Cos(th);
  694. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(curTr, x0, step, mass2);
  695. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  696. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  697. (*cov)(3,3) += angle2;
  698. Int_t iok = 0;
  699. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  700. curTr->SetWeight(*cov);
  701. trackBr[0] = *curTr;
  702. it = trackBr.begin();
  703. branchLay[nBr] = it->second;
  704. ++it;
  705. // Take branch with greater length
  706. for ( ; it != trackBr.end(); ++it) {
  707. if (it->second.GetLength() > branchLay[nBr].GetLength()) branchLay[nBr] = it->second;
  708. }
  709. ++nBr;
  710. }
  711. } // for (Int_t ibr = 0; ibr < nBr0;
  712. maxInLay = TMath::Max (maxInLay,nBr);
  713. if (nBr == 0) {
  714. if (ok == -1) {
  715. //cout << track->GetTrackID() << " " << lay << endl;
  716. //return ok; // stop track
  717. //cout << branch[0].GetNode() << ", " << branch[0].GetNodeNew() << endl;
  718. //cout << " Break: layer = " << lay << " node = " << branch[0].GetNode() << endl;
  719. break;
  720. }
  721. cout << " !!! MpdTrackFinderIts5spd::RunKalmanFilter: it can not happen |" << endl;
  722. exit(0);
  723. } else if (nBr <= maxBr) {
  724. for (Int_t i = 0; i < nBr; ++i) {
  725. branch[i].Reset();
  726. branch[i] = branchLay[i];
  727. }
  728. } else {
  729. // Too many branches - take the best ones
  730. for (Int_t i = 0; i < nBr; ++i) {
  731. Double_t c2 = TMath::Min (branchLay[i].GetChi2(),200.);
  732. quality[i] = branchLay[i].GetNofHits() + (0.999-c2/201);
  733. }
  734. TMath::Sort(nBr,quality,indx_lay);
  735. for (Int_t i = 0; i < maxBr; ++i) {
  736. branch[i].Reset();
  737. branch[i] = branchLay[indx_lay[i]];
  738. }
  739. nBr = maxBr;
  740. }
  741. } // for (Int_t lay = layOK-1; lay >= 0;
  742. // Select the best branch
  743. //cout << " Branches: " << nBr << " " << maxInLay << endl;
  744. Int_t ibest = 0;
  745. for (Int_t i = 1; i < nBr; ++i) {
  746. if (branch[i].GetNofHits() > branch[ibest].GetNofHits()) ibest = i;
  747. else if (branch[i].GetNofHits() == branch[ibest].GetNofHits() &&
  748. branch[i].GetChi2() < branch[ibest].GetChi2()) ibest = i;
  749. }
  750. track->Reset();
  751. *track = branch[ibest];
  752. return 0;
  753. }
  754. //__________________________________________________________________________
  755. Bool_t MpdTrackFinderIts5spd::NavigateToLayer(Int_t lay, MpdItsKalmanTrack *curTr, std::map<Int_t,MpdItsKalmanTrack> &trackBr,
  756. std::multimap<Int_t,Int_t>& hitsInWin)
  757. {
  758. ///< Navigate track to layer in modular geometry
  759. MpdItsKalmanTrack origTr(*curTr); // original track
  760. Double_t posNew = curTr->GetPosNew();
  761. TMatrixD parNew = *curTr->GetParamNew();
  762. TString nodeNew = curTr->GetNodeNew();
  763. // cout << " Nodes: " << nodeNew << " " << curTr->GetNode() << " Positions: " << posNew << " " << curTr->GetPos() << endl;
  764. // curTr->GetParamNew()->Print();
  765. // curTr->GetParam()->Print();
  766. TString curPath = curTr->GetNode();
  767. MpdKalmanHit hitTmp;
  768. hitTmp.SetType(MpdKalmanHit::kFixedR);
  769. hitTmp.SetPos(fRad[lay]);
  770. if (!MpdKalmanFilter::Instance()->PropagateParamR(curTr,&hitTmp,kFALSE)) {
  771. curTr->SetPosNew(posNew);
  772. curTr->SetParamNew(parNew);
  773. curTr->SetNodeNew(nodeNew);
  774. curTr->SetNode(curPath);
  775. //ok = -1;
  776. //break;
  777. return kFALSE;
  778. }
  779. // Find hits near the extrapolated track position
  780. //Double_t window = 1.0; // +- 1 cm window around extrapolated track position
  781. Double_t window = 2.0; // +- 1 cm window around extrapolated track position
  782. Double_t rphiTr = curTr->GetParamNew(0), zTr = curTr->GetParamNew(1);
  783. set<Int_t> indsRphi, indsZ;
  784. multimap<Double_t,Int_t>::iterator mitb = fHitMapRphi[lay].lower_bound(rphiTr-window), mit1;
  785. multimap<Double_t,Int_t>::iterator mite = fHitMapRphi[lay].upper_bound(rphiTr+window);
  786. for (mit1 = mitb; mit1 != mite; ++mit1) indsRphi.insert(mit1->second);
  787. mitb = fHitMapZ[lay].lower_bound(zTr-window);
  788. mite = fHitMapZ[lay].upper_bound(zTr+window);
  789. for (mit1 = mitb; mit1 != mite; ++mit1) indsZ.insert(mit1->second);
  790. vector<Int_t> indv(indsRphi.size()+indsZ.size());
  791. vector<Int_t>::iterator vit = set_intersection (indsRphi.begin(), indsRphi.end(), indsZ.begin(), indsZ.end(), indv.begin());
  792. indv.resize(vit-indv.begin());
  793. Int_t nWin = indv.size();
  794. //cout << " Hits in window: " << nWin << " " << curTr->GetTrackID() << endl;
  795. if (nWin == 0) {
  796. //cout << " !!! Navigate: No hits in window " << endl;
  797. //return kFALSE; // !!! FIXME - handle missing hit
  798. return kTRUE; // !!! FIXME - handle missing hit
  799. //exit(0);
  800. }
  801. for (Int_t j = 0; j < nWin; ++j) {
  802. MpdKalmanHit *h = (MpdKalmanHit*) fKHits->UncheckedAt(indv[j]);
  803. //cout << h->GetDetectorID() << " " << endl;
  804. if (trackBr.find(h->GetDetectorID()) != trackBr.end()) {
  805. // Sensor has already been found
  806. hitsInWin.insert(pair<Int_t,Int_t>(h->GetDetectorID(),indv[j]));
  807. continue;
  808. }
  809. MpdItsKalmanTrack tr(origTr);
  810. if (!MpdKalmanFilter::Instance()->PropagateToHit(&tr,h,kTRUE,kTRUE)) {
  811. //cout << " !!! Navigate: failed propagation " << endl;
  812. continue;
  813. }
  814. // Fork
  815. trackBr[h->GetDetectorID()] = tr;
  816. hitsInWin.insert(pair<Int_t,Int_t>(h->GetDetectorID(),indv[j]));
  817. }
  818. if (trackBr.size() > 14) {
  819. cout << "!!! Navigate: Too many hits in window " << endl;
  820. //exit(0);
  821. }
  822. return kTRUE;
  823. }
  824. //__________________________________________________________________________
  825. Int_t MpdTrackFinderIts5spd::TrackID(MpdKalmanHit *hit)
  826. {
  827. /// Return track ID of the hit
  828. FairHit *h = (FairHit*) fItsHits->UncheckedAt(hit->GetIndex());
  829. return ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  830. }
  831. //__________________________________________________________________________
  832. TVector2 MpdTrackFinderIts5spd::GetDistance(MpdKalmanTrack *track, MpdKalmanHit *hit)
  833. {
  834. /// Compute distance between track and hit
  835. Int_t lay = hit->GetLayer();
  836. Int_t lay2 = lay / 2;
  837. //Int_t side = lay % 2;
  838. Int_t module0 = ((MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex()))->Module();
  839. Double_t zTr = track->GetParamNew(1);
  840. Double_t zloc = zTr + fDz[lay2];
  841. Int_t modTr = Int_t (zloc / fZmod[lay2]);
  842. Int_t dMod = modTr - module0;
  843. Double_t dZ = 0;
  844. if (dMod != 0) {
  845. // Not in the same module0 - check Z-distance to module0 edge
  846. if (dMod > 0) dZ = zloc - (module0+1) * fZmod[lay2];
  847. else dZ = zloc - module0 * fZmod[lay2];
  848. if (TMath::Abs(dMod) > 2) return TVector2(TMath::Abs(dZ),999.); // not in neighbour modules
  849. }
  850. // Translate transverse track position to local system
  851. Double_t xloc = track->GetParamNew(0) * hit->GetCosSin(0) + zTr * hit->GetCosSin(1);
  852. Double_t phTr = xloc / track->GetPosNew();
  853. Double_t phHit = hit->GetMeas(0) / fRad[lay];
  854. //TVector2 dist(TMath::Abs(dZ), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit)-phTr));
  855. TVector2 dist(TMath::Abs(dZ),
  856. TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit,hit->GetCosSin(0))-phTr));
  857. //TVector2 dist(TMath::Abs(zTr-hit->GetZ()), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(track,hit)/hit->GetR()-ph));
  858. return dist;
  859. }
  860. //__________________________________________________________________________
  861. void MpdTrackFinderIts5spd::Write()
  862. {
  863. /// Write
  864. TFile histoFile("ItsRec.root","RECREATE");
  865. Writedir2current(fHistoDir);
  866. histoFile.Close();
  867. }
  868. //__________________________________________________________________________
  869. void MpdTrackFinderIts5spd::Writedir2current( TObject *obj )
  870. {
  871. /// Write
  872. if( !obj->IsFolder() ) obj->Write();
  873. else{
  874. TDirectory *cur = gDirectory;
  875. TDirectory *sub = cur->mkdir(obj->GetName());
  876. sub->cd();
  877. TList *listSub = ((TDirectory*)obj)->GetList();
  878. TIter it(listSub);
  879. while( TObject *obj1=it() ) Writedir2current(obj1);
  880. cur->cd();
  881. }
  882. }
  883. //__________________________________________________________________________
  884. void MpdTrackFinderIts5spd::AddHits()
  885. {
  886. /// Add hit objects to tracks and compute number of wrongly assigned hits
  887. /// (hits with ID different from ID of starting TPC track)
  888. Int_t nFound = fTracks->GetEntriesFast();
  889. for (Int_t i = 0; i < nFound; ++i) {
  890. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);
  891. // cout << track->GetNode() << " " << track->GetNodeNew() << endl;
  892. Double_t c2 = track->GetChi2();
  893. track->SetChi2(track->GetChi2Its());
  894. track->SetChi2Its(c2);
  895. Int_t nHits = track->GetNofHits();
  896. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  897. TClonesArray &trHits = *track->GetTrHits();
  898. // cout << nHits << " " << trHits.GetEntriesFast() << " " << track->GetTrackID() << endl;
  899. TObjArray *hits = track->GetHits();
  900. Int_t nWrong = 0, motherID = track->GetTrackID(); // nMirr = 0,
  901. // Get track mother ID
  902. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  903. while (mctrack->GetMotherId() >= 0) {
  904. motherID = mctrack->GetMotherId();
  905. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  906. }
  907. Int_t lastIndx = trHits.GetEntriesFast();
  908. for (Int_t j = 0; j < nHits; ++j) {
  909. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  910. hit->SetUniqueID(1); // flag ITS hits
  911. new (trHits[lastIndx+j]) MpdKalmanHit(*hit);
  912. // cout << " " << hit->GetLayer();
  913. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex());
  914. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  915. // cout << "-" << motherID1;
  916. // Get point mother ID
  917. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID1);
  918. while (mctrack->GetMotherId() >= 0) {
  919. motherID1 = mctrack->GetMotherId();
  920. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  921. }
  922. if (motherID1 != motherID) ++nWrong;
  923. }
  924. // if (nHits) cout << "\n" << nWrong << endl;
  925. track->SetNofWrong(nWrong);
  926. //MpdKalmanTrack *tpc = (MpdKalmanTrack*) fTpcTracks->UncheckedAt(track->GetUniqueID()-1);
  927. //track->SetChi2(track->GetChi2()+tpc->GetChi2());
  928. //track->SetLastLay();
  929. //track->GetParam()->Print();
  930. track->SetNofHits(track->GetNofTrHits()); // TPC and ITS hits
  931. track->SetNofIts(nHits);
  932. // cout << nHits << " " << track->GetNofTrHits() << " " << track->GetTrackID() << " " << track->GetChi2Its() << " " << track->GetChi2() << endl;
  933. // cout << " " << endl;
  934. }
  935. fTracks->Compress();
  936. }
  937. //__________________________________________________________________________
  938. Bool_t MpdTrackFinderIts5spd::Refit(MpdItsKalmanTrack *track, Double_t mass, Int_t charge)
  939. {
  940. /// Refit track in ITS+TPC using track hits (toward beam line) for some
  941. /// particle mass and charge hypothesis
  942. if (fTpcKF == nullptr) fTpcKF = (MpdTpcKalmanFilter*) FairRun::Instance()->GetTask("TPC Kalman filter");
  943. if (fTpcKF == nullptr) {
  944. cout << " !!! TPC Kalman Filter was not activated !!! " << endl;
  945. exit(1);
  946. }
  947. MpdKalmanGeoScheme *geoScheme = MpdKalmanFilter::Instance()->GetGeo();
  948. //MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(track->GetUniqueID()-1);
  949. track->SetUniqueID(0); // reset ITS flag
  950. //tr.GetWeightAtHit()->Print();
  951. //tr.GetParamAtHit()->Print();
  952. //AZ fTpcKF->Refit(track, mass, charge); // this mixes up hit order (ITS and TPC) !!!
  953. if (track->GetNofTrHits() > track->GetNofIts()) fTpcKF->Refit(track, mass, charge); // this mixes up hit order (ITS and TPC) !!!
  954. //tr.GetWeight()->Print();
  955. //tr.GetParam()->Print();
  956. if (track->GetNofIts() == 0) return kTRUE;
  957. MpdKalmanHit hTmp;
  958. hTmp.SetType(MpdKalmanHit::kFixedR);
  959. hTmp.SetPos(track->GetPos());
  960. track->SetParamNew(*track->GetParam());
  961. track->SetPos(track->GetPosNew());
  962. track->ReSetWeight();
  963. //TMatrixDSym w = *tr.GetWeight(); // save current weight matrix
  964. //MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  965. MpdKalmanFilter::Instance()->PropagateParamToHit(track,&hTmp,kFALSE);
  966. //tr.SetWeight(w); // restore original weight matrix (near TPC inner shell)
  967. track->SetDirection(MpdKalmanTrack::kInward);
  968. track->SetUniqueID(1);
  969. Int_t ntot = track->GetNofTrHits(), nhits = TMath::Min(15,ntot);
  970. Int_t indx0 = ntot - nhits;
  971. TMatrixD param(5,1);
  972. TMatrixDSym weight(5), pointWeight(5);
  973. TString mass2 = "";
  974. mass2 += mass * mass;
  975. Double_t vert[3] = {0.0,0.0,0.0};
  976. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07};
  977. for (Int_t ih = indx0; ih < ntot; ++ih) {
  978. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetTrHits()->UncheckedAt(ih);
  979. if (hit->GetUniqueID() == 0) continue; // skip TPC hits
  980. MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kFALSE, kTRUE); // do not adjust track length for the moment
  981. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hit, pointWeight, param);
  982. //cout << ih << " " << hit->GetPos() << " " << hit->GetUniqueID() << " " << dChi2 << "\n";
  983. track->SetChi2(track->GetChi2()+dChi2);
  984. weight = *track->GetWeight();
  985. weight += pointWeight;
  986. track->SetWeight(weight);
  987. track->SetParamNew(param);
  988. // Add multiple scattering in the sensor
  989. //*
  990. Double_t x0 = 9.36; //, step = 0.005 * 4.0; // rad. length
  991. TVector3 mom3 = track->Momentum3();
  992. mom3.SetMag(1.0);
  993. TVector3 norm = geoScheme->Normal(hit);
  994. Double_t cosPhi = norm.Dot(mom3);
  995. TMatrixDSym *cov = track->Weight2Cov();
  996. Double_t th = track->GetParamNew(3);
  997. Double_t cosTh = TMath::Cos(th);
  998. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, thick[lay]/cosPhistep, mass2);
  999. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, 4*thick[hit->GetLayer()]/cosPhi, mass2);
  1000. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1001. (*cov)(3,3) += angle2;
  1002. Int_t iok = 0;
  1003. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1004. track->SetWeight(*cov);
  1005. //*/
  1006. }
  1007. //------ Propagate track to the beam line
  1008. track->SetParam(*track->GetParamNew());
  1009. track->SetPos(track->GetPosNew());
  1010. Double_t pos = track->GetPos();
  1011. TMatrixD par = *track->GetParam();
  1012. TMatrixDSym cov = *track->Weight2Cov();
  1013. Double_t leng = track->GetLength();
  1014. TString nodeNew = track->GetNodeNew();
  1015. //cout << " 1: " << nodeNew << ", " << track->GetNode() << " " << track->GetHits()->GetEntriesFast() << endl;
  1016. // Go to beam pipe
  1017. hTmp.SetPos(fPipeR);
  1018. //hTmp.SetPos(2.9); // !!! FIXME
  1019. Int_t iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hTmp, kFALSE);
  1020. if (iok != 1) {
  1021. // Restore track
  1022. track->SetParam(par);
  1023. track->SetParamNew(par);
  1024. track->SetCovariance(cov);
  1025. track->ReSetWeight();
  1026. track->SetPos(pos);
  1027. track->SetPosNew(pos);
  1028. track->SetLength(leng);
  1029. //track->SetNode(node);
  1030. track->SetNodeNew(nodeNew);
  1031. } else {
  1032. // Add multiple scattering
  1033. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  1034. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  1035. TMatrixDSym* pcov = track->Weight2Cov();
  1036. Double_t th = track->GetParamNew(3);
  1037. Double_t cosTh = TMath::Cos(th);
  1038. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  1039. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  1040. (*pcov)(3,3) += angle2;
  1041. Int_t ok = 0;
  1042. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  1043. track->SetWeight(*pcov);
  1044. }
  1045. track->Weight2Cov(); // rebuild covar matrix (it was overwritten by weight above)
  1046. hTmp.SetPos(0.);
  1047. hTmp.SetMeas(0,track->GetParam(2)); // track Phi
  1048. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hTmp, kFALSE);
  1049. if (iok != 1) MpdKalmanFilter::Instance()->FindPca(track, vert);
  1050. else track->SetNodeNew(""); //??? 3-jan-2018
  1051. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  1052. return kTRUE;
  1053. }
  1054. //__________________________________________________________________________
  1055. ClassImp(MpdTrackFinderIts5spd);