MpdVectorFinder.cxx 97 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454
  1. // -------------------------------------------------------------------------
  2. // ----- source file -----
  3. // ----- Created 13/06/2014 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdVectorFinder.cxx
  6. *@author D.Zinchenko, A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD Inner Tracking System (ITS) using Vector Finder approach
  9. **/
  10. #include "MpdStsGeoPar.h"
  11. #include "MpdVectorFinder.h"
  12. #include "MpdVector.h"
  13. #include "MpdCodeTimer.h"
  14. #include "MpdConstField.h"
  15. #include "MpdItsHit5spd.h"
  16. #include "MpdItsKalmanTrack.h"
  17. #include "MpdKalmanFilter.h"
  18. #include "MpdKalmanGeoScheme.h"
  19. #include "MpdKalmanHit.h"
  20. #include "MpdKalmanTrack.h"
  21. #include "MpdMCTrack.h"
  22. #include "MpdStsPoint.h"
  23. //#include "MpdTpcKalmanTrack.h"
  24. #include "MpdTpcKalmanFilter.h"
  25. #include "FairGeoNode.h"
  26. #include "FairMCPoint.h"
  27. #include "FairRootManager.h"
  28. #include "FairRun.h"
  29. #include "FairRunAna.h"
  30. #include "FairRuntimeDb.h"
  31. #include "TGeoMatrix.h"
  32. #include "TGeoBBox.h"
  33. #include "TGeoNode.h"
  34. #include "TGeoTube.h"
  35. #include "TGeoManager.h"
  36. #include "TMath.h"
  37. #include "TFile.h"
  38. #include "TVector2.h"
  39. #include "TVector3.h"
  40. #include "TClonesArray.h"
  41. #include <TRandom.h>
  42. #include <TArrayI.h>
  43. #include <iostream>
  44. #include <map>
  45. #include <set>
  46. #include <string>
  47. using std::cout;
  48. using std::endl;
  49. //using std::map;
  50. const Double_t MpdVectorFinder::fgkChi2Cut = 10; //20; //100;
  51. static clock_t tStart = 0;
  52. static clock_t tFinish = 0;
  53. static clock_t tAll = 0;
  54. FILE *lunErr = nullptr; //fopen("error1.dat","w");
  55. //FILE* lunErr= fopen ("file.txt","w");
  56. std::ofstream fout_v("log.txt");
  57. //__________________________________________________________________________
  58. MpdVectorFinder::MpdVectorFinder(const char *name, Bool_t sa, Int_t iVerbose)
  59. :FairTask("MpdVectorFinder", iVerbose),
  60. fNPass(6),
  61. fExact(0), //(0),
  62. fNTracks(0),
  63. fTrackExact(nullptr),
  64. fGeo(0),
  65. fSa(sa), // by default its + tpc tracking is done
  66. fTpcKF(nullptr)
  67. {
  68. fKHits1 = new TClonesArray("MpdKalmanHit", 100);
  69. for (Int_t i = 0; i < 5; ++i) {
  70. fKHits[i] = new TClonesArray("MpdVector", 100);
  71. f2DHits[i] = new TClonesArray("MpdVector", 100);
  72. }
  73. //fKHits[0] = f2DHits[0]; //work
  74. fTracks = new TClonesArray("MpdItsKalmanTrack", 100);
  75. fTrackCand = new TClonesArray("MpdItsKalmanTrack", 100);
  76. if (fExact) fTrackExact = new TClonesArray("MpdItsKalmanTrack", 100);
  77. fHistoDir = 0x0;
  78. fhLays = new TH1F("hLaysITS","ITS layers",10,0,10);
  79. fLayPointers = 0x0;
  80. }
  81. //__________________________________________________________________________
  82. MpdVectorFinder::~MpdVectorFinder()
  83. {
  84. //delete fKHits;
  85. //delete fTracks;
  86. //delete fTrackCand;
  87. delete [] fLayPointers;
  88. delete fhLays;
  89. }
  90. //__________________________________________________________________________
  91. InitStatus MpdVectorFinder::Init()
  92. {
  93. //----- Get pipe radius
  94. fPipeR = ((TGeoTube*)gGeoManager->GetVolume("pipe1")->GetShape())->GetRmax();
  95. cout << "************** Pipe radius: " << fPipeR << " cm " << endl;
  96. //return ReInit();
  97. if (ReInit() != kSUCCESS) return kERROR;
  98. // Read database
  99. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  100. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  101. //cout << geoPar << endl;
  102. TString volName = "sts01 ", path = "";
  103. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  104. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  105. //----- Process modular geometry
  106. fGeo = 1;
  107. Int_t nLay = 5;
  108. //Double_t safety = 0.005;
  109. for (Int_t i = 0; i < nLay; ++i) {
  110. fNladders[i] = fNsectors[i] = 0;
  111. volName = "sts01ladder";
  112. volName += (i+1);
  113. path = "/cave_1/sts01_0/" + volName;
  114. path += "_";
  115. TString path0 = path;
  116. //----- Loop over all ladders to find the one with the smallest radius
  117. //fRad[i] = 999.9;
  118. fRad[i] = 0.;
  119. Int_t nladd = 0;
  120. for (Int_t jlad = 1; jlad < 99; ++jlad) {
  121. path = path0;
  122. path += jlad;
  123. if (!gGeoManager->CheckPath(path)) break;
  124. gGeoManager->cd(path);
  125. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  126. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  127. Double_t xyzL[3] = {0}, xyzM[3];
  128. gGeoManager->LocalToMaster(xyzL,xyzM);
  129. Double_t rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  130. fRad[i] += rad;
  131. ++nladd;
  132. /*AZ - 3.12.2017
  133. fRad[i] = TMath::Min (fRad[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[i] = TMath::Min (fRad[i],rad);
  138. xyzL[0] = -box->GetDX();
  139. gGeoManager->LocalToMaster(xyzL,xyzM);
  140. rad = TMath::Sqrt (xyzM[0] * xyzM[0] + xyzM[1] * xyzM[1]);
  141. fRad[i] = TMath::Min (fRad[i],rad);
  142. */
  143. }
  144. cout << " *** Layer # " << i << " min_R= " << fRad[i] << endl;
  145. TGeoVolume *ladd = gGeoManager->GetVolume(volName);
  146. if (!ladd) { nLay = i; break; }
  147. fRad[i] /= nladd; // mean radius
  148. /*
  149. TGeoBBox *box = (TGeoBBox*) ladd->GetShape();
  150. //safety = -box->GetDY();
  151. //safety = box->GetDY();
  152. safety = 2 * box->GetDY(); // new
  153. fRad[i] -= safety;
  154. */
  155. cout << " +++ Layer # " << i << " safety min_R= " << fRad[i] << endl;
  156. }
  157. FillGeoScheme();
  158. //----- Get cables
  159. TObjArray *vols = gGeoManager->GetListOfVolumes();
  160. Int_t nvols = vols->GetEntriesFast(), ncables = 0;
  161. for (Int_t i = 0; i < nvols; ++i) {
  162. TGeoVolume *vol = (TGeoVolume*) vols->UncheckedAt(i);
  163. TString cable = TString(vol->GetName());
  164. if (!(cable.Contains("sts") && cable.Contains("cable"))) continue;
  165. // cout << cable << endl;
  166. ++ncables;
  167. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  168. TString lad = cable;
  169. Int_t ip = lad.Index("cable");
  170. lad.Replace(ip,lad.Length()-ip+1,"");
  171. lad.Replace(lad.Length()-2,1,"");
  172. lad += "_1/";
  173. path = "/cave_1/sts01_0/" + lad + cable;
  174. path += "_1";
  175. gGeoManager->cd(path);
  176. Double_t xyzL[3] = {0}, xyzM[3];
  177. gGeoManager->LocalToMaster(xyzL,xyzM);
  178. // cout << xyzM[2] - box->GetDZ() << " " << xyzM[2] + box->GetDZ() << " " << box->GetDZ() << endl;
  179. if (xyzM[2] - box->GetDZ() > 0) {
  180. Int_t lay = TString(lad(lad.Length()-4,1)).Atoi();
  181. fCables[lay-1].insert(pair<Double_t,Double_t>(xyzM[2] - box->GetDZ(),xyzM[2] + box->GetDZ()));
  182. }
  183. }
  184. cout << "-I- MpdTrackFinderIts5spd: Intialization finished successfully" << endl;
  185. fout_v << "-Init - MpdTrackFinderIts5spd: Intialization finished successfully" << endl;
  186. return kSUCCESS;
  187. }
  188. //__________________________________________________________________________
  189. void MpdVectorFinder::FillGeoScheme()
  190. {
  191. /// Fill Kalman filter geometry manager info
  192. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  193. TGeoVolume *vSts = gGeoManager->GetVolume("sts01");
  194. for (Int_t layer = 1; layer < 999; ++layer) {
  195. // Loop over layers
  196. TString sladder = "sts01ladder";
  197. sladder += layer;
  198. TGeoVolume *vLay = gGeoManager->GetVolume(sladder);
  199. if (vLay == 0x0 && layer == 1) continue; // no first layer
  200. if (vLay == 0x0) break;
  201. sladder += "_";
  202. //----- Loop over ladders
  203. for (Int_t ladder = 1; ladder < 999; ++ladder) {
  204. TString sladder1 = sladder;
  205. sladder1 += ladder;
  206. TGeoNode *node = vSts->FindNode(sladder1);
  207. if (node == 0x0) break;
  208. ++fNladders[layer - 1];
  209. TGeoVolume *vLad = node->GetVolume();
  210. //cout << vLad->GetNodes()->GetEntriesFast() << " " << vLad->GetNdaughters() << endl;
  211. Int_t nDaught = vLad->GetNdaughters(), detID = -1, detIDsts = -1;
  212. TObjArray *daught = vLad->GetNodes();
  213. //----- Loop over ladder daughters
  214. Int_t iZ = 0;
  215. for (Int_t det = 0; det < nDaught; ++det) {
  216. TString sdet1 = ((TGeoNode*)(daught->UncheckedAt(det)))->GetName();
  217. ///if (!sdet1.Contains("sensor") && !sdet1.Contains("sector")) continue;
  218. if (ladder == 1) { ++fNsectors[layer-1];}
  219. Int_t det1 = TString(sdet1(sdet1.Index("_")+1,2)).Atoi();
  220. Int_t secType = -1;
  221. if (sdet1.Contains("sector")) secType = TString(sdet1(sdet1.Index("_")-2,1)).Atoi();
  222. ++iZ;
  223. Int_t side = 0.;
  224. detIDsts = fHitSts.SetDetId(layer, ladder, det1);
  225. detID = fHitSts.SetDetId(layer, ladder, iZ);
  226. fId2Id[layer-1].insert(pair<Int_t,Int_t>(detIDsts,detID)); // detIDsts == detID for pixel detectors (AZ)
  227. ///cout << "Detector ids " << detIDsts << " " << detID << endl;
  228. detID += 1000000 * (layer-1);
  229. TString detName = sladder1 + "/" + sdet1 + "#";
  230. detName += side;
  231. geo->SetDetId(detName, detID);
  232. TString path = "/cave_1/sts01_0/" + detName(0,detName.Length()-2);
  233. //cout << detName << " " << path << endl;
  234. gGeoManager->cd(path);
  235. geo->SetPath(detID, path);
  236. node = gGeoManager->GetCurrentNode();
  237. //cout << node->GetName() << " " << detID << endl;
  238. TGeoVolume *vol = node->GetVolume();
  239. TGeoBBox *box = (TGeoBBox*) vol->GetShape();
  240. TVector2 size(box->GetDX(), box->GetDZ());
  241. geo->SetSize(detID, size);
  242. Double_t xyzL[3] = {0}, xyzM[3], vecM[3];
  243. xyzL[1] = 1;
  244. gGeoManager->LocalToMasterVect(xyzL,vecM);
  245. xyzL[1] = 0;
  246. gGeoManager->LocalToMaster(xyzL,xyzM);
  247. Double_t sign = TMath::Sign (1.,xyzM[0]*vecM[0]+xyzM[1]*vecM[1]);
  248. TVector3 pos(xyzM[0], xyzM[1], xyzM[2]);
  249. geo->SetGlobalPos(detID, pos);
  250. TVector3 norm(sign*vecM[0], sign*vecM[1], sign*vecM[2]);
  251. geo->SetNormal(detID, norm);
  252. }
  253. }
  254. }
  255. fout_v << "- FillGeoScheme - done" << endl;
  256. }
  257. //__________________________________________________________________________
  258. InitStatus MpdVectorFinder::ReInit()
  259. {
  260. FairRootManager *manager = FairRootManager::Instance();
  261. fItsPoints = (TClonesArray *) manager->GetObject("StsPoint");
  262. fItsHits =(TClonesArray *) manager->GetObject("StsHit");
  263. if (fItsPoints == 0x0 || fItsHits == 0x0) return kERROR;
  264. fTpcTracks =(TClonesArray *) manager->GetObject("TpcKalmanTrack");
  265. fMCTracks =(TClonesArray *) manager->GetObject("MCTrack");
  266. manager->Register("ItsTrack", "Its", fTracks, kTRUE);
  267. manager->Register("ItsTrackCand", "ItsCand", fTrackCand, kTRUE);
  268. if (fTrackExact) manager->Register("ItsTrackExact", "ItsCand", fTrackExact, kTRUE);
  269. manager->Register("ItsKHits", "ItsKalHits", fKHits1, kFALSE);
  270. manager->Register("CellTrack0", "Its0", fKHits[0], kFALSE);
  271. manager->Register("CellTrack1", "Its1", fKHits[1], kFALSE);
  272. manager->Register("CellTrack2", "Its2", fKHits[2], kFALSE);
  273. manager->Register("CellTrack3", "Its3", fKHits[3], kFALSE);
  274. manager->Register("CellTrack4", "Its4", fKHits[4], kFALSE);
  275. manager->Register("CellTrack_0", "Its0_2D", f2DHits[0], kFALSE);
  276. manager->Register("CellTrack_1", "Its1_2D", f2DHits[1], kFALSE);
  277. manager->Register("CellTrack_2", "Its2_2D", f2DHits[2], kFALSE);
  278. manager->Register("CellTrack_3", "Its3_2D", f2DHits[3], kFALSE);
  279. manager->Register("CellTrack_4", "Its4_2D", f2DHits[4], kFALSE);
  280. //fNPass = 2; // 2 prochoda
  281. //fNPass = 3;
  282. fNPass = 6;
  283. fout_v << "- ReInit - done" << endl;
  284. return kSUCCESS;
  285. }
  286. //__________________________________________________________________________
  287. void MpdVectorFinder::Reset()
  288. {
  289. cout << " MpdVectorFinder::Reset " << endl;
  290. //fKHits->Clear();
  291. /// changed list of function calls to a loop
  292. for (Int_t i = 0; i < 5; ++i) { /// was 4
  293. fKHits[i]->Delete();
  294. f2DHits[i]->Delete();
  295. }
  296. fKHits1->Delete();
  297. fTracks->Delete();
  298. fTrackCand->Delete();
  299. if (fTrackExact) fTrackExact->Delete();
  300. delete [] fLayPointers;
  301. fLayPointers = nullptr;
  302. fCellMap.clear();
  303. fout_v << "- Reset - done" << endl;
  304. }
  305. //__________________________________________________________________________
  306. void MpdVectorFinder::SetParContainers()
  307. {
  308. // Get run and runtime database
  309. FairRunAna* run = FairRunAna::Instance();
  310. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  311. FairRuntimeDb* db = run->GetRuntimeDb();
  312. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  313. // Get STS geometry parameter container
  314. db->getContainer("MpdStsGeoPar");
  315. }
  316. //__________________________________________________________________________
  317. void MpdVectorFinder::Finish()
  318. {
  319. //Write();
  320. if (lunErr) fclose(lunErr);
  321. }
  322. //__________________________________________________________________________
  323. void MpdVectorFinder::Exec(Option_t * option)
  324. {
  325. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  326. static int eventCounter = 0;
  327. cout << " - - - - \n ItsRec event " << ++eventCounter << endl;
  328. Reset();
  329. tStart = clock();
  330. // Create Kalman hits
  331. if (fItsHits->GetEntriesFast() == 0) return;
  332. Build2DHits(); ///MakeKalmanHits stuff was added to Build2DHits
  333. /// if MpdKalmanHit->GetFlag == -1 then hit was used earlier (see ExcludeHits())
  334. /// fNPass is set in Reinit function
  335. /// i is for pass number:
  336. /// 0 - primary tracks, 5-layer - tight cuts
  337. // 1 - primary tracks, 5-layer - loose cuts
  338. /// 2 - primary tracks, 4-layer
  339. /// 3 - secondary tracks, 5-layer - tight cuts
  340. /// 4 - secondary tracks, 5-layer - loose cuts
  341. /// 5 - secondary tracks, 4-layer
  342. if (fTrackExact) fTrackExact->Delete();
  343. for (Int_t i = 0; i < fNPass; ++i) {
  344. /// resetting after each pass updated 19.2.2020
  345. for (Int_t j = 0; j < 5; ++j) {
  346. fKHits[j]->Delete();
  347. }
  348. fTrackCand->Delete(); //AZ-29.02.20
  349. fCellMap.clear();
  350. fVectorCands.clear();
  351. ///
  352. MakeTrackCandidates(i);// 21.02
  353. if (i < 3)
  354. ExtendCellTracks(i);// 21.02
  355. else
  356. ExtendSecondaryTracks(i);
  357. GetTrackSeeds(i);//attention!!
  358. DoTracking(i); //attention!!!
  359. //cout << " AZ before: " << fTrackCand->GetEntriesFast() << endl;
  360. RemoveDoubles();
  361. //cout << " AZ after: " << fTrackCand->GetEntriesFast() << endl;
  362. StoreTracks();
  363. fout_v << " Total number of found tracks: " << fTrackCand->GetEntriesFast() << endl;
  364. ExcludeHits(); // exclude used hits test
  365. Int_t usedHits = 0, unusedHits = 0;
  366. for (int j = 0; j < fKHits1->GetEntriesFast(); j++) {
  367. if (((MpdKalmanHit*)fKHits1->UncheckedAt(j))->GetFlag() < 0) {
  368. usedHits++;
  369. } else {
  370. unusedHits++;
  371. }
  372. }
  373. cout << "\n used/unused hits before making track candidates: " << usedHits << " " << unusedHits << endl;
  374. fout_v << "\n used/unused hits before making track candidates: " << usedHits << " " << unusedHits << endl;
  375. }
  376. GetShortTracks();
  377. AddHits(); // add hit objects to tracks
  378. tFinish = clock();
  379. tAll = tFinish - tStart;
  380. cout << " Total number of found tracks: " << fTrackCand->GetEntriesFast() << endl;
  381. cout << " Exec time: " << ((Float_t)tAll) / CLOCKS_PER_SEC << endl;
  382. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  383. fout_v.close();
  384. }
  385. //_________________________________________________________________________
  386. void MpdVectorFinder::Build2DHits()
  387. {
  388. /// Build track candidates from ITS 2-D hits.
  389. fItsHits->Sort();
  390. Int_t nHits = fItsHits->GetEntriesFast();//, layMax = 0, nKH = 0;
  391. Double_t errZ = 0.001, errX = 0.001; ///10um for Z and X// 250um in Z, 23um in R-Phi (local X)
  392. Double_t xloc;//, z;
  393. Double_t r;
  394. //Double_t dX = 0, dZ = 0;
  395. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  396. Int_t nKalmanHit = 0, nLayerHit = 0;
  397. for (Int_t ih = 0; ih < nHits; ++ih) {
  398. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(ih);
  399. if (h->GetFlag() < 0) continue; /// TODO check if necessary. Build2DHits is done only once
  400. Int_t lay = h->Layer() - 1;
  401. //AZ - for debug
  402. //Int_t id = ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  403. //if (id != 2346 && id != 2305 && id != 2876) continue;
  404. //AZ Get local Z
  405. //TString path = geo->Path(h->GetDetectorID()+1000000*lay);
  406. TString path = geo->Path(fId2Id[lay][h->GetDetectorID()]+1000000*lay);
  407. gGeoManager->cd(path);
  408. Double_t posLoc[3] = {0}, posMas[3] = {h->GetX(), h->GetY(), h->GetZ()};
  409. /// transform global coords to local
  410. gGeoManager->MasterToLocal(posMas,posLoc);
  411. //AZ
  412. r = TMath::Sqrt (h->GetX() * h->GetX() + h->GetY() * h->GetY());
  413. xloc = h->GetLocalX();
  414. //z = h->GetZ();
  415. //AZ Double_t meas[2] = {xloc + h->GetDx() * errX, z + dZ * errZ};
  416. Double_t meas[2] = {xloc + h->GetDx() * errX, posLoc[2] + h->GetDz() * errZ};
  417. //AZ Double_t err[2] = {errX, errZ};
  418. //AZ-161020 Double_t err[2] = {errX*1.2, errZ*1.2}; // effective measure
  419. Double_t err[2] = {errX*1.4, errZ*1.4}; // effective measure
  420. if (lay > 2) { err[0] *= 2; err[1] *= 2; } //AZ-161020 - thick layers
  421. Double_t cossin[2] = {TMath::Cos(fStereoA[0]), TMath::Sin(fStereoA[0])}; /// wtf is this
  422. MpdVector *hit = 0x0;//ms
  423. MpdItsHit5spd *hsts = h;
  424. Double_t coord[3] = {h->GetX(), h->GetY(), h->GetZ()};
  425. TVector3 point(coord);
  426. ///fout_v << "ih = " << ih << " layer = " << h->Layer()-1 << " detID = " << hsts->GetDetectorID() << " " << fId2Id[lay][hsts->GetDetectorID()] << endl;
  427. MpdKalmanHit *khit = 0x0;
  428. /// fId2Id[lay][hsts->GetDetectorID()] - encoded hit identifier (# of sensor, layer, ladder)
  429. khit = new ((*fKHits1)[nKalmanHit++]) MpdKalmanHit(lay * 1000000 + fId2Id[lay][hsts->GetDetectorID()], //detID
  430. 2, //nDim
  431. MpdKalmanHit::kFixedP, //HitType // different modules (planes) in local coordinates
  432. meas, // 2D measurement for kalman filter (x and z)
  433. err, // err
  434. cossin, //
  435. 0., // signal
  436. r, // dist
  437. ih); //index
  438. //AZ khit->SetUniqueID(0);
  439. khit->SetUniqueID(nKalmanHit); //AZ
  440. nLayerHit = f2DHits[lay]->GetEntriesFast();
  441. hit = new ((*f2DHits[lay])[nLayerHit]) MpdVector(lay * 1000000 + fId2Id[lay][hsts->GetDetectorID()], //detID
  442. 2, //nDim
  443. MpdVector::kFixedP, //HitType // different modules (planes) in local coordinates
  444. point, // TVector3 measurement
  445. err, // err
  446. cossin, //
  447. 0., // signal
  448. r, // dist
  449. ih, // index
  450. 0, // index1
  451. -1, // trackNo - deprecated?
  452. nullptr, // TrackPointer
  453. khit); // KalmanHit pointer
  454. hit->SetUniqueID(nLayerHit);
  455. ///cout << hit->GetKalmanHit() << " ";
  456. } // for (Int_t ih = 0; ih < nHits; ++ih)
  457. fout_v << "khits " << nKalmanHit << endl;
  458. cout << "- Build2DHits - done" << endl;
  459. }
  460. //_________________________________________________________________________
  461. void MpdVectorFinder::MakeTrackCandidates(Int_t iPass)
  462. {
  463. // Read 2D Hits and create trackCandidates from 2D Hits at first(last) layer
  464. //Int_t n = 5 - iPass; // last layer
  465. //n = TMath::Min (n, 4);
  466. const Int_t layNo[9] = {4, 4, 3, 4, 4, 3};
  467. Int_t n = layNo[iPass];
  468. Int_t nKH = 0, nHits = f2DHits[n]->GetEntriesFast();
  469. if (iPass < 3) {
  470. fout_v << "layer " << n << " hits " << nHits << endl;
  471. for (Int_t i = 0; i < nHits; ++i) {
  472. MpdVector *cellTr = (MpdVector*) f2DHits[n]->UncheckedAt(i);
  473. ///MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits-> UncheckedAt(cellTr->GetIndex(0));
  474. ///fout_v << h->GetTrackID() << " " << cellTr->GetMeas()[2] << endl;
  475. cellTr = new ((*fKHits[n])[nKH++]) MpdVector(*cellTr);
  476. ///fout_v << cellTr->GetKalmanHit() << " ";
  477. cellTr->SetCode(i);
  478. cellTr->SetCosSin(0, TMath::ATan2(cellTr->GetMeas()[1], cellTr->GetMeas()[0])); /// transverse angle for first layer hits
  479. Float_t r = TMath::Sqrt(cellTr->GetMeas()[0] * cellTr->GetMeas()[0] + cellTr->GetMeas()[1] * cellTr->GetMeas()[1]); /// distance from interaction point
  480. cellTr->SetCosSin(1, TMath::ATan2(r, cellTr->GetMeas()[2])); /// longitudinal angle
  481. /// added 17.1.2019
  482. cellTr->SetUniqueID(nKH - 1);
  483. ///fout_v << nKH - 1 << endl;
  484. ///fout_v << cellTr->GetKalmanHit() << " ";
  485. //AZ
  486. //MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(cellTr->GetIndex(0));
  487. //Int_t trID = ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  488. //cout << " AZ-TrackID: " << trID << endl;
  489. } // for (Int_t i = 0; ih < nHits; ++i)
  490. } else {/// added 29.1.2020
  491. // 4th and 5th pass are done for secondary particles
  492. //n = 5 - iPass + 2;
  493. nHits = f2DHits[n]->GetEntriesFast();
  494. /// First hit of track candidate is a hit from 1st layer
  495. for (Int_t i = 0; i < nHits; ++i) {
  496. MpdVector *cellTr = (MpdVector*) f2DHits[n]->UncheckedAt(i);
  497. if ((cellTr->GetKalmanHit())->GetFlag() < 0) {
  498. continue;
  499. }
  500. cellTr = new ((*fKHits[n])[nKH++]) MpdVector(*cellTr);
  501. cellTr->SetCode(i);
  502. cellTr->SetCosSin(0, TMath::ATan2(cellTr->GetMeas()[1], cellTr->GetMeas()[0])); /// transverse angle for first layer hits
  503. Float_t r = TMath::Sqrt(cellTr->GetMeas()[0] * cellTr->GetMeas()[0] + cellTr->GetMeas()[1] * cellTr->GetMeas()[1]); /// distance from interaction point
  504. cellTr->SetCosSin(1, TMath::ATan2(r, cellTr->GetMeas()[2])); //cellTr->GetMeas()[2] /// z coordinate instead of longitudinal angle? /// TODO 21.2.20
  505. cellTr->SetUniqueID(nKH - 1);
  506. //AZ
  507. cellTr->SetDeltaZ(-1/TMath::Tan(cellTr->GetCosSin(1)));
  508. //MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(cellTr->GetIndex(0));
  509. //Int_t trID = ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  510. //cout << " AZ1-TrackID: " << trID << " " << cellTr->GetMeas()[2] << endl;
  511. }
  512. cout << "Track Candidates on pass " << iPass + 1 << ": " << nKH << endl;
  513. }
  514. fout_v << "- MakeTrackCandidates - done" << endl;
  515. }
  516. //_________________________________________________________________________
  517. /// added 30.1.2020
  518. void MpdVectorFinder::ExtendSecondaryTracks(Int_t iPass)
  519. {
  520. //const Int_t nSigm = 3;
  521. int count = 0;
  522. //Double_t errZ = 0.001, errX = 0.001;/// 10um for Z and X // 250um in Z, 23um in R-Phi (local X)
  523. multimap <Float_t, MpdVector*> zMultimapHits[5], tMultimapHits[5];
  524. fout_v << "(transverse angle, longitudinal angle, z | MCtrackID)" << endl;
  525. //Int_t n = 4 - iPass + 2;
  526. const Int_t layNo[9] = {0, 0, 0, 3, 3, 2};
  527. //for (Int_t i234 = n; i234 >= 0; i234--) { /// loop over layers 3-0
  528. for (Int_t i234 = layNo[iPass]; i234 >= 0; i234--) { /// loop over layers 3-0
  529. Int_t nKH = 0;//, iprint = 0;
  530. cout << " AZsec: " << i234 << " " << f2DHits[i234]->GetEntriesFast() << endl; //AZ
  531. Double_t rmean = 0; //AZ
  532. for (Int_t i = 0; i < f2DHits[i234]->GetEntriesFast(); i++) {
  533. MpdVector *hit = (MpdVector*) f2DHits[i234]->UncheckedAt(i);
  534. if (hit->GetKalmanHit()->GetFlag() < 0) continue; //AZ
  535. TVector3 meas = hit->GetMeas();
  536. Double_t phi = TMath::ATan2(meas[1], meas[0]);
  537. zMultimapHits[i234].insert(pair<Float_t, MpdVector*>(meas[2], hit)); /// z coordinate
  538. tMultimapHits[i234].insert(pair<Float_t, MpdVector*>(phi, hit)); /// transverse, phi
  539. /// if transverse angle is close to pi, add it to other side of map with angle - 2pi
  540. if (TMath::Abs(TMath::Pi() - TMath::Abs(phi)) < TMath::Pi() / 9.0) {
  541. tMultimapHits[i234].insert(pair<Float_t, MpdVector*>(phi - TMath::Sign(1, phi) * 2 * TMath::Pi(), hit)); /// transverse, phi, duplicate hit
  542. }
  543. rmean += hit->GetDist(); //AZ
  544. }
  545. rmean /= zMultimapHits[i234].size(); //AZ
  546. int nTracks = fKHits[i234 + 1]->GetEntriesFast(); /// was -1
  547. cout << "secondary nTracks for layer " << i234+1 << ": " << nTracks << " " << rmean << endl;
  548. for (Int_t i = 0; i < nTracks; i++) {
  549. MpdVector *cand = (MpdVector*) fKHits[i234 + 1]->UncheckedAt(i); /// was -1
  550. TVector3 meas = cand->GetMeas();
  551. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(cand->GetIndex(0));
  552. Int_t trID0 = ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  553. Double_t pt = 0, epst_lower, epst_upper;
  554. //AZ Double_t epsz_lower = 5e-2, epsz_upper = 5e-2;
  555. //Double_t epsz_lower = 3, epsz_upper = 3;
  556. Double_t epsz_lower = 5, epsz_upper = 5, dphi = 0;
  557. MpdVector *prev = cand->GetPrevTrackPointer(), *preprev = nullptr;
  558. vector<pair<Double_t,Double_t> > epst(2,pair<Double_t,Double_t>(999,999)); // for fine structure (2 branches)
  559. if (prev != nullptr) {
  560. preprev = prev->GetPrevTrackPointer();
  561. Double_t circle[3] = {0};
  562. if (preprev == nullptr) {
  563. //if (1) {
  564. // Use primary vertex
  565. ///fout_v << cand->GetPrevTrackPointer() << endl;
  566. MpdVector zeroV;
  567. pt = EvalPt(prev, cand, &zeroV, circle);
  568. if (TMath::Abs(pt) < 0.02) continue;
  569. if (!PropagateHelix(rmean, meas, circle, dphi)) continue;
  570. Double_t ptabs = TMath::Abs(pt);
  571. if (i234 == 2) {
  572. // First region
  573. Double_t ptmin1 = TMath::Max (ptabs,0.1);
  574. ptmin1 = TMath::Min (ptmin1, 0.2);
  575. Double_t eps1 = 0.1 / TMath::Power (ptmin1-0.01, 0.3);
  576. Double_t ptmin2 = TMath::Max (ptabs,0.1);
  577. ptmin2 = TMath::Min (ptmin2, 0.3);
  578. Double_t eps2 = 0.04 / TMath::Power (ptmin2-0.01, 0.8);
  579. epst_lower = (pt > 0) ? eps1 : eps2;
  580. epst_upper = (pt > 0) ? eps2 : eps1;
  581. epst[0].first = -epst_lower;
  582. epst[0].second = epst_upper;
  583. //Second region
  584. ptmin1 = TMath::Max (ptabs, 0.02);
  585. Double_t phicor = 0.0569913 / TMath::Power (ptmin1-0.01+0.0627805, 0.894263);
  586. if (pt < 0) phicor = -phicor;
  587. eps1 = 0.01;
  588. eps2 = 0.06;
  589. epst_lower = (pt > 0) ? eps2 : eps1;
  590. epst_upper = (pt > 0) ? eps1 : eps2;
  591. epst[1].first = -epst_lower + phicor;
  592. epst[1].second = epst_upper + phicor;
  593. } else if (i234 == 1) {
  594. // First region
  595. Double_t eps1 = 0.004 / TMath::Power (ptabs-0.01, 1.2);
  596. eps1 = TMath::Max (eps1, 0.03);
  597. Double_t eps2 = 0.1;
  598. epst_lower = (pt > 0) ? eps2 : eps1;
  599. epst_upper = (pt > 0) ? eps1 : eps2;
  600. epst[0].first = -epst_lower;
  601. epst[0].second = epst_upper;
  602. // Second region
  603. Double_t ptmin1 = TMath::Max (ptabs, 0.02);
  604. Double_t phicor = 0.0276577 / TMath::Power (ptmin1-0.01, 0.827315);
  605. if (pt < 0) phicor = -phicor;
  606. Double_t ptmax1 = TMath::Min (ptmin1, 0.3);
  607. eps1 = 0.04 / TMath::Power (ptmax1-0.01, 0.5);
  608. eps1 = TMath::Min (eps1, 0.1);
  609. eps2 = 0.017 / TMath::Power (ptmin1-0.01, 0.4);
  610. eps2 = TMath::Min (eps2, 0.08);
  611. epst_lower = (pt > 0) ? eps2 : eps1;
  612. epst_upper = (pt > 0) ? eps1 : eps2;
  613. epst[1].first = -epst_lower + phicor;
  614. epst[1].second = epst_upper + phicor;
  615. }
  616. //epst_lower = epst_upper = 5 * epst_upper;
  617. //if (i234 == 1) epst_lower = epst_upper = 0.02 / TMath::Power (ptt-0.01, 2./3);
  618. //epsz_lower = epsz_upper = 1.;
  619. Double_t ptt1 = TMath::Min (0.1, ptabs);
  620. epsz_lower = epsz_upper = 0.7 / TMath::Power (ptt1-0.015, 0.5);
  621. if (i234 == 1) epsz_lower = epsz_upper = 1.5;
  622. //epsz_lower = epsz_upper = 5 * epsz_upper;
  623. } else {
  624. // Do not use primary vertex
  625. //Double_t circle[3] = {0};
  626. pt = EvalPt(preprev, prev, cand, circle);
  627. if (TMath::Abs(pt) < 0.02) continue;
  628. if (!PropagateHelix(rmean, meas, circle, dphi)) continue;
  629. Double_t ptt = TMath::Abs(pt);
  630. ptt = TMath::Min (0.1, ptt);
  631. //epst_lower = epst_upper = 0.1;
  632. epst_lower = epst_upper = 0.04 / TMath::Power (ptt-0.01, 0.2);
  633. //if (i234 == 0) epst_lower = epst_upper = 0.02 / TMath::Power (ptt-0.01, 2./3);
  634. epsz_lower = epsz_upper = 1.;
  635. if (iPass == 3) {
  636. epst_lower /= 2;
  637. epst_upper /= 2;
  638. epsz_lower /= 2;
  639. epsz_upper /= 2;
  640. }
  641. if (i234 == 0) {
  642. /*
  643. //epst_lower = epst_upper = 0.2;
  644. Double_t ptt1 = TMath::Abs(pt);
  645. ptt1 = TMath::Min (0.15, ptt1);
  646. ptt1 = TMath::Max (0.05, ptt1);
  647. //epst_lower = 0.025 / TMath::Power (ptt-0.01, 2./3);
  648. //epst_upper = 0.02 / TMath::Power (ptt-0.01, 2./3);
  649. epst_lower = 0.02 / TMath::Power (ptt-0.01, 2./3);
  650. epst_upper = 0.025 / TMath::Power (ptt-0.01, 2./3);
  651. */
  652. epst_lower = 0.1;
  653. epst_upper = 0.15;
  654. epsz_lower = epsz_upper = 1.5;
  655. }
  656. epst[0].first = -epst_lower;
  657. epst[0].second = epst_upper;
  658. }
  659. } else {
  660. // First window
  661. epst_lower = epst_upper = 0.3; //5e-1;
  662. if (i234 == 2) epst_lower = epst_upper = 0.4;
  663. //epst_lower = epst_upper = epst_upper*5;
  664. //epsz_lower = epsz_upper = 4.0; //1.;
  665. //if (i234 == 2) epsz_lower = epsz_upper = 6.0;
  666. epsz_lower = epsz_upper = 10.0; //1.;
  667. if (i234 == 2) epsz_lower = epsz_upper = 20.0;
  668. if (iPass == 3) {
  669. epst_lower /= 2;
  670. epst_upper /= 2;
  671. epsz_lower /= 2;
  672. epsz_upper /= 2;
  673. }
  674. //epsz_lower = epsz_upper = epsz_upper*5;
  675. epst[0].first = -epst_lower;
  676. epst[0].second = epst_upper;
  677. }
  678. Double_t lastTrackZ = meas[2];
  679. Double_t drad = TMath::Abs(cand->GetDist()-rmean); //AZ
  680. Double_t deltaZ = cand->GetDeltaZ() * drad; //AZ cand->GetDeltaZ() is tan(theta)
  681. // multimap<Float_t, MpdVector*>::iterator itlow = zMultimapHits[i234].lower_bound(lastTrackZ - epsz);
  682. // multimap<Float_t, MpdVector*>::iterator ittop = zMultimapHits[i234].upper_bound(lastTrackZ + epsz);
  683. //AZ multimap<Float_t, MpdVector*>::iterator itlow = zMultimapHits[i234].lower_bound(lastTrackZ - epsz_lower);
  684. //AZ multimap<Float_t, MpdVector*>::iterator ittop = zMultimapHits[i234].upper_bound(lastTrackZ + epsz_upper);
  685. multimap<Float_t, MpdVector*>::iterator itlow = zMultimapHits[i234].lower_bound(lastTrackZ - epsz_lower + deltaZ);
  686. multimap<Float_t, MpdVector*>::iterator ittop = zMultimapHits[i234].upper_bound(lastTrackZ + epsz_upper + deltaZ);
  687. Double_t lastTrackT = 0;
  688. //AZ if (preprev) lastTrackT = dphi; // circle from 3 hits
  689. if (prev) lastTrackT = dphi; //AZ-141020 - circle from 2 hits and (0,0) or from 3 hits
  690. else {
  691. lastTrackT = TMath::ATan2(meas[1], meas[0]);
  692. lastTrackT -= dphi;
  693. }
  694. //multimap<Float_t, MpdVector*>::iterator ittlow = tMultimapHits[i234].lower_bound(lastTrackT - epst_lower);
  695. //multimap<Float_t, MpdVector*>::iterator itttop = tMultimapHits[i234].upper_bound(lastTrackT + epst_upper);
  696. multimap<Float_t, MpdVector*>::iterator ittlow = tMultimapHits[i234].lower_bound(lastTrackT + epst[0].first);
  697. multimap<Float_t, MpdVector*>::iterator itttop = tMultimapHits[i234].upper_bound(lastTrackT + epst[0].second);
  698. set<MpdVector*> z, t, intersect;
  699. for (multimap<Float_t, MpdVector*>::iterator itr = itlow; itr != ittop; ++itr) {
  700. z.insert((*itr).second);
  701. }
  702. for (multimap<Float_t, MpdVector*>::iterator itr = ittlow; itr != itttop; ++itr) {
  703. t.insert((*itr).second);
  704. }
  705. if (epst[1].first < 900) {
  706. // Fine structure - second region
  707. //lastTrackT = epst[1][2];
  708. ittlow = tMultimapHits[i234].lower_bound(lastTrackT + epst[1].first);
  709. itttop = tMultimapHits[i234].upper_bound(lastTrackT + epst[1].second);
  710. for (multimap<Float_t, MpdVector*>::iterator itr = ittlow; itr != itttop; ++itr)
  711. t.insert((*itr).second);
  712. }
  713. set_intersection(z.begin(), z.end(), t.begin(), t.end(), std::inserter(intersect, intersect.begin()));
  714. for (set<MpdVector*>::iterator itr = intersect.begin(); itr != intersect.end(); ++itr) {
  715. MpdVector *track = (*itr);
  716. if (track->GetKalmanHit()->GetFlag() < 0) continue; //AZ
  717. MpdItsHit5spd *hh = (MpdItsHit5spd*) fItsHits->UncheckedAt(track->GetIndex(0));
  718. Int_t trID = ((MpdStsPoint*) fItsPoints->UncheckedAt(hh->GetTrackID()))->GetTrackID();
  719. if (fExact && trID != trID0) continue; // exact ID match
  720. MpdVector *trNew = nullptr;
  721. trNew = new ((*fKHits[i234])[nKH++]) MpdVector(*track); /// was h->Layer()
  722. trNew->SetUniqueID(nKH - 1);
  723. trNew->SetNofDim(2);
  724. trNew->SetCosSin(0,0);
  725. trNew->SetCosSin(1,0);
  726. trNew->SetPrevTrackPointer(cand);
  727. trNew->SetCode(cand->GetCode());
  728. trNew->SetCode(i);
  729. /// 21.2.20
  730. trNew->SetDeltaZ(trNew->GetMeas()[2] - trNew->GetPrevTrackPointer()->GetMeas()[2]);
  731. ///cout << "trNew " << trNew->GetMeas()[2] << " " << trNew->GetPrevTrackPointer()->GetMeas()[2] << endl;
  732. /// 25.2.20
  733. //AZ
  734. //trNew->SetDeltaZ(trNew->GetDeltaZ()/TMath::Abs(trNew->GetDist()-trNew->GetPrevTrackPointer()->GetDist()));
  735. trNew->SetDeltaZ(trNew->GetDeltaZ()/(trNew->GetMeas()-trNew->GetPrevTrackPointer()->GetMeas()).Pt());
  736. //cout << " AZ1-TrackID: " << trID << " " << i234 << endl;
  737. ///
  738. Int_t trackID = 0;
  739. Int_t MCtrackID = 0;
  740. if (i234 == 0) {
  741. Int_t f = 0;
  742. //AZ Int_t l = iPass - 3;
  743. //AZ for (int j = 0; j < 4 - l; j++) {
  744. for (int j = 0; j <= layNo[iPass]; j++) {
  745. trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(trNew->GetIndex(0)))->GetTrackID();
  746. Int_t tempTrackID = MCtrackID;
  747. MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  748. if ((tempTrackID != MCtrackID) && (j != 0)) {
  749. f = 1;
  750. }
  751. fout_v << "(" << TMath::ATan2(trNew->GetMeas()[1], trNew->GetMeas()[0])
  752. /*<< " " << TMath::ATan2(v, u)*/
  753. << " " << TMath::ATan2(TMath::Sqrt(trNew->GetMeas()[0] * trNew->GetMeas()[0] + trNew->GetMeas()[1] * trNew->GetMeas()[1]), trNew->GetMeas()[2])
  754. << " " << trNew->GetMeas()[2] << "|" //<< trNew->GetMeas()[1] << "|" << trNew->GetMeas()[2] << " " */
  755. //<< " " << trNew->GetKalmanHit()->GetMeas(0)
  756. << /*trackID <<*/ " " << MCtrackID << ")";
  757. //MpdVector *trOld = trNew;
  758. trNew = trNew->GetPrevTrackPointer();
  759. }
  760. trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(trNew->GetIndex(0)))->GetTrackID();
  761. Int_t tempTrackID = MCtrackID;
  762. MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  763. fout_v << "(" << TMath::ATan2(trNew->GetMeas()[1], trNew->GetMeas()[0])
  764. << " " /*<< TMath::ATan2(v, u) << " "*/
  765. << TMath::ATan2(TMath::Sqrt(trNew->GetMeas()[0] * trNew->GetMeas()[0] + trNew->GetMeas()[1] * trNew->GetMeas()[1]), trNew->GetMeas()[2])
  766. << " " << trNew->GetMeas()[2] << "|" //<< trNew->GetMeas()[1] << "|" << trNew->GetMeas()[2] << /*" " << trackID*/
  767. //<< " " << trNew->GetKalmanHit()->GetMeas(0)
  768. << " " << MCtrackID << ")" << endl;
  769. if (tempTrackID != MCtrackID) {
  770. f = 1;
  771. }
  772. if ( f == 1) {
  773. /// fout_v << "!!!" << endl;
  774. }
  775. count += f;
  776. }
  777. } // for (set<MpdVector*>::iterator itr = intersect.begin();
  778. } // for (Int_t i = 0; i < nTracks;
  779. } // for (Int_t i234 = n; i234 >= 0;
  780. fout_v << "Count of false reconstructed tracks " << count << endl;
  781. fout_v << "- ExtendSecondaryTracks - done" << endl;
  782. }
  783. //_________________________________________________________________________
  784. void MpdVectorFinder::ExtendCellTracks(Int_t iPass)
  785. {
  786. /// Extend cell tracks to layers 4-1
  787. /// pass 2 - start from the last but one layer
  788. //const Int_t nSigm = 3;
  789. int count = 0;
  790. //Double_t errZ = 0.001, errX = 0.001;/// 10um for Z and X // 250um in Z, 23um in R-Phi (local X)
  791. multimap <Float_t, MpdVector*> lMultimapHits[5], tMultimapHits[5];
  792. fout_v << "(transverse angle, longitudinal angle, z | MCtrackID)" << endl;
  793. Int_t n = 6 - 2 - iPass; // layer numbers 0,1,2,3,4
  794. n = TMath::Min (3, n);
  795. for (Int_t i234 = n; i234 >= 0; i234--) { /// loop over layers 3-0
  796. Int_t nKH = 0;//, iprint = 0;
  797. cout << " AZprim: " << i234 << " " << f2DHits[i234]->GetEntriesFast() << endl; //AZ
  798. ///nKH = fKHits[i234]->GetEntriesFast();
  799. // transverse and longitudinal hit maps
  800. for (Int_t i = 0; i < f2DHits[i234]->GetEntriesFast(); i++) {
  801. MpdVector *hit = (MpdVector*) f2DHits[i234]->UncheckedAt(i);
  802. TVector3 meas = hit->GetMeas();
  803. Double_t phi = TMath::ATan2(meas[1], meas[0]);
  804. ///cout << meas[0] << " " << meas[1] << " " << meas[2] << endl;
  805. lMultimapHits[i234].insert(pair<Float_t, MpdVector*>(TMath::ATan2(TMath::Sqrt(meas[0] * meas[0] + meas[1] * meas[1]), meas[2]), hit)); /// longitudinal, theta
  806. tMultimapHits[i234].insert(pair<Float_t, MpdVector*>(phi, hit)); /// transverse, phi
  807. /// possible solution: if transverse angle is close to pi, add it to other side of map with angle - 2pi
  808. if (TMath::Abs(TMath::Pi() - TMath::Abs(phi)) < TMath::Pi() / 9.0) {
  809. tMultimapHits[i234].insert(pair<Float_t, MpdVector*>(phi - TMath::Sign(1, phi) * 2 * TMath::Pi(), hit)); /// transverse, phi, duplicate hit
  810. }
  811. }
  812. int nTracks = fKHits[i234 + 1]->GetEntriesFast(); /// was -1
  813. cout << "primary nTracks for layer " << i234+1 << ": " << nTracks << endl;
  814. for (Int_t i = 0; i < nTracks; i++) {
  815. MpdVector *cand = (MpdVector*) fKHits[i234 + 1]->UncheckedAt(i); /// was -1
  816. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(cand->GetIndex(0));
  817. Int_t trID0 = ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  818. if (i234 == 3) {
  819. ///fout_v << h->GetTrackID() << endl;
  820. ///fout_v << cand->GetUniqueID() << endl;
  821. }
  822. Double_t epsl, epst, pt = 0, dphi = 0;
  823. if (cand->GetPrevTrackPointer() != 0) {
  824. ///fout_v << cand->GetPrevTrackPointer() << endl;
  825. pt = EvalPt(cand->GetPrevTrackPointer(), cand); /// 0-4-5 instead of 0-5-4
  826. if (TMath::Abs(pt) < 0.02) continue;
  827. //pt = TMath::Sign (TMath::Max(0.02,TMath::Abs(pt)), pt);
  828. ///fout_v << "pt " << pt << endl;
  829. //Double_t par0 = 0.00371767;
  830. Double_t par0 = 0.00601267; //AZ-131020: new geometry
  831. //if (i234 < 2) par0 = 0.00243389;
  832. if (i234 < 2) par0 = 0.0013629; //AZ-131020: new geometry
  833. Double_t ptt = TMath::Abs(pt);
  834. dphi = TMath::Sign(par0,pt) / TMath::Max(ptt,0.01);
  835. //if (iPass > 0) ptt = TMath::Min (0.1, ptt);
  836. if (iPass > 0) ptt = TMath::Min (0.2, ptt); //AZ-131020: new geometry
  837. //epst = 0.002 / TMath::Power (ptt-0.01, 2./3);
  838. epst = 0.002 / TMath::Power (ptt-0.015, 7./8); //AZ-131020: new geometry
  839. if (iPass == 0) epst = TMath::Min (epst, 0.01);
  840. //else if (ptt < 0.1) epst *= 2;
  841. //AZ epsl = 1e-2;//5e-3;//4, 15 /// 5e-3
  842. epsl = 0.003 / TMath::Power (ptt-0.01, 2./3);
  843. if (iPass == 0) epsl = TMath::Min (epsl, 0.01);
  844. //else if (ptt < 0.1) epsl *= 2;
  845. ///epst = 5e-2;
  846. ///epsl = 5e-2;
  847. } else {
  848. epst = 0.05; //0.2; // 5e-2;
  849. epsl = 0.01; //0.05; //2e-2;
  850. if (iPass > 0) {
  851. epst = 0.2; //0.2;
  852. if (iPass == 2) epst = 0.3; //AZ - 131020: new geo - larger distance to previous layer
  853. epsl = 0.05; //0.05;
  854. }
  855. }
  856. ///fout_v << "Cut parameters: transverse " << epst << "; longitudinal " << epsl << ";" << endl;
  857. TVector3 meas = cand->GetMeas();
  858. Double_t lastTrackL = TMath::ATan2(TMath::Sqrt(meas[0] * meas[0] + meas[1] * meas[1]), meas[2]);
  859. multimap<Float_t, MpdVector*>::iterator itlow = lMultimapHits[i234].lower_bound(lastTrackL - epsl);
  860. multimap<Float_t, MpdVector*>::iterator ittop = lMultimapHits[i234].upper_bound(lastTrackL + epsl);
  861. Double_t lastTrackT = TMath::ATan2(meas[1], meas[0]);
  862. lastTrackT -= dphi;
  863. multimap<Float_t, MpdVector*>::iterator ittlow = tMultimapHits[i234].lower_bound(lastTrackT - epst);
  864. multimap<Float_t, MpdVector*>::iterator itttop = tMultimapHits[i234].upper_bound(lastTrackT + epst);
  865. set<MpdVector*> l, t, intersect;
  866. for (multimap<Float_t, MpdVector*>::iterator itr = itlow; itr != ittop; ++itr) {
  867. l.insert((*itr).second);
  868. }
  869. for (multimap<Float_t, MpdVector*>::iterator itr = ittlow; itr != itttop; ++itr) {
  870. t.insert((*itr).second);
  871. }
  872. set_intersection(l.begin(), l.end(), t.begin(), t.end(), std::inserter(intersect, intersect.begin()));
  873. ///fout_v << "layer " << i234 << " long in window " << l.size() << " transverse in window " << t.size() << " intersection " << intersect.size() << " pt " << pt << endl;
  874. for (set<MpdVector*>::iterator itr = intersect.begin(); itr != intersect.end(); ++itr) {
  875. MpdVector *track = (*itr);
  876. MpdItsHit5spd *hh = (MpdItsHit5spd*) fItsHits->UncheckedAt(track->GetIndex(0));
  877. Int_t trID = ((MpdStsPoint*) fItsPoints->UncheckedAt(hh->GetTrackID()))->GetTrackID();
  878. if (fExact && trID != trID0) continue; // exact ID match
  879. MpdVector *trNew = nullptr;
  880. trNew = new ((*fKHits[i234])[nKH++]) MpdVector(*track); /// was h->Layer()
  881. trNew->SetUniqueID(nKH - 1);
  882. trNew->SetNofDim(2);
  883. trNew->SetCosSin(0,0);
  884. trNew->SetCosSin(1,0);
  885. trNew->SetPrevTrackPointer(cand);
  886. trNew->SetCode(cand->GetCode());
  887. trNew->SetCode(i);
  888. Int_t trackID = 0;
  889. Int_t MCtrackID = 0;
  890. if (i234 == 0) {
  891. Int_t f = 0;
  892. Int_t ll = iPass;
  893. if (iPass > 0) ll -= 1;
  894. for (int j = 0; j < 4 - ll; j++) {
  895. trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(trNew->GetIndex(0)))->GetTrackID();
  896. Int_t tempTrackID = MCtrackID;
  897. MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  898. if ((tempTrackID != MCtrackID) && (j != 0)) {
  899. f = 1;
  900. }
  901. fout_v << "(" << TMath::ATan2(trNew->GetMeas()[1], trNew->GetMeas()[0])
  902. /*<< " " << TMath::ATan2(v, u)*/
  903. << " " << TMath::ATan2(TMath::Sqrt(trNew->GetMeas()[0] * trNew->GetMeas()[0] + trNew->GetMeas()[1] * trNew->GetMeas()[1]), trNew->GetMeas()[2])
  904. << " " << trNew->GetMeas()[2] << "|" //<< trNew->GetMeas()[1] << "|" << trNew->GetMeas()[2] << " " */
  905. //<< " " << trNew->GetKalmanHit()->GetMeas(0)
  906. << /*trackID <<*/ " " << MCtrackID << ")";
  907. MpdVector *trOld = trNew;
  908. trNew = trNew->GetPrevTrackPointer();
  909. }
  910. trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(trNew->GetIndex(0)))->GetTrackID();
  911. Int_t tempTrackID = MCtrackID;
  912. MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  913. fout_v << "(" << TMath::ATan2(trNew->GetMeas()[1], trNew->GetMeas()[0])
  914. << " " /*<< TMath::ATan2(v, u) << " "*/
  915. << TMath::ATan2(TMath::Sqrt(trNew->GetMeas()[0] * trNew->GetMeas()[0] + trNew->GetMeas()[1] * trNew->GetMeas()[1]), trNew->GetMeas()[2])
  916. << " " << trNew->GetMeas()[2] << "|" //<< trNew->GetMeas()[1] << "|" << trNew->GetMeas()[2] << /*" " << trackID*/
  917. //<< " " << trNew->GetKalmanHit()->GetMeas(0)
  918. << " " << MCtrackID << ")" << endl;
  919. if (tempTrackID != MCtrackID) {
  920. f = 1;
  921. }
  922. if ( f == 1) {
  923. /// fout_v << "!!!" << endl;
  924. }
  925. count += f;
  926. }
  927. }
  928. }
  929. } // for (Int_t i234 = 1; i234 < 4;)
  930. fout_v << "Count of false reconstructed tracks " << count << endl;
  931. fout_v << "- ExtendCellTracks - done" << endl;
  932. }
  933. //__________________________________________________________________________
  934. void MpdVectorFinder::GetTrackSeeds(Int_t iPass)
  935. {
  936. /// Build ITS track seeds from Cell tracks
  937. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  938. Int_t nCand = 0;
  939. TVector3 vert(0.,0.,0.), pmom;
  940. MpdKalmanHit hit;
  941. hit.SetType(MpdKalmanHit::kFixedR);
  942. Int_t layEnd = 0;
  943. if (iPass == fNPass-2) layEnd = 1;
  944. set<MpdVector*> setVecs;
  945. // Loop over last layers
  946. for (Int_t lay = 0; lay <= layEnd; ++lay) {
  947. Int_t nTracks = fKHits[lay]->GetEntriesFast();
  948. fout_v << "seed ITS tracks:"<< nTracks << " fCellMap size " << fCellMap.size() << endl;
  949. for (Int_t itr = 0; itr < nTracks; ++itr) {
  950. map<Int_t,MpdVector*> mapTr;
  951. MpdVector *vec = (MpdVector*) fKHits[lay]->UncheckedAt(itr);
  952. if (lay && setVecs.count(vec)) continue; // hit has been already considered
  953. mapTr[vec->GetLayer()] = vec;
  954. // Check if this track has not been checked already during earlier passes
  955. /// if (/*iPass &&*/ fCellMap.find(track5->GetCode()) != fCellMap.end()) {
  956. /// cout << " Found: " << track5->GetCode() << " " << fCellMap[track5->GetCode()] << endl;
  957. /// continue;
  958. /// }
  959. fCellMap.insert(pair<TString,Int_t>(vec->GetCode(),1));
  960. ///fout_v << "track " << itr << "; track1 code " << track1->GetCode() << endl; /// TODO not useful
  961. for (Int_t jl = 0; jl < 4; ++jl) {
  962. vec = vec->GetPrevTrackPointer();
  963. if (vec == nullptr) break;
  964. mapTr[vec->GetLayer()] = vec;
  965. // Save vector from layer 1 not to double-use it
  966. if (layEnd > 0 && vec->GetLayer() == 1) setVecs.insert(vec);
  967. }
  968. //MpdVector *track2 = track1->GetPrevTrackPointer();
  969. //MpdVector *track3 = track2->GetPrevTrackPointer();
  970. //MpdVector *track4 = track3->GetPrevTrackPointer();
  971. //MpdVector *track5 = track4->GetPrevTrackPointer();
  972. ///fout_v << "x measure first - last: " << track1->GetMeas().X() << " " << track5->GetMeas().X() << endl;
  973. //AZ MpdItsKalmanTrack *track = new ((*fTrackCand)[nCand++]) MpdItsKalmanTrack(*track1, vert, iPass); ///TODO vert is not used?
  974. ///fout_v << "trcand noftrhits " << track->GetNofTrHits() << endl; /// track->GetNofTrHits mostly = 5
  975. //AZ Double_t pt = EvalPt(track2,track1); /// TODO maybe some other hits should be used for pt evaluation
  976. Double_t pt = 0;
  977. //if (track5) pt = EvalPt(track5,track4,track1); /// TODO maybe some other hits should be used for pt evaluation - AZ
  978. //else pt = EvalPt(track4,track3,track1); //AZ
  979. //AZ-241020 if (mapTr.count(4)) pt = EvalPt(mapTr[4],mapTr[3],mapTr.begin()->second); /// TODO maybe some other hits should be used for pt evaluation - AZ
  980. //AZ-241020 else pt = EvalPt(mapTr[3],mapTr[2],mapTr.begin()->second); //AZ
  981. MpdVector* vecs[3]; //uninitialized pointer array?
  982. Int_t iv = 2;
  983. map<Int_t,MpdVector*>::reverse_iterator rmit = mapTr.rbegin();
  984. if (mapTr.size() == 5) ++rmit; // skip one hit
  985. for ( ; rmit != mapTr.rend(); ++rmit) {
  986. vecs[iv--] = rmit->second;
  987. if (iv < 0) break;
  988. }
  989. pt = EvalPt(vecs[2],vecs[1],vecs[0]); //AZ
  990. if (TMath::Abs(pt) < 0.025) continue; //AZ
  991. //MpdItsKalmanTrack *track = new ((*fTrackCand)[nCand++]) MpdItsKalmanTrack(*track1, vert); ///TODO vert is not used?
  992. MpdItsKalmanTrack *track = new ((*fTrackCand)[nCand++]) MpdItsKalmanTrack(*(mapTr.begin()->second), vert); ///TODO vert is not used?
  993. fVectorCands.push_back(mapTr);
  994. fout_v << "GetTrackSeeds::pt " << pt << endl;
  995. //TVector3 v1 = track1->GetMeas();
  996. //TVector3 v2 = track2->GetMeas();
  997. //TVector3 v3 = track3->GetMeas();
  998. //TVector3 v4 = track4->GetMeas();
  999. Double_t phiOut, phiIn, rOut, rIn;
  1000. map<Int_t,MpdVector*>::reverse_iterator mit = mapTr.rbegin();
  1001. MpdVector *vecStart = mit->second;
  1002. ++mit;
  1003. MpdVector *vecNext = mit->second;
  1004. TVector3 vstart = vecStart->GetMeas();
  1005. TVector3 vnext = vecNext->GetMeas();
  1006. TVector3 dv = vstart - vnext;
  1007. phiOut = vnext.Phi();
  1008. phiIn = vstart.Phi();
  1009. rOut = vnext.Pt();
  1010. rIn = vstart.Pt();
  1011. ///fout_v << v1.Pt() << " " << v2.Pt() << endl; /// 17.1.2019
  1012. //AZ track->SetUniqueID(itr+1);
  1013. track->SetUniqueID(nCand);
  1014. ///cout << itr + 1 << endl;
  1015. ///< parametrs track in TPC:
  1016. ///< 0: RPhi - coordinate in R-Phi direction
  1017. ///< 1: Z - longitudinal coordinate
  1018. ///< 2: Phi - local azimuthal angle
  1019. ///< 3: Theta - local dip angle (angle w.r.t. the transverse plane)
  1020. ///< 4: q/Pt - signed inverse Pt
  1021. track->SetPos(rIn);//<- ispravleno!
  1022. track->SetParam (4, 1./pt); // q/Pt
  1023. // TODO don't quite understand this
  1024. //AZ track->SetParam (3, TMath::PiOver2() - track5->GetCosSin(1)); //longitudinal angL
  1025. //if (iPass < 3) track->SetParam (3, TMath::PiOver2() - track5->GetCosSin(1)); //AZ
  1026. if (iPass < 3) track->SetParam (3, TMath::PiOver2() - mapTr.rbegin()->second->GetCosSin(1)); //AZ
  1027. //if (0); //AZ-241020 (iPass < 3) track->SetParam (3, TMath::PiOver2() - vecStart->GetCosSin(1)); //AZ
  1028. else {
  1029. //if (-track1->GetDeltaZ() >= 0) track->SetParam (3, TMath::PiOver2() - TMath::ATan(-track1->GetDeltaZ())); //AZ
  1030. //else track->SetParam (3, TMath::PiOver2() - (TMath::Pi()+TMath::ATan(-track1->GetDeltaZ()))); //AZ
  1031. //AZ track->SetParam (3, TMath::PiOver2() - (v2-v1).Theta());
  1032. //track->SetParam (3, TMath::PiOver2() - (v5-v4).Theta());
  1033. track->SetParam (3, TMath::PiOver2() - dv.Theta());
  1034. }
  1035. track->SetParam (2, dv.Phi()); //Phi - rough estimate
  1036. // Adjust Phi
  1037. Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0., 0., 0.);
  1038. Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
  1039. Double_t rCirc = TMath::Abs (pt / factor);
  1040. Double_t ph = TMath::ASin (dv.Pt() / 2 / rCirc); /// v4 - v5 /// was v3 - v4 // v1 - v2
  1041. //AZ-241020 track->SetParam (2, track->GetParam(2) + TMath::Sign(ph,pt));// - TMath::Sign(ph,pt)
  1042. track->SetParam (2, track->GetParam(2) + TMath::Sign(ph,pt));// - TMath::Sign(ph,pt)
  1043. track->SetParam (1, vstart.Z()); // Z - coordinate
  1044. track->SetParam (0, phiIn*rIn); /// phiIn*rIn // length of circle arc angle phiIn
  1045. fout_v << " trackcand num: " << itr + 1 << " track parameters: "
  1046. << track->GetParam(0) << " " << track->GetParam(1) << " " << track->GetParam(2) << " " << track->GetParam(3) << " " << track->GetParam(4) << endl;
  1047. Double_t parOut[4] = {rOut,phiOut,0.,0.};
  1048. Double_t parIn[4] = {rIn,phiIn,0.,0.};
  1049. //EvalCovar(parOut,parIn,track,track5); /// was track4 /// asuming evalcovar works correctly
  1050. EvalCovar(parOut,parIn,track,vecStart); /// was track4 /// asuming evalcovar works correctly
  1051. track->SetPosNew(track->GetPos());
  1052. track->SetParamNew(*track->GetParam());
  1053. track->SetDirection(MpdKalmanTrack::kOutward);
  1054. hit.SetPos(rIn+5.0);///track5->GetMeas().Pt() /// rIn - 0.5 // point further than 5th layer
  1055. MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  1056. ///fout_v << "GetTrackSeeds " << hit.GetPos() << " " << track->GetPos() << endl;
  1057. ///fout_v << nCand-1 << " " << track->GetTrackID() << endl;
  1058. ///fout_v << track->GetHits()->GetEntriesFast() << " " << track->GetTrHits()->GetEntriesFast() << endl;
  1059. /// TODO 11.7.2019 where does Chi2 come from here
  1060. track->GetHits()->Clear();
  1061. /// track->GetChi2() equals 0 here
  1062. track->SetChi2Its(track->GetChi2()); // temporary storage
  1063. track->SetChi2(0.);
  1064. track->SetDirection(MpdKalmanTrack::kInward); /// was kOutward
  1065. Int_t trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(vecStart->GetIndex(0)))->GetTrackID();
  1066. track->SetTrackID(((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID());
  1067. if (fTrackExact) {
  1068. // For QA
  1069. MpdItsKalmanTrack *trE = new ((*fTrackExact)[fTrackExact->GetEntriesFast()]) MpdItsKalmanTrack(*track);
  1070. //Int_t nh = (Bool_t)track1 + (Bool_t)track2 + (Bool_t)track3 + (Bool_t)track4 + (Bool_t)track5;
  1071. Int_t nh = mapTr.size();
  1072. trE->SetNofHits(nh);
  1073. trE->SetNofIts(iPass);
  1074. trE->SetLastLay(mapTr.begin()->first);
  1075. }
  1076. } // for (Int_t itr = 0; itr < nTracks;
  1077. } // for (Int_t lay = 0; lay <= layEnd;
  1078. fout_v << " Number of ITS track candidates: " << nCand << endl;
  1079. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1080. fout_v << "- GetTrackSeeds - done" << endl;
  1081. }
  1082. //__________________________________________________________________________
  1083. //AZ Double_t MpdVectorFinder::EvalPt(const MpdVector *track1, const MpdVector *track2)
  1084. Double_t MpdVectorFinder::EvalPt(const MpdVector *track1, const MpdVector *track2, const MpdVector *track0, Double_t *circle)
  1085. {
  1086. /// Evaluate signed track Pt (curvature) using 3 points
  1087. TVector3 v1 = track1->GetMeas();///attention!!!
  1088. TVector3 v2 = track2->GetMeas();
  1089. TVector3 v0; //AZ
  1090. if (track0) v0 = track0->GetMeas(); //AZ - for secondaries
  1091. TVector2 vec1(v1.X()-v0.X(),v1.Y()-v0.Y());
  1092. TVector2 vec2(v2.X()-v0.X(),v2.Y()-v0.Y());
  1093. TVector2 vec21 = vec1 - vec2;
  1094. Double_t cosAlpha = vec2 * vec21 / vec2.Mod() / vec21.Mod();
  1095. Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  1096. Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
  1097. Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
  1098. Double_t phi1 = vec1.Phi();
  1099. Double_t phi2 = vec2.Phi();
  1100. //AZ Double_t charge = phi1 - MpdKalmanFilter::Instance()->Proxim(phi1,phi2);
  1101. Double_t charge = 0;
  1102. if (track0) charge = vec21.DeltaPhi(vec2); //AZ
  1103. else charge = phi1 - MpdKalmanFilter::Instance()->Proxim(phi1,phi2);
  1104. if (track1->GetLayer() > track2->GetLayer()) charge = -charge;
  1105. charge = TMath::Sign(1,charge);
  1106. if (track0) {
  1107. //cout << rad * TMath::Sin(TMath::ACos(cosAlpha)) << " " << vec1.Mod()/2 << endl;
  1108. Double_t perp = rad * cosAlpha;
  1109. TVector3 v31 = v1;
  1110. v31 += v0;
  1111. v31 *= 0.5;
  1112. Double_t phi31 = vec1.Phi();
  1113. phi31 -= (charge*TMath::PiOver2());
  1114. TVector2 perpV(TMath::Cos(phi31),TMath::Sin(phi31));
  1115. perpV *= perp;
  1116. TVector2 cent(v31.X(),v31.Y());
  1117. cent += perpV;
  1118. /*
  1119. TVector2 vv3(v1.X(),v1.Y()), vv2(v2.X(),v2.Y());
  1120. vv3 -= cent;
  1121. vv2 -= cent;
  1122. cout << rad << " " << vv3.Mod() << " " << vv2.Mod() << " " << " " << cent.Mod() << " " << iq << endl;
  1123. */
  1124. if (circle) {
  1125. circle[0] = cent.X();
  1126. circle[1] = cent.Y();
  1127. circle[2] = TMath::Abs(rad);
  1128. }
  1129. }
  1130. return -factor * TMath::Abs(rad) * charge;
  1131. }
  1132. ///-----------------------------------------------------------------------
  1133. Bool_t MpdVectorFinder::PropagateHelix(Double_t rlay, TVector3 pos1, Double_t *circle, Double_t &phi)
  1134. {
  1135. // Propagate track as a helix to a circular layer
  1136. // Transverse view - intersection of 2 circles (track and layer)
  1137. Double_t d2 = circle[0] * circle[0] + circle[1] * circle[1]; // dist**2 from (0,0) to track center
  1138. Double_t d = TMath::Sqrt(d2);
  1139. Double_t a = rlay * rlay - circle[2] * circle[2] + d2;
  1140. a /= (2 * d); // dist from (0,0) to cross
  1141. Double_t perp2 = rlay * rlay - a * a;
  1142. if (perp2 < 0) return kFALSE;
  1143. TVector3 crossp(circle[0],circle[1],0);
  1144. crossp *= (a/d);
  1145. Double_t phic = crossp.Phi();
  1146. if (crossp.Cross(pos1).Z() < 0) phic += TMath::PiOver2();
  1147. else phic -= TMath::PiOver2();
  1148. TVector3 perpV(TMath::Cos(phic),TMath::Sin(phic),0);
  1149. perpV *= TMath::Sqrt(perp2);
  1150. TVector3 point(crossp);
  1151. point -= perpV;
  1152. phi = point.Phi();
  1153. return kTRUE;
  1154. }
  1155. ///-----------------------------------------------------------------------
  1156. Double_t MpdVectorFinder::EvalPhiPt(const MpdVector *track1, const MpdVector *track2)
  1157. {
  1158. TVector3 v1 = track1->GetMeas();///attention!!!
  1159. TVector3 v2 = track2->GetMeas();
  1160. TVector2 vec1(v1.X(),v1.Y());
  1161. TVector2 vec2(v2.X(),v2.Y());
  1162. TVector2 vec21 = vec1 - vec2;
  1163. Double_t cosAlpha = vec2 * vec21 / vec2.Mod() / vec21.Mod();
  1164. Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  1165. TVector2 temp(v1.Y(), -v1.X());
  1166. TVector2 center = vec1 / 2. + temp * rad * cosAlpha / temp.Mod();
  1167. temp = vec2 - center;
  1168. TVector2 res(temp.Y(), -temp.X());
  1169. return TMath::ATan2(res.Y(), res.X());
  1170. }
  1171. //__________________________________________________________________________
  1172. void MpdVectorFinder::EvalCovar(Double_t *parOut, Double_t *parIn, MpdItsKalmanTrack *track, const MpdVector *track1)
  1173. {
  1174. /// Evaluate covariance matrix for track seed
  1175. Double_t rIn = parIn[0], phiIn = parIn[1];
  1176. Double_t rOut = parOut[0], phiOut = parOut[1];
  1177. Double_t xIn = rIn * TMath::Cos(phiIn);
  1178. Double_t yIn = rIn * TMath::Sin(phiIn);
  1179. Double_t xOut = rOut * TMath::Cos(phiOut);
  1180. Double_t yOut = rOut * TMath::Sin(phiOut);
  1181. parOut[2] = xOut;
  1182. parOut[3] = yOut;
  1183. parIn[2] = xIn;
  1184. parIn[3] = yIn;
  1185. TMatrixD ww(5,5);
  1186. ww(0,0) = track1->GetErr(0) * track1->GetErr(0); // <RphiRphi> //26.11 error x
  1187. ww(0,0) *= 4.; //extra factor of 4
  1188. ww(1,1) = track1->GetErr(1) * track1->GetErr(1); // <zz> //error z
  1189. ww(1,1) *= 4.; //AZ extra factor of 8
  1190. Double_t dx = parOut[2] - parIn[2], dy = parOut[3] - parIn[3];
  1191. Double_t dist2 = dx * dx + dy * dy;
  1192. Double_t sinPhi = TMath::Sin (track->GetParam(2));
  1193. Double_t cosPhi = TMath::Cos (track->GetParam(2));
  1194. Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
  1195. Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
  1196. ww(2,2) = (pOut * pOut + pIn * pIn) / dist2 * ww(0,0); // <PhiPhi>
  1197. //cout << ww(2,2) << " " << ww(0,0)*2/dist2 << endl;
  1198. ww(2,2) = ww(0,0)*2/dist2; //AZ-251020
  1199. //ww(2,2) *= 2.; // extra factor of 2
  1200. ww(2,2) *= 50.; //AZ
  1201. Double_t tanThe = TMath::Tan(track->GetParam(3));
  1202. Double_t dRad = parOut[0] - parIn[0];
  1203. Double_t denom = dRad * (1.+tanThe*tanThe);
  1204. ww(3,3) = ww(1,1) * 2. / denom / denom; // <TheThe>
  1205. ww(3,3) *= 50; //AZ
  1206. //AZ ww(1,1) *= 8.; //AZ extra factor of 8
  1207. ww(4,4) = (track->GetParam(4)*1.) * (track->GetParam(4)*1.); // error 100%
  1208. //ww(4,4) = 100; //AZ
  1209. Int_t iok = 0;
  1210. TMatrixD wwTmp = ww;
  1211. MpdKalmanFilter::Instance()->MnvertLocal(ww.GetMatrixArray(), 5, 5, 5, iok);
  1212. track->SetWeight(ww);
  1213. TVector3 pos, posOut, pmom;
  1214. MpdItsHit5spd *hit = (MpdItsHit5spd*) fItsHits->UncheckedAt(track1->GetIndex());
  1215. MpdStsPoint *p = (MpdStsPoint*)fItsPoints->UncheckedAt(hit->GetRefIndex());
  1216. p->Position(pos);
  1217. p->PositionOut(posOut);
  1218. p->Momentum(pmom);
  1219. Double_t phi = pmom.Phi();
  1220. //Double_t th = pmom.Theta();
  1221. Double_t r = (pos.Pt() + posOut.Pt()) / 2;//radius
  1222. Double_t phi1 = pos.Phi();
  1223. Double_t phi2 = posOut.Phi();
  1224. phi1 += MpdKalmanFilter::Instance()->Proxim(phi1, phi2);
  1225. phi1 /= 2;
  1226. Double_t zzz = (pos.Z() + posOut.Z()) / 2;
  1227. //track->SetParam(2,phi); // !!!!! exact value - just for test
  1228. //track->SetParam(3,TMath::PiOver2()-th); // !!!!! exact value - just for test
  1229. if (lunErr&&p->GetTrackID()>2000) fprintf(lunErr,"%12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e\n",
  1230. r*(MpdKalmanFilter::Instance()->Proxim(phi1,track->GetParam(0)/r)-phi1),
  1231. track->GetParam(1)-zzz,MpdKalmanFilter::Instance()->Proxim(phi,track->GetParam(2))-phi,
  1232. TMath::PiOver2()-pmom.Theta()-track->GetParam(3), 1/track->GetParam(4),
  1233. TMath::Sqrt(wwTmp(0,0)),TMath::Sqrt(wwTmp(1,1)),TMath::Sqrt(wwTmp(2,2)),TMath::Sqrt(wwTmp(3,3)));
  1234. }
  1235. //__________________________________________________________________________
  1236. void MpdVectorFinder::DoTracking(Int_t iPass)
  1237. {
  1238. /// Run Kalman tracking
  1239. fout_v << "- DoTracking -" << endl;
  1240. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1241. Int_t nCand = fTrackCand->GetEntriesFast(), iok = 0;// new 05.03.14
  1242. fout_v << "DoTracking candidates on pass " << iPass << ": " << nCand << endl;
  1243. cout << "DoTracking candidates on pass " << iPass << ": " << nCand << endl;
  1244. Int_t lay0 = ((MpdKalmanHit*)fKHits1->Last())->GetLayer();//ms 06.05 /// TODO should be last layer for mpdvector? was First()
  1245. for (Int_t i = 0; i < nCand; ++i) {
  1246. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);// new 05.03.14
  1247. //if (iPass == 2) fout_v << "fTrackCand GetNofTrHits and GetNofHits " << track->GetNofTrHits() << " " << track->GetNofHits() << endl; /// 17.1.2019
  1248. ///fout_v << " Track seed No. " << i << ", ID: " << track->GetTrackID() << ", Hits: " << track->GetNofTrHits() << endl;
  1249. //for (Int_t j = 0; j < track->GetNofTrHits(); ++j) {
  1250. //MpdKalmanHit *h = (MpdKalmanHit* )track->GetTrHits()->UncheckedAt(j);
  1251. //MpdStsHit *hh = (MpdStsHit*) fItsHits->UncheckedAt(h->GetIndex());
  1252. //Int_t id = ((FairMCPoint*) fItsPoints->UncheckedAt(hh->GetRefIndex()))->GetTrackID();
  1253. ///fout_v << j << " " << h->GetDist() << " " << h->GetLayer() << endl;
  1254. //}
  1255. /// RunKalmanFilterCell makes track->GetNofHits() from 0 to 3-5, adding them if they pass chi2 cut
  1256. iok = RunKalmanFilterCell(track, iPass); /// removed if (fGeo) or smth like that// from cell track
  1257. /// 26.2.20 was only (track->GetNofHits() < 3), without else
  1258. //AZ
  1259. //AZ if (iok == -1) {
  1260. if ((iok == -1) || (track->GetNofHits() < 3)) {
  1261. fTrackCand->RemoveAt(i);
  1262. continue;
  1263. //AZ-161020 } else if ((iPass < 3) && ((track->GetNofHits() < 4) || (track->GetNofHits() == 4) && (track->Pt() < 0.1))) {
  1264. //} else if ((iPass < 4) && ((track->GetNofHits() < 4) || (track->GetNofHits() == 4) && (track->Pt() < 0.1))) {
  1265. } else if ((iPass == fNPass-3) && (track->GetNofHits() < 5)) {
  1266. fTrackCand->RemoveAt(i);
  1267. continue;
  1268. } else if ((iPass < fNPass-1) && ((track->GetNofHits() < 4) || ((track->GetNofHits() == 4) && (track->Pt() < 0.1)))) {
  1269. fTrackCand->RemoveAt(i);
  1270. continue;
  1271. }
  1272. track->SetParam(*track->GetParamNew()); //AZ
  1273. track->SetPos(track->GetPosNew()); //AZ
  1274. /// added 20.2.20
  1275. /* temporary output for found secondary track hits with their layers, distance from IP and original MC track id
  1276. if (iPass >= 3) {
  1277. cout << "GetTrHits!";
  1278. /// get trackCand hits id for iPass 3 and 4
  1279. for (Int_t j = 0; j < track->GetNofTrHits(); ++j) { /// TODO maybe GetNofHits will be more representative
  1280. MpdKalmanHit *h = (MpdKalmanHit*) track->GetTrHits()->UncheckedAt(j);
  1281. Int_t trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(h->GetIndex(0)))->GetTrackID();
  1282. Int_t MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  1283. cout << "(" << h->GetLayer() << ", " << h->GetPos() << ", " << MCtrackID << ") ";
  1284. }
  1285. cout << " " << track->GetChi2() << endl;
  1286. // cout << "GetHits!";
  1287. // for (Int_t j = 0; j < track->GetNofHits(); ++j) { /// TODO maybe GetNofHits will be more representative
  1288. // MpdKalmanHit *h = (MpdKalmanHit*) track->GetHits()->UncheckedAt(j);
  1289. // Int_t trackID = ((MpdItsHit5spd*) fItsHits->UncheckedAt(h->GetIndex(0)))->GetTrackID();
  1290. // Int_t MCtrackID = ((MpdStsPoint*) fItsPoints->UncheckedAt(trackID))->GetTrackID();
  1291. // cout << "(" << h->GetLayer() << ", " << h->GetPos() << ", " << MCtrackID << ") ";
  1292. // }
  1293. // cout << endl << endl;
  1294. }
  1295. */
  1296. // Propagate track to the beam line
  1297. /// fSa stands for Stand Alone Flag. If it's value == 1 then no matching is done.
  1298. Bool_t okk = kTRUE;
  1299. if (fSa) okk = GoToBeamLine(track);
  1300. if (!okk) fTrackCand->RemoveAt(i);
  1301. } // for (Int_t i = 0; i < nCand;
  1302. fTrackCand->Compress();
  1303. fout_v << "DoTracking candidates after removal " << fTrackCand->GetEntriesFast() << endl;
  1304. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1305. }
  1306. //__________________________________________________________________________
  1307. Bool_t MpdVectorFinder::GoToBeamLine(MpdItsKalmanTrack *track)
  1308. {
  1309. // Go to beam line
  1310. Double_t vert[3] = {0.0, 0.0, 0.0};
  1311. TString mass2 = "0.0194797849"; // pion mass squared
  1312. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07};
  1313. MpdKalmanGeoScheme *geoScheme = MpdKalmanFilter::Instance()->GetGeo();
  1314. MpdKalmanHit *hitp;
  1315. TMatrixD param(5,1);
  1316. TMatrixDSym weight(5), pointWeight(5);
  1317. if (track->GetDirection() == MpdKalmanTrack::kOutward) {
  1318. // Going outward - refit inward
  1319. track->SetDirection(MpdKalmanTrack::kInward);
  1320. TMatrixDSym *w = track->GetWeight();
  1321. Int_t nHits = track->GetNofHits(), nHits2 = nHits * nHits;
  1322. track->SetParam(*track->GetParamNew());
  1323. track->SetPos(track->GetPosNew());
  1324. track->SetNode(track->GetNodeNew());
  1325. track->SetChi2(0);
  1326. for (Int_t i = 0; i < 5; ++i) {
  1327. for (Int_t j = i; j < 5; ++j) {
  1328. if (j == i) (*w)(i,j) /= nHits2;
  1329. else (*w)(i,j) = (*w)(j,i) = 0.;
  1330. }
  1331. }
  1332. //(*w)(4,4) = 1 / track->GetParamNew(4) / track->GetParamNew(4);
  1333. Int_t ibeg = track->GetHits()->GetEntriesFast() - 1, iend = -1, iDir = -1;
  1334. for (Int_t i = ibeg; i != iend; i+=iDir) {
  1335. hitp = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
  1336. if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hitp, kFALSE, kTRUE)) return kFALSE;
  1337. if (1./TMath::Abs(track->GetParamNew(4)) < 0.025) return kFALSE;
  1338. //MpdKalmanFilter::Instance()->PropagateToHit(track, hitp, kFALSE, kTRUE);
  1339. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hitp, pointWeight, param);
  1340. track->SetChi2(track->GetChi2()+dChi2);
  1341. weight = *track->GetWeight();
  1342. weight += pointWeight;
  1343. track->SetWeight(weight);
  1344. track->SetParamNew(param);
  1345. // Add scattering in Si
  1346. Int_t lay = hitp->GetLayer();
  1347. cout << " aaaa " << hitp->GetPos() << " " << lay << " " << dChi2 << " " << 1/param(4,0) << endl;
  1348. Double_t x0 = 9.36; // rad. length
  1349. TMatrixDSym *cov = track->Weight2Cov();
  1350. Double_t th = track->GetParamNew(3);
  1351. Double_t cosTh = TMath::Cos(th);
  1352. TVector3 norm = geoScheme->Normal(hitp);
  1353. //norm[2] = 0.0;
  1354. norm.SetMag(1.0);
  1355. MpdKalmanTrack trTmp(*track);
  1356. trTmp.SetParam(*trTmp.GetParamNew());
  1357. TVector3 mom3 = trTmp.Momentum3();
  1358. //mom3[2] = 0.0;
  1359. mom3.SetMag(1.0);
  1360. Double_t cosPhi = norm.Dot(mom3);
  1361. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick[lay]*3/x0/cosPhi, mass2);
  1362. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, 2*thick[lay]/cosPhi, mass2);
  1363. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1364. (*cov)(3,3) += angle2;
  1365. Int_t iok = 0;
  1366. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1367. track->SetWeight(*cov);
  1368. }
  1369. } // if (track->GetDirection() == MpdKalmanTrack::kOutward)
  1370. //AZ track->SetParam(*track->GetParamNew());
  1371. //AZ track->SetPos(track->GetPosNew());
  1372. Double_t pos = track->GetPosNew();
  1373. TMatrixD par = *track->GetParamNew();
  1374. TMatrixDSym cov = *track->Weight2Cov();
  1375. Double_t leng = track->GetLength();
  1376. TString nodeNew = track->GetNodeNew();
  1377. // Go to beam pipe
  1378. MpdKalmanHit hit;
  1379. hit.SetType(MpdKalmanHit::kFixedR);
  1380. hit.SetPos(fPipeR);
  1381. Int_t iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  1382. //if (1./TMath::Abs(track->GetParamNew(4)) < 0.025) return kFALSE;
  1383. if (iok != 1) {
  1384. // Restore track
  1385. track->SetParam(par);
  1386. track->SetParamNew(par);
  1387. track->SetCovariance(cov);
  1388. track->ReSetWeight();
  1389. track->SetPos(pos);
  1390. track->SetPosNew(pos);
  1391. track->SetLength(leng);
  1392. track->SetNodeNew(nodeNew);
  1393. } else {
  1394. // Add multiple scattering
  1395. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  1396. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  1397. TMatrixDSym* pcov = track->Weight2Cov();
  1398. Double_t th = track->GetParamNew(3);
  1399. Double_t cosTh = TMath::Cos(th);
  1400. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  1401. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  1402. (*pcov)(3,3) += angle2;
  1403. Int_t ok = 0;
  1404. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  1405. track->SetWeight(*pcov);
  1406. }
  1407. cov = *track->Weight2Cov();
  1408. hit.SetPos(0.);
  1409. hit.SetMeas(0,track->GetParam(2)); // track Phi
  1410. //AZ-030121 pos = track->GetPosNew(); //AZ-231020
  1411. //fout_v << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  1412. //AZ-030121 iok = 0; //AZ-231020 MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  1413. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  1414. if (iok != 1) {
  1415. track->SetNodeNew(""); // for FindPca
  1416. MpdKalmanFilter::Instance()->FindPca(track, vert);
  1417. }
  1418. //AZ-030121 track->SetNode(track->GetNodeNew()); //AZ-231020
  1419. //AZ-030121 track->SetPos(pos); //AZ-231020
  1420. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  1421. return kTRUE;
  1422. }
  1423. //__________________________________________________________________________
  1424. Int_t MpdVectorFinder::RunKalmanFilterCell(MpdItsKalmanTrack *track, Int_t iPass)
  1425. {
  1426. /// Run Kalman filter (fitter) for the hits from the cell track
  1427. /// (might not work when propagating outward!!!)
  1428. TString mass2 = "0.0194797849"; // pion mass squared
  1429. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07};
  1430. MpdKalmanGeoScheme *geoScheme = MpdKalmanFilter::Instance()->GetGeo();
  1431. Int_t layMax = ((MpdKalmanHit*)fKHits1->First())->GetLayer();
  1432. MpdKalmanHit *hitOK = 0x0;
  1433. MpdKalmanHit hitTmp;
  1434. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  1435. //Int_t layBeg = 0, layEnd = -1, dLay = -1, layOK = -1; /// old coment
  1436. //Int_t layEnd = -1, dLay = -1, layOK = -1;
  1437. //if (trackDir == MpdKalmanTrack::kInward) { /// was kOutward
  1438. // layEnd = layMax + 1;
  1439. // dLay = 1;
  1440. //}
  1441. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  1442. TMatrixD param(5,1), paramTmp(5,1);
  1443. Double_t posNew = 0.0; //saveZ = 0.0, saveLeng = 0.0, dChi2Min = 0.0,
  1444. Int_t ok = 0;
  1445. /// TODO what is happening here?
  1446. ///fout_v << track->GetUniqueID() << endl;
  1447. /*AZ-211020
  1448. MpdVector *cellTr = nullptr;
  1449. MpdVector *tracks[5];
  1450. tracks[0] = (MpdVector*) fKHits[0]->UncheckedAt(track->GetUniqueID() - 1); // layer 1
  1451. ///fout_v << track->GetUniqueID() - 1;
  1452. /// added 19.2.2020, iPass changed to l for pass #2 purposes
  1453. Int_t l = iPass;
  1454. //AZ if (iPass >= 2) l = iPass - 2;
  1455. if (iPass > 0) --l;
  1456. if (iPass >= 3) l = iPass - 3;
  1457. ///
  1458. for (Int_t i = 1; i <= 4 - l; i++) { // for layers 2-5
  1459. tracks[i] = (MpdVector*) fKHits[i]->UncheckedAt(tracks[i - 1]->GetPrevTrackPointer()->GetUniqueID());
  1460. }
  1461. */
  1462. //for (Int_t lay = 4 - l; lay >= 0; lay--) {
  1463. // Get CellTrack from the inner layer
  1464. //cellTr = (MpdVector*) fKHits[lay]->UncheckedAt(tracks[lay]->GetUniqueID()); // for layers 1-5
  1465. map<Int_t,MpdVector*> &mapTr = fVectorCands[track->GetUniqueID()-1];
  1466. for (map<Int_t,MpdVector*>::reverse_iterator mit = mapTr.rbegin(); mit != mapTr.rend(); ++mit) {
  1467. //for (map<Int_t,MpdVector*>::iterator mit = mapTr.begin(); mit != mapTr.end(); ++mit) {
  1468. ///if (lay == 0)
  1469. /// fout_v << "pt " << EvalPt(cellTr,cellTr->GetPrevTrackPointer()) << endl;
  1470. ///MpdKalmanHit *hit = (MpdKalmanHit*) fKHits1->UncheckedAt(cellTr->GetIndex(0));
  1471. //AZ MpdKalmanHit *hit = cellTr->GetKalmanHit();
  1472. MpdKalmanHit *hit = mit->second->GetKalmanHit();
  1473. ///fout_v << hit << " " << hit->GetUniqueID() << endl;
  1474. ///fout_v << track->GetPos() << " " << TMath::Sqrt(cellTr->GetMeas()[0] * cellTr->GetMeas()[0] + cellTr->GetMeas()[1] * cellTr->GetMeas()[1] + cellTr->GetMeas()[2] * cellTr->GetMeas()[2]) << endl;
  1475. // Propagate to hit (if it is not very close to the track)
  1476. if (TMath::Abs(hit->GetPos()-track->GetPosNew()) > 1.e-4) {
  1477. Double_t leng = track->GetLength();
  1478. posNew = track->GetPosNew();
  1479. TMatrixD parNew = *track->GetParamNew();
  1480. TString nodeNew = track->GetNodeNew();
  1481. TString curPath = track->GetNode();
  1482. ///fout_v << track->GetLength() << " " << track << " " << hit << endl;
  1483. /// PropagateToHit updates track->GetLength()
  1484. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit,kTRUE,kTRUE)) {
  1485. // Restore initial parameters for the failed track
  1486. ///cout << "restore" << endl; /// TODO this is happening too often
  1487. track->SetPosNew(posNew);
  1488. track->SetParamNew(parNew);
  1489. track->SetLength(leng);
  1490. track->SetNodeNew(nodeNew);
  1491. track->SetNode(curPath);
  1492. ok = -1;
  1493. break;
  1494. }
  1495. //Double_t step = track->GetLength() - leng;
  1496. /// fout_v << "!!!kalmanfiltercell" << " " << hit->GetLayer() << " " << track->GetLength() << " " << leng << " " << step << " " << " " << TMath::Sqrt(hit->GetMeas(0) * hit->GetMeas(0) + hit->GetMeas(1) * hit->GetMeas(1)) << " " << track->GetPosNew() << " " << track->GetUniqueID() << endl;
  1497. } // if (TMath::Abs(hit->GetPos()-track->GetPosNew()) > 1.e-4)
  1498. // Exclude used hits!
  1499. if (hit->GetFlag() < 0) continue; //new 04.03 work 15.04
  1500. // !!! Exact ID match
  1501. if (fExact && TrackID(hit) != track->GetTrackID()) continue;
  1502. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track,hit,pointWeight,param);
  1503. //fout_v << "dChi2 " << dChi2 << endl;
  1504. //AZ if ((TMath::Abs(dChi2) < fgkChi2Cut) || (iPass == 2)) {
  1505. if (TMath::Abs(dChi2) < fgkChi2Cut) { ///AZ
  1506. //if (iPass == 2) fout_v << "chi2 cut passed" << endl;
  1507. track->GetHits()->Add(hit); /// this changes GetNofHits() result
  1508. track->SetChi2(track->GetChi2()+dChi2);
  1509. TMatrixDSym w = *track->GetWeight();
  1510. w += pointWeight;
  1511. track->SetWeight(w);
  1512. track->SetParamNew(param);
  1513. // Save track params at last hit
  1514. track->SetLengAtHit(track->GetLength());
  1515. track->SetParamAtHit(param);
  1516. track->SetWeightAtHit(*track->GetWeight());
  1517. track->SetPosAtHit(track->GetPosNew());
  1518. }
  1519. // Add scattering in Si
  1520. Int_t lay = hit->GetLayer();
  1521. if (1) {
  1522. Double_t x0 = 9.36; // rad. length
  1523. //Double_t x0 = 30420; // rad. length of Air
  1524. TMatrixDSym *cov = track->Weight2Cov();
  1525. Double_t th = track->GetParamNew(3);
  1526. Double_t cosTh = TMath::Cos(th);
  1527. TVector3 norm = geoScheme->Normal(hit);
  1528. //norm[2] = 0.0;
  1529. norm.SetMag(1.0);
  1530. MpdKalmanTrack trTmp(*track);
  1531. trTmp.SetParam(*trTmp.GetParamNew());
  1532. TVector3 mom3 = trTmp.Momentum3();
  1533. //mom3[2] = 0.0;
  1534. mom3.SetMag(1.0);
  1535. Double_t cosPhi = norm.Dot(mom3);
  1536. //AZ Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, 0.005*2/x0, mass2);
  1537. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick[lay]*2/x0, mass2); //AZ-161020
  1538. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick[lay]*3/x0, mass2); //AZ-231020
  1539. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick[lay]*3/x0/cosPhi, mass2); //AZ-231020
  1540. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, 3*thick[lay]/cosPhi, mass2); //AZ-231020
  1541. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1542. (*cov)(3,3) += angle2;
  1543. Int_t iok = 0;
  1544. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1545. track->SetWeight(*cov);
  1546. }
  1547. /*
  1548. Double_t nCables = 0, x0 = 0.0116 * 2 / 9.36; // in rad. length - 116um cable per side
  1549. // Find number of cables crossed
  1550. TString path = gGeoManager->GetPath();
  1551. ///fout_v << path << endl;
  1552. if (!path.Contains("sensor") && !path.Contains("sector")) {
  1553. /// fout_v << " !!! MpdVector::RunKalmanFilter - Outside detector !!! " << endl;
  1554. exit(0);
  1555. }
  1556. Double_t v7[3] = {track->GetParamNew(0), track->GetPosNew(), track->GetParamNew(1)}, v77[3];
  1557. gGeoManager->LocalToMaster(v7,v77);
  1558. Double_t zTr = TMath::Abs (v77[2]); // global Z
  1559. map<Double_t,Double_t>::iterator it;
  1560. for (it = fCables[lay/2].begin(); it != fCables[lay/2].end(); ++it) {
  1561. if (zTr < it->first || zTr > it->second) continue;
  1562. ++nCables;
  1563. }
  1564. if (nCables) {
  1565. x0 *= nCables;
  1566. TMatrixDSym *cov = track->Weight2Cov();
  1567. Double_t th = track->GetParamNew(3);
  1568. Double_t cosTh = TMath::Cos(th);
  1569. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0);
  1570. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  1571. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1572. (*cov)(3,3) += angle2;
  1573. Int_t iok = 0;
  1574. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1575. track->SetWeight(*cov);
  1576. }
  1577. */
  1578. //} // for (Int_t lay = 0; lay < 5; ++ lay)
  1579. } // for (map<Int_t,MpdVector*>::reverse_iterator mit = mapTr.rbegin();
  1580. //AZ return 0;
  1581. return ok;
  1582. }
  1583. //__________________________________________________________________________
  1584. void MpdVectorFinder::RemoveDoubles()
  1585. {
  1586. /// Remove double tracks (keep the ones with better quality)
  1587. Int_t ntracks = fTrackCand->GetEntriesFast(); //new 05.03.14
  1588. fout_v << "- RemoveDoubles -" << endl;
  1589. fout_v << " Total tracks: " << ntracks << endl;
  1590. /// added new version 21.1.2019
  1591. // Fill map with track quality
  1592. /// MpdVector
  1593. multimap<Double_t,MpdItsKalmanTrack*> qualMap;
  1594. multimap<Double_t,MpdItsKalmanTrack*>::iterator mit1, mit2;
  1595. for (Int_t i = 0; i < ntracks; i++) {
  1596. MpdItsKalmanTrack *tr = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1597. //AZ Double_t quality = tr->GetNofHits() + TMath::Min(tr->GetChi2(),999.0) / 1000.0;
  1598. Double_t quality = tr->GetNofHits() + (1 - TMath::Min(tr->GetChi2(),999.0) / 1000.0);
  1599. qualMap.insert(pair<Double_t,MpdItsKalmanTrack*>(-quality,tr));
  1600. }
  1601. MpdItsKalmanTrack *tr1, *tr2;
  1602. for (mit1 = qualMap.begin(); mit1 != qualMap.end(); ++mit1) {
  1603. ///fout_v << mit1->first << " " << mit1->second << endl;
  1604. tr1 = mit1->second;
  1605. if (tr1 == nullptr) continue; // killed track
  1606. mit2 = mit1;
  1607. ++mit2;
  1608. for ( ; mit2 != qualMap.end(); ++mit2) {
  1609. tr2 = mit2->second;
  1610. if (tr2 == nullptr) continue; // killed track
  1611. bool x = AreTracksDoubles(tr1, tr2);
  1612. ///fout_v << x << endl;
  1613. if (!x) continue;
  1614. ///fout_v << "tracks doubles" << endl;
  1615. if (mit1->first < mit2->first) { fTrackCand->Remove(tr2); mit2->second = nullptr; }
  1616. else {
  1617. if (tr1->GetChi2() <= tr2->GetChi2()) { fTrackCand->Remove(tr2); mit2->second = nullptr; }
  1618. else { fTrackCand->Remove(tr1); mit1->second = nullptr; break; }
  1619. }
  1620. }
  1621. }
  1622. fTrackCand->Compress();
  1623. fNTracks = fTrackCand->GetEntriesFast();
  1624. }
  1625. //__________________________________________________________________________
  1626. Bool_t MpdVectorFinder::AreTracksDoubles(MpdItsKalmanTrack *tr1, MpdItsKalmanTrack *tr2)
  1627. {
  1628. /// Searching common hits in 2 tracks to determine doubles
  1629. // track1 contains fewer hits than track2
  1630. /// MpdVector?
  1631. MpdItsKalmanTrack *track1 = tr1, *track2 = tr2;
  1632. if (tr1->GetNofHits() > tr2->GetNofHits()) //12.02
  1633. track1 = tr2, track2 = tr1;
  1634. Int_t limCommonPoint = (track1->GetNofHits()+1) / 2; // at least so many common hits should be found
  1635. limCommonPoint = TMath::Min (limCommonPoint, 2); //AZ 2 common hits at max
  1636. TObjArray *hits1 = track1->GetHits(), *hits2 = track2->GetHits();
  1637. Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast(), nHitsCommon = 0, j = nh2 - 1;
  1638. for (Int_t i = nh1 - 1; i >= 0; i--){
  1639. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
  1640. for ( ; j >= 0; j--){
  1641. MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
  1642. // Is the hit common for two tracks?
  1643. if (hit1 == hit2) {
  1644. nHitsCommon++;
  1645. if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
  1646. break;
  1647. }
  1648. if (hit2->GetLayer() > hit1->GetLayer()) break; // already closer to beam /// was < /// what does this do?
  1649. //if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam /// was < /// what does this do?
  1650. }
  1651. if ((nh1 - i - 1) + limCommonPoint - nHitsCommon > nh1) return kFALSE; // there'll be not enough common hits already
  1652. }
  1653. if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
  1654. return kTRUE;
  1655. }
  1656. //__________________________________________________________________________
  1657. Int_t MpdVectorFinder::TrackID(MpdKalmanHit *hit, Int_t indx)
  1658. {
  1659. /// Return track ID of the hit
  1660. //AZ-141020 FairHit *h = (FairHit*) fItsHits->UncheckedAt(hit->GetIndex(indx));
  1661. //AZ-141020 return ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1662. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex(indx));
  1663. return ((MpdStsPoint*) fItsPoints->UncheckedAt(h->GetTrackID()))->GetTrackID();
  1664. }
  1665. //__________________________________________________________________________
  1666. TVector2 MpdVectorFinder::GetDistance(MpdKalmanTrack *track, MpdKalmanHit *hit)
  1667. {
  1668. /// Compute distance between track and hit
  1669. Int_t lay = hit->GetLayer();
  1670. Int_t lay2 = lay / 2;
  1671. //Int_t side = lay % 2;
  1672. ///Int_t module0 = ((MpdStsHit*) fItsHits->UncheckedAt(hit->GetIndex()))->Module();
  1673. Int_t module0 = ((MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex()))->Module();
  1674. Double_t zTr = track->GetParamNew(1);
  1675. Double_t zloc = zTr + fDz[lay2];
  1676. Int_t modTr = Int_t (zloc / fZmod[lay2]);
  1677. Int_t dMod = modTr - module0;
  1678. Double_t dZ = 0;
  1679. if (dMod != 0) {
  1680. // Not in the same module0 - check Z-distance to module0 edge
  1681. if (dMod > 0) dZ = zloc - (module0+1) * fZmod[lay2];
  1682. else dZ = zloc - module0 * fZmod[lay2];
  1683. if (TMath::Abs(dMod) > 2) return TVector2(TMath::Abs(dZ),999.); // not in neighbour modules
  1684. }
  1685. // Translate transverse track position to local system
  1686. Double_t xloc = track->GetParamNew(0) * hit->GetCosSin(0) + zTr * hit->GetCosSin(1);
  1687. Double_t phTr = xloc / track->GetPosNew();
  1688. Double_t phHit = hit->GetMeas(0) / fRad[lay];
  1689. //TVector2 dist(TMath::Abs(dZ), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit)-phTr));
  1690. TVector2 dist(TMath::Abs(dZ),
  1691. TMath::Abs(MpdKalmanFilter::Instance()->Proxim(phTr,phHit,hit->GetCosSin(0))-phTr));
  1692. //TVector2 dist(TMath::Abs(zTr-hit->GetZ()), TMath::Abs(MpdKalmanFilter::Instance()->Proxim(track,hit)/hit->GetR()-ph));
  1693. return dist;
  1694. }
  1695. //__________________________________________________________________________
  1696. void MpdVectorFinder::Write()
  1697. {
  1698. /// Write
  1699. TFile histoFile("ItsRec.root","RECREATE");
  1700. Writedir2current(fHistoDir);
  1701. histoFile.Close();
  1702. }
  1703. //__________________________________________________________________________
  1704. void MpdVectorFinder::Writedir2current( TObject *obj )
  1705. {
  1706. /// Write
  1707. if( !obj->IsFolder() ) obj->Write();
  1708. else{
  1709. TDirectory *cur = gDirectory;
  1710. TDirectory *sub = cur->mkdir(obj->GetName());
  1711. sub->cd();
  1712. TList *listSub = ((TDirectory*)obj)->GetList();
  1713. TIter it(listSub);
  1714. while( TObject *obj1=it() ) Writedir2current(obj1);
  1715. cur->cd();
  1716. }
  1717. }
  1718. //__________________________________________________________________________
  1719. void MpdVectorFinder::StoreTracks()
  1720. {
  1721. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1722. /// Transfer tracks from fTrackCand to fTracks
  1723. Int_t nFound = fTracks->GetEntriesFast();
  1724. ///cout << "fTrackCand entries " << fTrackCand->GetEntriesFast() << endl;
  1725. for (Int_t i = 0; i < fNTracks; ++i) {
  1726. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1727. track->Weight2Cov();
  1728. //AZ-211020 MpdVector *cellTr = (MpdVector*) fKHits[0]->UncheckedAt(track->GetUniqueID()-1);
  1729. MpdVector *cellTr = (MpdVector*) fVectorCands[track->GetUniqueID()-1].begin()->second;
  1730. /// added 16/9/2019, changing uniqueid
  1731. //AZ track->SetUniqueID(i + 1); /// uniqueid cannot be 0
  1732. track->SetUniqueID(nFound + 1); //AZ - 17.04.2020
  1733. new ((*fTracks)[nFound++]) MpdItsKalmanTrack(*track);
  1734. cout << track->GetNofTrHits() << " " << track->GetNofHits() << ", ";
  1735. TString code = cellTr->GetCode();
  1736. fCellMap[code] = -1; // used hit combination
  1737. //fTrackCand->RemoveAt(i);
  1738. }
  1739. fNTracks = 0;
  1740. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1741. }
  1742. //__________________________________________________________________________
  1743. void MpdVectorFinder::AddHits()
  1744. {
  1745. /// Add hit objects to tracks and compute number of wrongly assigned hits
  1746. /// (hits with ID different from ID of starting TPC track)
  1747. Int_t wrongCount = 0;
  1748. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1749. Int_t nFound = fTracks->GetEntriesFast();//work 07
  1750. //Int_t nFound = fTrackCand->GetEntriesFast();
  1751. fout_v << "nFound " << nFound << endl;
  1752. for (Int_t i = 0; i < nFound; ++i) {
  1753. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);
  1754. // MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1755. ///fout_v << track->GetNode() << " " << track->GetNodeNew() << endl;
  1756. Double_t c2 = track->GetChi2();
  1757. track->SetChi2(track->GetChi2Its());
  1758. track->SetChi2Its(c2);
  1759. Int_t nHits = track->GetNofHits();
  1760. ///if(nHits < 4) continue; // candidate not write track /// was 7
  1761. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1762. TClonesArray &trHits = *track->GetTrHits();
  1763. trHits.Delete(); // AZ GetNofTrHits() will be 0 after this
  1764. TObjArray *hits = track->GetHits();
  1765. SetTrackID(track); // set track ID as ID of majority of its hits
  1766. ///fout_v << nHits << " " << trHits.GetEntriesFast() << " " << track->GetTrackID() << " " << track->GetNofTrHits() << endl;
  1767. Int_t nWrong = 0, motherID = track->GetTrackID(); //nMirr = 0,
  1768. // Get track mother ID
  1769. ///cout << "motherID " << motherID << endl;
  1770. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  1771. while (mctrack->GetMotherId() >= 0) {
  1772. motherID = mctrack->GetMotherId();
  1773. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1774. }
  1775. //Int_t lastIndx = trHits.GetEntriesFast();
  1776. for (Int_t j = 0; j < nHits; ++j) {
  1777. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1778. //AZ hit->SetUniqueID(1); // flag ITS hits
  1779. new (trHits[j]) MpdKalmanHit(*hit);
  1780. /// fout_v << " " << hit->GetLayer();
  1781. Int_t nOver = hit->Index()->GetSize();
  1782. for (Int_t iov = 0; iov < nOver; ++iov) {
  1783. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex(iov));
  1784. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1785. /// fout_v << "-" << motherID1;
  1786. // Get point mother ID
  1787. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID1);
  1788. while (mctrack->GetMotherId() >= 0) {
  1789. motherID1 = mctrack->GetMotherId();
  1790. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1791. }
  1792. if (motherID1 != motherID) ++nWrong;
  1793. }
  1794. }
  1795. /// if (nHits) cout << "\n Wrongs: " << nWrong << endl;
  1796. if (nWrong > 0) {
  1797. wrongCount++;
  1798. /// cout << "!!!" << endl;
  1799. }
  1800. track->SetNofWrong(nWrong);
  1801. track->SetNofHits(track->GetNofTrHits()); // TPC and ITS hits
  1802. track->SetNofIts(nHits);
  1803. /// fout_v << nHits << " " << track->GetNofTrHits() << " " << track->GetTrackID() << " "
  1804. /// << track->GetChi2Its() << " " << track->GetChi2() << endl;
  1805. }
  1806. /// fout_v << "wrongCount: " << wrongCount << endl;
  1807. fTracks->Compress();
  1808. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1809. }
  1810. //__________________________________________________________________________
  1811. void MpdVectorFinder::SetTrackID(MpdItsKalmanTrack* track)
  1812. {
  1813. /// Set track ID as ID of majority of its hits
  1814. const Int_t idMax = 9999;
  1815. vector<Int_t> ids(idMax);
  1816. Int_t nHits = track->GetNofHits(), locMax = 0, size = idMax, id0 = -1;
  1817. TObjArray *hits = track->GetHits();
  1818. for (Int_t i = 0; i < nHits; ++i) {
  1819. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  1820. Int_t id = GetHitID(hit);
  1821. if (i == 0) id0 = id; // first hit
  1822. if (id >= size) { size = id + 1; ids.resize(size); }
  1823. ++ids[id];
  1824. if (ids[id] > ids[locMax]) locMax = id;
  1825. }
  1826. if (ids[locMax] > 1) track->SetTrackID(locMax);
  1827. else track->SetTrackID(id0);
  1828. }
  1829. //__________________________________________________________________________
  1830. Int_t MpdVectorFinder::GetHitID(MpdKalmanHit *hit)
  1831. {
  1832. /// Get hit ID from MCPoint ID
  1833. Int_t nOver = hit->Index()->GetSize();
  1834. for (Int_t iov = 0; iov < nOver; ++iov) {
  1835. MpdItsHit5spd *h = (MpdItsHit5spd*) fItsHits->UncheckedAt(hit->GetIndex(iov));
  1836. Int_t motherID1 = ((FairMCPoint*) fItsPoints->UncheckedAt(h->GetRefIndex()))->GetTrackID();
  1837. return motherID1;
  1838. }
  1839. return 0; // 0 points in hit?
  1840. }
  1841. //__________________________________________________________________________
  1842. void MpdVectorFinder::ExcludeHits()
  1843. {
  1844. /// Exclude hits, already used for tracking, from consideration during the next passes
  1845. fout_v << "- ExcludeHits -" << endl;
  1846. Int_t nReco = fTracks->GetEntriesFast();// work 07
  1847. //Int_t nReco = fTrackCand->GetEntriesFast();
  1848. fout_v << " nReco: " << nReco << endl;
  1849. for (Int_t i = 0; i < nReco; ++i) {
  1850. MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTracks->UncheckedAt(i);// work 07
  1851. //MpdItsKalmanTrack *track = (MpdItsKalmanTrack*) fTrackCand->UncheckedAt(i);
  1852. Int_t nhitsKF = track->GetNofHits();
  1853. TObjArray *hits = track->GetHits();
  1854. for (Int_t j = 0; j < nhitsKF; ++j) {
  1855. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1856. MpdItsHit5spd *h = (MpdItsHit5spd*)fItsHits->UncheckedAt(hit->GetIndex());
  1857. hit->SetFlag(-1);// work 15.04
  1858. h->SetFlag(-1); // work 15.04
  1859. }
  1860. }
  1861. }
  1862. //__________________________________________________________________________
  1863. Double_t MpdVectorFinder::Interp(Double_t angt, Int_t choice, Int_t lay)
  1864. {
  1865. // Perform parabolic interpolation for Angl (choice = 0) and dAngt (choice = 1)
  1866. const Int_t np = 7;
  1867. 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};
  1868. // const Double_t xp[np] = {0.05,0.125,0.2, 0.25, 0.50, 1.00, 2.00}; //3.09
  1869. //Double_t dpp05[np] = {0.18, 0.069, 0.043, 0.034, 0.017, 0.0086, 0.0044}; //mean angt
  1870. //Double_t err05[np] = {0.0018, 0.00032, 0.00018, 0.00014, 0.000092, 0.00007, 0.00006}; // error mean angt
  1871. //const Double_t angs[2][np] = {{0.042, 0.019, 0.017, 0.016, 0.017, 0.019, 0.018}, //sigma angl
  1872. // {0.054, 0.013, 0.0069, 0.006, 0.0036, 0.0034, 0.003}}; // sigma dangt Box(250)
  1873. //const Double_t dangs[2][np] = {{0.0015, 0.00067, 0.00069, 0.00065, 0.00073, 0.00065, 0.00068}, // error sigma angl
  1874. // {0.0027, 0.00045, 0.00027, 0.00018, 0.00014, 0.0001, 0.0001}}; //error sigma dangt
  1875. /*
  1876. const Double_t angls[3][np] = {{0.047, 0.027, 0.025, 0.023, 0.027, 0.026, 0.028},
  1877. {0.04, 0.022, 0.021, 0.02, 0.021, 0.021, 0.021},
  1878. {0.041, 0.014, 0.013, 0.012, 0.012, 0.012, 0.012}}; // sigma deltaz 2,3,4 layers
  1879. const Double_t angts[2][np] = {{0.062, 0.013, 0.0075, 0.0058, 0.0041, 0.0039, 0.0036},
  1880. {0.059, 0.013, 0.0074, 0.006, 0.0035, 0.0029, 0.0026}}; //sigma dangt 3,4 layers
  1881. */
  1882. const Double_t angls[3][np] = {{0.045, 0.026, 0.025, 0.023, 0.027, 0.027, 0.028},
  1883. {0.02, 0.011, 0.01, 0.01, 0.011, 0.011, 0.01},
  1884. {0.016, 0.0045, 0.0037, 0.0036, 0.0033, 0.0033, 0.0031}}; // sigma deltaz 2,3,4 layers
  1885. const Double_t angts[2][np] = {{0.031, 0.0068, 0.0037, 0.0029, 0.0021, 0.0019, 0.0018},
  1886. {0.019, 0.0039, 0.0021, 0.0016, 0.00091, 0.00059, 0.0005}}; //sigma dangt 3,4 layers
  1887. Double_t x0 = 0.1 / angt;
  1888. //Double_t x0 = angt;
  1889. Int_t ok = 0, inds[3] = {np - 3, np - 2, np - 1};
  1890. for (Int_t i = 0; i < np; ++i) {
  1891. if (xp[i] < x0) continue;
  1892. inds[1] = i;
  1893. inds[0] = i - 1;
  1894. inds[2] = i + 1;
  1895. ok = 1;
  1896. break;
  1897. }
  1898. if (ok) {
  1899. if (inds[0] < 0) for (Int_t i = 0; i < 3; ++i) ++inds[i];
  1900. if (inds[2] == np) for (Int_t i = 0; i < 3; ++i) --inds[i];
  1901. }
  1902. Double_t a{0.},b{0.},c{0.},dy01{0.},dy02{0.};
  1903. // Parabolic interpolation
  1904. Double_t dx01 = xp[inds[1]] - xp[inds[0]], dx02 = xp[inds[2]] - xp[inds[0]];
  1905. Double_t x01 = xp[inds[1]] + xp[inds[0]], x02 = xp[inds[2]] + xp[inds[0]];
  1906. // Double_t dy01 = angs[choice][inds[1]] - angs[choice][inds[0]];
  1907. // Double_t dy02 = angs[choice][inds[2]] - angs[choice][inds[0]];
  1908. if (choice == 0) {
  1909. dy01 = angls[lay][inds[1]] - angls[lay][inds[0]];
  1910. dy02 = angls[lay][inds[2]] - angls[lay][inds[0]];
  1911. }
  1912. if (choice == 1) {
  1913. dy01 = angts[lay][inds[1]] - angts[lay][inds[0]];
  1914. dy02 = angts[lay][inds[2]] - angts[lay][inds[0]];
  1915. }
  1916. Double_t slope = dy01 / dx01;
  1917. a = dy02 - slope * dx02;
  1918. a /= (dx02 * x02 - dx02 * x01);
  1919. b = slope - a * x01;
  1920. if (choice == 0) {
  1921. // c = angs[choice][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1922. c = angls[lay][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1923. }
  1924. if (choice == 1) {
  1925. // c = angs[choice][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1926. c = angts[lay][inds[0]] - a * xp[inds[0]] * xp[inds[0]] - b * xp[inds[0]];
  1927. }
  1928. //for (Int_t i = 0; i < 3; ++i) cout << xp[inds[i]] << " " << angs[choice][inds[i]] << " ";
  1929. //cout << endl << x0 << " " << a * x0 * x0 + b * x0 + c << endl;
  1930. return a * x0 * x0 + b * x0 + c;
  1931. }
  1932. //__________________________________________________________________________
  1933. void MpdVectorFinder::GetShortTracks()
  1934. {
  1935. /// Collect remaining tracks with small number of hits.
  1936. map<TString,Int_t>::iterator it;
  1937. const Int_t nLays = 5;
  1938. Int_t inds[nLays];
  1939. for (it = fCellMap.begin(); it != fCellMap.end(); ++it) {
  1940. //cout << it->first << " " << it->second <<endl;
  1941. if (it->second < 0) continue;
  1942. TString code = it->first, cc;
  1943. Int_t ibeg = 0, leng = -1;
  1944. //cout << code << endl;
  1945. for (Int_t i = 0; i < nLays; ++i) {
  1946. leng = code.Index("-",leng);
  1947. cc = code(ibeg,leng-ibeg);
  1948. inds[i] = cc.Atoi();
  1949. ibeg = leng + 1;
  1950. ++leng;
  1951. }
  1952. }
  1953. //Int_t left = 0;
  1954. /// fout_v << " Leftovers: " << endl;
  1955. for (it = fCellMap.begin(); it != fCellMap.end(); ++it) {
  1956. ///if (it->second > 0) fout_v << ++left << " " << it->first << endl;
  1957. }
  1958. }
  1959. //__________________________________________________________________________
  1960. Bool_t MpdVectorFinder::Refit(MpdItsKalmanTrack *track, Double_t mass, Int_t charge)
  1961. {
  1962. /// Refit track in ITS+TPC using track hits (toward beam line) for some
  1963. /// particle mass and charge hypothesis
  1964. if (fTpcKF == nullptr) fTpcKF = (MpdTpcKalmanFilter*) FairRun::Instance()->GetTask("TPC Kalman filter");
  1965. if (fTpcKF == nullptr) {
  1966. cout << " !!! TPC Kalman Filter was not activated !!! " << endl;
  1967. exit(1);
  1968. }
  1969. MpdKalmanGeoScheme *geoScheme = MpdKalmanFilter::Instance()->GetGeo();
  1970. MpdKalmanHit hTmp;
  1971. if (track->GetNofTrHits() > track->GetNofIts()) {
  1972. // Refit inside TPC
  1973. //MpdTpcKalmanTrack *tpc = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(track->GetUniqueID()-1);
  1974. track->SetUniqueID(0); // reset ITS flag
  1975. //tr.GetWeightAtHit()->Print();
  1976. //tr.GetParamAtHit()->Print();
  1977. fTpcKF->Refit(track, mass, charge); // this mixes up hit order (ITS and TPC) !!!
  1978. //tr.GetWeight()->Print();
  1979. //tr.GetParam()->Print();
  1980. if (track->GetNofIts() == 0) return kTRUE;
  1981. hTmp.SetType(MpdKalmanHit::kFixedR);
  1982. hTmp.SetPos(track->GetPos());
  1983. track->SetParamNew(*track->GetParam());
  1984. track->SetPos(track->GetPosNew());
  1985. track->ReSetWeight();
  1986. //TMatrixDSym w = *tr.GetWeight(); // save current weight matrix
  1987. //MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kFALSE);
  1988. MpdKalmanFilter::Instance()->PropagateParamToHit(track,&hTmp,kFALSE);
  1989. //tr.SetWeight(w); // restore original weight matrix (near TPC inner shell)
  1990. track->SetDirection(MpdKalmanTrack::kInward);
  1991. track->SetUniqueID(1);
  1992. } else {
  1993. track->SetDirection(MpdKalmanTrack::kOutward);
  1994. hTmp.SetType(MpdKalmanHit::kFixedR);
  1995. hTmp.SetPos(track->GetPos());
  1996. track->SetPos(track->GetPosNew());
  1997. track->SetParamNew(*track->GetParam());
  1998. Bool_t ok = MpdKalmanFilter::Instance()->PropagateParamToHit(track,&hTmp,kFALSE);
  1999. track->ReSetWeight();
  2000. ok = MpdKalmanFilter::Instance()->Refit(track, -1, 0);
  2001. if (!ok) return ok;
  2002. hTmp.SetPos(((MpdKalmanHit*)track->GetTrHits()->First())->GetDist()+1.0);
  2003. MpdKalmanFilter::Instance()->PropagateParamToHit(track,&hTmp,kFALSE);
  2004. track->SetDirection(MpdKalmanTrack::kInward);
  2005. track->SetPos(track->GetPosNew());
  2006. track->SetParam(*track->GetParamNew());
  2007. Int_t nHits = track->GetNofHits();
  2008. nHits *= nHits;
  2009. TMatrixDSym *www = track->GetWeight();
  2010. for (Int_t i = 0; i < 5; ++i) {
  2011. for (Int_t j = i; j < 5; ++j) {
  2012. if (j == i) (*www)(i,j) /= nHits;
  2013. else (*www)(i,j) = (*www)(j,i) = 0.;
  2014. }
  2015. }
  2016. }
  2017. Int_t ntot = track->GetNofTrHits(), nhits = TMath::Min(15,ntot);
  2018. Int_t indx0 = ntot - nhits;
  2019. TMatrixD param(5,1);
  2020. TMatrixDSym weight(5), pointWeight(5);
  2021. TString mass2 = "";
  2022. mass2 += mass * mass;
  2023. Double_t vert[3] = {0.0,0.0,0.0};
  2024. const Double_t thick[9] = {0.005, 0.005, 0.005, 0.07, 0.07};
  2025. for (Int_t ih = indx0; ih < ntot; ++ih) {
  2026. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetTrHits()->UncheckedAt(ih);
  2027. if (hit->GetUniqueID() == 0) continue; // skip TPC hits
  2028. MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kFALSE, kTRUE); // do not adjust track length for the moment
  2029. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hit, pointWeight, param);
  2030. //cout << ih << " " << hit->GetPos() << " " << hit->GetUniqueID() << " " << dChi2 << "\n";
  2031. track->SetChi2(track->GetChi2()+dChi2);
  2032. weight = *track->GetWeight();
  2033. weight += pointWeight;
  2034. track->SetWeight(weight);
  2035. track->SetParamNew(param);
  2036. // Add multiple scattering in the sensor
  2037. //*
  2038. Double_t x0 = 9.36; //, step = 0.005 * 4.0; // rad. length
  2039. TVector3 mom3 = track->Momentum3();
  2040. mom3.SetMag(1.0);
  2041. TVector3 norm = geoScheme->Normal(hit);
  2042. Double_t cosPhi = norm.Dot(mom3);
  2043. TMatrixDSym *cov = track->Weight2Cov();
  2044. Double_t th = track->GetParamNew(3);
  2045. Double_t cosTh = TMath::Cos(th);
  2046. //Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, thick[lay]/cosPhistep, mass2);
  2047. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, 4*thick[hit->GetLayer()]/cosPhi, mass2);
  2048. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  2049. (*cov)(3,3) += angle2;
  2050. Int_t iok = 0;
  2051. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  2052. track->SetWeight(*cov);
  2053. //*/
  2054. }
  2055. //------ Propagate track to the beam line
  2056. track->SetParam(*track->GetParamNew());
  2057. track->SetPos(track->GetPosNew());
  2058. Double_t pos = track->GetPos();
  2059. TMatrixD par = *track->GetParam();
  2060. TMatrixDSym cov = *track->Weight2Cov();
  2061. Double_t leng = track->GetLength();
  2062. TString nodeNew = track->GetNodeNew();
  2063. //cout << " 1: " << nodeNew << ", " << track->GetNode() << " " << track->GetHits()->GetEntriesFast() << endl;
  2064. // Go to beam pipe
  2065. hTmp.SetPos(fPipeR);
  2066. //hTmp.SetPos(2.9); // !!! FIXME
  2067. Int_t iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hTmp, kFALSE);
  2068. if (iok != 1) {
  2069. // Restore track
  2070. track->SetParam(par);
  2071. track->SetParamNew(par);
  2072. track->SetCovariance(cov);
  2073. track->ReSetWeight();
  2074. track->SetPos(pos);
  2075. track->SetPosNew(pos);
  2076. track->SetLength(leng);
  2077. //track->SetNode(node);
  2078. track->SetNodeNew(nodeNew);
  2079. } else {
  2080. // Add multiple scattering
  2081. //Double_t dX = 0.05 / 8.9; // 0.5 mm of Al
  2082. Double_t dX = 0.1 / 35.28; // 1. mm of Be
  2083. TMatrixDSym* pcov = track->Weight2Cov();
  2084. Double_t th = track->GetParamNew(3);
  2085. Double_t cosTh = TMath::Cos(th);
  2086. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dX);
  2087. (*pcov)(2,2) += (angle2 / cosTh / cosTh);
  2088. (*pcov)(3,3) += angle2;
  2089. Int_t ok = 0;
  2090. MpdKalmanFilter::Instance()->MnvertLocal(pcov->GetMatrixArray(), 5, 5, 5, ok);
  2091. track->SetWeight(*pcov);
  2092. }
  2093. track->Weight2Cov(); // rebuild covar matrix (it was overwritten by weight above)
  2094. hTmp.SetPos(0.);
  2095. hTmp.SetMeas(0,track->GetParam(2)); // track Phi
  2096. iok = MpdKalmanFilter::Instance()->PropagateToHit(track, &hTmp, kFALSE);
  2097. if (iok != 1) {
  2098. track->SetNodeNew(""); // for FindPca
  2099. MpdKalmanFilter::Instance()->FindPca(track, vert);
  2100. }
  2101. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  2102. return kTRUE;
  2103. }
  2104. //__________________________________________________________________________
  2105. ClassImp(MpdVectorFinder);