MpdCellAutomat.cxx 71 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768
  1. // -------------------------------------------------------------------------
  2. // ----- source file -----
  3. // ----- Created 13/06/2014 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdCellAutomat.cxx
  6. *@author M.Strelchenko, A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD Inner Tracking System (ITS) using cellular automaton (CA)
  9. **/
  10. #include "MpdStsGeoPar.h"
  11. //#include "MpdTrackFinderIts.h"
  12. #include "MpdCellAutomat.h"
  13. #include "MpdCellTrack.h"
  14. #include "MpdCodeTimer.h"
  15. #include "MpdConstField.h"
  16. #include "MpdKalmanFilter.h"
  17. #include "MpdKalmanGeoScheme.h"
  18. #include "MpdKalmanHit.h"
  19. #include "MpdKalmanTrack.h"
  20. #include "MpdStsHit.h"
  21. #include "MpdStsPoint.h"
  22. #include "MpdItsKalmanTrack.h"
  23. //#include "MpdTpcKalmanTrack.h"
  24. //#include "MpdTpcKalmanFilter.h"
  25. #include "FairGeoNode.h"
  26. #include "FairMCPoint.h"
  27. #include "MpdMCTrack.h"
  28. #include "FairRootManager.h"
  29. #include "FairRun.h"
  30. #include "FairRunAna.h"
  31. #include "FairRuntimeDb.h"
  32. #include "TGeoMatrix.h"
  33. #include "TGeoBBox.h"
  34. #include "TGeoNode.h"
  35. #include "TGeoTube.h"
  36. #include "TGeoManager.h"
  37. #include "TMath.h"
  38. //#include "TFile.h"
  39. #include "TVector2.h"
  40. #include "TVector3.h"
  41. #include "TClonesArray.h"
  42. #include <TRandom.h>
  43. #include <TArrayI.h>
  44. #include <iostream>
  45. #include <map>
  46. #include <set>
  47. using std::cout;
  48. using std::endl;
  49. //using std::map;
  50. const Double_t MpdCellAutomat::fgkChi2Cut = 20; //20; //100;
  51. FILE *lunErr = 0x0; //fopen("error1.dat","w");
  52. //FILE* lunErr= fopen ("file.txt","w");
  53. //__________________________________________________________________________
  54. MpdCellAutomat::MpdCellAutomat(const char *name, Int_t iVerbose )
  55. :FairTask(name, iVerbose),
  56. fExact(0),
  57. fNTracks(0),
  58. fNPass(2),
  59. fGeo(0)
  60. {
  61. fKHits1 = new TClonesArray("MpdKalmanHit", 100);
  62. for (Int_t i = 0; i < 4; ++i) {
  63. fKHits[i] = new TClonesArray("MpdCellTrack", 100);
  64. f2DHits[i] = new TClonesArray("MpdCellTrack", 100);
  65. }
  66. //fKHits[0] = f2DHits[0]; //work
  67. fTracks = new TClonesArray("MpdItsKalmanTrack", 100);
  68. fTrackCand = new TClonesArray("MpdItsKalmanTrack", 100);
  69. fHistoDir = 0x0;
  70. fhLays = new TH1F("hLaysITS","ITS layers",10,0,10);
  71. fLayPointers = 0x0;
  72. }
  73. //__________________________________________________________________________
  74. MpdCellAutomat::~MpdCellAutomat()
  75. {
  76. //delete fKHits;
  77. //delete fTracks;
  78. //delete fTrackCand;
  79. delete [] fLayPointers;
  80. delete fhLays;
  81. }
  82. //__________________________________________________________________________
  83. InitStatus MpdCellAutomat::Init()
  84. {
  85. //return ReInit();
  86. if (ReInit() != kSUCCESS) return kERROR;
  87. // Read database
  88. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  89. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  90. //cout << geoPar << endl;
  91. TString volName = "sts01 ", path = "";
  92. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  93. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  94. Int_t nLay = 4;
  95. Double_t size = 6.2;
  96. for (Int_t i = 0; i < nLay; ++i) {
  97. volName[5] = 97 + i; // 'a', 'b', ..
  98. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->FindObject(volName));
  99. if (sensVol == 0x0) {
  100. fGeo = 1; // modular geometry
  101. break;
  102. }
  103. TArrayD* params = sensVol->getParameters();
  104. fRad[2*i] = params->At(0);
  105. fRad[2*i+1] = params->At(1);
  106. //dR = (rMax-rMin) / 50;
  107. fDz[i] = params->At(2);
  108. Int_t nMods = Int_t (fDz[i] * 2. / size + 0.1);
  109. fZmod[i] = fDz[i] * 2. / nMods;
  110. cout << " *** STS sensitive volume: " << sensVol->GetName() << " " << params->At(0)
  111. << " " << fDz[i] << " " << fZmod[i] << endl;
  112. }
  113. //fStereoA[0] = -7.5;
  114. //fStereoA[1] = 7.5;
  115. fStereoA[0] = 7.5;
  116. fStereoA[1] = 0.0;
  117. //for (Int_t i = 0; i < 2; ++i) fStereoA[i] =
  118. //TMath::Tan (fStereoA[i]*TMath::DegToRad());
  119. for (Int_t i = 0; i < 2; ++i) fStereoA[i] *= TMath::DegToRad();
  120. if (fGeo) {
  121. // Process modular geometry
  122. Double_t safety = 0.03;
  123. for (Int_t i = 0; i < nLay; ++i) {
  124. fNladders[2*i] = fNladders[2*i+1] = fNsectors[2*i] = fNsectors[2*i+1] = 0;
  125. volName = "sts01ladder";
  126. volName += (i+1);
  127. path = "/cave_1/sts01_0/" + volName;
  128. path += "_";
  129. TString path0 = path;
  130. //volName = "/cave_1/sts01_0/sts01ladder";
  131. //volName += (i+1);
  132. // Loop over all ladders to find the one with the smallest radius
  133. fRad[2*i+1] = fRad[2*i] = 999.9;
  134. for (Int_t jlad = 1; jlad < 99; ++jlad) {
  135. path = path0;
  136. path += jlad;
  137. if (!gGeoManager->CheckPath(path)) break;
  138. gGeoManager->cd(path);
  139. Double_t xyzL[3] = {0}, xyzM[3];
  140. gGeoManager->LocalToMaster(xyzL,xyzM);
  141. Double_t rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  142. fRad[2*i] = TMath::Min (fRad[2*i],rad);
  143. }
  144. fRad[2*i+1] = fRad[2*i];
  145. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  146. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  147. //safety = -box->GetDY();
  148. safety = box->GetDY();
  149. fRad[2*i] -= safety;
  150. fRad[2*i+1] -= safety;
  151. if (i == 0) { fRad[2*i] -= safety; fRad[2*i+1] -= safety; }
  152. }
  153. FillGeoScheme();
  154. }
  155. // Get pipe radius
  156. fPipeR = ((TGeoTube*)gGeoManager->GetVolume("pipe1")->GetShape())->GetRmax();
  157. // Get cables
  158. TObjArray *vols = gGeoManager->GetListOfVolumes();
  159. Int_t nvols = vols->GetEntriesFast(), ncables = 0;
  160. for (Int_t i = 0; i < nvols; ++i) {
  161. TGeoVolume *vol = (TGeoVolume*) vols->UncheckedAt(i);
  162. TString cable = TString(vol->GetName());
  163. if (!(cable.Contains("sts") && cable.Contains("cable"))) continue;
  164. //cout << cable << endl;
  165. ++ncables;
  166. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  167. TString lad = cable;
  168. Int_t ip = lad.Index("cable");
  169. lad.Replace(ip,lad.Length()-ip+1,"");
  170. lad.Replace(lad.Length()-2,1,"");
  171. lad += "_1/";
  172. path = "/cave_1/sts01_0/" + lad + cable;
  173. path += "_1";
  174. gGeoManager->cd(path);
  175. Double_t xyzL[3] = {0}, xyzM[3];
  176. gGeoManager->LocalToMaster(xyzL,xyzM);
  177. //cout << xyzM[2] - box->GetDZ() << " " << xyzM[2] + box->GetDZ() << " " << box->GetDZ() << endl;
  178. if (xyzM[2] - box->GetDZ() > 0) {
  179. Int_t lay = TString(lad(lad.Length()-4,1)).Atoi();
  180. fCables[lay-1].insert(pair<Double_t,Double_t>(xyzM[2] - box->GetDZ(),xyzM[2] + box->GetDZ()));
  181. }
  182. }
  183. return kSUCCESS;
  184. }
  185. //__________________________________________________________________________
  186. void MpdCellAutomat::FillGeoScheme()
  187. {
  188. /// Fill Kalman filter geometry manager info
  189. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  190. TGeoVolume *vSts = gGeoManager->GetVolume("sts01");
  191. for (Int_t layer = 1; layer < 999; ++layer) {
  192. // Loop over layers
  193. TString sladder = "sts01ladder";
  194. sladder += layer;
  195. TGeoVolume *vLay = gGeoManager->GetVolume(sladder);
  196. if (vLay == 0x0 && layer == 1) continue; // no first layer
  197. if (vLay == 0x0) break;
  198. sladder += "_";
  199. //TString sdet = "sts01detector";
  200. TString sdet = "sts01sensor";
  201. sdet += layer;
  202. sdet += "_";
  203. //TString sensor = "sts01sensor";
  204. TString sensor = "";
  205. //sensor += layer;
  206. //sensor += "_0";
  207. for (Int_t ladder = 1; ladder < 999; ++ladder) {
  208. // Loop over ladders
  209. TString sladder1 = sladder;
  210. sladder1 += ladder;
  211. TGeoNode *node = vSts->FindNode(sladder1);
  212. if (node == 0x0) break;
  213. ++fNladders[2*(layer-1)];
  214. ++fNladders[2*(layer-1)+1];
  215. TGeoVolume *vLad = node->GetVolume();
  216. //cout << vLad->GetNodes()->GetEntriesFast() << " " << vLad->GetNdaughters() << endl;
  217. Int_t nDaught = vLad->GetNdaughters(), detID = -1, detIDsts = -1;
  218. TObjArray *daught = vLad->GetNodes();
  219. //for (Int_t j = 0; j < vLad->GetNdaughters(); ++j)
  220. //cout << ((TGeoNode*)(vLad->GetNodes()->UncheckedAt(j)))->GetName() << endl;
  221. Int_t iZ = 0;
  222. for (Int_t det = 0; det < nDaught; ++det) {
  223. // Loop over ladder daughters
  224. TString sdet1 = ((TGeoNode*)(daught->UncheckedAt(det)))->GetName();
  225. if (!sdet1.Contains("sensor") && !sdet1.Contains("sector")) continue;
  226. if (ladder == 1) { ++fNsectors[2*(layer-1)]; ++fNsectors[2*(layer-1)+1]; }
  227. Int_t det1 = TString(sdet1(sdet1.Index("_")+1,2)).Atoi();
  228. Int_t secType = -1;
  229. if (sdet1.Contains("sector")) secType = TString(sdet1(sdet1.Index("_")-2,1)).Atoi();
  230. ++iZ;
  231. for (Int_t side = 0; side < 2; ++side) {
  232. detIDsts = fHitSts.SetDetId(secType, layer, ladder, det1, side);
  233. detID = fHitSts.SetDetId(0, layer, ladder, iZ, side);
  234. //detIDsts += 1000000 * ((layer-1) * 2 + side);
  235. fId2Id[(layer-1) * 2 + side].insert(pair<Int_t,Int_t>(detIDsts,detID));
  236. detID += 1000000 * ((layer-1) * 2 + side);
  237. TString detName = sladder1 + "/" + sdet1 + "#";
  238. detName += side;
  239. //if (side) geo->SetDetId(detName, detID&((2<<26)-2)); // store odd side
  240. //if (!side) geo->SetDetId(detName, detID); // store even side
  241. geo->SetDetId(detName, detID);
  242. //TString path = "/cave_1/sts01_0/" + detName + "/" + sensor;
  243. TString path = "/cave_1/sts01_0/" + detName(0,detName.Length()-2);
  244. //cout << detName << " " << path << endl;
  245. gGeoManager->cd(path);
  246. //if (side) geo->SetPath(detID&((2<<26)-2), path); // store odd side
  247. //if (!side) geo->SetPath(detID, path); // store even side
  248. geo->SetPath(detID, path);
  249. node = gGeoManager->GetCurrentNode();
  250. //cout << node->GetName() << " " << detID << endl;
  251. TGeoVolume *vol = node->GetVolume();
  252. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  253. TVector2 size(box->GetDX(), box->GetDZ());
  254. geo->SetSize(detID, size);
  255. Double_t xyzL[3] = {0}, xyzM[3], vecM[3];
  256. xyzL[1] = 1;
  257. gGeoManager->LocalToMasterVect(xyzL,vecM);
  258. //cout << vecM[0] << " " << vecM[1] << " " << vecM[2] << endl;
  259. //TVector3 norm(vecM[0], vecM[1], vecM[2]);
  260. //geo->SetNormal(detID, norm);
  261. xyzL[1] = 0;
  262. gGeoManager->LocalToMaster(xyzL,xyzM);
  263. Double_t sign = TMath::Sign (1.,xyzM[0]*vecM[0]+xyzM[1]*vecM[1]);
  264. if (detID % 2) xyzL[1] = 0.015 * sign;
  265. else xyzL[1] = -0.015 * sign;
  266. gGeoManager->LocalToMaster(xyzL,xyzM);
  267. //cout << xyzM[0] << " " << xyzM[1] << " " << xyzM[2] << endl;
  268. TVector3 pos(xyzM[0], xyzM[1], xyzM[2]);
  269. geo->SetGlobalPos(detID, pos);
  270. //xyzL[1] = sign;
  271. //gGeoManager->LocalToMasterVect(xyzL,vecM);
  272. //cout << vecM[0] << " " << vecM[1] << " " << vecM[2] << endl;
  273. TVector3 norm(sign*vecM[0], sign*vecM[1], sign*vecM[2]);
  274. geo->SetNormal(detID, norm);
  275. }
  276. }
  277. }
  278. }
  279. }
  280. //__________________________________________________________________________
  281. InitStatus MpdCellAutomat::ReInit()
  282. {
  283. fItsPoints = (TClonesArray *) FairRootManager::Instance()->GetObject("StsPoint");
  284. fItsHits =(TClonesArray *) FairRootManager::Instance()->GetObject("StsHit");
  285. if (fItsPoints == 0x0 || fItsHits == 0x0) return kERROR;
  286. fTpcTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  287. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  288. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  289. //fPrimVtx = (FairVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
  290. FairRootManager::Instance()->Register("ItsTrack", "Its", fTracks, kTRUE);
  291. FairRootManager::Instance()->Register("ItsTrackCand", "ItsCand", fTrackCand, kTRUE);
  292. FairRootManager::Instance()->Register("CellTrack0", "Its0", fKHits[0], kTRUE);
  293. FairRootManager::Instance()->Register("CellTrack1", "Its1", fKHits[1], kTRUE);
  294. FairRootManager::Instance()->Register("CellTrack2", "Its2", fKHits[2], kTRUE);
  295. FairRootManager::Instance()->Register("CellTrack3", "Its3", fKHits[3], kTRUE);
  296. FairRootManager::Instance()->Register("CellTrack_0", "Its0_2D", f2DHits[0], kTRUE); // 21.04
  297. FairRootManager::Instance()->Register("CellTrack_1", "Its1_2D", f2DHits[1], kTRUE); // 21.04
  298. FairRootManager::Instance()->Register("CellTrack_2", "Its2_2D", f2DHits[2], kTRUE); // 21.04
  299. FairRootManager::Instance()->Register("CellTrack_3", "Its3_2D", f2DHits[3], kTRUE); // 21.04
  300. // fNPass = 2; // 2 prochoda
  301. //fNPass = 2;
  302. fNPass = 3;
  303. return kSUCCESS;
  304. }
  305. //__________________________________________________________________________
  306. void MpdCellAutomat::Reset()
  307. {
  308. cout << " MpdCellAutomat::Reset " << endl;
  309. //fKHits->Clear();
  310. for (Int_t i = 0; i < 4; ++i) {
  311. fKHits[i]->Delete();
  312. f2DHits[i]->Delete();
  313. }
  314. fKHits1->Delete();
  315. fTracks->Delete();
  316. fTrackCand->Delete();
  317. delete [] fLayPointers;
  318. fLayPointers = NULL;
  319. fCellMap.clear();
  320. }
  321. //__________________________________________________________________________
  322. void MpdCellAutomat::SetParContainers()
  323. {
  324. // Get run and runtime database
  325. FairRunAna* run = FairRunAna::Instance();
  326. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  327. FairRuntimeDb* db = run->GetRuntimeDb();
  328. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  329. // Get STS geometry parameter container
  330. db->getContainer("MpdStsGeoPar");
  331. }
  332. //__________________________________________________________________________
  333. void MpdCellAutomat::Finish()
  334. {
  335. //Write();
  336. }
  337. //__________________________________________________________________________
  338. void MpdCellAutomat::Exec(Option_t * option)
  339. {
  340. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  341. static int eventCounter = 0;
  342. cout << " - - - - \n ItsRec event " << ++eventCounter << endl;
  343. Reset();
  344. // Create Kalman hits
  345. if (fItsHits->GetEntriesFast() == 0) return;
  346. //MakeTrackCandidates(); // 21.02
  347. //ExtendCellTracks(); // 21.02
  348. //return;
  349. Build2DHits();
  350. MakeKalmanHits();
  351. for (Int_t i = 0; i < fNPass; ++i) {
  352. fTrackCand->Delete(); // 05.03
  353. fKHits[0]->Delete();// MpdCellTrack
  354. fKHits[1]->Delete();
  355. fKHits[2]->Delete();
  356. fKHits[3]->Delete();
  357. MakeTrackCandidates(i);// 21.02
  358. ExtendCellTracks(i);// 21.02
  359. GetTrackSeeds(i);//attention!!
  360. /*
  361. Int_t nHitsOk = 0, nHits = fKHits->GetEntriesFast();
  362. for (Int_t j = 0; j < nHits; ++j){
  363. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->Unchecked(j);
  364. if (hit->GetFlag() >= 0) ++nHitsOk;
  365. }
  366. */
  367. // cout << " Total number of hits for tracking: " << fKHits->GetEntriesFast() << endl;
  368. cout << " Total number of track seeds: " << fTracks->GetEntriesFast() << endl;
  369. cout << " Total number of track seeds fTrackCand: " << fTrackCand->GetEntriesFast() << endl;
  370. DoTracking(i); //attention!!!
  371. RemoveDoubles();
  372. StoreTracks();
  373. cout << " Total number of found tracks: " << fTrackCand->GetEntriesFast() << endl;
  374. //if (i != fNPass - 1) ExcludeHits(); // exclude used hits work
  375. ExcludeHits(); // exclude used hits test
  376. }
  377. GetShortTracks();
  378. AddHits(); // add hit objects to tracks
  379. cout << " Total number of found tracks: " << fTrackCand->GetEntriesFast() << endl;
  380. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  381. }
  382. //__________________________________________________________________________
  383. void MpdCellAutomat::MakeKalmanHits()
  384. {
  385. /// Create Kalman hits from ITS hits.
  386. Int_t nHits = fItsHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0;
  387. Double_t r, z, xloc, errZ = 0.012, errX = 0.0023; // 120um in Z, 23um in R-Phi (local X)
  388. //Double_t r, z, xloc, errZ = 0.12, errX = 0.0023; // 1.2mm in Z, 23um in R-Phi (local X)
  389. //Double_t r, z, xloc, errZ = 50.0, errX = 0.01; // 50cm in Z, 100um in R-Phi (local X)
  390. for (Int_t ih = 0; ih < nHits; ++ih) {
  391. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(ih);
  392. r = TMath::Sqrt (h->GetX() * h->GetX() + h->GetY() * h->GetY());
  393. z = h->GetZ();
  394. xloc = h->GetLocalX();
  395. //cout << ih << " " << h->Layer()-1 << endl;
  396. //lay = h->Layer() * 2 + h->Side() + 1;
  397. lay = (h->Layer() - 1) * 2 + h->Side();
  398. // Add error
  399. Double_t dX = 0, dZ = 0;
  400. //gRandom->Rannor(dX,dZ);
  401. //if (errZ > 2) dZ = 0.0; // 1-D case
  402. dZ = 0.0; // 1-D case
  403. //Double_t meas[2] = {xloc+dX*errX, z+dZ*errZ};
  404. Double_t meas[2] = {xloc+h->GetDx()*errX, z+dZ*errZ};
  405. Double_t err[2] = {errX, errZ};
  406. Double_t cossin[2] = {TMath::Cos(fStereoA[h->Side()]), TMath::Sin(fStereoA[h->Side()])};
  407. //(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)
  408. //MpdKalmanStripHit *hit = new ((*fKHits)[nKH++]) MpdKalmanStripHit(r, fStereoA[h->Side()],
  409. // xloc+dX*errX, z+dZ*errZ, errX, errZ, lay, ih);
  410. //hit->SetType(MpdKalmanHit::kFixedR);
  411. MpdKalmanHit *hit = 0x0;
  412. //cout << h->GetDetectorID() << " " << fId2Id[lay][h->GetDetectorID()] << endl;
  413. //if (fGeo && h->GetUniqueID()) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+h->GetDetectorID(), 1,
  414. if (fGeo && h->GetUniqueID()) hit =
  415. new ((*fKHits1)[nKH++]) MpdKalmanHit(lay*1000000+fId2Id[lay][h->GetDetectorID()], 1,
  416. MpdKalmanHit::kFixedP, meas, err, cossin, 0., r, ih);
  417. // Mask out sector number - sensor layout
  418. else if (fGeo) hit = new ((*fKHits1)[nKH++]) MpdKalmanHit(lay*1000000+(h->GetDetectorID()&((2<<12)-1)), 1,
  419. MpdKalmanHit::kFixedP, meas, err, cossin, 0., r, ih);
  420. else hit = new ((*fKHits1)[nKH++]) MpdKalmanHit(lay*1000000+nKH-1, 1, MpdKalmanHit::kFixedR,
  421. meas, err, cossin, 0., r, ih);
  422. hit->SetUniqueID(0);
  423. // Add second measurement - just for test at the moment
  424. //!!!
  425. //hit->SetNofDim(2);
  426. //!!!
  427. layMax = TMath::Max (lay, layMax);
  428. }
  429. cout << " Max layer = " << layMax << " " << fKHits1->GetEntriesFast() << endl;
  430. //fKHits1->Sort(); // in descending order in R
  431. //cout << ((MpdKalmanHit*)fKHits->UncheckedAt(0))->GetPos() << endl;
  432. cout << ((MpdKalmanHit*)fKHits1->UncheckedAt(0))->GetDist() << endl;
  433. }
  434. //_________________________________________________________________________
  435. void MpdCellAutomat::Build2DHits()
  436. {
  437. /// Build track candidates from ITS 2-D hits.
  438. fItsHits->Sort();
  439. Int_t nHits = fItsHits->GetEntriesFast(), layMax = 0, lay = 0, nKH = 0;
  440. //Double_t errZ = 0.012, errX = 0.0023; // 120um in Z, 23um in R-Phi (local X)
  441. Double_t errZ = 0.025, errX = 0.0023; // 250um in Z, 23um in R-Phi (local X)
  442. //Double_t r, z, xloc, errZ = 0.12, errX = 0.0023; // 1.2mm in Z, 23um in R-Phi (local X)
  443. //Double_t r, z, xloc, errZ = 50.0, errX = 0.01; // 50cm in Z, 100um in R-Phi (local X)
  444. Double_t xloc[2];
  445. Double_t meas[2];
  446. Double_t r;
  447. Double_t dX = 0, dZ = 0;
  448. //gRandom->Rannor(dX,dZ);
  449. Int_t lay0 = 1; // layer 1
  450. fLayBeg[0] = 0; // index for layer 1
  451. for (Int_t i1 = 0; i1 < nHits; ++i1) {
  452. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(i1);
  453. Int_t lay = h->Layer();
  454. if (lay0 != lay) {
  455. fLayBeg[lay - 1] = i1;
  456. lay0 = lay;
  457. }
  458. }
  459. for (Int_t i234 = 0; i234 < 4; ++i234) { // Loop over layers 1-4
  460. Int_t nKH = 0, iprint = 0;
  461. printf("i234 =%d\n", i234);
  462. printf("h->lay=%d\n",lay);
  463. Int_t hEnd = nHits;
  464. if (i234 != 3) hEnd = fLayBeg[i234 + 1];
  465. for (Int_t ih = fLayBeg[i234]; ih < hEnd; ++ih) {
  466. // for (Int_t ih = fLayBeg[0]; ih < fLayBeg[1]; ++ih) { // work 28.04
  467. // for (Int_t ih = 0; ih < nHits; ++ih) { // old
  468. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(ih);
  469. if (h->GetFlag() < 0) continue;
  470. lay = (h->Layer() - 1) * 2 + h->Side();
  471. r = TMath::Sqrt (h->GetX() * h->GetX() + h->GetY() * h->GetY());
  472. // printf( "r=%d\n",r);
  473. //gRandom->Rannor(dX,dZ);
  474. xloc[0]=h->GetLocalX()+h->GetDx()*errX;
  475. //cout << ih << " " << h->Layer()-1 << endl;
  476. // Double_t dX = 0, dZ = 0;
  477. // gRandom->Rannor(dX,dZ);
  478. // dZ = 0.0;
  479. Double_t err[2] = {errX, errZ};
  480. Double_t cossin[2] = {TMath::Cos(fStereoA[h->Side()]), TMath::Sin(fStereoA[h->Side()])}; //23.08
  481. MpdCellTrack *hit = 0x0;//ms
  482. MpdStsHit *hsts = h;
  483. if (h->Side()) hsts = NULL;
  484. for (Int_t ih1 = ih + 1; ih1 < hEnd; ++ih1) {
  485. // printf("Dim:xloc[0]=%f,xloc[1]=%f,nHits=%d\n",xloc[0],xloc[1],nHits);
  486. MpdStsHit *h1 = (MpdStsHit*) fItsHits->UncheckedAt(ih1);
  487. if (h1->GetFlag() < 0) continue;
  488. if (h->Side() == h1->Side()) continue;
  489. if (h->Layer() != h1->Layer()) continue;
  490. if (h->Ladder() != h1->Ladder()) continue;
  491. if (h->Detector() != h1->Detector()) continue;
  492. //gRandom->Rannor(dX,dZ);
  493. xloc[1]=h1->GetLocalX()+h1->GetDx()*errX;
  494. Double_t cossinn[2] = {TMath::Cos(fStereoA[h1->Side()]), TMath::Sin(fStereoA[h1->Side()])};
  495. if (hsts == NULL) hsts = h1;
  496. //Make 2-dim coordinate
  497. // Double_t x1=((xloc1[1]*sterSin[0]-xloc[0]*sterSin[1])/(sterCos[1]*sterSin[0]-sterCos[0]*sterSin[1]));
  498. Double_t x1=((xloc[1]*cossin[1]-xloc[0]*cossinn[1])/(cossinn[0]*cossin[1]-cossin[0]*cossinn[1]));
  499. //Double_t z1=((xloc[0]-(x1)*sterCos[0])/sterSin[0]);
  500. Double_t z1 = TMath::Abs(cossin[1]) > 1.e-6 ? (xloc[0] - x1 * cossin[0]) / cossin[1] :
  501. (xloc[1] - x1 * cossinn[0]) / cossinn[1];
  502. MpdKalmanHit hitTmp;
  503. //Int_t lay1 = (h1->Layer() - 1) * 2 + hsts->Side();
  504. Int_t lay1 = (hsts->Layer() - 1) * 2 + hsts->Side();
  505. hitTmp.SetDetectorID(lay1*1000000+fId2Id[lay1][hsts->GetDetectorID()]);
  506. Double_t sizeZ = MpdKalmanFilter::Instance()->GetGeo()->Size(&hitTmp).Y();
  507. if (TMath::Abs(z1) > sizeZ) continue;
  508. // printf("sizeZ=%f %d %d %d %d %d %d %d %d %f %f\n",sizeZ,h->Layer(),h1->Layer(),h->Ladder(),h1->Ladder(),h->Detector(),h1->Detector(),h->GetTrackID(),h1->GetTrackID(),x1,z1);
  509. TString path = MpdKalmanFilter::Instance()->GetGeo()->Path(hitTmp.GetDetectorID());
  510. meas[0] = x1;
  511. meas[1] = z1;
  512. gGeoManager->cd(path);
  513. Double_t v3[3] = {meas[0], 0.0, meas[1]}, v31[3] = {0.0};
  514. gGeoManager->LocalToMaster(v3, v31);
  515. if (h->Layer() == 1 && TMath::Abs(v31[2]) > 18.0) continue; // Eta < 1.99
  516. else if ( h->Layer() == 2 && TMath::Abs(v31[2]) > 40.0) continue; // rad 11
  517. else if ( h->Layer() == 3 && TMath::Abs(v31[2]) > 61.0) continue; // rad 17
  518. else if ( h->Layer() == 4 && TMath::Abs(v31[2]) > 83.0) continue; // rad 23
  519. TVector3 point1(v31);
  520. Float_t x1y1[3] = {v31[0],v31[1],v31[2]};
  521. Float_t r1 = TMath::Sqrt(x1y1[1] * x1y1[1] + x1y1[0] * x1y1[0]);
  522. Float_t angL = TMath::ATan(x1y1[2] / r1); //absoliyt angle
  523. /*
  524. ////printf("angL=%f\n",angL);
  525. //(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)
  526. //MpdKalmanStripHit *hit = new ((*fKHits)[nKH++]) MpdKalmanStripHit(r, fStereoA[h->Side()],
  527. // xloc+dX*errX, z+dZ*errZ, errX, errZ, lay, ih);
  528. //hit->SetType(MpdKalmanHit::kFixedR);
  529. //cout << h->GetDetectorID() << " " << fId2Id[lay][h->GetDetectorID()] << endl;
  530. //if (fGeo && h->GetUniqueID()) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+h->GetDetectorID(), 1,
  531. // if (fGeo && h->GetUniqueID()) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay1*1000000+fId2Id[lay1][hsts->GetDetectorID()], 1, MpdKalmanHit::kFixedP,meas, err, cossin, 0., r, ih);//ms
  532. */
  533. if (fGeo && h->GetUniqueID()) hit = new ((*f2DHits[h->Layer()-1])[nKH++]) MpdCellTrack(lay1*1000000+fId2Id[lay1][hsts->GetDetectorID()], 1, MpdCellTrack::kFixedP,point1, err,
  534. cossin, 0., r, ih, ih1, -1); //ms 17.06
  535. // Mask out sector number - sensor layout
  536. // else if (fGeo) hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+(h->GetDetectorID()&((2<<12)-1)), 1, MpdKalmanHit::kFixedP,meas, err, cossin, 0., r, ih);//ms
  537. // else hit = new ((*fKHits)[nKH++]) MpdKalmanHit(lay*1000000+nKH-1, 1, MpdKalmanHit::kFixedR,meas, err, cossin, 0., r, ih); //ms
  538. else if (fGeo) hit = new ((*f2DHits[h->Layer()-1])[nKH++]) MpdCellTrack(lay*1000000+(h->GetDetectorID()&((2<<12)-1)), 1, MpdCellTrack::kFixedP,point1, err, cossin, 0., r, ih, ih1, -1); //ms 17.06
  539. else hit = new ((*f2DHits[h->Layer()-1])[nKH++]) MpdCellTrack(lay*1000000+nKH-1, 1, MpdCellTrack::kFixedR,point1, err, cossin, 0., r, ih, ih1, -1); //ms 17.06
  540. hit->SetUniqueID(0);
  541. // Add second measurement - just for test at the moment
  542. hit->SetNofDim(2);
  543. // hit->SetIndex(ih1);
  544. hit->SetCosSin(1,angL);
  545. layMax = TMath::Max (lay1, layMax);
  546. } // for (Int_t ih1 = 0; ih1 < nHits; ++ih1)
  547. } // for (Int_t ih = 0; ih < nHits; ++ih)
  548. printf("Schetchik:nKH=%d\n",nKH);
  549. } // for (Int_t i234 = 0; i234 < 4; ++i234)
  550. //printf("fLayBeg[1]=%d\n",fLayBeg[0]);
  551. printf("fLayBeg[1]=%d\n",fLayBeg[1]);
  552. printf("fLayBeg[2]=%d\n",fLayBeg[2]);
  553. printf("fLayBeg[3]=%d\n",fLayBeg[3]);
  554. cout << "------------| |--------------" << endl;
  555. cout << "------------------ STOP Build2DHits ---------------------" << endl;
  556. //cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << endl;//ms 06.05
  557. /*
  558. for (Int_t i = 0; i < 4; ++i) {
  559. fKHits[i]->Sort();
  560. }
  561. */
  562. //fKHits->Sort(); // in descending order in R
  563. //cout << ((MpdKalmanHit*)fKHits->UncheckedAt(0))->GetPos() << endl;
  564. //cout << ((MpdCellTrack*)fKHits->UncheckedAt(0))->GetDist() << endl; //ms 06.05
  565. }
  566. //_________________________________________________________________________
  567. void MpdCellAutomat::MakeTrackCandidates(Int_t iPass)
  568. {
  569. // Read 2D Hits and create trackCandidates from 2D Hits at first layer
  570. Int_t nKH = 0, nHits = f2DHits[0]->GetEntriesFast();
  571. for (Int_t i = 0; i < nHits; ++i) {
  572. MpdCellTrack *cellTr = (MpdCellTrack*) f2DHits[0]->UncheckedAt(i);
  573. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(cellTr->GetIndex(0));
  574. if (h->GetFlag() < 0) continue;
  575. h = (MpdStsHit*) fItsHits->UncheckedAt(cellTr->GetIndex(1));
  576. if (h->GetFlag() < 0)continue;
  577. if (fGeo && h->GetUniqueID()) cellTr = new ((*fKHits[0])[nKH++]) MpdCellTrack(*cellTr); /// work!!
  578. else if (fGeo) cellTr = new ((*fKHits[0])[nKH++]) MpdCellTrack(*cellTr); //
  579. else cellTr = new ((*fKHits[0])[nKH++]) MpdCellTrack(*cellTr); //
  580. cellTr->SetCode(i);
  581. } // for (Int_t i = 0; ih < nHits; ++i)
  582. printf("Schetchik:nKH=%d\n",nKH);
  583. cout << "-------STOP MakeTrackCand---------" << endl;
  584. }
  585. //_________________________________________________________________________
  586. // relizz
  587. void MpdCellAutomat::ExtendCellTracks(Int_t iPass)
  588. {
  589. /// Extend cell tracks to layers 2-4
  590. const Int_t nSigm = 3;
  591. Double_t errZ = 0.025, errX = 0.0023; // 250um in Z, 23um in R-Phi (local X)
  592. for (Int_t i234 = 1; i234 < 4; ++i234) { // Loop over layers 1-3
  593. Int_t nKH = 0, iprint = 0;
  594. printf("i234 =%d\n", i234);
  595. Int_t nHits2 = f2DHits[i234]->GetEntriesFast(); //26.05
  596. for (Int_t i = 0; i < nHits2; ++i) {
  597. MpdCellTrack *hit = (MpdCellTrack*) f2DHits[i234]->UncheckedAt(i);
  598. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex(0));
  599. if (h->GetFlag()<0 ) continue; //work
  600. h = (MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex(1));
  601. if (h->GetFlag()<0 )continue;
  602. TVector3 v311 = hit->GetMeas();
  603. Int_t nTracks = fKHits[i234-1]->GetEntriesFast(); //prediduchii sloi fkHits ->massiv ykazateley// work version 6.06.2014
  604. if (iprint == 0) { ++iprint; cout << "nTracks= " << nTracks << endl; }
  605. for (Int_t itr = 0; itr < nTracks; ++itr){
  606. Float_t angT, angL, deltaZ, dangL; //ang1->angT(transverse-"poperechnii" ); ang2->angL(longitudinal-"prodolnii");
  607. MpdCellTrack *track = (MpdCellTrack*) fKHits[i234-1]->UncheckedAt(itr); // work version 6.06.2014
  608. /// printf ("CosSin(0)1111=%f\n",track->GetCosSin(0));
  609. if ( track->GetPrevTrack() < 0) {
  610. // Layer 2
  611. TVector3 v1 = track->GetMeas();
  612. Float_t x1y1[3] = {v311[0] - v1[0], v311[1] - v1[1], v311[2] - v1[2]}; //vector from layer 1 to lay 2
  613. Float_t r1 = TMath::Sqrt(x1y1[1] * x1y1[1] + x1y1[0] * x1y1[0]);
  614. angL = TMath::ATan (x1y1[2] / r1);
  615. dangL = angL - track->GetCosSin(1);
  616. //cout << "dangL: " << dangL << endl;
  617. // if (TMath::Abs(dangL) > 0.047 * nSigm) continue;
  618. if (TMath::Abs(dangL) > 0.045 * nSigm) continue; // sigma lay == 2 pt 0.05 Gev
  619. Float_t phi1 = TMath::ATan2(x1y1[1],x1y1[0]), phi0 = TMath::ATan2(v1[1],v1[0]);
  620. angT = MpdKalmanFilter::Instance()->Proxim(phi0, phi1) - phi0;
  621. //cout << "angT: " << r1 << " " << angT << endl;
  622. if ( iPass == 0 && TMath::Abs(angT) > 0.023) continue; //Pt > 0.4 GeV/c
  623. else if ( iPass == 1 && TMath::Abs(angT) > 0.043) continue; // Pt > 0.2 GeV/c
  624. else if ( iPass == 2 && TMath::Abs(angT) > 0.18) continue; // Pt > 0.05 GeV/c
  625. if (TMath::Abs(dangL) > nSigm * Interp(TMath::Abs(angT),0,i234-1)) continue;
  626. } else {
  627. MpdCellTrack *track1 = (MpdCellTrack*) fKHits[i234-2]->UncheckedAt(track->GetPrevTrack());// trek s previos layer // work version 6.06.2014
  628. TVector3 v1 = track->GetMeas();
  629. TVector3 v5 = (track->GetMeas() - track1->GetMeas()); // track to previous layer
  630. Float_t x1y1[3] = {v311[0] - v1.X(), v311[1] - v1.Y(), v311[2]- v1.Z()}; //vector to this layer
  631. Float_t r1 = TMath::Sqrt(x1y1[1] * x1y1[1] + x1y1[0] * x1y1[0]);
  632. angL = TMath::ATan (x1y1[2] / r1);
  633. //AZ - rough estimate
  634. if (TMath::Abs(angL - track->GetCosSin(1)) > 0.045 * nSigm) continue; // sigma lay == 2 pt 0.05 Gev
  635. Double_t angLmean = angL;
  636. if (h->Layer() == 3) {
  637. // Take average of 2 values
  638. angLmean = (angLmean + track->GetCosSin(1)) / 2;
  639. } else {
  640. // Take average of 3 values
  641. angLmean = (angLmean + 2 * track->GetCosSin(1)) / 3;
  642. }
  643. //dangL = angL - track->GetCosSin(1);
  644. dangL = angLmean - track->GetCosSin(1);
  645. //cout << "dangL: " << dangL << endl;
  646. if (h->Layer() == 3 && TMath::Abs(dangL) > 0.02 * nSigm) continue; // sigma lay == 3 pt 0.05 GeV
  647. else if (h->Layer() == 4 && TMath::Abs(dangL) > 0.016 * nSigm) continue; // sigma lay == 4 pt 0.05 GeV
  648. Float_t phi1 = TMath::ATan2(x1y1[1],x1y1[0]), phi0 = TMath::ATan2(v5[1],v5[0]);
  649. angT = MpdKalmanFilter::Instance()->Proxim(phi0, phi1) - phi0;
  650. //AZ - rough estimate
  651. if ( iPass == 0 && TMath::Abs(angT) > 0.023) continue; //Pt > 0.4 GeV/c
  652. else if ( iPass == 1 && TMath::Abs(angT) > 0.043) continue; // Pt > 0.2 GeV/c
  653. else if ( iPass == 2 && TMath::Abs(angT) > 0.18) continue; // Pt > 0.05 GeV/c
  654. Double_t angTmean = angT;
  655. if (h->Layer() == 3) {
  656. // Take average of 2 values
  657. angTmean = (angTmean + track->GetCosSin(0)) / 2;
  658. } else {
  659. // Take average of 3 values
  660. angTmean = (angTmean + 2 * track->GetCosSin(0)) / 3;
  661. }
  662. //cout << "dangT: " << angTmean - track->GetCosSin(0) << endl;
  663. //if (TMath::Abs(angT - track->GetCosSin(0)) > nSigm * Interp(TMath::Abs(angTmean),1,h->Layer()-3)) continue;///dangt
  664. if (TMath::Abs(angTmean - track->GetCosSin(0)) > nSigm * Interp(TMath::Abs(angTmean),1,i234-2)) continue;///dangt
  665. if (TMath::Abs(dangL) > nSigm * Interp(TMath::Abs(angTmean),0,i234-1)) continue;
  666. angL = angLmean;
  667. angT = angTmean;
  668. }
  669. MpdCellTrack *trNew = NULL;
  670. if (fGeo && h->GetUniqueID()) trNew = new ((*fKHits[h->Layer()-1])[nKH++]) MpdCellTrack(*hit);
  671. else if (fGeo) trNew = new ((*fKHits[h->Layer()-1])[nKH++]) MpdCellTrack(*hit);
  672. else trNew = new ((*fKHits[h->Layer()-1])[nKH++]) MpdCellTrack(*hit);
  673. //printf("nomer treka=%d\n",i);
  674. trNew->SetUniqueID(0);
  675. trNew->SetNofDim(2);
  676. trNew->SetCosSin(0,angT);
  677. trNew->SetCosSin(1,angL);
  678. trNew->SetPrevTrack(itr);
  679. trNew->SetCode(track->GetCode());
  680. trNew->SetCode(i);
  681. //cout << i234 << " " << trNew->GetCode() << endl;
  682. } // for (Int_t itr = 0; itr < nTracks;
  683. } // for (Int_t i = 0; i < nHits2;
  684. printf("nKH=%d\n", nKH);
  685. } // for (Int_t i234 = 1; i234 < 4;)
  686. }
  687. //__________________________________________________________________________
  688. void MpdCellAutomat::GetTrackSeeds(Int_t iPass)
  689. {
  690. /// Build ITS track seeds from Cell tracks
  691. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  692. Int_t nCand = 0;
  693. TVector3 vert(0.,0.,0.),pmom;
  694. Int_t nTracks = fKHits[3]->GetEntriesFast();
  695. cout << "seed ITS tracks:"<< nTracks << " " << fCellMap.size() << endl;
  696. printf ("nTracks56789=%d\n",nTracks);
  697. MpdKalmanHit hit;
  698. hit.SetType(MpdKalmanHit::kFixedR);
  699. for (Int_t itr = 0; itr < nTracks; ++itr) {
  700. MpdCellTrack *track1 = (MpdCellTrack*) fKHits[3]->UncheckedAt(itr); // 4 layer
  701. // Check if this track has not been checked already during earlier passes
  702. if (iPass && fCellMap.find(track1->GetCode()) != fCellMap.end()) {
  703. cout << " Found: " << track1->GetCode() << " " << fCellMap[track1->GetCode()] << endl;
  704. continue;
  705. }
  706. fCellMap.insert(pair<TString,Int_t>(track1->GetCode(),1));
  707. //cout << itr << " " << track1->GetCode() << endl;
  708. MpdCellTrack *track2 = (MpdCellTrack*) fKHits[2]->UncheckedAt(track1->GetPrevTrack()); // index treka 3 layer
  709. MpdCellTrack *track3 = (MpdCellTrack*) fKHits[1]->UncheckedAt(track2->GetPrevTrack()); // index treka 2 layer
  710. MpdCellTrack *track4 = (MpdCellTrack*) fKHits[0]->UncheckedAt(track3->GetPrevTrack()); // 1 layer
  711. ///tpc->GetParam()->Print();
  712. ///MpdEctKalmanTrack *track = new ((*fTracks)[nCand++]) MpdEctKalmanTrack(itr, *tpc);
  713. // MpdItsKalmanTrack *track = new ((*fTracks)[nCand++]) MpdItsKalmanTrack(*track1, vert); //sozdatia constructor iz CellTrack vvodit v ITS track constructor
  714. MpdItsKalmanTrack *track = new ((*fTrackCand)[nCand++]) MpdItsKalmanTrack(*track1, vert); //new 05.03.14
  715. // is Celltrack cozdavalsia constrycror v ITS Track!!!! ATTENTION
  716. ///MpdITSKalmanTrack!!!
  717. Double_t pt = EvalPt(track1,track2);
  718. TVector3 v1 = track1->GetMeas();
  719. TVector3 v2 = track2->GetMeas();
  720. TVector3 v3 = track3->GetMeas();
  721. TVector3 v4 = track4->GetMeas();
  722. /*
  723. Double_t phiOut = v1.Phi(); // get azimuth angle
  724. Double_t phiIn = v2.Phi();
  725. Double_t rOut = v1.Pt();// 26.11 rastoynie ot tochki vzaimod-ya
  726. // Double_t phiOut = hitOut->GetMeas(0) / rOut;
  727. Double_t rIn = v2.Pt(); //26.11 rastoyanie
  728. */
  729. //< blok ot 25.12.2013
  730. Double_t phiOut = v3.Phi();
  731. Double_t phiIn = v4.Phi();
  732. Double_t rOut = v3.Pt();
  733. Double_t rIn = v4.Pt();
  734. //<end blok
  735. track->SetUniqueID(itr+1); ///
  736. ///< parametrs track in TPC:
  737. ///< 0: RPhi - coordinate in R-Phi direction
  738. ///< 1: Z - longitudinal coordinate
  739. ///< 2: Phi - local azimuthal angle
  740. ///< 3: Theta - local dip angle (angle w.r.t. the transverse plane)
  741. ///< 4: q/Pt - signed inverse Pt
  742. /*
  743. track->SetPos(rOut);
  744. track->SetParam (4, 1./pt); // q/Pt
  745. track->SetParam (3, track1->GetCosSin(1)); //longitudinal angL
  746. //track->SetParam (2, track1->GetCosSin(0)); //transverse angT
  747. track->SetParam (2, (v1-v2).Phi()); //transverse angT
  748. track->SetParam (1, v1.Z()); // Z - coordinate
  749. track->SetParam( 0,phiOut*rOut);
  750. */
  751. //< blok ot 25.12.2013
  752. track->SetPos(rIn);//<- ispravleno!
  753. track->SetParam (4, 1./pt); // q/Pt
  754. track->SetParam (3, track1->GetCosSin(1)); //longitudinal angL
  755. //track->SetParam (2, track1->GetCosSin(0)); //transverse angT
  756. track->SetParam (2, (v3-v4).Phi()); //Phi - rough estimate
  757. // Adjust Phi
  758. Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
  759. Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
  760. Double_t rCirc = TMath::Abs (pt / factor);
  761. Double_t ph = TMath::ASin ((v3-v4).Pt() / 2 / rCirc);
  762. track->SetParam (2, track->GetParam(2) - TMath::Sign(ph,pt));
  763. track->SetParam (1, v4.Z()); // Z - coordinate
  764. track->SetParam( 0,phiIn*rIn);
  765. //end blok
  766. //old code!
  767. /*
  768. track->SetParam(*track->GetParamAtHit());
  769. track->SetParamNew(*track->GetParamAtHit());
  770. track->SetPosNew(track->GetPos());
  771. track->SetWeight(*track->GetWeightAtHit());
  772. track->SetLength(track->GetLengAtHit());
  773. if (track.fHits == 0x0) return;
  774. //const MpdCellTrack *hitOut;// ukazatel nelzya nichego bratiat 02/12
  775. // const MpdCellTrack *hitIn;// ukazatel
  776. // Double_t rOut = hitOut->GetPos();// GetPos() rasstoyanie ot tochki vzaimodeistviya? old
  777. */
  778. Double_t parOut[4] = {rOut,phiOut,0.,0.};
  779. Double_t parIn[4] = {rIn,phiIn,0.,0.};
  780. //EvalCovar(parOut,parIn,track,track1); //coment ot 25.12 work version
  781. EvalCovar(parOut,parIn,track,track4); // new ot 25.12 track4 -> 1 layer old
  782. track->SetPosNew(track->GetPos());
  783. track->SetParamNew(*track->GetParam());
  784. //track->ReSetWeight();
  785. //TMatrixDSym w = *track->GetWeight(); // save current weight matrix
  786. hit.SetPos(rIn-0.5);// radius tochki 4 layer
  787. MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  788. //track->SetWeight(w); // restore original weight matrix (near TPC inner shell)
  789. ///cout << nCand-1 << " " << track->GetTrackID() << endl;
  790. ///cout << track->GetHits()->GetEntriesFast() << " " << track->GetTrHits()->GetEntriesFast() << endl;
  791. track->GetHits()->Clear();
  792. track->SetChi2Its(track->GetChi2()); // temporary storage
  793. track->SetChi2(0.);
  794. track->SetDirection(MpdKalmanTrack::kOutward);
  795. }
  796. cout << " Number of ITS track candidates: " << nCand << endl;
  797. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  798. }
  799. //__________________________________________________________________________
  800. Double_t MpdCellAutomat::EvalPt(const MpdCellTrack *track1, const MpdCellTrack *track2)
  801. {
  802. /// Evaluate signed track Pt (curvature) assuming the track coming from the ///eval pt daet odin parametr a tam4
  803. /// primary vertex
  804. TVector3 v1 = track1->GetMeas();///attention!!!
  805. TVector3 v2 = track2->GetMeas();
  806. TVector2 vec1(v1.X(),v1.Y());
  807. TVector2 vec2(v2.X(),v2.Y());
  808. TVector2 vec21 = vec1 - vec2;
  809. Double_t cosAlpha = vec2 * vec21 / vec2.Mod() / vec21.Mod();
  810. Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  811. Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
  812. Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
  813. Double_t phi1 = vec1.Phi();
  814. Double_t phi2 = vec2.Phi();
  815. Double_t charge = phi1 - MpdKalmanFilter::Instance()->Proxim(phi1,phi2);
  816. if (track1->GetLayer() > track2->GetLayer()) charge = -charge;
  817. return factor * TMath::Abs(rad) * TMath::Sign(1., -charge);
  818. }
  819. //__________________________________________________________________________
  820. void MpdCellAutomat::EvalCovar(Double_t *parOut, Double_t *parIn, MpdItsKalmanTrack *track,const MpdCellTrack *track1)
  821. {
  822. /// Evaluate covariance matrix for track seed
  823. Double_t rIn = parIn[0], phiIn = parIn[1];
  824. Double_t rOut = parOut[0], phiOut = parOut[1];
  825. Double_t xIn = rIn * TMath::Cos(phiIn);
  826. Double_t yIn = rIn * TMath::Sin(phiIn);
  827. Double_t xOut = rOut * TMath::Cos(phiOut);
  828. Double_t yOut = rOut * TMath::Sin(phiOut);
  829. parOut[2] = xOut;
  830. parOut[3] = yOut;
  831. parIn[2] = xIn;
  832. parIn[3] = yIn;
  833. TVector3 v4 = track1->GetMeas();
  834. TMatrixD ww(5,5);
  835. ww(0,0) = track1->GetErr(0) * track1->GetErr(0); // <RphiRphi> //26.11 error x
  836. ww(0,0) *= 4.; //extra factor of 4
  837. ww(1,1) = track1->GetErr(1) * track1->GetErr(1); // <zz> //error z
  838. Double_t dx = parOut[2] - parIn[2], dy = parOut[3] - parIn[3];
  839. Double_t dist2 = dx * dx + dy * dy;
  840. Double_t sinPhi = TMath::Sin (track->GetParam(2));
  841. Double_t cosPhi = TMath::Cos (track->GetParam(2));
  842. Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
  843. Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
  844. ww(2,2) = (pOut * pOut + pIn * pIn) / dist2 * ww(0,0); // <PhiPhi>
  845. ww(2,2) *= 2.; // extra factor of 2
  846. Double_t tanThe = TMath::Tan(track->GetParam(3));
  847. Double_t dRad = parOut[0] - parIn[0];
  848. Double_t denom = dRad * (1.+tanThe*tanThe);
  849. ww(3,3) = ww(1,1) * 2. / denom / denom; // <TheThe>
  850. ww(3,3) *= 1;
  851. ww(1,1) *= 8.; //AZ extra factor of 8
  852. //ww(4,4) = (track->GetParam(4)*0.5) * (track->GetParam(4)*0.5); // error 50%
  853. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
  854. ww(4,4) = (track->GetParam(4)*1.) * (track->GetParam(4)*1.); // error 100%
  855. //fWeight->Print();
  856. //fWeight->Invert(); // weight matrix
  857. Int_t iok = 0;
  858. TMatrixD wwTmp = ww;
  859. MpdKalmanFilter::Instance()->MnvertLocal(ww.GetMatrixArray(), 5, 5, 5, iok);
  860. track->SetWeight(ww);
  861. //fWeight->Print();
  862. // Obtain errors
  863. //if (lunErr && idEct == tofP->GetTrackID()) fprintf(lunErr,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e \n",track->GetParam(2),MpdKalmanFilter::Instance()->Proxim(track->GetParam(2),pmom.Phi()),track->GetParam(3),MpdKalmanFilter::Instance()->Proxim(track->GetParam(3),TMath::PiOver2()-pmom.Theta()),1./TMath::Sqrt((*track->GetWeight())(2,2)),1./TMath::Sqrt((*track->GetWeight())(3,3)),pmom.Mag(),rTof,pmom.Pt(),pt);
  864. TVector3 pos, posOut, pmom;
  865. // if (lunErr) fprintf(lunErr,"%-12.8s %-12.5s %-10.5s %-10.8s %-10.9s\n","|r*phi|","| Y |","| Z |","| Phi |","| Theta |");
  866. MpdStsHit *hit = (MpdStsHit*) fItsHits->UncheckedAt(track1->GetIndex());
  867. MpdStsPoint *p = (MpdStsPoint*)fItsPoints->UncheckedAt(hit->GetRefIndex());
  868. //hit->Position(pos);
  869. p->Position(pos);
  870. p->PositionOut(posOut);
  871. p->Momentum(pmom);
  872. Double_t phi = pmom.Phi();
  873. Double_t th = pmom.Theta();
  874. Double_t r = (pos.Pt() + posOut.Pt()) / 2;//radius
  875. Double_t phi1 = pos.Phi();
  876. Double_t phi2 = posOut.Phi();
  877. phi1 += MpdKalmanFilter::Instance()->Proxim(phi1, phi2);
  878. phi1 /= 2;
  879. Double_t zzz = (pos.Z() + posOut.Z()) / 2;
  880. /*
  881. if (lunErr) fprintf(lunErr,"%10.3f %10.3f %10.3f %10.3f %10.3f\n",r*phi1,TMath::PiOver2()-pmom.Theta(),
  882. pos.Z(),phi,th);
  883. */
  884. /*
  885. if (lunErr) fprintf(lunErr,"%10.15s %10.15s %10.15s %10.15s \n","Param(0)-r*phi","Param(1)-z",
  886. "Param(2)-Phi","Param(3)-Theta");
  887. if (lunErr&&p->GetTrackID()<1) fprintf(lunErr,"%12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e\n",
  888. r*(MpdKalmanFilter::Instance()->Proxim(phi1,track->GetParam(0)/r)-phi1),
  889. track->GetParam(1)-zzz,MpdKalmanFilter::Instance()->Proxim(phi,track->GetParam(2))-phi,
  890. TMath::PiOver2()-pmom.Theta()-track->GetParam(3),
  891. TMath::Sqrt(wwTmp(0,0)),TMath::Sqrt(wwTmp(1,1)),TMath::Sqrt(wwTmp(2,2)),TMath::Sqrt(wwTmp(3,3)), phi); //9 parametr "phi" writing in "file.txt"
  892. */
  893. if (lunErr&&p->GetTrackID()<1) fprintf(lunErr,"%12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e\n",
  894. r*(MpdKalmanFilter::Instance()->Proxim(phi1,track->GetParam(0)/r)-phi1),
  895. track->GetParam(1)-zzz,MpdKalmanFilter::Instance()->Proxim(phi,track->GetParam(2))-phi,
  896. TMath::PiOver2()-pmom.Theta()-track->GetParam(3),
  897. TMath::Sqrt(wwTmp(0,0)),TMath::Sqrt(wwTmp(1,1)),TMath::Sqrt(wwTmp(2,2)),TMath::Sqrt(wwTmp(3,3)));
  898. /*
  899. Int_t nHits = fItsHits->GetEntriesFast();
  900. for (Int_t i = 0; i < nHits; ++i){
  901. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(i);
  902. MpdStsPoint *h1 = (MpdStsPoint*) fItsPoints->UncheckedAt(h->GetRefIndex());
  903. ///primer coda
  904. }
  905. //if (lunErr) fclose(lunErr);
  906. //exit(0);
  907. */
  908. }
  909. //__________________________________________________________________________
  910. void MpdCellAutomat::DoTracking(Int_t iPass)
  911. {
  912. /// Run Kalman tracking
  913. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  914. Double_t vert[3] = {0.0,0.0,0.0};
  915. //Int_t nCand = fTracks->GetEntriesFast(), iok = 0;// old
  916. Int_t nCand = fTrackCand->GetEntriesFast(), iok = 0;// new 05.03.14
  917. Int_t lay0 = ((MpdKalmanHit*)fKHits1->First())->GetLayer();//ms 06.05
  918. for (Int_t i = 0; i < nCand; ++i) {
  919. // MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);// old
  920. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);// new 05.03.14
  921. //cout << " Track seed No. " << i << ", ID: " << track->GetTrackID() << ", Hits: " << track->GetNofTrHits() << endl;
  922. for (Int_t j = 0; j < track->GetNofTrHits(); ++j) {
  923. MpdKalmanHit *h = (MpdKalmanHit* )track->GetTrHits()->UncheckedAt(j);
  924. //MpdStsHit *hh = (MpdStsHit*) fItsHits->UncheckedAt(h->GetIndex());
  925. //Int_t id = ((FairMCPoint*) fItsPoints->UncheckedAt(hh->GetRefIndex()))->GetTrackID();
  926. //cout << j << " " << h->GetDist() << " " << h->GetLayer() << endl;
  927. }
  928. if (fGeo) iok = RunKalmanFilterCell(track); // from cell track
  929. //if (fGeo) iok = RunKalmanFilterCyl(track, lay0); // modular geometry
  930. // else iok = RunKalmanFilterCyl(track, lay0); // cylindrical geometry
  931. if (iok == -1) {
  932. // fTracks->RemoveAt(i); //old
  933. fTrackCand->RemoveAt(i); //new 05.03.14
  934. continue;
  935. }
  936. // Mark hits as being used
  937. /*
  938. TObjArray *hits = track->GetHits();// track(MPdItsKalmantrack)
  939. for (Int_t j = 0; j < nHits; ++j) {
  940. hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  941. hit->SetFlag(-1);
  942. }
  943. */
  944. //if (track->GetNofHits() == 0) continue; // no hits added
  945. // Propagate track to the beam line
  946. track->SetParam(*track->GetParamNew());
  947. track->SetPos(track->GetPosNew());
  948. Double_t pos = track->GetPos();
  949. TMatrixD par = *track->GetParam();
  950. TMatrixDSym cov = *track->Weight2Cov();
  951. Double_t leng = track->GetLength();
  952. TString nodeNew = track->GetNodeNew();
  953. //cout << " 1: " << nodeNew << ", " << track->GetNode() << endl;
  954. // Go to beam pipe
  955. MpdKalmanHit hit;
  956. hit.SetType(MpdKalmanHit::kFixedR);
  957. hit.SetPos(fPipeR);
  958. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  959. if (iok != 1) {
  960. // Restore track
  961. track->SetParam(par);
  962. track->SetParamNew(par);
  963. track->SetCovariance(cov);
  964. track->ReSetWeight();
  965. track->SetPos(pos);
  966. track->SetPosNew(pos);
  967. track->SetLength(leng);
  968. //track->SetNode(node);
  969. //cout << " 2: " << nodeNew << ", " << track->GetNode() << endl;
  970. track->SetNodeNew(nodeNew);
  971. } else {
  972. // Add multiple scattering
  973. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  974. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  975. TMatrixDSym* pcov = track->Weight2Cov();
  976. Double_t th = track->GetParamNew(3);
  977. Double_t cosTh = TMath::Cos(th);
  978. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  979. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  980. (*pcov)(3,3) += angle2;
  981. Int_t ok = 0;
  982. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  983. track->SetWeight(*pcov);
  984. }
  985. cov = *track->Weight2Cov();
  986. hit.SetPos(0.);
  987. hit.SetMeas(0,track->GetParam(2)); // track Phi
  988. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  989. //Double_t pos = ((MpdKalmanHit*)track->GetHits()->Last())->GetPos();
  990. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  991. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  992. if (iok != 1) MpdKalmanFilter::Instance()->FindPca(track, vert);
  993. //track->SetPos(pos); // restore position
  994. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  995. //track->GetCovariance()->Print();
  996. } // for (Int_t i = 0; i < nCand;
  997. //fTracks->Compress(); //old
  998. fTrackCand->Compress(); // new 05.03.14
  999. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__); //04.03
  1000. }
  1001. //__________________________________________________________________________
  1002. Int_t MpdCellAutomat::RunKalmanFilterCell(MpdItsKalmanTrack *track)
  1003. {
  1004. /// Run Kalman filter (fitter) for the hits from the cell track
  1005. /// (might not work when propagating outward!!!)
  1006. Int_t layMax = ((MpdKalmanHit*)fKHits1->First())->GetLayer();
  1007. MpdKalmanHit *hitOK = 0x0;
  1008. MpdKalmanHit hitTmp;
  1009. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  1010. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1; /// old coment
  1011. Int_t layEnd = -1, dLay = -1, layOK = -1;
  1012. if (trackDir == MpdKalmanTrack::kOutward) {
  1013. layEnd = layMax + 1;
  1014. dLay = 1;
  1015. }
  1016. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  1017. TMatrixD param(5,1), paramTmp(5,1);
  1018. Double_t saveZ = 0.0, saveLeng = 0.0, dChi2Min = 0.0, posNew = 0.0;
  1019. Int_t ok = 0;
  1020. MpdCellTrack *cellTr = NULL;
  1021. MpdCellTrack *track1 = (MpdCellTrack*) fKHits[3]->UncheckedAt(track->GetUniqueID()-1); //4
  1022. MpdCellTrack *track2 = (MpdCellTrack*) fKHits[2]->UncheckedAt(track1->GetPrevTrack()); //3
  1023. MpdCellTrack *track3 = (MpdCellTrack*) fKHits[1]->UncheckedAt(track2->GetPrevTrack()); //2
  1024. MpdCellTrack *track4 = (MpdCellTrack*) fKHits[0]->UncheckedAt(track3->GetPrevTrack()); //1
  1025. for (Int_t lay = 0; lay < 4; ++lay) {
  1026. // Get CellTrack from the outermost layer
  1027. if (lay == 0) cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(track3->GetPrevTrack()); // 1 layer
  1028. else if (lay == 1) cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(track2->GetPrevTrack()); // 2 layer
  1029. if (lay == 2) cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(track1->GetPrevTrack()); // 3 layer
  1030. else if (lay ==3) cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(track->GetUniqueID()-1); // 4 layer
  1031. /*
  1032. for (Int_t lay = 3; lay > -1; --lay) {
  1033. // Get CellTrack from the outermost layer
  1034. if (lay == 3) cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(track->GetUniqueID()-1); // 4 layer
  1035. else cellTr = (MpdCellTrack*) fKHits[lay]->UncheckedAt(cellTr->GetPrevTrack()); // previous layer
  1036. */
  1037. for (Int_t ihit = 0; ihit < 2; ++ihit) {
  1038. // Loop over 2 sides of one layer
  1039. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits1->UncheckedAt(cellTr->GetIndex(ihit));
  1040. // Exclude used hits
  1041. // if (hit->GetFlag() < 0) continue; //04.03
  1042. // Propagate to hit (if it is not very close to the track)
  1043. if (TMath::Abs(hit->GetPos()-track->GetPosNew()) > 1.e-4) {
  1044. Double_t leng = track->GetLength();
  1045. Double_t posNew = track->GetPosNew();
  1046. TMatrixD parNew = *track->GetParamNew();
  1047. TString nodeNew = track->GetNodeNew();
  1048. TString curPath = track->GetNode();
  1049. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit,kTRUE,kTRUE)) {
  1050. // Restore initial parameters for the failed track ///vostanovlenie ischodnich parametrov for bad track
  1051. track->SetPosNew(posNew);
  1052. track->SetParamNew(parNew);
  1053. track->SetLength(leng);
  1054. track->SetNodeNew(nodeNew);
  1055. track->SetNode(curPath);
  1056. ok = -1;
  1057. break;
  1058. }
  1059. Double_t step = track->GetLength() - leng;
  1060. Int_t lay1 = hit->GetLayer();
  1061. //if (lay1 % 2 == 0 && step > 1.e-4) {
  1062. if (lay1 % 2 != 0 && step > 1.e-4) {
  1063. // Crossing silicon layer - add mult. scat. in the sensor
  1064. Double_t x0 = 9.36; // rad. length
  1065. TMatrixDSym *cov = track->Weight2Cov();
  1066. Double_t th = track->GetParamNew(3);
  1067. Double_t cosTh = TMath::Cos(th);
  1068. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, step);
  1069. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl; ///old coment
  1070. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1071. (*cov)(3,3) += angle2;
  1072. Int_t iok = 0;
  1073. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1074. track->SetWeight(*cov);
  1075. } else if (0 && lay1 % 2 != 0 && step > 1.e-4 && fCables[lay1/2].size() > 0) {
  1076. // Crossing silicon layer - add mult. scat. in the cable
  1077. Double_t nCables = 0, x0 = 0.0116 * 2 / 9.36; // in rad. length - 116um cable per side
  1078. // Find number of cables crossed //poisk nomera cables skrechenich
  1079. TString path = gGeoManager->GetPath();
  1080. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1081. cout << " !!! MpdCellAutomat::RunKalmanFilter - Outside detector !!! " << endl;
  1082. exit(0);
  1083. }
  1084. Double_t v7[3] = {track->GetParamNew(0), track->GetPosNew(), track->GetParamNew(1)}, v77[3];
  1085. gGeoManager->LocalToMaster(v7,v77);
  1086. Double_t zTr = TMath::Abs (v77[2]); // global Z
  1087. //cout << zTr << endl;
  1088. map<Double_t,Double_t>::iterator it;
  1089. for (it = fCables[lay1/2].begin(); it != fCables[lay1/2].end(); ++it) {
  1090. if (zTr < it->first || zTr > it->second) continue;
  1091. ++nCables;
  1092. }
  1093. //cout << " Cables: " << nCables << endl;
  1094. if (nCables) {
  1095. x0 *= nCables;
  1096. TMatrixDSym *cov = track->Weight2Cov();
  1097. Double_t th = track->GetParamNew(3);
  1098. Double_t cosTh = TMath::Cos(th);
  1099. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0);
  1100. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  1101. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1102. (*cov)(3,3) += angle2;
  1103. Int_t iok = 0;
  1104. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1105. track->SetWeight(*cov);
  1106. }
  1107. }
  1108. }
  1109. //cout << hit->GetDetectorID() << endl;
  1110. // Exclude used hits!
  1111. //if (hit->GetFlag() != 1) continue;// ? 04.03 coment work version
  1112. if (hit->GetFlag() < 0) continue; //new 04.03 work 15.04
  1113. // !!! Exact ID match
  1114. if (fExact && TrackID(hit) != track->GetTrackID()) continue;
  1115. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterStripLocal(track,hit,pointWeight,param,posNew);
  1116. // Add Z-contribution (if track is outside the detector)
  1117. Double_t sizeZ = MpdKalmanFilter::Instance()->GetGeo()->Size(hit).Y();
  1118. //if (TMath::Abs(branchTr->GetParamNew(1)) > sizeZ) {
  1119. if (TMath::Abs(param(1,0)) > sizeZ) {
  1120. // Outside detector
  1121. //Double_t dChi2z = (TMath::Abs(branchTr->GetParamNew(1)) - sizeZ) / hit->GetErr(1);
  1122. Double_t dChi2z = (TMath::Abs(param(1,0)) - sizeZ) / hit->GetErr(1);
  1123. dChi2 += dChi2z * dChi2z;
  1124. }
  1125. //if (hit->GetNofDim() == 1) cout << " lay, c2, id1: " << lay << " " << dChi2 << " " << TrackID(hit) << endl;
  1126. //else cout << " lay, c2, id1, id2: " << lay << " " << dChi2 << " " << TrackID(hit) << " " << TrackID(hit,1) << endl;
  1127. if (TMath::Abs(dChi2) < fgkChi2Cut) {
  1128. //if (dChi2 < fgkChi2Cut && lay % 2 == 0) {
  1129. track->GetHits()->Add(hit);
  1130. track->SetChi2(track->GetChi2()+dChi2);
  1131. TMatrixDSym w = *track->GetWeight();
  1132. w += pointWeight;
  1133. track->SetWeight(w);
  1134. track->SetParamNew(param);
  1135. // Save track params at last hit
  1136. track->SetLengAtHit(track->GetLength());
  1137. track->SetParamAtHit(param);
  1138. track->SetWeightAtHit(*track->GetWeight());
  1139. }
  1140. Int_t lay1 = hit->GetLayer();
  1141. if (lay1 % 2 != 0 && fCables[lay1/2].size() > 0) {
  1142. // Crossing silicon layer - add mult. scat. in the cable
  1143. Double_t nCables = 0, x0 = 0.0116 * 2 / 9.36; // in rad. length - 116um cable per side
  1144. // Find number of cables crossed //poisk nomera cables skrechenich
  1145. TString path = gGeoManager->GetPath();
  1146. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1147. cout << " !!! MpdCellAutomat::RunKalmanFilter - Outside detector !!! " << endl;
  1148. exit(0);
  1149. }
  1150. Double_t v7[3] = {track->GetParamNew(0), track->GetPosNew(), track->GetParamNew(1)}, v77[3];
  1151. gGeoManager->LocalToMaster(v7,v77);
  1152. Double_t zTr = TMath::Abs (v77[2]); // global Z
  1153. //cout << zTr << endl;
  1154. map<Double_t,Double_t>::iterator it;
  1155. for (it = fCables[lay1/2].begin(); it != fCables[lay1/2].end(); ++it) {
  1156. if (zTr < it->first || zTr > it->second) continue;
  1157. ++nCables;
  1158. }
  1159. //cout << " Cables: " << nCables << endl;
  1160. if (nCables) {
  1161. x0 *= nCables;
  1162. TMatrixDSym *cov = track->Weight2Cov();
  1163. Double_t th = track->GetParamNew(3);
  1164. Double_t cosTh = TMath::Cos(th);
  1165. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0);
  1166. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  1167. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1168. (*cov)(3,3) += angle2;
  1169. Int_t iok = 0;
  1170. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1171. track->SetWeight(*cov);
  1172. }
  1173. }
  1174. } // for (Int_t ihit = 0; ihit < 2;
  1175. } // for (Int_t lay = 3; lay > -1;
  1176. return 0;
  1177. }
  1178. //__________________________________________________________________________
  1179. void MpdCellAutomat::RemoveDoubles()
  1180. {
  1181. /// Remove double tracks (keep the ones with better quality)
  1182. Int_t ntracks = fTrackCand->GetEntriesFast(); //new 05.03.14
  1183. cout << " Total tracks: " << ntracks << endl;
  1184. MpdItsKalmanTrack *tr1, *tr2;
  1185. MpdCellTrack *cellTr = NULL;
  1186. TString code;
  1187. for (Int_t i = 0; i < ntracks; i++) {
  1188. tr1 = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1189. if (tr1 == 0x0) continue;
  1190. for (Int_t j = i+1; j < ntracks; j++) { // j = 0 -> j = i+1
  1191. //if (j == i) continue; // add comment
  1192. tr2 = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(j);
  1193. if (tr2 == 0x0) continue;
  1194. //Int_t nHitsCommon = GetNofCommonHits(tr1, tr2);
  1195. //if ((float)nHitsCommon / TMath::Min(tr1->GetNofHits(),tr2->GetNofHits()) < 0.5) continue;
  1196. if (!AreTracksDoubles(tr1, tr2)) continue;
  1197. if (tr2->GetNofHits() < tr1->GetNofHits()) {
  1198. fTrackCand->RemoveAt(j);
  1199. cellTr = (MpdCellTrack*) fKHits[3]->UncheckedAt(tr2->GetUniqueID()-1);
  1200. code = cellTr->GetCode();
  1201. fCellMap[code] = -1; // used hit combination
  1202. //cout << " Removed1: " << code << " " << tr2->GetNofHits() << endl;
  1203. } else {
  1204. if ((tr2->GetNofHits() > tr1->GetNofHits()) || (tr2->GetChi2() < tr1->GetChi2())) {
  1205. fTrackCand->RemoveAt(i);
  1206. cellTr = (MpdCellTrack*) fKHits[3]->UncheckedAt(tr1->GetUniqueID()-1);
  1207. code = cellTr->GetCode();
  1208. fCellMap[code] = -1; // used hit combination
  1209. //cout << " Removed2: " << code << " " << tr1->GetNofHits() << endl;
  1210. break;
  1211. } else {
  1212. fTrackCand->RemoveAt(j);
  1213. cellTr = (MpdCellTrack*) fKHits[3]->UncheckedAt(tr2->GetUniqueID()-1);
  1214. code = cellTr->GetCode();
  1215. fCellMap[code] = -1; // used hit combination
  1216. //cout << " Removed3: " << code << " " << tr2->GetNofHits() << endl;
  1217. }
  1218. }
  1219. } // for j
  1220. } // for i
  1221. fTrackCand->Compress();
  1222. fNTracks = fTrackCand->GetEntriesFast();
  1223. }
  1224. //__________________________________________________________________________
  1225. Bool_t MpdCellAutomat::AreTracksDoubles(MpdItsKalmanTrack *tr1, MpdItsKalmanTrack *tr2)
  1226. {
  1227. /// Searching common hits in 2 tracks to determine doubles
  1228. // track1 contains fewer hits than track2
  1229. MpdItsKalmanTrack *track1, *track2;
  1230. if (tr1->GetNofHits() > tr2->GetNofHits()) //12.02
  1231. track1 = tr2, track2 = tr1;
  1232. else
  1233. track1 = tr1, track2 = tr2;
  1234. //Int_t limCommonPoint = (track1->GetNofHits()+1) / 2; // at least many common hits should be found //12.02
  1235. Int_t limCommonPoint = 3; // new 25.03
  1236. TObjArray *hits1 = track1->GetHits(), *hits2 = track2->GetHits();//12.02
  1237. Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast(), nHitsCommon = 0, j = 0;
  1238. for (Int_t i = 0; i < nh1; i++){
  1239. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
  1240. for ( ; j < nh2; j++){
  1241. MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
  1242. //is hit common for two tracks compared
  1243. if (hit1 == hit2) {
  1244. nHitsCommon++;
  1245. break;
  1246. }
  1247. //if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam
  1248. if (hit2->GetLayer() > hit1->GetLayer()) break; // already farther from beam
  1249. }
  1250. if (i+limCommonPoint-nHitsCommon > nh1) return kFALSE; // there'll be not enough common hits already
  1251. }
  1252. //if count of common hits is greater limit
  1253. if (nHitsCommon < limCommonPoint) return kFALSE;
  1254. // Test
  1255. //MpdCellTrack *cell1 = (MpdCellTrack*) fKHits[3]->UncheckedAt(tr1->GetUniqueID()-1);
  1256. //MpdCellTrack *cell2 = (MpdCellTrack*) fKHits[3]->UncheckedAt(tr2->GetUniqueID()-1);
  1257. //cout << nHitsCommon << " " << cell1->GetCode() << " " << cell2->GetCode() << endl;
  1258. return kTRUE;
  1259. }
  1260. //__________________________________________________________________________
  1261. Int_t MpdCellAutomat::TrackID(MpdKalmanHit *hit, Int_t indx)
  1262. {
  1263. /// Return track ID of the hit
  1264. FairHit *h = (FairHit*) fItsHits->UncheckedAt(hit->GetIndex(indx));
  1265. return ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1266. }
  1267. //__________________________________________________________________________
  1268. TVector2 MpdCellAutomat::GetDistance(MpdKalmanTrack *track, MpdKalmanHit *hit)
  1269. {
  1270. /// Compute distance between track and hit
  1271. Int_t lay = hit->GetLayer();
  1272. Int_t lay2 = lay / 2;
  1273. Int_t side = lay % 2;
  1274. Int_t module = ((MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex()))->Module();
  1275. Double_t zTr = track->GetParamNew(1);
  1276. Double_t zloc = zTr + fDz[lay2];
  1277. Int_t modTr = Int_t (zloc / fZmod[lay2]);
  1278. Int_t dMod = modTr - module;
  1279. Double_t dZ = 0;
  1280. if (dMod != 0) {
  1281. // Not in the same module - check Z-distance to module edge
  1282. if (dMod > 0) dZ = zloc - (module+1) * fZmod[lay2];
  1283. else dZ = zloc - module * fZmod[lay2];
  1284. if (TMath::Abs(dMod) > 2) return TVector2(TMath::Abs(dZ),999.); // not in neighbour modules
  1285. }
  1286. // Translate transverse track position to local system
  1287. Double_t xloc = track->GetParamNew(0) * hit->GetCosSin(0) + zTr * hit->GetCosSin(1);
  1288. Double_t phTr = xloc / track->GetPosNew();
  1289. Double_t phHit = hit->GetMeas(0) / fRad[lay];
  1290. //TVector2 dist(TMath::Abs(dZ), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit)-phTr));
  1291. TVector2 dist(TMath::Abs(dZ),
  1292. TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit,hit->GetCosSin(0))-phTr));
  1293. //TVector2 dist(TMath::Abs(zTr-hit->GetZ()), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(track,hit)/hit->GetR()-ph));
  1294. return dist;
  1295. }
  1296. //__________________________________________________________________________
  1297. void MpdCellAutomat::Write()
  1298. {
  1299. /// Write
  1300. TFile histoFile("ItsRec.root","RECREATE");
  1301. Writedir2current(fHistoDir);
  1302. histoFile.Close();
  1303. }
  1304. //__________________________________________________________________________
  1305. void MpdCellAutomat::Writedir2current( TObject *obj )
  1306. {
  1307. /// Write
  1308. if( !obj->IsFolder() ) obj->Write();
  1309. else{
  1310. TDirectory *cur = gDirectory;
  1311. TDirectory *sub = cur->mkdir(obj->GetName());
  1312. sub->cd();
  1313. TList *listSub = ((TDirectory*)obj)->GetList();
  1314. TIter it(listSub);
  1315. while( TObject *obj1=it() ) Writedir2current(obj1);
  1316. cur->cd();
  1317. }
  1318. }
  1319. //__________________________________________________________________________
  1320. void MpdCellAutomat::StoreTracks()
  1321. {
  1322. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1323. /// Transfer tracks from fTrackCand to fTracks
  1324. Int_t nFound = fTracks->GetEntriesFast();
  1325. for (Int_t i = 0; i < fNTracks; ++i) {
  1326. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1327. //
  1328. if (track->GetNofHits() < 7) continue; // !!! store only long tracks
  1329. //
  1330. track->Weight2Cov();
  1331. new ((*fTracks)[nFound++]) MpdItsKalmanTrack(*track);
  1332. MpdCellTrack *cellTr = (MpdCellTrack*) fKHits[3]->UncheckedAt(track->GetUniqueID()-1);
  1333. TString code = cellTr->GetCode();
  1334. fCellMap[code] = -1; // used hit combination
  1335. //fTrackCand->RemoveAt(i);
  1336. }
  1337. fNTracks = 0;
  1338. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1339. }
  1340. //__________________________________________________________________________
  1341. void MpdCellAutomat::AddHits()
  1342. {
  1343. /// Add hit objects to tracks and compute number of wrongly assigned hits
  1344. /// (hits with ID different from ID of starting TPC track)
  1345. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1346. Int_t nFound = fTracks->GetEntriesFast();//work 07
  1347. //Int_t nFound = fTrackCand->GetEntriesFast();
  1348. for (Int_t i = 0; i < nFound; ++i) {
  1349. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);//work 07
  1350. // MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1351. cout << track->GetNode() << " " << track->GetNodeNew() << endl;
  1352. Double_t c2 = track->GetChi2();
  1353. track->SetChi2(track->GetChi2Its());
  1354. track->SetChi2Its(c2);
  1355. Int_t nHits = track->GetNofHits();
  1356. if(nHits < 7) continue; // candidate not write track
  1357. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1358. TClonesArray &trHits = *track->GetTrHits();
  1359. TObjArray *hits = track->GetHits();
  1360. SetTrackID(track); // set track ID as ID of majority of its hits
  1361. cout << nHits << " " << trHits.GetEntriesFast() << " " << track->GetTrackID() << endl;
  1362. Int_t nWrong = 0, nMirr = 0, motherID = track->GetTrackID();
  1363. // Get track mother ID
  1364. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  1365. while (mctrack->GetMotherId() >= 0) {
  1366. motherID = mctrack->GetMotherId();
  1367. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1368. }
  1369. Int_t lastIndx = trHits.GetEntriesFast();
  1370. for (Int_t j = 0; j < nHits; ++j) {
  1371. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1372. hit->SetUniqueID(1); // flag ITS hits
  1373. new (trHits[lastIndx+j]) MpdKalmanHit(*hit);
  1374. cout << " " << hit->GetLayer();
  1375. Int_t nOver = hit->Index()->GetSize();
  1376. for (Int_t iov = 0; iov < nOver; ++iov) {
  1377. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex(iov));
  1378. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1379. cout << "-" << motherID1;
  1380. // Get point mother ID
  1381. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID1);
  1382. while (mctrack->GetMotherId() >= 0) {
  1383. motherID1 = mctrack->GetMotherId();
  1384. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1385. }
  1386. if (motherID1 != motherID) ++nWrong;
  1387. }
  1388. }
  1389. if (nHits) cout << "\n Wrongs: " << nWrong << endl;
  1390. track->SetNofWrong(nWrong);
  1391. //MpdKalmanTrack *tpc = (MpdKalmanTrack*) fTpcTracks->UncheckedAt(track->GetUniqueID()-1);
  1392. //track->SetChi2(track->GetChi2()+tpc->GetChi2());
  1393. //track->SetLastLay();
  1394. //track->GetParam()->Print();
  1395. track->SetNofHits(track->GetNofTrHits()); // TPC and ITS hits
  1396. track->SetNofIts(nHits);
  1397. cout << nHits << " " << track->GetNofTrHits() << " " << track->GetTrackID() << " "
  1398. << track->GetChi2Its() << " " << track->GetChi2() << endl;
  1399. }
  1400. fTracks->Compress();//work 07
  1401. //fTrackCand->Compress();//??
  1402. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1403. }
  1404. //__________________________________________________________________________
  1405. void MpdCellAutomat::SetTrackID(MpdItsKalmanTrack* track)
  1406. {
  1407. /// Set track ID as ID of majority of its hits
  1408. const Int_t idMax = 9999;
  1409. vector<Int_t> ids(idMax);
  1410. Int_t nHits = track->GetNofHits(), locMax = 0, size = idMax, id0 = -1;
  1411. TObjArray *hits = track->GetHits();
  1412. for (Int_t i = 0; i < nHits; ++i) {
  1413. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  1414. Int_t id = GetHitID(hit);
  1415. if (i == 0) id0 = id; // first hit
  1416. if (id >= size) { size = id + 1; ids.resize(size); }
  1417. ++ids[id];
  1418. if (ids[id] > ids[locMax]) locMax = id;
  1419. }
  1420. if (ids[locMax] > 1) track->SetTrackID(locMax);
  1421. else track->SetTrackID(id0);
  1422. }
  1423. //__________________________________________________________________________
  1424. Int_t MpdCellAutomat::GetHitID(MpdKalmanHit *hit)
  1425. {
  1426. /// Get hit ID from MCPoint ID
  1427. Int_t nOver = hit->Index()->GetSize();
  1428. for (Int_t iov = 0; iov < nOver; ++iov) {
  1429. MpdStsHit *h = (MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex(iov));
  1430. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1431. return motherID1;
  1432. }
  1433. }
  1434. //__________________________________________________________________________
  1435. void MpdCellAutomat::ExcludeHits()
  1436. {
  1437. /// Exclude hits, already used for tracking, from consideration during the next passes
  1438. Int_t nReco = fTracks->GetEntriesFast();// work 07
  1439. //Int_t nReco = fTrackCand->GetEntriesFast();
  1440. cout << " nReco: " << nReco << endl;
  1441. for (Int_t i = 0; i < nReco; ++i) {
  1442. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);// work 07
  1443. //MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1444. Int_t nhitsKF = track->GetNofHits();
  1445. TObjArray *hits = track->GetHits();
  1446. for (Int_t j = 0; j < nhitsKF; ++j) {
  1447. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1448. MpdStsHit *h = (MpdStsHit*)fItsHits->UncheckedAt(hit->GetIndex());
  1449. hit->SetFlag(-1);// work 15.04
  1450. h->SetFlag(-1); // work 15.04
  1451. }
  1452. }
  1453. }
  1454. //__________________________________________________________________________
  1455. Double_t MpdCellAutomat::Interp(Double_t angt, Int_t choice, Int_t lay)
  1456. {
  1457. // Perform parabolic interpolation for Angl (choice = 0) and dAngt (choice = 1)
  1458. const Int_t np = 7;
  1459. const Double_t xp[np] = {0.1/0.18, 0.1/0.069, 0.1/0.043, 0.1/0.034, 0.1/0.017, 0.1/0.0086, 0.1/0.0044};
  1460. // const Double_t xp[np] = {0.05,0.125,0.2, 0.25, 0.50, 1.00, 2.00}; //3.09
  1461. Double_t dpp05[np] = {0.18, 0.069, 0.043, 0.034, 0.017, 0.0086, 0.0044}; //mean angt
  1462. Double_t err05[np] = {0.0018, 0.00032, 0.00018, 0.00014, 0.000092, 0.00007, 0.00006}; // error mean angt
  1463. const Double_t angs[2][np] = {{0.042, 0.019, 0.017, 0.016, 0.017, 0.019, 0.018}, //sigma angl
  1464. {0.054, 0.013, 0.0069, 0.006, 0.0036, 0.0034, 0.003}}; // sigma dangt Box(250)
  1465. const Double_t dangs[2][np] = {{0.0015, 0.00067, 0.00069, 0.00065, 0.00073, 0.00065, 0.00068}, // error sigma angl
  1466. {0.0027, 0.00045, 0.00027, 0.00018, 0.00014, 0.0001, 0.0001}}; //error sigma dangt
  1467. /*
  1468. const Double_t angls[3][np] = {{0.047, 0.027, 0.025, 0.023, 0.027, 0.026, 0.028},
  1469. {0.04, 0.022, 0.021, 0.02, 0.021, 0.021, 0.021},
  1470. {0.041, 0.014, 0.013, 0.012, 0.012, 0.012, 0.012}}; // sigma deltaz 2,3,4 layers
  1471. const Double_t angts[2][np] = {{0.062, 0.013, 0.0075, 0.0058, 0.0041, 0.0039, 0.0036},
  1472. {0.059, 0.013, 0.0074, 0.006, 0.0035, 0.0029, 0.0026}}; //sigma dangt 3,4 layers
  1473. */
  1474. const Double_t angls[3][np] = {{0.045, 0.026, 0.025, 0.023, 0.027, 0.027, 0.028},
  1475. {0.02, 0.011, 0.01, 0.01, 0.011, 0.011, 0.01},
  1476. {0.016, 0.0045, 0.0037, 0.0036, 0.0033, 0.0033, 0.0031}}; // sigma deltaz 2,3,4 layers
  1477. const Double_t angts[2][np] = {{0.031, 0.0068, 0.0037, 0.0029, 0.0021, 0.0019, 0.0018},
  1478. {0.019, 0.0039, 0.0021, 0.0016, 0.00091, 0.00059, 0.0005}}; //sigma dangt 3,4 layers
  1479. Double_t x0 = 0.1 / angt;
  1480. //Double_t x0 = angt;
  1481. Int_t ok = 0, inds[3] = {np - 3, np - 2, np - 1};
  1482. for (Int_t i = 0; i < np; ++i) {
  1483. if (xp[i] < x0) continue;
  1484. inds[1] = i;
  1485. inds[0] = i - 1;
  1486. inds[2] = i + 1;
  1487. ok = 1;
  1488. break;
  1489. }
  1490. if (ok) {
  1491. if (inds[0] < 0) for (Int_t i = 0; i < 3; ++i) ++inds[i];
  1492. if (inds[2] == np) for (Int_t i = 0; i < 3; ++i) --inds[i];
  1493. }
  1494. Double_t a,b,c,dy01,dy02;
  1495. // Parabolic interpolation
  1496. Double_t dx01 = xp[inds[1]] - xp[inds[0]], dx02 = xp[inds[2]] - xp[inds[0]];
  1497. Double_t x01 = xp[inds[1]] + xp[inds[0]], x02 = xp[inds[2]] + xp[inds[0]];
  1498. // Double_t dy01 = angs[choice][inds[1]] - angs[choice][inds[0]];
  1499. // Double_t dy02 = angs[choice][inds[2]] - angs[choice][inds[0]];
  1500. if (choice == 0) {
  1501. dy01 = angls[lay][inds[1]] - angls[lay][inds[0]];
  1502. dy02 = angls[lay][inds[2]] - angls[lay][inds[0]];
  1503. }
  1504. if (choice == 1) {
  1505. dy01 = angts[lay][inds[1]] - angts[lay][inds[0]];
  1506. dy02 = angts[lay][inds[2]] - angts[lay][inds[0]];
  1507. }
  1508. Double_t slope = dy01 / dx01;
  1509. a = dy02 - slope * dx02;
  1510. a /= (dx02 * x02 - dx02 * x01);
  1511. b = slope - a * x01;
  1512. if (choice == 0) {
  1513. // c = angs[choice][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1514. c = angls[lay][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1515. }
  1516. if (choice == 1) {
  1517. // c = angs[choice][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1518. c = angts[lay][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1519. }
  1520. //for (Int_t i = 0; i < 3; ++i) cout << xp[inds[i]] << " " << angs[choice][inds[i]] << " ";
  1521. //cout << endl << x0 << " " << a * x0 * x0 + b * x0 + c << endl;
  1522. return a * x0 * x0 + b * x0 + c;
  1523. }
  1524. //__________________________________________________________________________
  1525. void MpdCellAutomat::GetShortTracks()
  1526. {
  1527. /// Collect remaining tracks with small number of hits.
  1528. map<TString,Int_t>::iterator it;
  1529. const Int_t nLays = 4;
  1530. Int_t inds[nLays];
  1531. for (it = fCellMap.begin(); it != fCellMap.end(); ++it) {
  1532. //cout << it->first << " " << it->second <<endl;
  1533. if (it->second < 0) continue;
  1534. TString code = it->first, cc;
  1535. Int_t ibeg = 0, leng = -1;
  1536. //cout << code << endl;
  1537. for (Int_t i = 0; i < nLays; ++i) {
  1538. leng = code.Index("-",leng);
  1539. cc = code(ibeg,leng-ibeg);
  1540. //cout << ibeg << " " << leng << endl;
  1541. inds[i] = cc.Atoi();
  1542. ibeg = leng + 1;
  1543. ++leng;
  1544. //cout << i << " " << cc << endl;
  1545. }
  1546. for (Int_t i = 0; i < nLays; ++i) {
  1547. MpdCellTrack *tr = (MpdCellTrack*) f2DHits[i]->UncheckedAt(inds[i]);
  1548. MpdStsHit *hit = (MpdStsHit*) fItsHits->UncheckedAt(tr->GetIndex());
  1549. if (hit->GetFlag() < 0) { fCellMap[code] = -1; break; }
  1550. }
  1551. }
  1552. Int_t left = 0;
  1553. cout << " Leftovers: " << endl;
  1554. for (it = fCellMap.begin(); it != fCellMap.end(); ++it) {
  1555. if (it->second > 0) cout << ++left << " " << it->first << endl;
  1556. }
  1557. }
  1558. ClassImp(MpdCellAutomat);