MpdTrackFinderIts.cxx 57 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530
  1. // -------------------------------------------------------------------------
  2. // ----- MpdTrackFinderIts source file -----
  3. // ----- Created 21/07/09 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdTrackFinderIts.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD Inner Tracking System (ITS) using seeds from TPC
  9. **/
  10. #include "MpdStsGeoPar.h"
  11. #include "MpdTrackFinderIts.h"
  12. #include "MpdKalmanFilter.h"
  13. #include "MpdKalmanGeoScheme.h"
  14. #include "MpdKalmanHit.h"
  15. #include "MpdKalmanTrack.h"
  16. //#include "MpdKalmanStripHit.h"
  17. #include "MpdStsHit.h"
  18. #include "MpdStsPoint.h"
  19. #include "MpdItsKalmanTrack.h"
  20. #include "MpdTpcKalmanTrack.h"
  21. #include "FairGeoNode.h"
  22. #include "FairMCPoint.h"
  23. #include "MpdMCTrack.h"
  24. #include "FairRootManager.h"
  25. #include "FairRun.h"
  26. #include "FairRunAna.h"
  27. #include "FairRuntimeDb.h"
  28. #include "TGeoMatrix.h"
  29. #include "TGeoBBox.h"
  30. #include "TGeoNode.h"
  31. #include "TGeoTube.h"
  32. #include "TGeoManager.h"
  33. #include "TMath.h"
  34. //#include "TFile.h"
  35. //#include "TLorentzVector.h"
  36. #include "TVector2.h"
  37. #include "TClonesArray.h"
  38. #include <TRandom.h>
  39. #include <iostream>
  40. #include <map>
  41. using std::cout;
  42. using std::endl;
  43. //using std::map;
  44. const Double_t MpdTrackFinderIts::fgkChi2Cut = 20; //20; //100;
  45. //const Double_t MpdTrackFinderIts::fgkChi2Cut = 50; // for track fitting !!!
  46. //__________________________________________________________________________
  47. MpdTrackFinderIts::MpdTrackFinderIts(const char *name, Int_t iVerbose )
  48. :FairTask(name, iVerbose),
  49. fGeo(0)
  50. {
  51. //fKHits = new TClonesArray("MpdKalmanStripHit", 100);
  52. fKHits = new TClonesArray("MpdKalmanHit", 100);
  53. fTracks = new TClonesArray("MpdItsKalmanTrack", 100);
  54. fHistoDir = 0x0;
  55. fhLays = new TH1F("hLaysITS","ITS layers",10,0,10);
  56. fLayPointers = 0x0;
  57. }
  58. //__________________________________________________________________________
  59. MpdTrackFinderIts::~MpdTrackFinderIts()
  60. {
  61. //delete fKHits;
  62. //delete fTracks;
  63. //delete fTrackCand;
  64. delete [] fLayPointers;
  65. delete fhLays;
  66. }
  67. //__________________________________________________________________________
  68. InitStatus MpdTrackFinderIts::Init()
  69. {
  70. //return ReInit();
  71. if (ReInit() != kSUCCESS) return kERROR;
  72. // Read database
  73. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  74. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  75. //cout << geoPar << endl;
  76. TString volName = "sts01 ", path = "";
  77. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  78. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  79. Int_t nLay = 5;
  80. Double_t size = 6.2;
  81. for (Int_t i = 0; i < nLay; ++i) {
  82. volName[5] = 97 + i; // 'a', 'b', ..
  83. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->FindObject(volName));
  84. if (sensVol == 0x0) {
  85. if (i == 0) fGeo = 1; // modular geometry
  86. break;
  87. }
  88. TArrayD* params = sensVol->getParameters();
  89. fRad[2*i] = params->At(0);
  90. fRad[2*i+1] = params->At(1);
  91. //dR = (rMax-rMin) / 50;
  92. fDz[i] = params->At(2);
  93. Int_t nMods = Int_t (fDz[i] * 2. / size + 0.1);
  94. fZmod[i] = fDz[i] * 2. / nMods;
  95. cout << " *** STS sensitive volume: " << sensVol->GetName() << " " << params->At(0)
  96. << " " << fDz[i] << " " << fZmod[i] << endl;
  97. }
  98. //fStereoA[0] = -7.5;
  99. //fStereoA[1] = 7.5;
  100. fStereoA[0] = 7.5;
  101. fStereoA[1] = 0.0;
  102. //for (Int_t i = 0; i < 2; ++i) fStereoA[i] =
  103. //TMath::Tan (fStereoA[i]*TMath::DegToRad());
  104. for (Int_t i = 0; i < 2; ++i) fStereoA[i] *= TMath::DegToRad();
  105. if (fGeo) {
  106. // Process modular geometry
  107. Double_t safety = 0.03;
  108. for (Int_t i = 0; i < nLay; ++i) {
  109. fNladders[2*i] = fNladders[2*i+1] = fNsectors[2*i] = fNsectors[2*i+1] = 0;
  110. volName = "sts01ladder";
  111. volName += (i+1);
  112. path = "/cave_1/sts01_0/" + volName;
  113. path += "_";
  114. TString path0 = path;
  115. //volName = "/cave_1/sts01_0/sts01ladder";
  116. //volName += (i+1);
  117. // Loop over all ladders to find the one with the smallest radius
  118. fRad[2*i+1] = fRad[2*i] = 999.9;
  119. for (Int_t jlad = 1; jlad < 99; ++jlad) {
  120. path = path0;
  121. path += jlad;
  122. if (!gGeoManager->CheckPath(path)) break;
  123. gGeoManager->cd(path);
  124. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  125. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  126. Double_t xyzL[3] = {0}, xyzM[3];
  127. gGeoManager->LocalToMaster(xyzL,xyzM);
  128. Double_t rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  129. fRad[2*i] = TMath::Min (fRad[2*i],rad);
  130. xyzL[0] = box->GetDX();
  131. gGeoManager->LocalToMaster(xyzL,xyzM);
  132. rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  133. fRad[2*i] = TMath::Min (fRad[2*i],rad);
  134. xyzL[0] = -box->GetDX();
  135. gGeoManager->LocalToMaster(xyzL,xyzM);
  136. rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  137. fRad[2*i] = TMath::Min (fRad[2*i],rad);
  138. }
  139. fRad[2*i+1] = fRad[2*i];
  140. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  141. if (ladd == NULL) { nLay = i; break; }
  142. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  143. //safety = -box->GetDY();
  144. //safety = box->GetDY();
  145. safety = 2 * box->GetDY(); // new
  146. fRad[2*i] -= safety;
  147. fRad[2*i+1] -= safety;
  148. //new if (i == 0) { fRad[2*i] -= safety; fRad[2*i+1] -= safety; }
  149. }
  150. FillGeoScheme();
  151. }
  152. // Get pipe radius
  153. fPipeR = ((TGeoTube*)gGeoManager->GetVolume("pipe1")->GetShape())->GetRmax();
  154. // Get cables
  155. TObjArray *vols = gGeoManager->GetListOfVolumes();
  156. Int_t nvols = vols->GetEntriesFast(), ncables = 0;
  157. for (Int_t i = 0; i < nvols; ++i) {
  158. TGeoVolume *vol = (TGeoVolume*) vols->UncheckedAt(i);
  159. TString cable = TString(vol->GetName());
  160. if (!(cable.Contains("sts") && cable.Contains("cable"))) continue;
  161. //cout << cable << endl;
  162. ++ncables;
  163. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  164. TString lad = cable;
  165. Int_t ip = lad.Index("cable");
  166. lad.Replace(ip,lad.Length()-ip+1,"");
  167. lad.Replace(lad.Length()-2,1,"");
  168. lad += "_1/";
  169. path = "/cave_1/sts01_0/" + lad + cable;
  170. path += "_1";
  171. gGeoManager->cd(path);
  172. Double_t xyzL[3] = {0}, xyzM[3];
  173. gGeoManager->LocalToMaster(xyzL,xyzM);
  174. //cout << xyzM[2] - box->GetDZ() << " " << xyzM[2] + box->GetDZ() << " " << box->GetDZ() << endl;
  175. if (xyzM[2] - box->GetDZ() > 0) {
  176. Int_t lay = TString(lad(lad.Length()-4,1)).Atoi();
  177. fCables[lay-1].insert(pair<Double_t,Double_t>(xyzM[2] - box->GetDZ(),xyzM[2] + box->GetDZ()));
  178. }
  179. }
  180. return kSUCCESS;
  181. }
  182. //__________________________________________________________________________
  183. void MpdTrackFinderIts::FillGeoScheme()
  184. {
  185. /// Fill Kalman filter geometry manager info
  186. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  187. TGeoVolume *vSts = gGeoManager->GetVolume("sts01");
  188. for (Int_t layer = 1; layer < 999; ++layer) {
  189. // Loop over layers
  190. TString sladder = "sts01ladder";
  191. sladder += layer;
  192. TGeoVolume *vLay = gGeoManager->GetVolume(sladder);
  193. if (vLay == 0x0 && layer == 1) continue; // no first layer
  194. if (vLay == 0x0) break;
  195. sladder += "_";
  196. //TString sdet = "sts01detector";
  197. TString sdet = "sts01sensor";
  198. sdet += layer;
  199. sdet += "_";
  200. //TString sensor = "sts01sensor";
  201. TString sensor = "";
  202. //sensor += layer;
  203. //sensor += "_0";
  204. for (Int_t ladder = 1; ladder < 999; ++ladder) {
  205. // Loop over ladders
  206. TString sladder1 = sladder;
  207. sladder1 += ladder;
  208. TGeoNode *node = vSts->FindNode(sladder1);
  209. if (node == 0x0) break;
  210. ++fNladders[2*(layer-1)];
  211. ++fNladders[2*(layer-1)+1];
  212. TGeoVolume *vLad = node->GetVolume();
  213. //cout << vLad->GetNodes()->GetEntriesFast() << " " << vLad->GetNdaughters() << endl;
  214. Int_t nDaught = vLad->GetNdaughters(), detID = -1, detIDsts = -1;
  215. TObjArray *daught = vLad->GetNodes();
  216. //for (Int_t j = 0; j < vLad->GetNdaughters(); ++j)
  217. //cout << ((TGeoNode*)(vLad->GetNodes()->UncheckedAt(j)))->GetName() << endl;
  218. Int_t iZ = 0;
  219. for (Int_t det = 0; det < nDaught; ++det) {
  220. // Loop over ladder daughters
  221. TString sdet1 = ((TGeoNode*)(daught->UncheckedAt(det)))->GetName();
  222. if (!sdet1.Contains("sensor") && !sdet1.Contains("sector")) continue;
  223. if (ladder == 1) { ++fNsectors[2*(layer-1)]; ++fNsectors[2*(layer-1)+1]; }
  224. Int_t det1 = TString(sdet1(sdet1.Index("_")+1,2)).Atoi();
  225. Int_t secType = -1;
  226. if (sdet1.Contains("sector")) secType = TString(sdet1(sdet1.Index("_")-2,1)).Atoi();
  227. ++iZ;
  228. for (Int_t side = 0; side < 2; ++side) {
  229. detIDsts = fHitSts.SetDetId(secType, layer, ladder, det1, side);
  230. detID = fHitSts.SetDetId(0, layer, ladder, iZ, side);
  231. //detIDsts += 1000000 * ((layer-1) * 2 + side);
  232. fId2Id[(layer-1) * 2 + side].insert(pair<Int_t,Int_t>(detIDsts,detID));
  233. detID += 1000000 * ((layer-1) * 2 + side);
  234. TString detName = sladder1 + "/" + sdet1 + "#";
  235. detName += side;
  236. //if (side) geo->SetDetId(detName, detID&((2<<26)-2)); // store odd side
  237. //if (!side) geo->SetDetId(detName, detID); // store even side
  238. geo->SetDetId(detName, detID);
  239. //TString path = "/cave_1/sts01_0/" + detName + "/" + sensor;
  240. TString path = "/cave_1/sts01_0/" + detName(0,detName.Length()-2);
  241. //cout << detName << " " << path << endl;
  242. gGeoManager->cd(path);
  243. //if (side) geo->SetPath(detID&((2<<26)-2), path); // store odd side
  244. //if (!side) geo->SetPath(detID, path); // store even side
  245. geo->SetPath(detID, path);
  246. node = gGeoManager->GetCurrentNode();
  247. //cout << node->GetName() << " " << detID << endl;
  248. TGeoVolume *vol = node->GetVolume();
  249. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  250. TVector2 size(box->GetDX(), box->GetDZ());
  251. geo->SetSize(detID, size);
  252. Double_t xyzL[3] = {0}, xyzM[3], vecM[3];
  253. xyzL[1] = 1;
  254. gGeoManager->LocalToMasterVect(xyzL,vecM);
  255. //cout << vecM[0] << " " << vecM[1] << " " << vecM[2] << endl;
  256. //TVector3 norm(vecM[0], vecM[1], vecM[2]);
  257. //geo->SetNormal(detID, norm);
  258. xyzL[1] = 0;
  259. gGeoManager->LocalToMaster(xyzL,xyzM);
  260. Double_t sign = TMath::Sign (1.,xyzM[0]*vecM[0]+xyzM[1]*vecM[1]);
  261. if (detID % 2) xyzL[1] = 0.015 * sign;
  262. else xyzL[1] = -0.015 * sign;
  263. gGeoManager->LocalToMaster(xyzL,xyzM);
  264. //cout << xyzM[0] << " " << xyzM[1] << " " << xyzM[2] << endl;
  265. TVector3 pos(xyzM[0], xyzM[1], xyzM[2]);
  266. geo->SetGlobalPos(detID, pos);
  267. //xyzL[1] = sign;
  268. //gGeoManager->LocalToMasterVect(xyzL,vecM);
  269. //cout << vecM[0] << " " << vecM[1] << " " << vecM[2] << endl;
  270. TVector3 norm(sign*vecM[0], sign*vecM[1], sign*vecM[2]);
  271. geo->SetNormal(detID, norm);
  272. }
  273. }
  274. }
  275. }
  276. }
  277. //__________________________________________________________________________
  278. InitStatus MpdTrackFinderIts::ReInit()
  279. {
  280. fItsPoints = (TClonesArray *) FairRootManager::Instance()->GetObject("StsPoint");
  281. fItsHits =(TClonesArray *) FairRootManager::Instance()->GetObject("StsHit");
  282. if (fItsPoints == 0x0 || fItsHits == 0x0) return kERROR;
  283. fTpcTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  284. fEctTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("EctTrack");
  285. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  286. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  287. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  288. FairRootManager::Instance()->Register("ItsTrack", "Its", fTracks, kTRUE);
  289. fNPass = 1;
  290. return kSUCCESS;
  291. }
  292. //__________________________________________________________________________
  293. void MpdTrackFinderIts::Reset()
  294. {
  295. ///
  296. cout << " MpdTrackFinderIts::Reset " << endl;
  297. //fKHits->Clear();
  298. fKHits->Delete();
  299. fTracks->Delete();
  300. delete [] fLayPointers;
  301. fLayPointers = NULL;
  302. }
  303. //__________________________________________________________________________
  304. void MpdTrackFinderIts::SetParContainers()
  305. {
  306. // Get run and runtime database
  307. FairRunAna* run = FairRunAna::Instance();
  308. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  309. FairRuntimeDb* db = run->GetRuntimeDb();
  310. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  311. // Get STS geometry parameter container
  312. db->getContainer("MpdStsGeoPar");
  313. }
  314. //__________________________________________________________________________
  315. void MpdTrackFinderIts::Finish()
  316. {
  317. //Write();
  318. }
  319. //__________________________________________________________________________
  320. void MpdTrackFinderIts::Exec(Option_t * option)
  321. {
  322. static int eventCounter = 0;
  323. cout << " - - - - \n ItsRec event " << ++eventCounter << endl;
  324. Reset();
  325. // Create Kalman hits
  326. if (fItsHits->GetEntriesFast() == 0) return;
  327. MakeKalmanHits();
  328. fTracks->Clear();
  329. for (Int_t i = 0; i < fNPass; ++i) {
  330. //fTracks->Clear();
  331. GetTrackSeeds(0); // TPC tracks
  332. GetTrackSeeds(1); // ECT tracks
  333. cout << " Total number of hits for tracking: " << fKHits->GetEntriesFast() << endl;
  334. cout << " Total number of track seeds: " << fTracks->GetEntriesFast() << endl;
  335. DoTracking(i);
  336. //StoreTracks();
  337. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  338. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  339. }
  340. AddHits(); // add hit objects to tracks
  341. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  342. }
  343. //__________________________________________________________________________
  344. void MpdTrackFinderIts::MakeKalmanHits()
  345. {
  346. /// Create Kalman hits from ITS hits.
  347. fhLays->Reset();
  348. Int_t nHits = fItsHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0;
  349. Double_t r, z, xloc, errZ = 0.012, errX = 0.0023; // 120um in Z, 23um in R-Phi (local X)
  350. //Double_t r, z, xloc, errZ = 0.12, errX = 0.0023; // 1.2mm in Z, 23um in R-Phi (local X)
  351. //Double_t r, z, xloc, errZ = 50.0, errX = 0.01; // 50cm in Z, 100um in R-Phi (local X)
  352. for (Int_t ih = 0; ih < nHits; ++ih) {
  353. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(ih);
  354. r = TMath::Sqrt (h->GetX() * h->GetX() + h->GetY() * h->GetY());
  355. z = h->GetZ();
  356. xloc = h->GetLocalX();
  357. //cout << ih << " " << h->Layer()-1 << endl;
  358. //lay = h->Layer() * 2 + h->Side() + 1;
  359. lay = (h->Layer() - 1) * 2 + h->Side();
  360. // Add error
  361. Double_t dX = 0, dZ = 0;
  362. gRandom->Rannor(dX,dZ);
  363. //if (errZ > 2) dZ = 0.0; // 1-D case
  364. dZ = 0.0; // 1-D case
  365. Double_t meas[2] = {xloc+dX*errX, z+dZ*errZ};
  366. Double_t err[2] = {errX, errZ};
  367. Double_t cossin[2] = {TMath::Cos(fStereoA[h->Side()]), TMath::Sin(fStereoA[h->Side()])};
  368. //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
  369. //MpdKalmanStripHit *hit = new ((*fKHits)[nKH++]) MpdKalmanStripHit(r, fStereoA[h->Side()],
  370. // xloc+dX*errX, z+dZ*errZ, errX, errZ, lay, ih);
  371. //hit->SetType(MpdKalmanHit::kFixedR);
  372. MpdKalmanHit *hit = 0x0;
  373. //cout << h->GetDetectorID() << " " << fId2Id[lay][h->GetDetectorID()] << endl;
  374. //if (fGeo && h->GetUniqueID()) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+h->GetDetectorID(), 1,
  375. if (fGeo && h->GetUniqueID()) hit =
  376. new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+fId2Id[lay][h->GetDetectorID()], 1,
  377. MpdKalmanHit::kFixedP, meas, err, cossin, 0., r, ih);
  378. // Mask out sector number - sensor layout
  379. else if (fGeo) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+(h->GetDetectorID()&((2<<12)-1)), 1,
  380. MpdKalmanHit::kFixedP, meas, err, cossin, 0., r, ih);
  381. else hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+nKH-1, 1, MpdKalmanHit::kFixedR,
  382. meas, err, cossin, 0., r, ih);
  383. hit->SetUniqueID(0);
  384. // Add second measurement - just for test at the moment
  385. //!!!
  386. //hit->SetNofDim(2);
  387. //!!!
  388. layMax = TMath::Max (lay, layMax);
  389. fhLays->Fill(lay+0.1);
  390. }
  391. cout << " Max layer = " << layMax << ", hits: " << fKHits->GetEntriesFast() << endl;
  392. fKHits->Sort(); // in descending order in R
  393. //cout << ((MpdKalmanHit*)fKHits->UncheckedAt(0))->GetPos() << endl;
  394. cout << ((MpdKalmanHit*)fKHits->UncheckedAt(0))->GetDist() << endl;
  395. fLayPointers = new Int_t [layMax+1];
  396. Int_t ipos = 0;
  397. for (Int_t i = layMax; i > -1; --i) {
  398. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  399. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  400. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  401. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  402. if (cont == 0) {
  403. fLayPointers[i] = -1;
  404. continue;
  405. }
  406. fLayPointers[i] = ipos;
  407. ipos += cont;
  408. }
  409. }
  410. //__________________________________________________________________________
  411. void MpdTrackFinderIts::GetTrackSeeds(Int_t iPass)
  412. {
  413. /// Build ITS track seeds from TPC tracks
  414. TClonesArray *seeds = (iPass == 0) ? fTpcTracks : fEctTracks;
  415. if (seeds == NULL) return;
  416. Int_t nSeeds = seeds->GetEntriesFast();
  417. cout << " Seed tracks: " << nSeeds << endl;
  418. Int_t nCand = fTracks->GetEntriesFast();
  419. MpdKalmanHit hit;
  420. hit.SetType(MpdKalmanHit::kFixedR);
  421. for (Int_t itr = 0; itr < nSeeds; ++itr) {
  422. MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) seeds->UncheckedAt(itr);
  423. if (tpc->GetType() != MpdKalmanTrack::kBarrel) continue;
  424. //tpc->GetParam()->Print();
  425. //MpdEctKalmanTrack *track = new ((*fTracks)[nCand++]) MpdEctKalmanTrack(itr, *tpc);
  426. MpdItsKalmanTrack *track = NULL;
  427. if (TString(tpc->ClassName()).Contains("Tpc")) track = new ((*fTracks)[nCand++]) MpdItsKalmanTrack(*tpc);
  428. else track = new ((*fTracks)[nCand++]) MpdItsKalmanTrack(*((MpdEctKalmanTrack*)tpc));
  429. track->SetUniqueID(itr+1);
  430. /*
  431. track->SetParam(*track->GetParamAtHit());
  432. track->SetParamNew(*track->GetParamAtHit());
  433. track->SetPosNew(track->GetPos());
  434. track->SetWeight(*track->GetWeightAtHit());
  435. track->SetLength(track->GetLengAtHit());
  436. */
  437. hit.SetPos(track->GetPos());
  438. track->SetParamNew(*track->GetParam());
  439. track->SetPos(track->GetPosNew());
  440. track->ReSetWeight();
  441. TMatrixDSym w = *track->GetWeight(); // save current weight matrix
  442. MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  443. track->SetWeight(w); // restore original weight matrix (near TPC inner shell)
  444. //cout << nCand-1 << " " << track->GetTrackID() << endl;
  445. //cout << track->GetHits()->GetEntriesFast() << " " << track->GetTrHits()->GetEntriesFast() << endl;
  446. track->GetHits()->Clear();
  447. track->SetChi2Its(track->GetChi2()); // temporary storage
  448. track->SetChi2(0.);
  449. }
  450. cout << " Number of ITS track candidates: " << nCand << endl;
  451. }
  452. //__________________________________________________________________________
  453. void MpdTrackFinderIts::DoTracking(Int_t iPass)
  454. {
  455. /// Run Kalman tracking
  456. Double_t vert[3] = {0.0,0.0,0.0};
  457. Int_t nCand = fTracks->GetEntriesFast(), iok = 0;
  458. Int_t lay0 = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  459. for (Int_t i = 0; i < nCand; ++i) {
  460. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);
  461. //cout << " Track seed No. " << i << ", ID: " << track->GetTrackID() << ", Hits: " << track->GetNofTrHits() << endl;
  462. for (Int_t j = 0; j < track->GetNofTrHits(); ++j) {
  463. MpdKalmanHit *h = (MpdKalmanHit* )track->GetTrHits()->UncheckedAt(j);
  464. //MpdStsHit *hh = (MpdStsHit*) fItsHits->UncheckedAt(h->GetIndex());
  465. //Int_t id = ((FairMCPoint*) fItsPoints->UncheckedAt(hh->GetRefIndex()))->GetTrackID();
  466. //cout << j << " " << h->GetDist() << " " << h->GetLayer() << endl;
  467. }
  468. // Reset weight matrix - for debug
  469. /*
  470. Int_t nHits = track->GetNofTrHits();
  471. nHits *= nHits;
  472. for (Int_t ii = 0; ii < 5; ++ii) {
  473. for (Int_t j = i; j < 5; ++j) {
  474. if (j == ii) (*track->GetWeight())(ii,j) /= nHits;
  475. else (*track->GetWeight())(ii,j) = (*track->GetWeight())(j,ii) = 0.;
  476. }
  477. }
  478. */
  479. if (fGeo) iok = RunKalmanFilterMod(track, lay0); // modular geometry
  480. else iok = RunKalmanFilterCyl(track, lay0); // cylindrical geometry
  481. if (iok == -1) {
  482. fTracks->RemoveAt(i);
  483. continue;
  484. }
  485. //if (track->GetNofHits() == 0) continue; // no hits added
  486. // Propagate track to the beam line
  487. track->SetParam(*track->GetParamNew());
  488. track->SetPos(track->GetPosNew());
  489. Double_t pos = track->GetPos();
  490. TMatrixD par = *track->GetParam();
  491. TMatrixDSym cov = *track->Weight2Cov();
  492. Double_t leng = track->GetLength();
  493. TString nodeNew = track->GetNodeNew();
  494. //cout << " 1: " << nodeNew << ", " << track->GetNode() << endl;
  495. // Go to beam pipe
  496. MpdKalmanHit hit;
  497. hit.SetType(MpdKalmanHit::kFixedR);
  498. hit.SetPos(fPipeR);
  499. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  500. if (iok != 1) {
  501. // Restore track
  502. track->SetParam(par);
  503. track->SetParamNew(par);
  504. track->SetCovariance(cov);
  505. track->ReSetWeight();
  506. track->SetPos(pos);
  507. track->SetPosNew(pos);
  508. track->SetLength(leng);
  509. //track->SetNode(node);
  510. //cout << " 2: " << nodeNew << ", " << track->GetNode() << endl;
  511. track->SetNodeNew(nodeNew);
  512. } else {
  513. // Add multiple scattering
  514. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  515. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  516. TMatrixDSym* pcov = track->Weight2Cov();
  517. Double_t th = track->GetParamNew(3);
  518. Double_t cosTh = TMath::Cos(th);
  519. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  520. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  521. (*pcov)(3,3) += angle2;
  522. Int_t ok = 0;
  523. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  524. track->SetWeight(*pcov);
  525. }
  526. cov = *track->Weight2Cov();
  527. hit.SetPos(0.);
  528. hit.SetMeas(0,track->GetParam(2)); // track Phi
  529. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  530. //Double_t pos = ((MpdKalmanHit*)track->GetHits()->Last())->GetPos();
  531. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  532. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  533. if (iok != 1) MpdKalmanFilter::Instance()->FindPca(track, vert);
  534. //track->SetPos(pos); // restore position
  535. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  536. //track->GetCovariance()->Print();
  537. } // for (Int_t i = 0; i < nCand;
  538. fTracks->Compress();
  539. }
  540. //__________________________________________________________________________
  541. Int_t MpdTrackFinderIts::RunKalmanFilterCyl(MpdItsKalmanTrack *track, Int_t layBeg)
  542. {
  543. /// Run Kalman filter in the cylindrical geometry (might not work when propagating outward!!!)
  544. const Int_t maxBr = 10, maxBrLay = 1000; // max number of branches
  545. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  546. MpdKalmanHit *hitOK = 0x0;
  547. MpdKalmanHit hitTmp;
  548. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  549. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1;
  550. Int_t layEnd = -1, dLay = -1, layOK = -1;
  551. if (trackDir == MpdKalmanTrack::kOutward) {
  552. layEnd = layMax + 1;
  553. dLay = 1;
  554. }
  555. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  556. TMatrixD param(5,1), paramTmp(5,1);
  557. Double_t saveZ = 0.0, saveLeng = 0.0, dChi2Min = 0.0, quality[maxBrLay] = {0.0};
  558. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  559. MpdItsKalmanTrack branch[maxBr];
  560. Int_t nBr = 1, nBr0 = 0, indx_lay[maxBrLay] = {0}, maxInLay = 0, ok = 0, nAdd = 0;
  561. branch[0] = *track;
  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: " << lay << " " << nLay << " " << indx0 << endl;
  567. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  568. MpdItsKalmanTrack branchLay[maxBrLay];
  569. nBr0 = nBr;
  570. nBr = 0;
  571. for (Int_t ibr = 0; ibr < nBr0; ++ibr) {
  572. MpdItsKalmanTrack *curTr = &branch[ibr];
  573. // Check for missing layer (2 sides)
  574. Int_t lastLay = 8, frozen = 0;
  575. MpdKalmanHit *h = (MpdKalmanHit*) curTr->GetHits()->Last();
  576. if (h) lastLay = h->GetLayer();
  577. if (TMath::Abs(lay - lastLay) > 2) frozen = 1;
  578. ok = nAdd = 0;
  579. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  580. if (frozen) break;
  581. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  582. Double_t posNew = 0;
  583. // Propagate to hit (if it is not very close to the track)
  584. if (TMath::Abs(hit->GetPos()-curTr->GetPosNew()) > 1.e-4) {
  585. Double_t leng = curTr->GetLength();
  586. posNew = curTr->GetPosNew();
  587. TMatrixD parNew = *curTr->GetParamNew();
  588. TString nodeNew = curTr->GetNodeNew();
  589. TString curPath = curTr->GetNode();
  590. if (!MpdKalmanFilter::Instance()->PropagateToHit(curTr,hit,kTRUE,kFALSE)) {
  591. // Restore initial parameters for the failed track
  592. curTr->SetPosNew(posNew);
  593. curTr->SetParamNew(parNew);
  594. curTr->SetLength(leng);
  595. curTr->SetNodeNew(nodeNew);
  596. curTr->SetNode(curPath);
  597. ok = -1;
  598. break;
  599. }
  600. Double_t step = curTr->GetLength() - leng;
  601. //*
  602. //if (hit->GetLayer() % 2 == 10 && step > 1.e-4) {
  603. if (hit->GetLayer() % 2 == 0 && step > 1.e-4) {
  604. // Crossing silicon layer - add mult. scat. in the sensor
  605. Double_t x0 = 9.36; // rad. length
  606. TMatrixDSym *cov = curTr->Weight2Cov();
  607. Double_t th = curTr->GetParamNew(3);
  608. Double_t cosTh = TMath::Cos(th);
  609. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(curTr, x0, step);
  610. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  611. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  612. (*cov)(3,3) += angle2;
  613. Int_t iok = 0;
  614. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  615. curTr->SetWeight(*cov);
  616. } else if (hit->GetLayer() % 2 != 0 && step > 1.e-4 && fCables[lay/2].size() > 0) {
  617. // Crossing silicon layer - add mult. scat. in the cable
  618. Double_t nCables = 0, x0 = 0.0116 * 2 / 9.36; // in rad. length - 116um cable per side
  619. // Find number of cables crossed
  620. TString path = gGeoManager->GetPath();
  621. if (!path.Contains("sensor") && !path.Contains("sector")) {
  622. cout << " !!! MpdTrackFinderIts::RunKalmanFilter - Outside detector !!! " << endl;
  623. exit(0);
  624. }
  625. Double_t v7[3] = {curTr->GetParamNew(0), curTr->GetPosNew(), curTr->GetParamNew(1)}, v77[3];
  626. gGeoManager->LocalToMaster(v7,v77);
  627. Double_t zTr = TMath::Abs (v77[2]); // global Z
  628. //cout << zTr << endl;
  629. map<Double_t,Double_t>::iterator it;
  630. for (it = fCables[lay/2].begin(); it != fCables[lay/2].end(); ++it) {
  631. if (zTr < it->first || zTr > it->second) continue;
  632. ++nCables;
  633. }
  634. //cout << " Cables: " << nCables << endl;
  635. if (nCables) {
  636. x0 *= nCables;
  637. TMatrixDSym *cov = curTr->Weight2Cov();
  638. Double_t th = curTr->GetParamNew(3);
  639. Double_t cosTh = TMath::Cos(th);
  640. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(curTr, x0);
  641. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  642. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  643. (*cov)(3,3) += angle2;
  644. Int_t iok = 0;
  645. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  646. curTr->SetWeight(*cov);
  647. }
  648. }
  649. //curTr->GetWeight()->Print();
  650. }
  651. //cout << hit->GetDetectorID() << endl;
  652. // Exclude used hits
  653. if (hit->GetFlag() != 1) continue;
  654. // !!! Exact ID match
  655. if (TrackID(hit) != track->GetTrackID()) continue;
  656. Double_t dChi2 = 0;
  657. // Simple (old) geometry
  658. TVector2 dist = GetDistance(curTr, hit);
  659. //if (track->GetTrackID() != -1335) cout << /*((TpcLheKalmanTrack*)fTpcTracks->UncheckedAt(track->GetTpcIndex()))->GetNofTrHits() <<*/ " " << dist.X() << " " << dist.Y() << " " << hit->GetLayer() << endl;
  660. if (dist.X() > 15.) continue; // distance in Z
  661. if (dist.Y() > 0.25) continue; // distance in Phi
  662. hitTmp = *hit;
  663. hitTmp.SetMeas(0,MpdKalmanFilter::Instance()->Proxim(curTr,hit,hit->GetCosSin(0)));
  664. /*cout << hit->GetPos() << " " << hitTmp.GetMeas(0) << " " << hitTmp.GetMeas(1) << endl;
  665. cout << TrackID(hit) << " " << track->GetTrackID() << " " << curTr->GetPosNew() << " "
  666. << curTr->GetParamNew(0) << " " << track->GetParamNew(1) << endl; */
  667. dChi2 = MpdKalmanFilter::Instance()->FilterStrip(curTr,&hitTmp,pointWeight,param);
  668. MpdItsKalmanTrack trTmp = *curTr;
  669. trTmp.SetParamNew(param);
  670. dist = GetDistance(&trTmp, hit);
  671. //cout << " Chi2: " << dChi2 << " " << dist.X() << " " << dist.Y() << " " << curTr->GetParamNew(0) << " " << hit->GetLayer() << endl;
  672. Double_t dChi2Z = dist.X() / hit->GetErr(1); // contrib. due to Z-residual
  673. dChi2 += dChi2Z * dChi2Z;
  674. //cout << dChi2 << " " << TrackID(hit) << endl;
  675. if (dChi2 < 0) cout << " !!! Negative Chi2: !!! " << dChi2 << endl;
  676. if (TMath::Abs(dChi2) < fgkChi2Cut) {
  677. //if (dChi2 < fgkChi2Cut && lay % 2 == 0) {
  678. if (nBr < maxBrLay) {
  679. branchLay[nBr] = *curTr;
  680. branchLay[nBr].GetHits()->Add(hit);
  681. branchLay[nBr].SetChi2(curTr->GetChi2()+dChi2);
  682. TMatrixDSym w = *curTr->GetWeight();
  683. w += pointWeight;
  684. branchLay[nBr].SetWeight(w);
  685. branchLay[nBr].SetPosNew(curTr->GetPosNew());
  686. //if (hit->GetType() != MpdKalmanHit::kFixedP) branchLay[nBr].SetPosNew(curTr->GetPosNew());
  687. //else branchLay[nBr].SetPosNew(posNew); // modular geometry
  688. branchLay[nBr].SetLength(curTr->GetLength());
  689. branchLay[nBr].SetParamNew(param);
  690. // Save track params at last hit
  691. branchLay[nBr].SetLengAtHit(curTr->GetLength());
  692. branchLay[nBr].SetParamAtHit(param);
  693. branchLay[nBr].SetWeightAtHit(*branchLay[nBr].GetWeight());
  694. ++nBr;
  695. ++nAdd;
  696. }
  697. }
  698. } // for (Int_t indx = indxBeg; indx != indxEnd;
  699. // Add branch with missing layer
  700. if (nAdd == 0 && nBr < maxBrLay && ok != -1) {
  701. branchLay[nBr] = *curTr;
  702. ++nBr;
  703. }
  704. } // for (Int_t ibr = 0; ibr < nBr0;
  705. maxInLay = TMath::Max (maxInLay,nBr);
  706. if (nBr == 0) {
  707. if (ok == -1) {
  708. //cout << track->GetTrackID() << " " << lay << endl;
  709. //return ok; // stop track
  710. cout << branch[0].GetNode() << ", " << branch[0].GetNodeNew() << endl;
  711. break;
  712. }
  713. cout << " !!! MpdTrackFinderIts::RunKalmanFilter: it can not happen |" << endl;
  714. exit(0);
  715. } else if (nBr <= maxBr) {
  716. for (Int_t i = 0; i < nBr; ++i) {
  717. branch[i].Reset();
  718. branch[i] = branchLay[i];
  719. }
  720. } else {
  721. // Too many branches - take the best ones
  722. for (Int_t i = 0; i < nBr; ++i) {
  723. Double_t c2 = TMath::Min (branchLay[i].GetChi2(),200.);
  724. quality[i] = branchLay[i].GetNofHits() + (0.999-c2/201);
  725. }
  726. TMath::Sort(nBr,quality,indx_lay);
  727. for (Int_t i = 0; i < maxBr; ++i) {
  728. branch[i].Reset();
  729. branch[i] = branchLay[indx_lay[i]];
  730. }
  731. nBr = maxBr;
  732. }
  733. } // for (Int_t lay = layOK-1; lay >= 0;
  734. // Select the best branch
  735. //cout << " Branches: " << nBr << " " << maxInLay << endl;
  736. Int_t ibest = 0;
  737. for (Int_t i = 1; i < nBr; ++i) {
  738. if (branch[i].GetNofHits() > branch[ibest].GetNofHits()) ibest = i;
  739. else if (branch[i].GetNofHits() == branch[ibest].GetNofHits() &&
  740. branch[i].GetChi2() < branch[ibest].GetChi2()) ibest = i;
  741. }
  742. track->Reset();
  743. *track = branch[ibest];
  744. return 0;
  745. }
  746. //__________________________________________________________________________
  747. Int_t MpdTrackFinderIts::RunKalmanFilterMod(MpdItsKalmanTrack *track, Int_t layBeg)
  748. {
  749. /// Run Kalman filter in the modular geometry (might not work when propagating outward!!!)
  750. const Int_t maxBr = 20, maxBrLay = 1000; // max number of branches
  751. //const Int_t maxBr = 100, maxBrLay = 1000; // max number of branches
  752. //cout << fHits->GetEntriesFast() << endl;
  753. //Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  754. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  755. MpdKalmanHit *hitOK = 0x0;
  756. MpdKalmanHit hitTmp;
  757. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  758. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1;
  759. Int_t layEnd = -1, dLay = -1, layOK = -1;
  760. if (trackDir == MpdKalmanTrack::kOutward) {
  761. layEnd = layMax + 1;
  762. dLay = 1;
  763. }
  764. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  765. TMatrixD param(5,1), paramTmp(5,1);
  766. Double_t saveZ = 0.0, saveLeng = 0.0, dChi2Min = 0.0, quality[maxBrLay] = {0.0}, posNew = 0.0;
  767. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  768. MpdItsKalmanTrack branch[maxBr];
  769. Int_t nBr = 1, nBr0 = 0, indx_lay[maxBrLay] = {0}, maxInLay = 0, ok = 0, nAdd = 0;
  770. branch[0] = *track;
  771. TString mass2 = "0.0194797849"; // pion mass squared
  772. if (fMCTracks) {
  773. // Get particle mass - ideal PID
  774. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(track->GetTrackID());
  775. TParticlePDG *pdgP = TDatabasePDG::Instance()->GetParticle(mctrack->GetPdgCode());
  776. if (pdgP) {
  777. Double_t mass = pdgP->Mass();
  778. if (mass < 0.1 || mass > 0.8) {
  779. // Electrons or protons (or heavier than protons)
  780. mass2 = "";
  781. mass2 += mass*mass;
  782. }
  783. }
  784. }
  785. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  786. Int_t nLay = GetNofHitsInLayer(lay);
  787. Int_t indx0 = GetHitsInLayer(lay);
  788. MpdKalmanHit *hitMin = 0x0;
  789. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  790. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  791. MpdItsKalmanTrack branchLay[maxBrLay];
  792. nBr0 = nBr;
  793. nBr = 0;
  794. for (Int_t ibr = 0; ibr < nBr0; ++ibr) {
  795. MpdItsKalmanTrack *curTr = &branch[ibr], *branchTr = NULL;
  796. // Check for missing layer (2 sides)
  797. Int_t lastLay = 8, frozen = 0;
  798. //Int_t lastLay = 10, frozen = 0;
  799. MpdKalmanHit *h = (MpdKalmanHit*) curTr->GetHits()->Last();
  800. if (h) lastLay = h->GetLayer();
  801. if (TMath::Abs(lay - lastLay) > 2) frozen = 1;
  802. ok = nAdd = 0;
  803. Int_t firstHit = -1, skipHit = -1;
  804. MpdItsKalmanTrack trackBr[3] = {*curTr, *curTr, *curTr};
  805. map<Int_t,Int_t> trackBrM;
  806. Bool_t navOK = NavigateToLayer(lay, curTr, trackBr, trackBrM); // navigate to layer
  807. if (!navOK) { ok = -1; continue; }
  808. // Debug
  809. //cout << " Overlaps: " << trackBrM.size() << endl;
  810. map<Int_t,Int_t>::iterator it;
  811. /*
  812. for (it = trackBrM.begin(); it != trackBrM.end(); ++it) {
  813. cout << " detID " << it->first << " " << fHitSts.GetLadder(it->first%1000000) << " "
  814. << fHitSts.GetSensor(it->first%1000000) << " " << it->second << endl;
  815. }
  816. */
  817. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  818. //if (frozen) break;
  819. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  820. if (hit->GetDetectorID() != firstHit) {
  821. firstHit = hit->GetDetectorID();
  822. if (trackBrM.find(hit->GetDetectorID()) == trackBrM.end()) {
  823. // Track and hit in different detectors - can happen near sensor edges
  824. // (hits in the same layer but in different detectors)
  825. skipHit = firstHit;
  826. continue; // next hit
  827. }
  828. Double_t leng = curTr->GetLength();
  829. branchTr = &trackBr[trackBrM[hit->GetDetectorID()]];
  830. //cout << branchTr->GetNodeNew() << endl;
  831. Double_t step = branchTr->GetLength() - leng;
  832. //cout << " Step " << step << " " << hit->GetLayer() << endl;
  833. //*
  834. //if (hit->GetLayer() % 2 == 10 && step > 1.e-4) {
  835. if (hit->GetLayer() % 2 == 0 && step > 1.e-4) {
  836. // Crossing silicon layer - add mult. scat. in the sensor
  837. Double_t x0 = 9.36; // rad. length
  838. TMatrixDSym *cov = branchTr->Weight2Cov();
  839. Double_t th = branchTr->GetParamNew(3);
  840. Double_t cosTh = TMath::Cos(th);
  841. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(branchTr, x0, step, mass2);
  842. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(branchTr, x0, step*4, mass2);
  843. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  844. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  845. (*cov)(3,3) += angle2;
  846. Int_t iok = 0;
  847. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  848. branchTr->SetWeight(*cov);
  849. } else if (hit->GetLayer() % 2 != 0 && step > 1.e-4 && fCables[lay/2].size() > 0) {
  850. /*
  851. // Crossing silicon layer - add mult. scat. in the cable
  852. //Double_t nCables = 0, x0 = 0.0116 * 2 / 9.36; // in rad. length - 116um cable per side
  853. Double_t nCables = 0, x0 = 0.0116 * 2; // in rad. length - 116um cable per side
  854. // Find number of cables crossed
  855. TString path = branchTr->GetNodeNew();
  856. if (!path.Contains("sensor") && !path.Contains("sector")) {
  857. cout << " !!! MpdTrackFinderIts::RunKalmanFilter - Outside detector !!! " << endl;
  858. exit(0);
  859. }
  860. gGeoManager->cd(path);
  861. Double_t v7[3] = {branchTr->GetParamNew(0), branchTr->GetPosNew(), branchTr->GetParamNew(1)}, v77[3];
  862. gGeoManager->LocalToMaster(v7,v77);
  863. Double_t zTr = TMath::Abs (v77[2]); // global Z
  864. cout << " zTr " << v77[0] << " " << v77[1] << " " << v77[2] << endl;
  865. cout << gGeoManager->GetPath() << " " << branchTr->GetNodeNew() << endl;
  866. map<Double_t,Double_t>::iterator it1;
  867. for (it1 = fCables[lay/2].begin(); it1 != fCables[lay/2].end(); ++it1) {
  868. if (zTr < it1->first || zTr > it1->second) continue;
  869. ++nCables;
  870. }
  871. cout << " Cables: " << nCables << endl;
  872. if (nCables) {
  873. x0 *= nCables;
  874. TMatrixDSym *cov = branchTr->Weight2Cov();
  875. Double_t th = branchTr->GetParamNew(3);
  876. Double_t cosTh = TMath::Cos(th);
  877. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(branchTr, 9.36, x0, mass2);
  878. cout << " Scat: " << hit->GetLayer() << " " << cosTh << " " << TMath::Sqrt(angle2) << endl;
  879. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  880. (*cov)(3,3) += angle2;
  881. Int_t iok = 0;
  882. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  883. branchTr->SetWeight(*cov);
  884. }
  885. */
  886. }
  887. } //if (hit->GetDetectorID() != firstHit)
  888. if (hit->GetDetectorID() == skipHit) {
  889. // Track and hit in different detectors
  890. continue; // next hit
  891. }
  892. //cout << hit->GetDetectorID() << endl;
  893. // Exclude used hits
  894. if (hit->GetFlag() != 1) continue;
  895. if (frozen) continue;
  896. // !!! Exact ID match
  897. //if (TrackID(hit) != track->GetTrackID()) continue;
  898. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterStripLocal(branchTr,hit,pointWeight,param,posNew);
  899. //cout << " Chi2: " << lay << " " << dChi2 << endl;
  900. // Add Z-contribution (if track is outside the detector)
  901. Double_t sizeZ = MpdKalmanFilter::Instance()->GetGeo()->Size(hit).Y();
  902. //if (TMath::Abs(branchTr->GetParamNew(1)) > sizeZ) {
  903. if (TMath::Abs(param(1,0)) > sizeZ) {
  904. // Outside detector
  905. //Double_t dChi2z = (TMath::Abs(branchTr->GetParamNew(1)) - sizeZ) / hit->GetErr(1);
  906. Double_t dChi2z = (TMath::Abs(param(1,0)) - sizeZ) / hit->GetErr(1);
  907. //AZ!!! dChi2 += dChi2z * dChi2z;
  908. }
  909. //cout << dChi2 << " " << TrackID(hit) << endl;
  910. if (TMath::Abs(dChi2) < fgkChi2Cut) {
  911. //if (dChi2 < fgkChi2Cut && lay % 2 == 0) {
  912. if (nBr < maxBrLay) {
  913. branchLay[nBr] = *branchTr;
  914. branchLay[nBr].GetHits()->Add(hit);
  915. branchLay[nBr].SetChi2(branchTr->GetChi2()+dChi2);
  916. TMatrixDSym w = *branchTr->GetWeight();
  917. w += pointWeight;
  918. branchLay[nBr].SetWeight(w);
  919. branchLay[nBr].SetPosNew(branchTr->GetPosNew());
  920. //if (hit->GetType() != MpdKalmanHit::kFixedP) branchLay[nBr].SetPosNew(curTr->GetPosNew());
  921. //else branchLay[nBr].SetPosNew(posNew); // modular geometry
  922. branchLay[nBr].SetLength(branchTr->GetLength());
  923. branchLay[nBr].SetParamNew(param);
  924. // Save track params at last hit
  925. branchLay[nBr].SetLengAtHit(branchTr->GetLength());
  926. branchLay[nBr].SetParamAtHit(param);
  927. branchLay[nBr].SetWeightAtHit(*branchLay[nBr].GetWeight());
  928. ++nBr;
  929. ++nAdd;
  930. }
  931. }
  932. } // for (Int_t indx = indxBeg; indx != indxEnd;
  933. // Add branch with missing layer
  934. if (nAdd == 0 && nBr < maxBrLay && ok != -1) {
  935. it = trackBrM.begin();
  936. branchLay[nBr] = trackBr[it->second];
  937. ++it;
  938. // Take branch with greater length
  939. for ( ; it != trackBrM.end(); ++it) {
  940. if (trackBr[it->second].GetLength() > branchLay[nBr].GetLength()) branchLay[nBr] = trackBr[it->second];
  941. }
  942. ++nBr;
  943. }
  944. } // for (Int_t ibr = 0; ibr < nBr0;
  945. maxInLay = TMath::Max (maxInLay,nBr);
  946. if (nBr == 0) {
  947. if (ok == -1) {
  948. //cout << track->GetTrackID() << " " << lay << endl;
  949. //return ok; // stop track
  950. cout << branch[0].GetNode() << ", " << branch[0].GetNodeNew() << endl;
  951. break;
  952. }
  953. cout << " !!! MpdTrackFinderIts::RunKalmanFilter: it can not happen |" << endl;
  954. exit(0);
  955. } else if (nBr <= maxBr) {
  956. for (Int_t i = 0; i < nBr; ++i) {
  957. branch[i].Reset();
  958. branch[i] = branchLay[i];
  959. }
  960. } else {
  961. // Too many branches - take the best ones
  962. for (Int_t i = 0; i < nBr; ++i) {
  963. Double_t c2 = TMath::Min (branchLay[i].GetChi2(),200.);
  964. quality[i] = branchLay[i].GetNofHits() + (0.999-c2/201);
  965. }
  966. TMath::Sort(nBr,quality,indx_lay);
  967. for (Int_t i = 0; i < maxBr; ++i) {
  968. branch[i].Reset();
  969. branch[i] = branchLay[indx_lay[i]];
  970. }
  971. nBr = maxBr;
  972. }
  973. } // for (Int_t lay = layOK-1; lay >= 0;
  974. // Select the best branch
  975. cout << " Branches: " << nBr << " " << maxInLay << endl;
  976. Int_t ibest = 0;
  977. for (Int_t i = 1; i < nBr; ++i) {
  978. cout << i << " " << branch[i].GetNofHits() << " " << branch[i].GetChi2() << endl;
  979. if (branch[i].GetNofHits() > branch[ibest].GetNofHits()) ibest = i;
  980. else if (branch[i].GetNofHits() == branch[ibest].GetNofHits() &&
  981. branch[i].GetChi2() < branch[ibest].GetChi2()) ibest = i;
  982. }
  983. track->Reset();
  984. *track = branch[ibest];
  985. return 0;
  986. }
  987. //__________________________________________________________________________
  988. Bool_t MpdTrackFinderIts::NavigateToLayer(Int_t lay, MpdItsKalmanTrack *curTr, MpdItsKalmanTrack *trackBr,
  989. std::map<Int_t,Int_t>& trackBrM)
  990. {
  991. ///< Navigate track to layer in modular geometry
  992. const Double_t posCable[4] = {0.08, 0.12, 0.24, 0.40}; // average cable positions (1-4) w.r.t. sensor
  993. Int_t inode = 0, detID[3] = {-1,-1,-1}, nCables = 0;
  994. TString curPath, nodeNew, mass2 = "0.0194797849"; // pion mass squared
  995. Double_t zTr = 0.0, x0 = 0.0, posNew = 0.0;
  996. MpdKalmanHit hitTmp;
  997. TGeoNode *node = NULL;
  998. TMatrixD parNew = *curTr->GetParamNew();
  999. //*
  1000. if (lay % 2 == 0) {
  1001. TString path = curTr->GetNodeNew();
  1002. Int_t last = path.Length() - 1;
  1003. Int_t ladd = path.Index("ladder") - 5;
  1004. TString name = path(ladd, last-ladd+1) + "#";
  1005. name += (lay % 2);
  1006. detID[0] = MpdKalmanFilter::Instance()->GetGeo()->DetId(name);
  1007. goto skip;
  1008. }
  1009. //*/
  1010. posNew = curTr->GetPosNew();
  1011. nodeNew = curTr->GetNodeNew();
  1012. //cout << " Nodes: " << nodeNew << " " << curTr->GetNode() << " " << posNew << " " << curTr->GetPos() << endl;
  1013. //curTr->GetParamNew()->Print();
  1014. //curTr->GetParam()->Print();
  1015. curPath = curTr->GetNode();
  1016. hitTmp.SetType(MpdKalmanHit::kFixedR);
  1017. hitTmp.SetPos(fRad[lay]);
  1018. if (!MpdKalmanFilter::Instance()->PropagateParamR(curTr,&hitTmp,kFALSE)) {
  1019. curTr->SetPosNew(posNew);
  1020. curTr->SetParamNew(parNew);
  1021. curTr->SetNodeNew(nodeNew);
  1022. curTr->SetNode(curPath);
  1023. //ok = -1;
  1024. //break;
  1025. return kFALSE;
  1026. }
  1027. // Find sensor (sector)
  1028. Double_t v7[7], v77[7];
  1029. v77[6] = -1;
  1030. curTr->SetNode(curTr->GetNodeNew());
  1031. MpdKalmanFilter::Instance()->SetGeantParamB(curTr, v7, 1);
  1032. //MpdKalmanFilter::Instance()->SetGeantParamB(curTr, v7, -1);
  1033. node = gGeoManager->FindNode(v7[0], v7[1], v7[2]);
  1034. //if (!TString(gGeoManager->GetPath()).Contains("sts01")) break; // outside ITS coverage
  1035. if (!TString(gGeoManager->GetPath()).Contains("sts01")) return kFALSE; // outside ITS coverage
  1036. zTr = v7[2];
  1037. for (Int_t iover = 0; iover < 2; ++iover) {
  1038. // Loop to find first detector
  1039. if (iover == 0) {
  1040. for (Int_t i2 = 0; i2 < 2; ++i2) {
  1041. TString path = gGeoManager->GetPath();
  1042. //cout << path << endl;
  1043. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1044. if (!path.Contains("ladder")) {
  1045. gGeoManager->SetStep(100);
  1046. gGeoManager->FindNextDaughterBoundary(v7,&v7[3],inode);
  1047. if (inode >= 0) { gGeoManager->CdDown(inode);
  1048. /*cout << " Node: " << inode << " " << gGeoManager->GetPath() << endl;*/ }
  1049. else break;
  1050. } else {
  1051. gGeoManager->MasterToLocal(v7,v77);
  1052. gGeoManager->MasterToLocalVect(&v7[3],&v77[3]);
  1053. gGeoManager->SetStep(100);
  1054. // Step in normal direction (in transverse plane)
  1055. Double_t r2 = v77[3] * v77[3] + v77[4] * v77[4];
  1056. v77[3] = 0;
  1057. v77[4] = TMath::Sqrt (r2) * TMath::Sign(1.,v77[4]);
  1058. //v77[4] = TMath::Sqrt (r2) * TMath::Sign(1.,-v77[4]*v77[0]); // for rotated ladders
  1059. gGeoManager->FindNextDaughterBoundary(v77,&v77[3],inode);
  1060. if (inode >= 0) {
  1061. gGeoManager->CdDown(inode);
  1062. for (Int_t i3 = 0; i3 < 3; ++i3) v77[i3] += v77[i3+3] * (gGeoManager->GetStep() + 0.04);
  1063. v77[6] = 1;
  1064. }
  1065. //cout << " ok " << gGeoManager->GetCurrentPoint()[0] << " " << gGeoManager->GetCurrentPoint()[1] << " " << gGeoManager->GetCurrentPoint()[2] << " " << gGeoManager->GetStep() << endl;
  1066. zTr += (gGeoManager->GetStep() * v77[5]);
  1067. //cout << zTr << endl;
  1068. }
  1069. } else break;
  1070. } // for (Int_t i2 = 0; i2 < 2;
  1071. } else { // if (iover == 0)
  1072. // Look for overlap (at first in Z)
  1073. TString path = gGeoManager->GetPath();
  1074. while (path.Contains("sensor") || path.Contains("sector")) {
  1075. // find ladder
  1076. gGeoManager->CdUp();
  1077. path = gGeoManager->GetPath();
  1078. }
  1079. if (v77[6] < 0) {
  1080. gGeoManager->MasterToLocal(v7,v77);
  1081. gGeoManager->MasterToLocalVect(&v7[3],&v77[3]);
  1082. Double_t r2 = v77[3] * v77[3] + v77[4] * v77[4];
  1083. v77[3] = 0;
  1084. v77[4] = TMath::Sqrt (r2) * TMath::Sign(1.,v77[4]);
  1085. for (Int_t i3 = 0; i3 < 3; ++i3) v77[i3] += v77[i3+3] * 0.04;
  1086. }
  1087. gGeoManager->FindNextDaughterBoundary(v77,&v77[3],inode);
  1088. if (inode >= 0) gGeoManager->CdDown(inode);
  1089. else {
  1090. // Look for overlap in adjacent ladders
  1091. Double_t secPos[3];
  1092. gGeoManager->LocalToMaster(v77,secPos);
  1093. gGeoManager->CdUp();
  1094. path = gGeoManager->GetPath();
  1095. Double_t r2 = TMath::Sqrt (v7[3] * v7[3] + v7[4] * v7[4]);
  1096. //v7[5] = 0.0;
  1097. //cout << secPos[0] << " " << secPos[1] << endl;
  1098. for (Int_t i3 = 0; i3 < 3; ++i3) secPos[i3] += v7[i3+3] * 0.1 / r2;
  1099. //cout << secPos[0] << " " << secPos[1] << " " << gGeoManager->GetStep() << endl;
  1100. gGeoManager->SetStep(4);
  1101. gGeoManager->FindNextDaughterBoundary(secPos,&v7[3],inode);
  1102. if (inode >= 0) {
  1103. // Found ladder - find sector (sensor)
  1104. for (Int_t i3 = 0; i3 < 3; ++i3) secPos[i3] += v7[i3+3] * (gGeoManager->GetStep() + 0.001);
  1105. gGeoManager->CdDown(inode);
  1106. gGeoManager->MasterToLocal(secPos,v77);
  1107. gGeoManager->MasterToLocalVect(&v7[3],&v77[3]);
  1108. //cout << path << " " << gGeoManager->GetPath() << " " << gGeoManager->GetStep() << endl;
  1109. gGeoManager->SetStep(1);
  1110. gGeoManager->FindNextDaughterBoundary(v77,&v77[3],inode);
  1111. //cout << inode << endl;
  1112. //if (inode >= 0) exit(0);
  1113. if (inode >= 0) gGeoManager->CdDown(inode);
  1114. }
  1115. }
  1116. } // if (iover == 0)
  1117. TString path = gGeoManager->GetPath();
  1118. if (!path.Contains("ladder")) continue;
  1119. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1120. //if (iover) continue;
  1121. continue;
  1122. //Fatal(" RunKalmanFilter ", path);
  1123. }
  1124. //cout << path << endl;
  1125. Int_t last = path.Length() - 1;
  1126. //if (path.Contains("sensor") || path.Contains("frame")) last = path.Last('/') - 1;
  1127. Int_t ladd = path.Index("ladder") - 5;
  1128. //TString name = path(ladd, last-ladd+1);
  1129. TString name = path(ladd, last-ladd+1) + "#";
  1130. //name += (!track->GetDirection());
  1131. name += (lay % 2);
  1132. //cout << name << endl;
  1133. detID[iover] = MpdKalmanFilter::Instance()->GetGeo()->DetId(name);
  1134. // Debugging
  1135. //TVector3 v3tmp;
  1136. //MpdStsHit hitDb1(detID[iover]%1000000,v3tmp,v3tmp,0.0,-1);
  1137. //MpdStsHit hitDb2(hit->GetDetectorID()%1000000,v3tmp,v3tmp,0.0,-1);
  1138. //cout << iover << " " << name << " " << detID[iover] << " " << hit->GetDetectorID() << " " << hitDb1.Ladder() << " " << hitDb1.Sector() << " " << hitDb2.Ladder() << " " << hitDb2.Sector() << " " << hit->GetMeas(0) << endl;
  1139. break;
  1140. } // for (Int_t iover = 0; iover < 2;
  1141. //if (detID[0] < 0 && detID[1] < 0) break;
  1142. if (detID[0] < 0 && detID[1] < 0) return kFALSE;
  1143. if (detID[0] < 0) { detID[0] = detID[1]; detID[1] = -1; }
  1144. // if (detID[0] >= 0 && detID[1] >= 0) exit(0); // Debugging
  1145. curTr->SetPosNew(posNew);
  1146. curTr->SetParamNew(parNew);
  1147. curTr->SetNodeNew(nodeNew);
  1148. curTr->SetNode(curPath);
  1149. //firstDet = 0;
  1150. if (lay % 2 != 0 && fCables[lay/2].size() > 0) {
  1151. // Find number of cables crossed
  1152. TString path = gGeoManager->GetPath();
  1153. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1154. cout << " !!! MpdTrackFinderIts::RunKalmanFilterMod - Outside detector !!! " << endl;
  1155. exit(0);
  1156. }
  1157. zTr = TMath::Abs (zTr); // global Z
  1158. map<Double_t,Double_t>::iterator it1;
  1159. for (it1 = fCables[lay/2].begin(); it1 != fCables[lay/2].end(); ++it1) {
  1160. if (zTr < it1->first || zTr > it1->second) continue;
  1161. ++nCables;
  1162. }
  1163. }
  1164. //cout << " zTr " << zTr << " " << nCables << endl;
  1165. //nCables = 0;
  1166. // Crossing silicon layer - add mult. scat. in the cable
  1167. //Double_t x0 = 0.0116 * 2 / 9.36 * nCables; // in rad. length - 116um cable per side
  1168. x0 = 0.0116 * 2 * nCables; // in rad. length - 116um cable per side
  1169. if (fMCTracks) {
  1170. // Get particle mass - ideal PID
  1171. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(curTr->GetTrackID());
  1172. TParticlePDG *pdgP = TDatabasePDG::Instance()->GetParticle(mctrack->GetPdgCode());
  1173. if (pdgP) {
  1174. Double_t mass = pdgP->Mass();
  1175. if (mass < 0.1 || mass > 0.8) {
  1176. // Electrons or protons (or heavier than protons)
  1177. mass2 = "";
  1178. mass2 += mass*mass;
  1179. }
  1180. }
  1181. }
  1182. skip:
  1183. Int_t firstBranch = 1;
  1184. // Propagate track to the found detector and look for the overlapping sensors and ladders
  1185. hitTmp.SetType(MpdKalmanHit::kFixedP);
  1186. for (Int_t iover = 0; iover < 3; ++iover) {
  1187. if (detID[iover] >= 0) {
  1188. // Already found detector
  1189. hitTmp.SetDetectorID(detID[iover]);
  1190. if (nCables) {
  1191. // Propagate to average cable position
  1192. MpdKalmanHit hitCable;
  1193. hitCable.SetDetectorID(-2); // just some fake detector ID
  1194. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  1195. TVector3 pos = geo->GlobalPos(&hitTmp);
  1196. //pos.Print();
  1197. pos += geo->Normal(&hitTmp) * posCable[nCables-1];
  1198. //pos.Print();
  1199. geo->SetGlobalPos(-2, pos, kTRUE);
  1200. geo->SetNormal(-2, geo->Normal(&hitTmp), kTRUE);
  1201. geo->SetPath(-2, geo->Path(detID[iover]), kTRUE);
  1202. if (!MpdKalmanFilter::Instance()->PropagateToHit(&trackBr[iover],&hitCable,kTRUE,kTRUE)) {
  1203. detID[iover] = -1;
  1204. continue;
  1205. }
  1206. TMatrixDSym *cov = trackBr[iover].Weight2Cov();
  1207. Double_t th = trackBr[iover].GetParamNew(3);
  1208. Double_t cosTh = TMath::Cos(th);
  1209. x0 /= cosTh;
  1210. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(&trackBr[iover], 9.36, x0*2, mass2);
  1211. //cout << TMath::Sqrt(angle2) << " " << cosTh << endl;
  1212. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1213. (*cov)(3,3) += angle2;
  1214. Int_t iok = 0;
  1215. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1216. trackBr[iover].SetWeight(*cov);
  1217. }
  1218. if (!MpdKalmanFilter::Instance()->PropagateToHit(&trackBr[iover],&hitTmp,kTRUE,kTRUE)) {
  1219. //cout << detID[iover] << " " << MpdKalmanFilter::Instance()->GetGeo()->Path(detID[iover]) << endl;
  1220. //cout << iover << " " << detID[iover] << " " << trackBr[iover].GetNode() << " " << trackBr[iover].GetNodeNew() << " " << trackBr[iover].GetNofHits() << endl;
  1221. detID[iover] = -1;
  1222. continue;
  1223. }
  1224. if (firstBranch) {
  1225. // Try to get the closest neighbours
  1226. // Check if the track is close to the detector edges
  1227. Double_t xloc = trackBr[iover].GetParamNew(0);
  1228. Double_t zloc = trackBr[iover].GetParamNew(1);
  1229. hitTmp.SetDetectorID(detID[iover]);
  1230. TVector2 size = MpdKalmanFilter::Instance()->GetGeo()->Size(&hitTmp);
  1231. Int_t idSTS = detID[iover] % 1000000;
  1232. //cout << xloc << " " << zloc << " " << size.X() << " " << size.Y() << endl;
  1233. if (TMath::Abs(TMath::Abs(xloc)-size.X()) < 0.3) {
  1234. Int_t ladder = fHitSts.GetLadder(idSTS);
  1235. if (xloc < 0) {
  1236. ++ladder;
  1237. if (ladder > fNladders[lay]) ladder = 1;
  1238. } else {
  1239. --ladder;
  1240. if (ladder == 0) ladder = fNladders[lay];
  1241. }
  1242. if (iover < 2 && detID[iover+1] < 0) detID[iover+1] = 1000000*lay + fHitSts.SetLadder(idSTS,ladder);
  1243. else if (iover < 1 && detID[iover+2] < 0) detID[iover+2] = 1000000*lay + fHitSts.SetLadder(idSTS,ladder);
  1244. }
  1245. if (TMath::Abs(TMath::Abs(zloc)-size.Y()) < 0.3) {
  1246. Int_t sector = fHitSts.GetSensor(idSTS);
  1247. if (zloc < 0) --sector;
  1248. else ++sector;
  1249. if (sector > 0 && sector <= fNsectors[lay]) {
  1250. if (iover < 2 && detID[iover+1] < 0) detID[iover+1] = 1000000*lay + fHitSts.SetSensor(idSTS,sector);
  1251. else if (iover < 1 && detID[iover+2] < 0) detID[iover+2] = 1000000*lay + fHitSts.SetSensor(idSTS,sector);
  1252. }
  1253. }
  1254. firstBranch = 0;
  1255. }
  1256. }
  1257. }
  1258. if (firstBranch) {
  1259. // Failed propagation
  1260. cout << " !!! Failed navigation to layer " << lay << endl;
  1261. //exit(0);
  1262. return kFALSE;
  1263. }
  1264. for (Int_t iover = 0; iover < 3; ++iover) {
  1265. if (detID[iover] < 0) continue;
  1266. trackBrM.insert(pair<Int_t,Int_t>(detID[iover],iover));
  1267. }
  1268. return kTRUE;
  1269. }
  1270. //__________________________________________________________________________
  1271. Int_t MpdTrackFinderIts::TrackID(MpdKalmanHit *hit)
  1272. {
  1273. /// Return track ID of the hit
  1274. FairHit *h = (FairHit*) fItsHits->UncheckedAt(hit->GetIndex());
  1275. return ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1276. }
  1277. //__________________________________________________________________________
  1278. TVector2 MpdTrackFinderIts::GetDistance(MpdKalmanTrack *track, MpdKalmanHit *hit)
  1279. {
  1280. /// Compute distance between track and hit
  1281. Int_t lay = hit->GetLayer();
  1282. Int_t lay2 = lay / 2;
  1283. Int_t side = lay % 2;
  1284. Int_t module = ((MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex()))->Module();
  1285. Double_t zTr = track->GetParamNew(1);
  1286. Double_t zloc = zTr + fDz[lay2];
  1287. Int_t modTr = Int_t (zloc / fZmod[lay2]);
  1288. Int_t dMod = modTr - module;
  1289. Double_t dZ = 0;
  1290. if (dMod != 0) {
  1291. // Not in the same module - check Z-distance to module edge
  1292. if (dMod > 0) dZ = zloc - (module+1) * fZmod[lay2];
  1293. else dZ = zloc - module * fZmod[lay2];
  1294. if (TMath::Abs(dMod) > 2) return TVector2(TMath::Abs(dZ),999.); // not in neighbour modules
  1295. }
  1296. // Translate transverse track position to local system
  1297. Double_t xloc = track->GetParamNew(0) * hit->GetCosSin(0) + zTr * hit->GetCosSin(1);
  1298. Double_t phTr = xloc / track->GetPosNew();
  1299. Double_t phHit = hit->GetMeas(0) / fRad[lay];
  1300. //TVector2 dist(TMath::Abs(dZ), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit)-phTr));
  1301. TVector2 dist(TMath::Abs(dZ),
  1302. TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit,hit->GetCosSin(0))-phTr));
  1303. //TVector2 dist(TMath::Abs(zTr-hit->GetZ()), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(track,hit)/hit->GetR()-ph));
  1304. return dist;
  1305. }
  1306. //__________________________________________________________________________
  1307. void MpdTrackFinderIts::Write()
  1308. {
  1309. /// Write
  1310. TFile histoFile("ItsRec.root","RECREATE");
  1311. Writedir2current(fHistoDir);
  1312. histoFile.Close();
  1313. }
  1314. //__________________________________________________________________________
  1315. void MpdTrackFinderIts::Writedir2current( TObject *obj )
  1316. {
  1317. /// Write
  1318. if( !obj->IsFolder() ) obj->Write();
  1319. else{
  1320. TDirectory *cur = gDirectory;
  1321. TDirectory *sub = cur->mkdir(obj->GetName());
  1322. sub->cd();
  1323. TList *listSub = ((TDirectory*)obj)->GetList();
  1324. TIter it(listSub);
  1325. while( TObject *obj1=it() ) Writedir2current(obj1);
  1326. cur->cd();
  1327. }
  1328. }
  1329. //__________________________________________________________________________
  1330. void MpdTrackFinderIts::AddHits()
  1331. {
  1332. /// Add hit objects to tracks and compute number of wrongly assigned hits
  1333. /// (hits with ID different from ID of starting TPC track)
  1334. Int_t nFound = fTracks->GetEntriesFast();
  1335. for (Int_t i = 0; i < nFound; ++i) {
  1336. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);
  1337. cout << track->GetNode() << " " << track->GetNodeNew() << endl;
  1338. Double_t c2 = track->GetChi2();
  1339. track->SetChi2(track->GetChi2Its());
  1340. track->SetChi2Its(c2);
  1341. Int_t nHits = track->GetNofHits();
  1342. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1343. TClonesArray &trHits = *track->GetTrHits();
  1344. cout << nHits << " " << trHits.GetEntriesFast() << " " << track->GetTrackID() << endl;
  1345. TObjArray *hits = track->GetHits();
  1346. Int_t nWrong = 0, nMirr = 0, motherID = track->GetTrackID();
  1347. // Get track mother ID
  1348. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  1349. while (mctrack->GetMotherId() >= 0) {
  1350. motherID = mctrack->GetMotherId();
  1351. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1352. }
  1353. Int_t lastIndx = trHits.GetEntriesFast();
  1354. for (Int_t j = 0; j < nHits; ++j) {
  1355. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1356. hit->SetUniqueID(1); // flag ITS hits
  1357. new (trHits[lastIndx+j]) MpdKalmanHit(*hit);
  1358. cout << " " << hit->GetLayer();
  1359. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex());
  1360. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1361. cout << "-" << motherID1;
  1362. // Get point mother ID
  1363. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID1);
  1364. while (mctrack->GetMotherId() >= 0) {
  1365. motherID1 = mctrack->GetMotherId();
  1366. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1367. }
  1368. if (motherID1 != motherID) ++nWrong;
  1369. }
  1370. if (nHits) cout << "\n" << nWrong << endl;
  1371. track->SetNofWrong(nWrong);
  1372. MpdKalmanTrack *tpc = (MpdKalmanTrack*) fTpcTracks->UncheckedAt(track->GetUniqueID()-1);
  1373. //track->SetChi2(track->GetChi2()+tpc->GetChi2());
  1374. //track->SetLastLay();
  1375. //track->GetParam()->Print();
  1376. track->SetNofHits(track->GetNofTrHits()); // TPC and ITS hits
  1377. track->SetNofIts(nHits);
  1378. cout << nHits << " " << track->GetNofTrHits() << " " << track->GetTrackID() << " "
  1379. << track->GetChi2Its() << " " << track->GetChi2() << endl;
  1380. }
  1381. fTracks->Compress();
  1382. }
  1383. ClassImp(MpdTrackFinderIts);