MpdEctTrackFinderTof.cxx 87 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256
  1. // -------------------------------------------------------------------------
  2. // ----- MpdEctTrackFinderTof source file -----
  3. // ----- Created 28/07/08 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdEctTrackFinderTof.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Track finder in MPD End-Cap Tracker (ECT) using seeds built from ECT,
  9. ** ETOF and vertex points
  10. **/
  11. #include "MpdEctTrackFinderTof.h"
  12. #include "MpdKalmanFilter.h"
  13. #include "MpdKalmanGeoScheme.h"
  14. #include "MpdKalmanHit.h"
  15. #include "MpdEctKalmanTrack.h"
  16. #include "MpdTpcKalmanTrack.h"
  17. #include "MpdStrawendcapGeoPar.h"
  18. #include "MpdStrawendcapPoint.h"
  19. #include "MpdTofHit.h"
  20. #include "MpdTofPoint.h"
  21. #include "MpdTofUtils.h"
  22. #include "MpdVertex.h"
  23. #include "MpdKfPrimaryVertexFinder.h"
  24. #include "FairMCEventHeader.h"
  25. #include "FairGeoNode.h"
  26. #include "MpdMCTrack.h"
  27. #include "FairRootManager.h"
  28. #include "FairRunAna.h"
  29. #include "FairRuntimeDb.h"
  30. #include "TMath.h"
  31. #include "TFile.h"
  32. //#include "TLorentzVector.h"
  33. #include "TVector2.h"
  34. //#include "TClonesArray.h"
  35. #include <TRandom.h>
  36. #include <iostream>
  37. #include <set>
  38. using std::cout;
  39. using std::endl;
  40. //using std::vector;
  41. const Double_t MpdEctTrackFinderTof::fgkChi2Cut = 10; //10; //20; //100;
  42. FILE *lun1 = 0x0; //fopen("error1.dat","w");
  43. FILE *lun2 = 0x0; //fopen("match.dat","w");
  44. //__________________________________________________________________________
  45. MpdEctTrackFinderTof::MpdEctTrackFinderTof(const char *name, Int_t iVerbose )
  46. :FairTask(name, iVerbose),
  47. fNPass(1),
  48. fTpc(kFALSE),
  49. fHistoDir(0x0),
  50. fExact(0)
  51. {
  52. //fKHits = new TClonesArray("MpdKalmanHit", 100);
  53. fTrackCand = new TClonesArray("MpdEctKalmanTrack", 100);
  54. fhLays = new TH1F("hLays1","ECT layers",fgkNlays+1,0,fgkNlays+1);
  55. fLayPointers = 0x0;
  56. }
  57. //__________________________________________________________________________
  58. MpdEctTrackFinderTof::~MpdEctTrackFinderTof()
  59. {
  60. //delete fKHits;
  61. //delete fTracks;
  62. //delete fTrackCand;
  63. delete [] fLayPointers;
  64. delete fhLays;
  65. }
  66. //__________________________________________________________________________
  67. InitStatus MpdEctTrackFinderTof::Init()
  68. {
  69. ReInit();
  70. InitGeo();
  71. return kSUCCESS;
  72. }
  73. //__________________________________________________________________________
  74. InitStatus MpdEctTrackFinderTof::ReInit()
  75. {
  76. if (fTpc) {
  77. fKHits =(TClonesArray *) FairRootManager::Instance()->GetObject("EctKalmanHit");
  78. fTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("EctTrack");
  79. //fTpcTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  80. } else {
  81. fKHits = new TClonesArray("MpdKalmanHit", 100);
  82. fTracks = new TClonesArray("MpdEctKalmanTrack", 100);
  83. FairRootManager::Instance()->Register("EctTrack", "Ect", fTracks, kTRUE);
  84. }
  85. fTpcTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("TpcKalmanTrack");
  86. fEctHits =(TClonesArray *) FairRootManager::Instance()->GetObject("STRAWPoint");
  87. fTofHits =(TClonesArray *) FairRootManager::Instance()->GetObject("ETOFHit");
  88. fTofPoints =(TClonesArray *) FairRootManager::Instance()->GetObject("ETOFPoint");
  89. fMCTracks =(TClonesArray *) FairRootManager::Instance()->GetObject("MCTrack");
  90. //fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("STSTrackMatch");
  91. fPrimVtx = (TClonesArray *) FairRootManager::Instance() ->GetObject("Vertex");
  92. fNPass = 1;
  93. return kSUCCESS;
  94. }
  95. //__________________________________________________________________________
  96. void MpdEctTrackFinderTof::Reset()
  97. {
  98. ///
  99. cout << " MpdEctTrackFinderTof::Reset " << endl;
  100. //fKHits->Clear();
  101. //fTracks->Clear("C");
  102. if (!fTpc) {
  103. fKHits->Clear();
  104. fTracks->Delete();
  105. }
  106. fTrackCand->Delete();
  107. delete [] fLayPointers;
  108. }
  109. //__________________________________________________________________________
  110. void MpdEctTrackFinderTof::SetParContainers()
  111. {
  112. // Get run and runtime database
  113. FairRunAna* run = FairRunAna::Instance();
  114. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  115. FairRuntimeDb* db = run->GetRuntimeDb();
  116. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  117. // Get ECT geometry parameter container
  118. db->getContainer("MpdStrawendcapGeoPar");
  119. }
  120. //__________________________________________________________________________
  121. void MpdEctTrackFinderTof::Finish()
  122. {
  123. //Write();
  124. }
  125. //__________________________________________________________________________
  126. void MpdEctTrackFinderTof::Exec(Option_t * option)
  127. {
  128. static int eventCounter = 0;
  129. cout << " EctRec event " << ++eventCounter << endl;
  130. Reset();
  131. // Create Kalman hits
  132. MakeKalmanHits();
  133. if (fKHits->GetEntriesFast() == 0) return;
  134. // Match with ETOF hits for tracks extrapolated from TPC
  135. if (fTpc) MatchEtof();
  136. for (Int_t i = 0; i < fNPass; ++i) {
  137. fTrackCand->Delete();
  138. GetTrackSeeds(i);
  139. Int_t nHitsOk = 0, nHits = fKHits->GetEntriesFast();
  140. for (Int_t j = 0; j < nHits; ++j) {
  141. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(j);
  142. if (hit->GetFlag() >= 0) ++nHitsOk;
  143. }
  144. cout << " Total number of hits for tracking: " << nHits << ", good: " << nHitsOk << endl;
  145. cout << " Total number of track seeds: " << fTrackCand->GetEntriesFast() << endl;
  146. DoTracking(i);
  147. RemoveDoubles(); // remove double tracks
  148. StoreTracks(i);
  149. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  150. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  151. }
  152. SelectTracks(0); // do track selection and compute shared hit multiplicities
  153. AddHits(); // add hit objects to tracks
  154. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  155. if (fTracks->GetEntriesFast() < 1) return;
  156. if (fTpcTracks) MatchTpc();
  157. GoToBeamLine();
  158. Smooth(); // apply constraints due to primary vertex
  159. GoOutward(); // for matching with outer detectors
  160. }
  161. //__________________________________________________________________________
  162. void MpdEctTrackFinderTof::InitGeo()
  163. {
  164. // Initialize detector geometry
  165. // Create configuration
  166. cout << " Creating configuration ... " << endl;
  167. map<TString, FairGeoNode*> layers;
  168. map<TString, FairGeoNode*> tubes;
  169. map<TString, FairGeoNode*> modules;
  170. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  171. // Create mother volume
  172. FairGeoNode *cave = new FairGeoNode();
  173. cave->setVolumeType(kFairGeoTopNode);
  174. MpdStrawendcapGeoPar *geoPar = (MpdStrawendcapGeoPar*) rtdb->getContainer("MpdStrawendcapGeoPar");
  175. TObjArray* passNodes = geoPar->GetGeoPassiveNodes();
  176. Int_t nVol = passNodes->GetEntriesFast();
  177. // Get volumes
  178. Int_t nTripl = 0, nMod = 0, nTube = 0;
  179. frMinMax[0] = frMinMax[2] = -1.;
  180. for (Int_t i = 0; i < nVol; ++i) {
  181. FairGeoNode* passVol = (FairGeoNode*) (passNodes->UncheckedAt(i));
  182. TString name = passVol->GetName();
  183. TArrayD* params = passVol->getParameters();
  184. //cout << i << " " << name << " " << (*params)[0] << " " << (*params)[1] << " "
  185. // << (*params)[2] << " " << passVol->GetListOfDaughters()->GetEntriesFast() << endl;
  186. if (name.Contains("layer")) {
  187. layers.insert(pair<TString, FairGeoNode*>(name,passVol));
  188. if (frMinMax[0] < 0) {
  189. frMinMax[0] = (*params)[0];
  190. frMinMax[1] = (*params)[1];
  191. }
  192. Ssiz_t pos = name.First("#") + 1;
  193. nTripl = TMath::Max (nTripl, (TString(name(pos,name.Length()))).Atoi());
  194. } else if (name.Contains("tube")) {
  195. tubes.insert(pair<TString, FairGeoNode*>(name,passVol));
  196. if (frMinMax[2] < 0) frMinMax[2] = (*params)[1] * 10.; // mm
  197. Ssiz_t pos = name.First("#") + 1;
  198. nTube = TMath::Max (nTube, (TString(name(pos,name.Length()))).Atoi());
  199. } else if (name.Contains("mod")) {
  200. modules.insert(pair<TString, FairGeoNode*>(name,passVol));
  201. Ssiz_t pos = name.First("#") + 1;
  202. nMod = TMath::Max (nMod, (TString(name(pos,name.Length()))).Atoi());
  203. }
  204. /*
  205. CbmGeoTransform *transf = passVol->getLabTransform();
  206. CbmGeoVector trans = transf->getTranslation();
  207. //cout << trans << endl;
  208. CbmGeoRotation rot = transf->getRotMatrix();
  209. //rot.print();
  210. passVol->setMother(cave);
  211. passVol->calcLabTransform();
  212. */
  213. }
  214. //fNofLays = nTripl * 3;
  215. //fNofLays *= 2;
  216. fNofTubes = nTube;
  217. // Set mother volumes and find out views ordering
  218. map<TString,FairGeoNode*>::iterator it;
  219. Int_t ishft[3] = {0};
  220. Double_t zzz[3] = {0.};
  221. for (it = layers.begin(); it != layers.end(); ++it) {
  222. FairGeoNode *lay = (*it).second;
  223. lay->setMother(cave);
  224. TString name = lay->GetName();
  225. if (name == "stt01layerradial#1") zzz[0] = lay->getLabTransform()->getTranslation().getZ();
  226. else if (name == "stt01layerleft#1") zzz[1] = lay->getLabTransform()->getTranslation().getZ();
  227. else if (name == "stt01layerright#1") zzz[2] = lay->getLabTransform()->getTranslation().getZ();
  228. //cout << lay->getLabTransform()->getTranslation() << endl;
  229. //lay->getLabTransform()->getRotMatrix().print();
  230. }
  231. TMath::Sort(3,zzz,ishft,kFALSE);
  232. // Get transforms
  233. for (it = tubes.begin(); it != tubes.end(); ++it) {
  234. FairGeoNode *tube = (*it).second;
  235. TString layName = tube->GetName();
  236. //if (layName.Contains("stt03") || layName.Contains("stt04")) continue; // skip z<0 for the moment !!!
  237. layName.ReplaceAll("tube","layer");
  238. Ssiz_t pos = layName.First("#") + 1;
  239. Int_t group, itube, ind = 0;
  240. sscanf(&layName[4],"%d",&group);
  241. --group;
  242. itube = (TString(layName(pos,layName.Length()))).Atoi() - 1;
  243. //cout << lay << " " << module << " " << itube << endl;
  244. if (itube >= fgkNtube) {
  245. cout << " Too many tubes !!! " << itube << endl;
  246. exit(0);
  247. }
  248. if (layName.Contains("left")) ind = 1;
  249. else if (layName.Contains("right")) ind = 2;
  250. // Loop over triplets
  251. for (Int_t i = 0; i < nTripl; ++i) {
  252. layName.Remove(pos);
  253. layName += (i+1);
  254. tube->setMother(layers[layName]);
  255. Int_t ilay = i * 3 + group * nTripl * 3 + ishft[ind];
  256. if (ilay >= fgkNlays) {
  257. cout << " Too many layers !!! " << ilay << endl;
  258. exit(0);
  259. }
  260. //cout << tube->GetName() << " " << lay << endl;
  261. FairGeoTransform *tran = tube->calcLabTransform();
  262. //cout << tran->getTranslation() << endl;
  263. //tran->getRotMatrix().print();
  264. fTransf[ilay][itube] = tran;
  265. //transf[ilay][itube] = layers[lay]->getLabTransform();
  266. if (itube == 0) fZplanes[ilay] = layers[layName]->getLabTransform()->getTranslation().getZ();
  267. }
  268. } // for (it = tubes.begin();
  269. }
  270. //__________________________________________________________________________
  271. void MpdEctTrackFinderTof::MakeKalmanHits()
  272. {
  273. /// Create Kalman hits from ECT points.
  274. if (!fTpc) {
  275. // Without TPC tracking
  276. Bool_t mirror = kTRUE; //kFALSE; //kTRUE; // add mirror hits if TRUE
  277. if (fExact) mirror = kFALSE;
  278. Int_t nHits = fEctHits->GetEntriesFast(), lay = 0, nKH = 0, itube = 0;
  279. //Double_t phi, r, coord, errR = 0.2, errRphi = 0.02; // 2mm in R, 200um in R-Phi
  280. Double_t phi, coord, errR = 0.02, errRphi = 0.02; //r, 200um in R, 200um in R-Phi
  281. /*Double_t firstHit[100][1000] = {{0},{0}};
  282. Int_t indxTube[100][1000];
  283. for (Int_t i = 0; i < 100; ++i) {
  284. for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
  285. }*/
  286. Double_t firstHit[fgkNlays+1][fgkNtube] = {{0},{0}};
  287. Int_t indxTube[fgkNlays+1][fgkNtube];
  288. Int_t iend = fgkNlays + 1;
  289. for (Int_t i = 0; i < iend; ++i) {
  290. for (Int_t j = 0; j < fgkNtube; ++j) indxTube[i][j] = -1;
  291. }
  292. for (Int_t ih = 0; ih < nHits; ++ih) {
  293. MpdStrawendcapPoint *h = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(ih);
  294. //if (h->GetZ() < 0) continue; // !!! exclude one side for now
  295. phi = h->GetPhi(); // tube Phi
  296. lay = h->GetDetectorID() / 1000;
  297. if (h->GetZ() < 0) lay += fgkNlays2;;
  298. itube = h->GetDetectorID() % 1000; // tube number
  299. // Extrapolate track to Z = Ztube
  300. Double_t dZ = h->GetZ() - h->GetTrackZ();
  301. Double_t dt = 0.; // dZ / h->GetPz();
  302. if (TMath::Abs(h->GetPz()) > 1.e-6 && h->GetPz() * h->GetZ() > 0) dt = dZ / h->GetPz();
  303. Double_t xNew = h->GetTrackX() + dt * h->GetPx();
  304. Double_t yNew = h->GetTrackY() + dt * h->GetPy();
  305. //cout << dZ << " " << h->GetTrackX() << " " << xNew << " " << h->GetTrackY() << " " << yNew << " " << lay << endl;
  306. //Double_t zNew = h->GetTrackZ() + dt * h->GetPz(); // just for cross-check
  307. // Transform to the tube local coordinate system
  308. Double_t cosSin[2] = {TMath::Cos(phi), TMath::Sin(phi)}; // phi - tube direction
  309. //Double_t xLoc = h->GetX() * cosPhi + h->GetY() * sinPhi; // cross-check
  310. //Double_t yLoc = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  311. //Double_t xLoc = xNew * cosSin[0] + yNew * cosSin[1];
  312. Double_t yLoc = -xNew * cosSin[1] + yNew * cosSin[0];
  313. //Double_t xLoc = (xNew - h->GetX()) * cosPhi + (yNew - h->GetY()) * sinPhi;
  314. //Double_t yLoc = -(xNew - h->GetX()) * sinPhi + (yNew - h->GetY()) * cosPhi;
  315. //r = xNew * xNew + yNew * yNew;
  316. //r = TMath::Sqrt (r);
  317. //r = TMath::Abs(xLoc);
  318. //r = xLoc;
  319. //cout << xLoc << " " << yLoc << " " << r << " " << h->GetPz() << endl;
  320. coord = yLoc;
  321. Double_t yWire = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  322. Double_t dY = TMath::Abs (yWire - coord);
  323. // Add error
  324. Double_t dRphi = 0, dR = 0;
  325. gRandom->Rannor(dRphi,dR); // add errors
  326. Double_t meas[2] = {coord+dRphi*errRphi, 0};
  327. Double_t err[2] = {errRphi, errR};
  328. MpdKalmanHit *hit = new ((*fKHits)[nKH]) MpdKalmanHit(lay*1000000, 1, MpdKalmanHit::kFixedZ,
  329. //MpdKalmanHit *hit = new ((*fKHits)[nKH]) MpdKalmanHit(h->GetDetectorID()*1000, 1, MpdKalmanHit::kFixedZ,
  330. meas, err, cosSin, 9., h->GetZ(), ih);
  331. //hit->SetXY(h->GetX(), h->GetY()); // wire coordinate
  332. // Add second measurement - just for test at the moment
  333. //!!!
  334. //hit->SetNofDim(2);
  335. //!!!
  336. //lay = hit.GetLayer();
  337. //layMax = TMath::Max (lay, layMax);
  338. //fhLays->Fill(lay+0.1);
  339. // Check if multiple hits per tube
  340. if (1 && indxTube[lay][itube] >= 0) {
  341. // Multiple hits
  342. MpdKalmanHit *hitOld = (MpdKalmanHit*) fKHits->UncheckedAt(indxTube[lay][itube]); // previous hit
  343. if (dY < firstHit[lay][itube]) {
  344. // Closer to anode wire - remove previous hit
  345. firstHit[lay][itube] = dY;
  346. const TArrayI *inds = hitOld->Index();
  347. Int_t nOld = inds->GetSize();
  348. for (Int_t io = 0; io < nOld; ++io) hit->SetIndex((*inds)[io]);
  349. fKHits->RemoveAt(indxTube[lay][itube]);
  350. if (mirror) fKHits->RemoveAt(indxTube[lay][itube]+1);
  351. indxTube[lay][itube] = nKH;
  352. } else {
  353. hitOld->SetIndex(hit->GetIndex());
  354. fKHits->RemoveAt(nKH); // remove last hit
  355. }
  356. } else {
  357. firstHit[lay][itube] = dY;
  358. indxTube[lay][itube] = nKH;
  359. //fhLays->Fill(lay+0.1);
  360. }
  361. if (mirror && fKHits->UncheckedAt(nKH)) {
  362. // Add mirror hit
  363. // wire point transverse coordinate in rotated system:
  364. //yLoc = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  365. //Double_t dy = hit->GetRphi() - yRot;
  366. Double_t dy = hit->GetMeas(0) - yWire;
  367. MpdKalmanHit *hitM = new ((*fKHits)[++nKH]) MpdKalmanHit(*hit);
  368. //hitM->SetRphi(yRot-dy);
  369. hitM->SetMeas(0, yWire-dy);
  370. hitM->SetMirror();
  371. //fhLays->Fill(lay+0.1);
  372. }
  373. ++nKH;
  374. }
  375. fKHits->Compress();
  376. fKHits->Sort(); // in descending order in abs(Z)
  377. } // if (!fTpc)
  378. if (fKHits->GetEntriesFast() == 0) return;
  379. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  380. fLayPointers = new Int_t [layMax+1];
  381. fhLays->Reset();
  382. /*
  383. Int_t ipos = 0;
  384. for (Int_t i = 0; i <= layMax; ++i) {
  385. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  386. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  387. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  388. Int_t cont = TMath::Nint (fhLays->GetCellContent(i+1,0));
  389. if (cont == 0) {
  390. fLayPointers[i] = -1;
  391. continue;
  392. }
  393. fLayPointers[i] = ipos;
  394. ipos += cont;
  395. }
  396. */
  397. for (Int_t i = 0; i <= layMax; ++i) fLayPointers[i] = -1;
  398. Int_t nKHits = fKHits->GetEntriesFast(), nOK = 0;
  399. for (Int_t i = 0; i < nKHits; ++i) {
  400. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(i);
  401. Int_t lay = hit->GetLayer();
  402. if (hit->GetFlag() > -1) { hit->SetFlag(1+hit->IsMirror()*1000); ++nOK; }
  403. fhLays->Fill(lay+0.1);
  404. if (fLayPointers[lay] < 0) fLayPointers[lay] = i;
  405. }
  406. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << " " << nOK << endl;
  407. //exit(0);
  408. }
  409. //__________________________________________________________________________
  410. void MpdEctTrackFinderTof::MatchEtof()
  411. {
  412. /// Match tracks, extended from TPC, with ETOF hits.
  413. static const Double_t rMax = 2.;
  414. MpdKalmanHit hTmp;
  415. hTmp.SetType(MpdKalmanHit::kFixedZ);
  416. hTmp.SetDist(((MpdTofHit*) fTofHits->First())->GetZ());
  417. Int_t nTof = fTofHits->GetEntriesFast();
  418. Int_t nTpc = fTracks->GetEntriesFast(), nMatch = 0;
  419. for (Int_t i = 0; i < nTpc; ++i) {
  420. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  421. MpdEctKalmanTrack track = *tr;
  422. //MpdKalmanFilter::Instance()->PropagateToHit(&track,&hTmp);
  423. track.SetParam(*track.GetParamNew());
  424. track.SetPos(track.GetPosNew());
  425. hTmp.SetDist(TMath::Abs(hTmp.GetDist())*TMath::Sign(1.,track.GetPos()));
  426. //MpdKalmanFilter::Instance()->PropagateParamZ(&track,&hTmp,kFALSE);
  427. MpdKalmanFilter::Instance()->PropagateParamZ(&track,&hTmp,kTRUE);
  428. Int_t iMin = 0; //id = track.GetTrackID(),
  429. Double_t rMin = 999999., xTr = track.GetParamNew(0), yTr = track.GetParamNew(1); //rID = 999,
  430. for (Int_t itof = 0; itof < nTof; ++itof) {
  431. MpdTofHit *tof = (MpdTofHit*) fTofHits->UncheckedAt(itof);
  432. //MpdTofPoint *tofP = (MpdTofPoint*) fTofPoints->UncheckedAt(tof->GetRefIndex());
  433. //if (tof->GetZ() < 0) continue; // !!! exclude one side for now
  434. Double_t dx = tof->GetX() - xTr, dy = tof->GetY() - yTr;
  435. Double_t r = dx * dx + dy * dy;
  436. //if (id == tofP->GetTrackID()) rID = r;
  437. if (r < rMin) { rMin = r; iMin = itof; }
  438. }
  439. //if (lun2) fprintf(lun2,"%d %10.3e %10.3e\n", id, TMath::Sqrt(rID), TMath::Sqrt(rMin));
  440. if (TMath::Sqrt(rMin) < rMax) {
  441. // Exclude ETOF hit
  442. //if (((MpdTofHit*) fTofHits->UncheckedAt(iMin))->GetFlag() != -1) {
  443. ++nMatch;
  444. MpdTofHit *tof = (MpdTofHit*) fTofHits->UncheckedAt(iMin);
  445. tof->SetFlag(tof->GetFlag() | MpdTofUtils::IsSelected);
  446. //3-may-12 tr->SetTofIndex(iMin);
  447. tr->SetMatchEtof();
  448. //}
  449. }
  450. // Modify track length
  451. tr->SetLength(track.GetLength());
  452. }
  453. cout << " *** Found " << nMatch << " matches with ETOF hits *** " << endl;
  454. }
  455. //__________________________________________________________________________
  456. void MpdEctTrackFinderTof::GetTrackSeeds(Int_t iPass)
  457. {
  458. /// Build track seeds from TOF, ECT hits and vertex
  459. Int_t tmp[99999] = {0}, idMax = 0, nID = 0;
  460. Int_t nTof = fTofHits->GetEntriesFast();
  461. Int_t nCand = 0, layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  462. const Double_t dphiMax = 0.1; //0.2; // cut on Phi between ECT and ETOF
  463. TVector3 vert(0.,0.,0.), pmom;
  464. cout << " Number of ETOF hits: " << nTof << endl;
  465. TVector3 posTof;
  466. for (Int_t itof = 0; itof < nTof; ++itof) {
  467. MpdTofHit *tof = (MpdTofHit*) fTofHits->UncheckedAt(itof);
  468. if (tof->GetFlag() & MpdTofUtils::IsSelected) continue;
  469. MpdTofPoint *tofP = (MpdTofPoint*) fTofPoints->UncheckedAt(tof->GetRefIndex());
  470. //if (tof->GetZ() < 0) continue; // !!! exclude one side for now
  471. tof->Position(posTof);
  472. Double_t rTof = posTof.Pt();
  473. Double_t drdz = rTof / tof->GetZ();
  474. Double_t phTof = posTof.Phi(), dph = 0, phEct = 0;
  475. tofP->Momentum(pmom);
  476. Int_t layBeg = layMax;
  477. if (tof->GetZ() > 0 && layMax > fgkNlays2) layBeg = layMax - fgkNlays2;
  478. Int_t layEnd = layBeg - 10;
  479. //for (Int_t lay = layMax; lay > 50; --lay) {
  480. for (Int_t lay = layBeg; lay > layEnd; --lay) {
  481. Int_t nLay = GetNofHitsInLayer(lay);
  482. Int_t indx0 = GetHitsInLayer(lay);
  483. Int_t view = (lay - 1) % 3; // view
  484. for (Int_t indx = 0; indx < nLay; ++indx) {
  485. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  486. // Exclude used hits
  487. if (hit->GetFlag() < 0) continue;
  488. FairMCPoint *ectP = (FairMCPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  489. Int_t idEct = ectP->GetTrackID();
  490. // !!! Exact match of TOF and ECT
  491. //if (fExact) {
  492. //Int_t itmp = 0;
  493. //if (fExact && idEct != tofP->GetTrackID()) continue;
  494. //if (HitMotherId(hit,tofP->GetTrackID(),itmp) != tofP->GetTrackID()) continue;
  495. //}
  496. Int_t detID = ectP->GetDetectorID();
  497. if (fWireMap.find(detID) == fWireMap.end()) {
  498. // Add position of wire point
  499. fWireMap.insert(pair<Int_t,pair<Double_t,Double_t> >(detID,
  500. pair<Double_t,Double_t>(ectP->GetX(),ectP->GetY())));
  501. }
  502. // Check Phi between ECT and ETOF
  503. Double_t rEct = drdz * hit->GetPos();
  504. if (view == 0) {
  505. // radial view
  506. phEct = hit->GetPhi();
  507. } else {
  508. // stereo view - find PHI of wire point at radius, extrapolated from ETOF
  509. //Double_t b = hit->GetXY(0) * hit->GetCos() + hit->GetXY(1) * hit->GetSin();
  510. //Double_t d = hit->GetXY(0) * hit->GetXY(0) + hit->GetXY(1) * hit->GetXY(1)
  511. Double_t b = fWireMap[detID].first * hit->GetCosSin(0) + fWireMap[detID].second * hit->GetCosSin(1);
  512. Double_t d = fWireMap[detID].first * fWireMap[detID].first + fWireMap[detID].second * fWireMap[detID].second
  513. - rEct * rEct;
  514. Double_t l1 = -b + TMath::Sqrt (b * b - d);
  515. Double_t l2 = -b - TMath::Sqrt (b * b - d);
  516. if (TMath::Abs(l1) > TMath::Abs(l2)) l1 = l2;
  517. //phEct = TMath::ATan2 (hit->GetXY(1) + l1 * hit->GetSin(),
  518. // hit->GetXY(0) + l1 * hit->GetCos());
  519. phEct = TMath::ATan2 (fWireMap[detID].second + l1 * hit->GetCosSin(1),
  520. fWireMap[detID].first + l1 * hit->GetCosSin(0));
  521. }
  522. dph = MpdKalmanFilter::Instance()->Proxim(phTof, phEct) - phTof;
  523. if (TMath::Abs(dph) > dphiMax) continue;
  524. // Add track seed
  525. //if (lun1 && idTpc == tof->GetTrackID()) fprintf(lun1,"%10.3e %10.3e %10.3e %10.3e\n",rTpc,rTof,dphi,slope-slope0);
  526. if (tofP->GetTrackID() < 99999) { tmp[tofP->GetTrackID()]++; idMax = TMath::Max(idMax,tofP->GetTrackID()); }
  527. MpdEctKalmanTrack *track =
  528. //new ((*fTracks)[nCand++]) MpdEctKalmanTrack(itof, indx0+indx, tof, hit, vert);
  529. new ((*fTrackCand)[nCand++]) MpdEctKalmanTrack(itof, idEct, tof, hit, vert);
  530. track->SetTrackID(idEct);
  531. EvalParams(tof, hit, track, rEct, phEct);
  532. // Eval. track Phi - linear extrapolation of 2 estimates
  533. //track->SetParam(2, phTof); // phi
  534. Double_t dx = tof->GetX() - rEct * TMath::Cos(phEct);
  535. Double_t dy = tof->GetY() - rEct * TMath::Sin(phEct);
  536. Double_t dr = TMath::Sqrt (dx * dx + dy * dy);
  537. Double_t phTofEct = TMath::ATan2 (dy, dx);
  538. Double_t ph = phTofEct;
  539. //Double_t ph = phEct - 1.0 * (MpdKalmanFilter::Instance()->Proxim(phTofEct,phTof)-phTofEct)/(rTof-dr) * rEct;
  540. track->SetParam(2, ph); // phi
  541. Double_t th = TMath::ATan2 (tof->GetZ()-hit->GetPos(), dr);
  542. track->SetParam(3, th); // dip angle
  543. Double_t pt = 1. / TMath::Abs(track->GetParam(4));
  544. /*
  545. track->SetParam(2, pmom.Phi()); // exact value - for test
  546. track->SetParam(3, TMath::PiOver2()-pmom.Theta()); // exact value - for test
  547. track->SetParam(4, TMath::Sign(1.,track->GetParam(4))/pmom.Pt()); // exact value - for test
  548. */
  549. /*
  550. if (track->GetParam(3) < 0) {
  551. // Going in backward direction
  552. --nCand;
  553. fTrackCand->RemoveAt(nCand);
  554. continue;
  555. }
  556. */
  557. track->SetParamNew(*track->GetParam());
  558. //track->GetParam()->Print();
  559. EvalCovar(tof, hit, track, rEct, phEct);
  560. //if (lun1 && idEct == tofP->GetTrackID()) fprintf(lun1,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e \n",track->GetParam(2),MpdKalmanFilter::Instance()->Proxim(track->GetParam(2),ph),track->GetParam(3),MpdKalmanFilter::Instance()->Proxim(track->GetParam(3),th),1./TMath::Sqrt((*track->GetWeight())(2,2)),1./TMath::Sqrt((*track->GetWeight())(3,3)),pmom.Mag(),rTof,pmom.Pt(),pt);
  561. if (lun1 && idEct == tofP->GetTrackID()) fprintf(lun1,"%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e \n",track->GetParam(2),MpdKalmanFilter::Instance()->Proxim(track->GetParam(2),pmom.Phi()),track->GetParam(3),MpdKalmanFilter::Instance()->Proxim(track->GetParam(3),TMath::PiOver2()-pmom.Theta()),1./TMath::Sqrt((*track->GetWeight())(2,2)),1./TMath::Sqrt((*track->GetWeight())(3,3)),pmom.Mag(),rTof,pmom.Pt(),pt);
  562. //cout << nCand-1 << " " << lay << " " << track->GetTrackID() << " " << tofP->GetTrackID() << endl;
  563. } // for (Int_t indx = 0; indx < nLay;
  564. } // for (Int_t lay = layMax; lay > 50;
  565. } // for (Int_t itof = 0; itof < nTof;
  566. for (Int_t j = 0; j <= idMax; ++j) if (tmp[j] > 0) ++nID;
  567. cout << " Number of ECT track candidates: " << nCand << " " << nID << endl;
  568. //if (lun1) fclose(lun1);
  569. //exit(0);
  570. }
  571. //__________________________________________________________________________
  572. void MpdEctTrackFinderTof::EvalParams(const MpdTofHit *tof, const MpdKalmanHit *ect,
  573. MpdEctKalmanTrack *track, Double_t rEct, Double_t phEct)
  574. {
  575. /// Evaluate track parameters
  576. // Evaluate signed track Pt (curvature) assuming the track coming from the
  577. // primary vertex
  578. //Double_t rTof = TMath::Sqrt (tof->GetX() * tof->GetX() + tof->GetY() * tof->GetY());
  579. Double_t phTof = TMath::ATan2 (tof->GetY(), tof->GetX());
  580. //TVector2 vecTof(rTof*TMath::Cos(phTof)-0.,rTof*TMath::Sin(phTof)-0.);
  581. TVector2 vecTof(tof->GetX()-0.,tof->GetY()-0.);
  582. TVector2 vecEct(rEct*TMath::Cos(phEct)-0.,rEct*TMath::Sin(phEct)-0.);
  583. TVector2 dvec = vecTof - vecEct;
  584. Double_t cosAlpha = vecEct * dvec / vecEct.Mod() / dvec.Mod();
  585. Double_t rad = vecTof.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  586. Double_t factor = 0.0015; // 0.3 * 0.5T * 0.01
  587. //Double_t charge = phTof - MpdKalmanFilter::Instance()->Proxim(phTof,phTpc);
  588. Double_t charge = MpdKalmanFilter::Instance()->Proxim(phTof,phEct) - phTof;
  589. Double_t pt = factor * TMath::Abs(rad);
  590. pt = TMath::Min (pt, 2.5) * TMath::Sign(1., -charge);
  591. track->SetParam(0, vecTof.X()); // X
  592. track->SetParam(1, vecTof.Y()); // Y
  593. track->SetParam(4, 1./pt); // 1/pt
  594. }
  595. //__________________________________________________________________________
  596. void MpdEctTrackFinderTof::EvalCovar(const MpdTofHit *tof, const MpdKalmanHit *ect,
  597. MpdEctKalmanTrack *track, Double_t rEct, Double_t phEct)
  598. {
  599. /// Evaluate covariance matrix for the track seed
  600. static const Double_t xErrTof = 2. / TMath::Sqrt(12.), yErrTof = xErrTof;
  601. //static const Double_t xErrTof = 0.1, yErrTof = xErrTof;
  602. TMatrixDSym weight(5);
  603. weight(0,0) = xErrTof * xErrTof; // <xx>
  604. //weight(0,0) *= 4.; // extra factor of 2
  605. weight(1,1) = yErrTof * yErrTof; // <yy>
  606. //weight(1,1) *= 4.; // extra factor of 2
  607. Double_t xEct = rEct * TMath::Cos(phEct);
  608. Double_t yEct = rEct * TMath::Sin(phEct);
  609. Double_t xTof = tof->GetX(), yTof = tof->GetY();
  610. Double_t dx = xTof - xEct, dy = yTof - yEct;
  611. Double_t dist2 = dx * dx + dy * dy;
  612. weight(2,2) = xErrTof * xErrTof / dist2; // <PhiPhi>
  613. //weight(2,2) *= 2.; // extra factor of sqrt(2)
  614. //Double_t tanThe = TMath::Tan (track->GetParam(3));
  615. Double_t rTof = TMath::Sqrt (tof->GetX() * tof->GetX() + tof->GetY() * tof->GetY());
  616. //Double_t dR = rTof - rEct;
  617. //Double_t dz = tof->GetZ() - ect->GetPos();
  618. Double_t dR = rTof;
  619. Double_t dz = tof->GetZ();
  620. Double_t dThdR = dz / (dR * dR + dz * dz);
  621. weight(3,3) = dThdR * dThdR * xErrTof * xErrTof; // <TheThe>
  622. //weight(3,3) /= 4.; // extra factor of 0.5
  623. weight(3,3) *= 4.; // extra factor of 2
  624. //weight(4,4) = (track->GetParam(4)*0.5) * (track->GetParam(4)*0.5); // error 50%
  625. weight(4,4) = (track->GetParam(4)*1.) * (track->GetParam(4)*1.); // error 100%
  626. //weight(4,4) = TMath::Max (1.5 * 1.5, weight(4,4));
  627. //weight(4,4) = (track->GetParam(4)*2.) * (track->GetParam(4)*2.); // error 200%
  628. //weight(4,4) = (track->GetParam(4)*4.) * (track->GetParam(4)*4.); // error 400%
  629. //weight.Print();
  630. //fWeight->Invert(); // weight matrix
  631. Int_t iok = 0;
  632. MpdKalmanFilter::Instance()->MnvertLocal(weight.GetMatrixArray(), 5, 5, 5, iok);
  633. //fWeight->Print();
  634. track->SetWeight(weight);
  635. }
  636. //__________________________________________________________________________
  637. void MpdEctTrackFinderTof::DoTracking(Int_t iPass)
  638. {
  639. /// Run Kalman tracking
  640. Int_t nCand = fTrackCand->GetEntriesFast(), iok = -1;
  641. TMatrixD param(5,1);
  642. TMatrixDSym weight(5), pointWeight(5);
  643. for (Int_t i = 0; i < nCand; ++i) {
  644. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  645. //cout << " Track seed No. " << i << " " << track->GetTrackID() << " " ;//<< endl;
  646. //if (track->GetTrackID() != 117) continue;
  647. //if (track->GetTrackID() != 10) continue;
  648. //if (track->GetTrackID() != 1105) continue;
  649. //if (track->GetTrackID() != 77) continue;
  650. //(*(track->GetParamNew()))(4,0) = -0.5; // just for check
  651. /*
  652. for (Int_t k = 0; k < 5; ++k) {
  653. for (Int_t j = i; j < 5; ++j) {
  654. if (j == i) (*track->GetWeight())(i,j) /= 100.;
  655. else (*track->GetWeight())(i,j) = (*track->GetWeight())(j,i) = 0.;
  656. }
  657. }
  658. */
  659. // Add first ECT hit
  660. //track->Weight2Cov()->Print();
  661. //track->GetParamNew()->Print();
  662. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(0);
  663. MpdKalmanFilter::Instance()->PropagateToHit(track, hit);
  664. //cout << track->GetLength() << endl;
  665. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(track,hit,pointWeight,param);
  666. track->SetChi2(track->GetChi2()+dChi2);
  667. weight = *track->GetWeight();
  668. weight += pointWeight;
  669. track->SetWeight(weight);
  670. track->SetParamNew(param);
  671. //if (i == 0) track->SetLength(0.); // for correct track length
  672. //track->Weight2Cov()->Print();
  673. //track->GetParamNew()->Print();
  674. // Start ECT tracking from different layers to account for missing hits
  675. Int_t laySkip = 2, lay0 = hit->GetLayer() - 1;
  676. //Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  677. MpdEctKalmanTrack tracks[10];
  678. for (Int_t j = 0; j < laySkip; ++j) {
  679. tracks[j] = *track;
  680. //cout << track->GetParamNew(0) << endl;
  681. //cout << i << " " << lay << " " << tracks[lay].GetNofHits() << " " << tracks[lay].GetChi2() << " " << tracks[lay].GetParam(0) << endl;
  682. if (j > 0) {
  683. lay0 = -1;
  684. // Find first layer after the missing one
  685. Int_t nh = tracks[j-1].GetNofHits(), lastLay = hit->GetLayer();
  686. TObjArray *hits = tracks[j-1].GetHits();
  687. for (Int_t jh = 1; jh < nh; ++jh) {
  688. hit = (MpdKalmanHit*) hits->UncheckedAt(jh);
  689. if (TMath::Abs(hit->GetLayer()-lastLay) == 1) {
  690. lastLay = hit->GetLayer();
  691. tracks[j].GetHits()->Add(hit);
  692. continue;
  693. }
  694. lay0 = hit->GetLayer() - 1;
  695. break;
  696. }
  697. if (lay0 >= 0) {
  698. // Refit track to the last kept hit
  699. nh = tracks[j].GetNofHits();
  700. for (Int_t jh = 1; jh < nh; ++jh) {
  701. hit = (MpdKalmanHit*) tracks[j].GetHits()->UncheckedAt(jh);
  702. MpdKalmanFilter::Instance()->PropagateToHit(&tracks[j], hit);
  703. dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(&tracks[j],hit,pointWeight,param);
  704. tracks[j].SetChi2(tracks[j].GetChi2()+dChi2);
  705. weight = *tracks[j].GetWeight();
  706. weight += pointWeight;
  707. tracks[j].SetWeight(weight);
  708. tracks[j].SetParamNew(param);
  709. if (jh == nh - 1) {
  710. tracks[j].SetWeightAtHit(weight);
  711. tracks[j].SetParamAtHit(param);
  712. tracks[j].SetLengAtHit(tracks[j].GetLength());
  713. }
  714. }
  715. }
  716. } // if (j > 0)
  717. //if (lay0 >= 0) iok = RunKalmanFilter(&tracks[j], lay0);
  718. if (lay0 >= 0) iok &= RunKalmanFilter(&tracks[j], lay0);
  719. else --laySkip;
  720. //cout << iok << " " << tracks[j].GetMisses() << " ";//<< endl;
  721. }
  722. //cout << endl;
  723. // Select the best track (with max number of hits)
  724. Int_t nHitsMax = tracks[0].GetNofHits(), iMax = 0;
  725. //if (iok > -9) {
  726. if (iok >= 0) {
  727. for (Int_t j = 1; j < laySkip; ++j) {
  728. Int_t nhits = tracks[j].GetNofHits();
  729. if (nhits > nHitsMax) {
  730. nHitsMax = nhits;
  731. iMax = j;
  732. } else if (nhits == nHitsMax) {
  733. if (tracks[j].GetChi2() < tracks[iMax].GetChi2()) {
  734. iMax = j;
  735. nHitsMax = nhits;
  736. }
  737. }
  738. }
  739. }
  740. fTrackCand->RemoveAt(i);
  741. if (iok < 0) continue;
  742. track = new ((*fTrackCand)[i]) MpdEctKalmanTrack(tracks[iMax]);
  743. // Refit
  744. if (1) {
  745. //track->Weight2Cov()->Print();
  746. //track->GetParamNew()->Print();
  747. Double_t pos = track->GetPosNew();
  748. track->SetDirection(MpdKalmanTrack::kOutward);
  749. MpdKalmanFilter::Instance()->Refit(track, -1);
  750. //track->GetParamNew()->Print();
  751. track->SetDirection(MpdKalmanTrack::kInward);
  752. MpdKalmanFilter::Instance()->Refit(track, 1);
  753. track->SetParamAtHit(*track->GetParamNew());
  754. track->SetWeightAtHit(*track->GetWeight());
  755. //track->SetLengAtHit(track->GetLength());
  756. MpdKalmanHit hitTmp;
  757. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  758. hitTmp.SetPos(pos);
  759. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp,kFALSE);
  760. }
  761. //iok = RunKalmanFilter(track, hit->GetLayer()-1);
  762. //cout << i << " " << track->GetNofHits() << endl;
  763. if (0 && iok == 0 && track->GetNofHits() < 0) {
  764. //*
  765. hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(0);
  766. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  767. track->SetPos(hit->GetPos()-dir);
  768. track->SetPosNew(track->GetPos());
  769. //for (Int_t j = 0; j < 5; ++j) track->SetParam(j,track->GetParam1()[j]);
  770. //track->SetParamNew(*track->GetParam());
  771. //MpdKalmanHitZ hitTmp(*hit);
  772. //hitTmp.SetZ(track->GetPosNew()-dir);
  773. //MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  774. MpdKalmanFilter::Instance()->Refit(track,1);
  775. //*/
  776. /*
  777. track->Print("");
  778. Double_t dir = TMath::Sign(1.,track->GetPosNew());
  779. MpdKalmanHitZ *hit = (MpdKalmanHitZ*) track->GetHits()->UncheckedAt(0);
  780. MpdKalmanHitZ hitTmp(*hit);
  781. hitTmp.SetZ(track->GetPosNew()+dir);
  782. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  783. track->Print("");
  784. MpdKalmanFilter::Instance()->Refit(track,-1);
  785. track->Print("");
  786. hitTmp.SetZ(track->GetPosNew()-dir);
  787. MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp);
  788. track->Print("");
  789. MpdKalmanFilter::Instance()->Refit(track,1);
  790. track->Print("");
  791. */
  792. }
  793. }
  794. }
  795. //__________________________________________________________________________
  796. Int_t MpdEctTrackFinderTof::RunKalmanFilter(MpdEctKalmanTrack *track, Int_t layBeg)
  797. {
  798. /// Run Kalman filter
  799. const Double_t rMin = 49.2, rMax = 110.5; // min and max radial size of ECT - to be taken elsewhere
  800. //cout << fHits->GetEntriesFast() << endl;
  801. Int_t layMax = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  802. MpdKalmanHit *hitOK = 0x0;
  803. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  804. //Int_t layBeg = layMax, layEnd = -1, dLay = -1, layOK = -1;
  805. Int_t layEnd = -1, dLay = -1;//, layOK = -1;
  806. if (track->GetPosNew() < 0) layEnd = fgkNlays2;
  807. if (trackDir == MpdKalmanTrack::kOutward) {
  808. layEnd = layMax + 1;
  809. dLay = 1;
  810. if (track->GetPosNew() > 0 && layMax > fgkNlays2) layEnd -= fgkNlays2;
  811. }
  812. //Int_t indxOK = hits->IndexOf(hitOK);
  813. //Int_t nHits = hits->GetEntriesFast();
  814. Int_t miss = 0, missSec = 0;
  815. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5);
  816. TMatrixD param(5,1), paramTmp(5,1);
  817. Double_t saveZ = 0;
  818. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  819. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  820. Int_t nLay = GetNofHitsInLayer(lay);
  821. Int_t indx0 = GetHitsInLayer(lay);
  822. Double_t dChi2Min = 1.e+6;
  823. MpdKalmanHit *hitMin = 0x0;
  824. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  825. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1;
  826. /*
  827. if (trackDir == MpdKalmanTrack::kOutward) {
  828. // !!! This part is due to the simplified hit merging (digitization) -
  829. // different hit position inside one layer - should be removed later
  830. indxBeg = nLay - 1;
  831. indxEnd = -1;
  832. dIndx = -1;
  833. }
  834. */
  835. Double_t radNew = -1;
  836. //for (Int_t indx = 0; indx < nLay; ++indx) {
  837. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  838. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  839. // !!! Exact ID match
  840. Int_t itmp = 0;
  841. //if (fExact && ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() != track->GetTrackID()) continue;
  842. if (fExact && HitMotherId(hit,track->GetTrackID(),itmp) != track->GetTrackID()) continue;
  843. // Exclude used hits
  844. if (hit->GetFlag() < 0) continue;
  845. //cout << " window:" << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << hit->GetZ() << " " << track->GetParamNew(1) << endl;
  846. // Check if the hit within some window (15x15cm for the moment - check!!!)
  847. //if (TMath::Abs(hit->GetRphi()-track->GetParamNew(0)) > 9) continue;
  848. //if (TMath::Abs(Proxim(track,hit)-track->GetParamNew(0)) > 15) continue;
  849. TVector2 dist = GetDistance(track, hit);
  850. if (dist.X() > 15.) continue; // distance in transverse to the tube direction
  851. if (hit->GetNofDim() > 1 && dist.Y() > 25.) continue; // distance in R
  852. //*if (hit->GetTrackID() == 186)*/ cout << " Hit: " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() << " " << hit->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetDetectorID() << " " << hit->GetRphi() << " " << hit->GetR() << " " << hit->GetZ() << " " << dist.X() << " " << dist.Y() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  853. //track->Print("");
  854. //hit->Print("");
  855. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit)) { track->SetMisses(miss); return -1; }
  856. //*if (hit->GetTrackID() == -607)*/ cout << /*hit->GetTrackID() <<*/ " " << hit->GetRphi() << " " << track->GetParamNew(0) << " " << track->GetParamNew(1) << " " << hit->GetZ() << " " << track->GetPosNew() << endl;
  857. //
  858. // For debugging
  859. /*
  860. TVector2 dist0 = GetDistance(track, hit);
  861. cout << dist0.X() << " ";
  862. MpdKalmanHitZ hitDbg = *hit;
  863. Double_t xDbg = hit->GetXY(0) * TMath::Cos(hit->GetPhi()) + hit->GetXY(1) * TMath::Sin(hit->GetPhi());
  864. Double_t yDbg = -hit->GetXY(0) * TMath::Sin(hit->GetPhi()) + hit->GetXY(1) * TMath::Cos(hit->GetPhi());
  865. hitDbg.SetRphi(yDbg);
  866. hitDbg.SetR(xDbg);
  867. dist = GetDistance(track, &hitDbg);
  868. cout << dist.X() << endl;
  869. */
  870. //Double_t radNew = track->GetRadNew();
  871. if (radNew < 0) radNew = track->GetRadNew();
  872. //if (radNew < rMin || radNew > rMax) return -1;
  873. // Going to interaction point
  874. if (radNew > rMax) break; // next layer
  875. //else if (radNew < rMin) return -1;
  876. else if (radNew < rMin) { track->SetMisses(miss); return 0; }
  877. //Double_t err = hit->GetRphiErr();
  878. //if (track->GetNofHits() == 0) hit->SetRphiErr(0.04);
  879. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHitZ(track,hit,pointWeight,param);
  880. //if (track->GetNofHits() == 0) hit->SetRphiErr(err);
  881. //if (param(3,0) < 0) { cout << " Check: " << param(3,0) << " " << dChi2 << " " << (*fParamNew)(3,0) << " " << hit->GetRphi() << " " << hit->GetZ() << endl; fParamNew->Print();}
  882. if (dChi2 < dChi2Min) {
  883. //if (dChi2 < dChi2Min && dist0.X() <= dist.X()) {
  884. //if (dChi2 < dChi2Min && dist.X() <= 0.2) {
  885. dChi2Min = dChi2;
  886. hitMin = hit;
  887. saveWeight = *track->GetWeight();
  888. saveZ = track->GetPosNew();
  889. // temporary storage for the current track
  890. paramTmp = param;
  891. pointWeightTmp = pointWeight;
  892. //cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetR() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hit->GetIndex()))->GetTrackID() << endl;
  893. //cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  894. //cout << hit->GetRphi() << " " << hit->GetZ() << endl;
  895. //cout << param(0,0) << " " << param(1,0) << endl;
  896. //paramTmp.Print();
  897. }
  898. } // for (Int_t indx = indxBeg; indx != indxEnd;
  899. Double_t cut = fgkChi2Cut;
  900. //if (track->GetNofHits() == 0) cut /= 2.;
  901. //if (dChi2Min < fgkChi2Cut) {
  902. if (dChi2Min < cut) {
  903. //layOK = lay;
  904. hitOK = hitMin;
  905. track->GetHits()->Add(hitOK);
  906. //miss = 0;
  907. // Restore track params at the best hit
  908. track->SetChi2(track->GetChi2()+dChi2Min);
  909. saveWeight += pointWeightTmp;
  910. track->SetWeight(saveWeight);
  911. track->SetPosNew(saveZ);
  912. track->SetParamNew(paramTmp);
  913. //cout << " *** Adding hit: " << hitOK->GetLayer() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetTrackID() << " " << ((FairMCPoint*)fEctHits->UncheckedAt(hitOK->GetIndex()))->GetDetectorID() << " " << dChi2Min << " " << track->GetChi2() << endl;
  914. //paramTmp.Print();
  915. // Check if the accepted hit is the same as the seeded hit
  916. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  917. if (track->GetNofHits() == 1) track->SetParam1();
  918. if (track->GetDirection() == MpdKalmanTrack::kInward) {
  919. // Save track params at last hit
  920. track->SetLengAtHit(track->GetLength());
  921. track->SetParamAtHit(paramTmp);
  922. track->SetWeightAtHit(saveWeight);
  923. track->SetPosAtHit(saveZ);
  924. }
  925. missSec = 0;
  926. } else if (lay > 0 && radNew < rMax && radNew > rMin) {
  927. // Check whether tracks goes through a tube
  928. Int_t layM1 = lay - 1;
  929. if (nLay == 0) {
  930. // No hits in the layer - extrapolate track
  931. MpdKalmanHit hitTmp;
  932. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  933. hitTmp.SetDist(fZplanes[layM1]);
  934. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,&hitTmp)) { track->SetMisses(miss); return -1; }
  935. }
  936. FairGeoVector vec(track->GetParamNew(0)*10,track->GetParamNew(1)*10,track->GetPosNew()*10); // mm
  937. //cout << vec.getX() << " " << vec.getY() << " " << vec.getZ() << endl;
  938. Int_t iTube1 = 0, iter = 0, idTubes = 0;//, iTubeOk = 0;
  939. Int_t iTube = Int_t (TMath::ATan2 (vec.getY(),vec.getX()) / TMath::TwoPi() * fNofTubes);
  940. if (iTube < 0) iTube += fNofTubes;
  941. iTube = iTube % fNofTubes;
  942. FairGeoVector vecLoc = fTransf[layM1][iTube]->transTo(vec);
  943. //cout << vecLoc.getX() << " " << vecLoc.getY() << " " << vecLoc.getZ() << " " << iTube << endl;
  944. //cout << transf[lay][0]->getTranslation() << endl;
  945. //transf[lay][0]->getRotMatrix().print();
  946. /*if (TMath::Abs(vecLoc.getZ()) > z2) cout << vecLoc << " " << ph << endl;
  947. if (TMath::Abs(vecLoc.getZ()) > z2) {
  948. iTube = nTubes / 2;
  949. vecLoc = transf[lay][iTube]->transTo(vec);
  950. cout << vecLoc;
  951. }*/
  952. Double_t dyMin = 0;
  953. while (1) {
  954. FairGeoVector vecLoc1 = fTransf[layM1][(iTube+1)%fNofTubes]->transTo(vec);
  955. Double_t dy = vecLoc.getY() - vecLoc1.getY();
  956. idTubes = TMath::Nint (vecLoc.getY() / dy);
  957. iTube1 = iTube + idTubes;
  958. if (iTube1 < 0) iTube1 += fNofTubes;
  959. iTube1 = iTube1 % fNofTubes;
  960. vecLoc1 = fTransf[layM1][iTube1]->transTo(vec);
  961. dyMin = TMath::Abs(vecLoc1.getY());
  962. //iTubeOk = iTube1;
  963. //cout << vecLoc1.getX() << " " << vecLoc1.getY() << " " << vecLoc1.getZ() << " " << iTubeOk << " " << dyMin << " " << iter << endl;
  964. if (dyMin <= frMinMax[2] || (iter && TMath::Abs(idTubes) < 2 && TMath::Abs(dyMin/dy) < 0.51)) break;
  965. iTube = iTube1;
  966. vecLoc = vecLoc1;
  967. ++iter;
  968. }
  969. //if (dyMin > 4.) ++miss; // too far from the tube
  970. //if (dyMin < frMinMax[2]+0.5) ++miss; // too far from the tube
  971. if (dyMin < frMinMax[2]-0.2) { ++miss; ++missSec; } // going thru the tube
  972. //cout << dyMin << " " << dChi2Min << endl;
  973. }
  974. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  975. if (missSec > 5) { track->SetMisses(miss); return -9; }
  976. } // for (Int_t lay = layBeg; lay != layEnd;
  977. track->SetMisses(miss);
  978. return 0;
  979. }
  980. //__________________________________________________________________________
  981. TVector2 MpdEctTrackFinderTof::GetDistance(MpdEctKalmanTrack *track, MpdKalmanHit *hit)
  982. {
  983. /// Compute distance between track and hit
  984. Double_t xTr, yTr;
  985. if (track->GetType() == MpdKalmanTrack::kEndcap) {
  986. xTr = track->GetParamNew(0);
  987. yTr = track->GetParamNew(1);
  988. } else {
  989. Double_t rPhi = track->GetParamNew(0);
  990. Double_t r = track->GetPosNew();
  991. Double_t ph = rPhi / r;
  992. xTr = r * TMath::Cos(ph);
  993. yTr = r * TMath::Sin(ph);
  994. }
  995. Double_t phi = hit->GetPhi();
  996. Double_t cosPhi = TMath::Cos(phi);
  997. Double_t sinPhi = TMath::Sin(phi);
  998. // Transform track coordinates to local tube coordinates
  999. Double_t xLoc = xTr * cosPhi + yTr * sinPhi;
  1000. Double_t yLoc = -xTr * sinPhi + yTr * cosPhi;
  1001. //Double_t xLoc = (xTr - hit->GetXY(0)) * cosPhi + (yTr - hit->GetXY(1)) * sinPhi;
  1002. //Double_t yLoc = -(xTr - hit->GetXY(0)) * sinPhi + (yTr - hit->GetXY(1)) * cosPhi;
  1003. TVector2 dist(TMath::Abs(yLoc-hit->GetMeas(0)), TMath::Abs(xLoc-hit->GetMeas(1)));
  1004. return dist;
  1005. }
  1006. //__________________________________________________________________________
  1007. Double_t MpdEctTrackFinderTof::Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit)
  1008. {
  1009. /// Adjust hit coord. R-Phi to be "around" track R-Phi - to avoid
  1010. /// discontinuity around +- Pi
  1011. Fatal("Proxim", " !!! Not implemented !!!");
  1012. /*
  1013. Double_t hitPhi = hit->GetRphi() / hit->GetR();
  1014. Double_t phi0 = track->GetParamNew(0) / track->GetPosNew();
  1015. return hit->GetR() * MpdKalmanFilter::Instance()->Proxim(phi0, hitPhi);
  1016. */
  1017. return 0;
  1018. }
  1019. //__________________________________________________________________________
  1020. void MpdEctTrackFinderTof::Write()
  1021. {
  1022. /// Write
  1023. TFile histoFile("EctRec.root","RECREATE");
  1024. Writedir2current(fHistoDir);
  1025. histoFile.Close();
  1026. }
  1027. //__________________________________________________________________________
  1028. void MpdEctTrackFinderTof::Writedir2current( TObject *obj )
  1029. {
  1030. /// Write
  1031. if( !obj->IsFolder() ) obj->Write();
  1032. else{
  1033. TDirectory *cur = gDirectory;
  1034. TDirectory *sub = cur->mkdir(obj->GetName());
  1035. sub->cd();
  1036. TList *listSub = ((TDirectory*)obj)->GetList();
  1037. TIter it(listSub);
  1038. while( TObject *obj1=it() ) Writedir2current(obj1);
  1039. cur->cd();
  1040. }
  1041. }
  1042. //__________________________________________________________________________
  1043. void MpdEctTrackFinderTof::RemoveDoubles()
  1044. {
  1045. /// Remove double tracks (if number of common hits greater than 50% of hits on track)
  1046. Int_t nFound = fTrackCand->GetEntriesFast(), nOK = 0;
  1047. for (Int_t i = 0; i < nFound; ++i) {
  1048. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  1049. if (track == 0x0) continue;
  1050. Int_t nh = track->GetNofHits();
  1051. //cout << i << " " << nh << " " << ++nOK << endl;
  1052. for (Int_t j = i + 1; j < nFound; ++j) {
  1053. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(j);
  1054. if (track1 == 0x0) continue;
  1055. Int_t nh1 = track1->GetNofHits();
  1056. Int_t nc = NofCommonHits(track, track1);
  1057. if (1.*nc/TMath::Min(nh,nh1) < 0.5) continue;
  1058. if (nh > nh1) fTrackCand->RemoveAt(j);
  1059. //if (nh > nh1+5) fTrackCand->RemoveAt(j);
  1060. else if (nh < nh1) {
  1061. //else if (nh < nh1-5) {
  1062. fTrackCand->RemoveAt(i);
  1063. --nOK;
  1064. break;
  1065. } else {
  1066. if (track->GetChi2() > track1->GetChi2()) {
  1067. //if (track->GetChi2()/nh > track1->GetChi2()/nh1) {
  1068. fTrackCand->RemoveAt(i);
  1069. --nOK;
  1070. break;
  1071. }
  1072. fTrackCand->RemoveAt(j);
  1073. }
  1074. }
  1075. }
  1076. fTrackCand->Compress();
  1077. // Remove double tracks (originated from the same ETOF hit)
  1078. //*
  1079. nFound = fTrackCand->GetEntriesFast();
  1080. for (Int_t i = 0; i < nFound; ++i) {
  1081. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  1082. if (track == 0x0) continue;
  1083. Int_t iTof = track->GetTofIndex();
  1084. Int_t nh = track->GetNofHits();
  1085. for (Int_t j = i + 1; j < nFound; ++j) {
  1086. MpdEctKalmanTrack *track1 = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(j);
  1087. if (track1 == 0x0) continue;
  1088. Int_t iTof1 = track1->GetTofIndex();
  1089. if (iTof1 != iTof) continue;
  1090. Int_t nh1 = track1->GetNofHits();
  1091. if (nh > nh1) fTrackCand->RemoveAt(j);
  1092. else if (nh < nh1) {
  1093. fTrackCand->RemoveAt(i);
  1094. break;
  1095. } else {
  1096. if (track->GetChi2() > track1->GetChi2()) {
  1097. fTrackCand->RemoveAt(i);
  1098. break;
  1099. }
  1100. fTrackCand->RemoveAt(j);
  1101. }
  1102. }
  1103. }
  1104. fTrackCand->Compress();
  1105. //*/
  1106. }
  1107. //__________________________________________________________________________
  1108. Int_t MpdEctTrackFinderTof::NofCommonHits(MpdEctKalmanTrack* track, MpdEctKalmanTrack* track1)
  1109. {
  1110. /// Compute number of common hits on 2 tracks
  1111. TObjArray *hits = track->GetHits(), *hits1 = track1->GetHits();
  1112. MpdKalmanHit *hit = (MpdKalmanHit*) hits->First();
  1113. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->First();
  1114. Int_t dir = 1;
  1115. if (hit->GetLayer() < ((MpdKalmanHit*) hits->Last())->GetLayer()) dir = -1;
  1116. Int_t nCom = 0;
  1117. while (hit && hit1) {
  1118. if (dir * hit->GetLayer() < dir * hit1->GetLayer()) {
  1119. hit1 = (MpdKalmanHit*) hits1->After(hit1);
  1120. continue;
  1121. }
  1122. if (hit == hit1) ++nCom;
  1123. hit = (MpdKalmanHit*) hits->After(hit);
  1124. }
  1125. return nCom;
  1126. }
  1127. //__________________________________________________________________________
  1128. void MpdEctTrackFinderTof::StoreTracks(Int_t iPass)
  1129. {
  1130. /// Transfer tracks from fTrackCand to fTracks
  1131. static const Int_t nHitMin = 10;
  1132. //static const Double_t etaMin = 1.3, etaMax = 1.5;
  1133. static const Double_t etaMin = 1.6, etaMax = 1.5;
  1134. Int_t nFound = fTracks->GetEntriesFast(), nCand = fTrackCand->GetEntriesFast();
  1135. for (Int_t i = 0; i < nCand; ++i) {
  1136. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTrackCand->UncheckedAt(i);
  1137. Double_t eta = TMath::Abs(track->Momentum3().Eta());
  1138. //if (iPass == 0 && track->GetChi2() / (track->GetNofHits()-5) > 3.0) continue;
  1139. //if (iPass == 0 && eta < etaMin) continue;
  1140. if (eta < etaMin) continue;
  1141. if (eta > etaMax && track->GetNofHits() < nHitMin) continue;
  1142. track->Weight2Cov();
  1143. //SetTrackID(track);
  1144. new ((*fTracks)[nFound++]) MpdEctKalmanTrack(*track);
  1145. //fTpcTracks->RemoveAt(track->GetTpcIndex());
  1146. //fTrackCand->RemoveAt(i);
  1147. }
  1148. }
  1149. //__________________________________________________________________________
  1150. void MpdEctTrackFinderTof::AddHits()
  1151. {
  1152. /// Add hit objects to tracks and compute number of wrongly assigned hits
  1153. /// (hits with ID different from ID of starting track)
  1154. Int_t nFound = fTracks->GetEntriesFast();
  1155. for (Int_t i = 0; i < nFound; ++i) {
  1156. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1157. if (track->GetTofIndex() < 0) continue; // track from TPC
  1158. Int_t nHits = track->GetNofHits();
  1159. if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1160. track->SetNofHits(nHits);
  1161. TClonesArray &trHits = *track->GetTrHits();
  1162. SetTrackID(track); // set track ID as ID of majority of its hits
  1163. TObjArray *hits = track->GetHits();
  1164. Int_t nWrong = 0, motherID = track->GetTrackID(), motherID1 = 0, nMirror = 0;
  1165. cout << i << " " << nHits << " " << track->GetTrackID() << " " << track->GetChi2() << " " << track->Phi()*TMath::RadToDeg()
  1166. << " " << track->Momentum3().Eta() << " " << 1./track->GetParam(4) << " "
  1167. << track->Momentum() << " " << track->GetMisses() << endl;
  1168. // Get track mother ID
  1169. MpdMCTrack *mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(motherID);
  1170. while (mctrack->GetMotherId() >= 0) {
  1171. motherID = mctrack->GetMotherId();
  1172. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1173. }
  1174. for (Int_t j = 0; j < nHits; ++j) {
  1175. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1176. new (trHits[j]) MpdKalmanHit(*hit);
  1177. Int_t iproj = (hit->GetLayer() - 1) % 3;
  1178. if (iproj == 0) cout << " R";
  1179. else if (iproj == 1) cout << " U";
  1180. else cout << " V";
  1181. cout << hit->GetLayer();
  1182. Int_t hitId = HitMotherId(hit,motherID,motherID1);
  1183. cout << "-" << hitId;
  1184. if (hit->Index()->GetSize() > 1) cout << "-" << hit->Index()->GetSize();
  1185. if (hit->IsMirror()) cout << "(M)";
  1186. if (motherID1 != motherID) nWrong++;
  1187. else if (hit->IsMirror()) ++nMirror;
  1188. }
  1189. cout << "\n" << nWrong << " " << track->GetTrackID() << " " << motherID << " " << nMirror << endl;
  1190. track->SetNofWrong(nWrong);
  1191. track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
  1192. }
  1193. fTracks->Compress();
  1194. }
  1195. //__________________________________________________________________________
  1196. Int_t MpdEctTrackFinderTof::HitMotherId(MpdKalmanHit* hit, Int_t idM, Int_t &id1)
  1197. {
  1198. /// Check if hit has the same mother ID as idM
  1199. Int_t nOver = hit->Index()->GetSize(), idHit = 0;
  1200. for (Int_t i = 0; i < nOver; ++i) {
  1201. FairMCPoint *h = (FairMCPoint*) fEctHits->UncheckedAt(hit->GetIndex(i));
  1202. id1 = idHit = h->GetTrackID();
  1203. // Get point mother ID
  1204. MpdMCTrack* mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(id1);
  1205. while (mctrack->GetMotherId() >= 0) {
  1206. id1 = mctrack->GetMotherId();
  1207. mctrack = (MpdMCTrack*) fMCTracks->UncheckedAt(mctrack->GetMotherId());
  1208. }
  1209. if (id1 == idM) return idHit;
  1210. }
  1211. return idHit;
  1212. }
  1213. //__________________________________________________________________________
  1214. void MpdEctTrackFinderTof::SetTrackID(MpdEctKalmanTrack* track)
  1215. {
  1216. /// Set track ID as ID of majority of its hits
  1217. multiset<Int_t> ids;
  1218. multiset<Int_t>::iterator it;
  1219. Int_t nHits = track->GetNofHits(), nMax = 0, idMax = -1;
  1220. TObjArray *hits = track->GetHits();
  1221. for (Int_t i = 0; i < nHits; ++i) {
  1222. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  1223. MpdStrawendcapPoint *p = (MpdStrawendcapPoint*) fEctHits->UncheckedAt(hit->GetIndex());
  1224. Int_t id = p->GetTrackID();
  1225. ids.insert(id);
  1226. }
  1227. for (it = ids.begin(); it != ids.end(); ++it) {
  1228. Int_t nC = ids.count(*it);
  1229. if (nC > nMax) {
  1230. nMax = nC;
  1231. idMax = *it;
  1232. }
  1233. }
  1234. if (nMax > 1) track->SetTrackID(idMax);
  1235. }
  1236. //__________________________________________________________________________
  1237. void MpdEctTrackFinderTof::SelectTracks(Int_t iPass)
  1238. {
  1239. /// Do track selection and compute shared hit multiplicities
  1240. //if (iPass) return;
  1241. Int_t nFound = fTracks->GetEntriesFast(), nHitsTot = fKHits->GetEntriesFast();
  1242. Int_t *index = new Int_t [nHitsTot];
  1243. for (Int_t i = 0; i < nHitsTot; ++i) index[i] = 0;
  1244. for (Int_t i = 0; i < nFound; ++i) {
  1245. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1246. if (track->GetTofIndex() < 0) continue; // track from TPC
  1247. Int_t nHits = track->GetNofHits();
  1248. //if (nHits == 0) { fTracks->RemoveAt(i); continue; }
  1249. if (nHits == 0) continue;
  1250. TObjArray *hits = track->GetHits();
  1251. for (Int_t j = 0; j < nHits; ++j) {
  1252. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1253. Int_t ind = fKHits->IndexOf(hit);
  1254. if (index[ind]) hit->SetFlag(TMath::Abs(hit->GetFlag())+10);
  1255. ++index[ind];
  1256. }
  1257. }
  1258. //fTracks->Compress();
  1259. delete [] index;
  1260. // Remove ghost tracks (with too many shared hits)
  1261. nFound = fTracks->GetEntriesFast();
  1262. Int_t *nh = new Int_t [nFound];
  1263. Double_t *etas = new Double_t [nFound];
  1264. index = new Int_t [nFound];
  1265. for (Int_t i = 0; i < nFound; ++i) {
  1266. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1267. if (track->GetTofIndex() < 0) continue; // track from TPC
  1268. nh[i] = track->GetNofHits();
  1269. etas[i] = TMath::Abs (track->Momentum3().Eta());
  1270. }
  1271. TMath::Sort(nFound,nh,index,kFALSE); // in ascending order
  1272. //TMath::Sort(nFound,etas,index,kTRUE); // in descending order
  1273. for (Int_t i = 0; i < nFound; ++i) {
  1274. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(index[i]);
  1275. if (track->GetTofIndex() < 0) continue; // track from TPC
  1276. Int_t nHits = track->GetNofHits();
  1277. TObjArray *hits = track->GetHits();
  1278. Int_t nShared = 0;
  1279. for (Int_t j = 0; j < nHits; ++j) {
  1280. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1281. if (hit->GetFlag()%1000 >= 10) ++nShared;
  1282. }
  1283. if (1.*nShared / nHits >= 0.5) {
  1284. // Remove track
  1285. // Update hit sharing
  1286. for (Int_t j = 0; j < nHits; ++j) {
  1287. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1288. if (hit->GetFlag()%1000 >= 10) hit->SetFlag(hit->GetFlag()-10);
  1289. }
  1290. fTracks->RemoveAt(index[i]);
  1291. cout << " *** Removing track (too many shared hits): " << i << " "
  1292. << track->GetTrackID() << " " << nHits << " " << nShared << endl;
  1293. }
  1294. }
  1295. fTracks->Compress();
  1296. // Apply rejection based on the number of missing hits
  1297. nFound = fTracks->GetEntriesFast();
  1298. for (Int_t i = 0; i < nFound; ++i) {
  1299. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1300. if (track->GetTofIndex() < 0) continue; // track from TPC
  1301. Int_t nHits = track->GetNofHits();
  1302. if (track->GetMisses() >= nHits) {
  1303. fTracks->RemoveAt(i);
  1304. cout << " *** Removing track (too many misses): " << i << " "
  1305. << track->GetTrackID() << " " << nHits << " " << track->GetMisses() << endl;
  1306. }
  1307. }
  1308. fTracks->Compress();
  1309. delete [] nh;
  1310. delete [] etas;
  1311. delete [] index;
  1312. }
  1313. //__________________________________________________________________________
  1314. void MpdEctTrackFinderTof::MatchTpc()
  1315. {
  1316. /// Match ETOF-ECT track with tracks from TPC.
  1317. Int_t nTpc = fTpcTracks->GetEntriesFast();//, nMatch = 0;
  1318. map<Int_t,pair<Double_t,Double_t> > tpcM, ectM;
  1319. set<Int_t> tpcEctS;
  1320. cout << " TPC tracks: " << nTpc << endl;
  1321. for (Int_t i = 0; i < nTpc; ++i) {
  1322. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(i);
  1323. // Exclude TPC tracks used for ECT tracking
  1324. //if (tr == 0x0 || tr->GetUniqueID() == 1) continue;
  1325. if (tr == 0x0) continue;
  1326. if (tr->GetUniqueID() == 1) tpcEctS.insert(i);
  1327. //if (tr->GetParam(3) < 0) continue; // !!! going in backward hemisphere
  1328. if (TMath::Abs(tr->GetParam(3)) < 0.3 ||
  1329. tr->Momentum() < 0.07) continue; // !!! track in central region or with low momentum
  1330. MpdEctKalmanTrack tpc = MpdEctKalmanTrack(i, *tr);
  1331. /*2-may-12
  1332. Double_t distFirst = ((MpdKalmanHit*) tpc.GetHits()->First())->GetDist();
  1333. Double_t distLast = ((MpdKalmanHit*) tpc.GetHits()->Last())->GetDist();
  1334. tpc.SetPos(TMath::Min(distFirst,distLast));
  1335. tpc.SetPosNew(tpc.GetPos());
  1336. Bool_t ok = MpdKalmanFilter::Instance()->Refit(&tpc, -1); // from last to first hit
  1337. //cout << i << " " << tpc.GetHits()->GetEntriesFast() << " " << tpc.GetParamNew(3) << " " << tpc.GetParamNew(1) << " " << tpc.GetPosNew() << endl;
  1338. //if (!ok || tpc.GetParamNew(3) < 0) {
  1339. if (!ok) {
  1340. // propagation failure or going in backward hemisphere
  1341. //fTrackCand->RemoveLast();
  1342. //--nCand;
  1343. continue;
  1344. }
  1345. */
  1346. //cout << i << " " << tpc.GetHits()->GetEntriesFast() << " " << tpc.GetType() << " " << tpc.GetParamNew(1)
  1347. // << " " << tpc.GetParamNew(3) << " " << tpc.GetPosNew() << endl;
  1348. // Propagate to TPC end-plate
  1349. MpdKalmanHit hit;
  1350. hit.SetType(MpdKalmanHit::kFixedZ);
  1351. Double_t zEnd = 150.; // get it from geometry !!!
  1352. hit.SetPos(TMath::Sign(zEnd,tpc.GetParam(3)));
  1353. MpdKalmanFilter::Instance()->PropagateToHit(&tpc, &hit, kFALSE);
  1354. Double_t rad = tpc.GetParamNew(0) * tpc.GetParamNew(0) + tpc.GetParamNew(1) * tpc.GetParamNew(1);
  1355. rad = TMath::Sqrt(rad);
  1356. //cout << tpc.GetTrackID() << " " << tpc.GetType() << " " << tpc.GetPosNew() << " " << tpc.GetParamNew(3) << " " << rad << endl;
  1357. if (rad > 100.) continue; // outer radius of TPC
  1358. tpcM.insert(pair<Int_t,pair<Double_t,Double_t> > (i,pair<Double_t,Double_t>(tpc.GetParamNew(0),tpc.GetParamNew(1))));
  1359. }
  1360. Int_t nFound = fTracks->GetEntriesFast();
  1361. cout << " ECT tracks: " << nFound << endl;
  1362. for (Int_t i = 0; i < nFound; ++i) {
  1363. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1364. if (tr->IsFromTpc()) continue; // track from TPC
  1365. // Propagate to TPC end-plate
  1366. MpdKalmanHit hit;
  1367. hit.SetType(MpdKalmanHit::kFixedZ);
  1368. Double_t zEnd = 150.; // get it from geometry !!!
  1369. //hit.SetPos(TMath::Sign(zEnd,tr->GetPosNew()));
  1370. hit.SetPos(TMath::Sign(zEnd,tr->GetParam(3)));
  1371. MpdKalmanFilter::Instance()->PropagateToHit(tr, &hit, kTRUE);
  1372. Double_t rad = tr->GetParamNew(0) * tr->GetParamNew(0) + tr->GetParamNew(1) * tr->GetParamNew(1);
  1373. rad = TMath::Sqrt(rad);
  1374. //cout << tr->GetTrackID() << " " << tr->GetPosNew() << " " << tr->GetParamNew(3) << " " << rad << endl;
  1375. //if (rad < 34.) continue; // inner radius of TPC
  1376. Double_t thick = 0.06; // TPC end-plate material thickness in rad. lengths
  1377. PassWall(tr, thick);
  1378. ectM.insert(pair<Int_t,pair<Double_t,Double_t> > (i,pair<Double_t,Double_t>(tr->GetParamNew(0),tr->GetParamNew(1))));
  1379. }
  1380. map<Int_t,pair<Double_t,Double_t> >::iterator it, it1;
  1381. multimap<Int_t,Int_t> matchM;
  1382. Double_t rMax = 0.5; // 5 mm matching radius
  1383. for (it = ectM.begin(); it != ectM.end(); ++it) {
  1384. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fTracks->UncheckedAt(it->first);
  1385. Int_t id = tr->GetTrackID(), dir = Int_t(tr->GetPosNew());
  1386. for (it1 = tpcM.begin(); it1 != tpcM.end(); ++it1) {
  1387. MpdTpcKalmanTrack *tr1 = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(it1->first);
  1388. Int_t id1 = tr1->GetTrackID(), dir1 = Int_t(tr1->GetParam(3));
  1389. if (dir * dir1 < 0) continue; // different arms
  1390. Double_t dx = it->second.first - it1->second.first;
  1391. Double_t dy = it->second.second - it1->second.second;
  1392. Double_t dr = dx * dx + dy * dy;
  1393. dr = TMath::Sqrt (dr);
  1394. if (id == id1) {
  1395. // Exact ID matching
  1396. Double_t r = it->second.first * it->second.first + it->second.second * it->second.second;
  1397. if (lun2) fprintf(lun2,"%6d %10.3e %10.3e %1d %3d %3d\n", id, TMath::Sqrt(r), dr, tr1->GetUniqueID(),
  1398. tr->GetNofHits(), tr1->GetNofHits());
  1399. }
  1400. if (dr < rMax) matchM.insert(pair<Int_t,Int_t>(it->first,it1->first));
  1401. }
  1402. }
  1403. cout << " Number of ETOF tracks to TPC: " << ectM.size() << ", potential matches: " << matchM.size() << endl;
  1404. // Find matches
  1405. multimap<Int_t,Int_t>::iterator itm, itm1;
  1406. pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
  1407. map<Int_t,MpdEctKalmanTrack> trackM;
  1408. TMatrixDSym pointWeight(5);
  1409. TMatrixD param(5,1);
  1410. Int_t indx0 = -1;
  1411. for (itm = matchM.begin(); itm != matchM.end(); ++itm) {
  1412. Int_t indx = itm->first;
  1413. if (indx == indx0) continue;
  1414. indx0 = indx;
  1415. MpdEctKalmanTrack *tr = (MpdEctKalmanTrack*) fTracks->UncheckedAt(indx);
  1416. MpdEctKalmanTrack ect(*tr);
  1417. ret = matchM.equal_range(indx);
  1418. trackM.clear();
  1419. for (itm1 = ret.first; itm1 != ret.second; ++itm1) {
  1420. MpdTpcKalmanTrack *tr1 = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(itm1->second);
  1421. TClonesArray *hits = tr1->GetTrHits();
  1422. Int_t nHits = hits->GetEntriesFast(), lay = -1, lay0 = -1;
  1423. MpdEctKalmanTrack trTMP = ect;
  1424. trTMP.GetHits()->Clear();
  1425. // Follow track thru TPC
  1426. map<Int_t,Double_t>& stepMap = tr1->GetSteps();
  1427. for (Int_t j = 0; j < nHits; ++j) {
  1428. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1429. lay0 = lay;
  1430. lay = hit->GetLayer();
  1431. Double_t stepBack = -1.0;
  1432. if (stepMap.find(lay) != stepMap.end() && stepMap.find(lay0) != stepMap.end())
  1433. stepBack = TMath::Abs (stepMap[lay] - stepMap[lay0]);
  1434. if (!MpdKalmanFilter::Instance()->PropagateToHit(&trTMP,hit,kTRUE,kTRUE,stepBack)) {
  1435. trTMP.SetNodeNew(trTMP.GetNode()); // restore node
  1436. continue;
  1437. }
  1438. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(&trTMP,hit,pointWeight,param);
  1439. //cout << j << " " << dChi2 << " " << hit->GetPos() << endl;
  1440. if (dChi2 > fgkChi2Cut * 4) continue;
  1441. trTMP.SetParam(param);
  1442. trTMP.SetParamNew(param);
  1443. (*trTMP.GetWeight()) += pointWeight;
  1444. trTMP.GetHits()->Add(hit);
  1445. trTMP.SetChi2(trTMP.GetChi2()+dChi2);
  1446. // Set params at hit
  1447. trTMP.SetParamAtHit(*trTMP.GetParamNew());
  1448. trTMP.SetWeightAtHit(*trTMP.GetWeight());
  1449. trTMP.SetLengAtHit(trTMP.GetLength());
  1450. trTMP.SetPosAtHit(trTMP.GetPosNew());
  1451. }
  1452. trTMP.SetPartID(tr1->GetPartID());
  1453. trackM.insert(pair<Int_t,MpdEctKalmanTrack>(itm1->second,trTMP));
  1454. }
  1455. // Select the best track (with max number of hits)
  1456. map<Int_t,MpdEctKalmanTrack>::iterator iter;
  1457. Int_t nHitsMax = 0, iMax = -1;
  1458. for (iter = trackM.begin(); iter != trackM.end(); ++iter) {
  1459. Int_t nhits = iter->second.GetNofHits();
  1460. if (nhits == 0) continue;
  1461. if (nhits > nHitsMax) {
  1462. nHitsMax = nhits;
  1463. iMax = iter->first;
  1464. } else if (nhits == nHitsMax) {
  1465. if (iter->second.GetChi2() < trackM[iMax].GetChi2()) {
  1466. iMax = iter->first;
  1467. nHitsMax = nhits;
  1468. }
  1469. }
  1470. }
  1471. if (iMax >= 0) {
  1472. cout << " ECT-TPC match found: " << tr->GetTrackID() << " " << trackM[iMax].GetTrackID() << " " << iMax << endl;
  1473. trackM[iMax].SetMatchTpc();
  1474. trackM[iMax].SetTpcIndex(iMax);
  1475. // Check if there is a TPC-ECT track (clone)
  1476. Bool_t removeThis = kFALSE;
  1477. if (tpcEctS.find(iMax) != tpcEctS.end()) {
  1478. for (Int_t j = 0; j < nFound; ++j) {
  1479. MpdEctKalmanTrack *tr1 = (MpdEctKalmanTrack*) fTracks->UncheckedAt(j);
  1480. if (tr1 == 0x0) continue;
  1481. if (!tr1->IsFromTpc()) continue; // track not from TPC
  1482. if (tr1->GetTpcIndex() != iMax) continue;
  1483. cout << tr->GetNofHits() << " " << tr1->GetNofHits() << " " << tr1->GetChi2() << endl;
  1484. if (tr->GetNofHits() > tr1->GetNofHits()) fTracks->RemoveAt(j);
  1485. else if (tr->GetNofHits() < tr1->GetNofHits()) removeThis = kTRUE;
  1486. else if (tr->GetChi2() > tr1->GetChi2()) removeThis = kTRUE;
  1487. break;
  1488. }
  1489. }
  1490. cout << tr->GetChi2() << " " << tr->GetNofTrHits() << endl;
  1491. fTracks->RemoveAt(indx);
  1492. if (!removeThis) {
  1493. // Add TPC hits
  1494. TClonesArray *oldHits = trackM[iMax].GetTrHits();
  1495. TObjArray *newHits = trackM[iMax].GetHits();
  1496. Int_t nAdd = newHits->GetEntriesFast(), nOld = oldHits->GetEntriesFast();
  1497. MpdEctKalmanTrack *newTrack = new ((*fTracks)[itm->first]) MpdEctKalmanTrack(trackM[iMax]);
  1498. oldHits = newTrack->GetTrHits();
  1499. for (Int_t ih = 0; ih < nAdd; ++ih) {
  1500. new ((*oldHits)[nOld++]) MpdKalmanHit (*((MpdKalmanHit*)newHits->UncheckedAt(ih)));
  1501. }
  1502. cout << trackM[iMax].GetChi2() << " " << trackM[iMax].GetNofHits() << " " << trackM[iMax].GetNofTrHits() << endl;
  1503. newHits = newTrack->GetHits();
  1504. newHits->Clear();
  1505. for (Int_t ih = 0; ih < nOld; ++ih) newHits->Add(oldHits->UncheckedAt(ih));
  1506. newTrack->SetNofHits(newTrack->GetNofTrHits());
  1507. newTrack->SetParamAtHit(*trackM[iMax].GetParamAtHit());
  1508. newTrack->SetWeightAtHit(*trackM[iMax].GetWeightAtHit());
  1509. newTrack->SetLengAtHit(trackM[iMax].GetLengAtHit());
  1510. newTrack->SetPosAtHit(trackM[iMax].GetPosAtHit());
  1511. }
  1512. //fTracks->RemoveAt(indx);
  1513. }
  1514. } // for (itm = matchM.begin(); ...
  1515. trackM.clear();
  1516. fTracks->Compress();
  1517. }
  1518. //__________________________________________________________________________
  1519. void MpdEctTrackFinderTof::PassWall(MpdEctKalmanTrack *track, Double_t thick)
  1520. {
  1521. /// Propagate track thru TPC end plate (include MS)
  1522. // Add multiple scattering
  1523. TMatrixDSym *cov = track->Weight2Cov();
  1524. Double_t th = track->GetParamNew(3);
  1525. Double_t cosTh = TMath::Cos(th);
  1526. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, thick);
  1527. (*cov)(2,2) += (angle2 / cosTh / cosTh );
  1528. (*cov)(3,3) += angle2;
  1529. Int_t iok = 0;
  1530. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1531. track->SetWeight(*cov);
  1532. }
  1533. //__________________________________________________________________________
  1534. void MpdEctTrackFinderTof::GoToBeamLine()
  1535. {
  1536. /// Propagate tracks to the beam line (mult. scat. and energy loss corrections)
  1537. const Int_t nR = 7;
  1538. const Double_t rad[nR] = {27.035, 27.57, 28.105, 30.64, 33.15, 33.665, 34.178};
  1539. const Double_t dx[nR] = {0.0031, 0.00085, 0.0031, 0.0003, 0.001, 0.00085, 0.001};
  1540. Bool_t ok = 0;
  1541. Int_t nReco = fTracks->GetEntriesFast();
  1542. for (Int_t i = 0; i < nReco; ++i) {
  1543. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1544. //if (track->IsFromTpc()) continue; // track from TPC
  1545. if (track->IsFromTpc()) MergeWithTpc(track); // merge track with TPC track
  1546. // Check track radial position
  1547. Double_t rTr = track->GetPosNew();
  1548. if (track->GetType() == MpdKalmanTrack::kEndcap) {
  1549. rTr = track->GetParamNew(0) * track->GetParamNew(0) + track->GetParamNew(1) * track->GetParamNew(1);
  1550. rTr = TMath::Sqrt (rTr);
  1551. }
  1552. MpdKalmanHit hit;
  1553. hit.SetDetectorID(-999);
  1554. hit.SetType(MpdKalmanHit::kFixedR);
  1555. // Check if goes inward
  1556. MpdEctKalmanTrack trTmp = *track;
  1557. hit.SetPos(0.0);
  1558. Double_t leng = trTmp.GetLength();
  1559. //ok = MpdKalmanFilter::Instance()->PropagateToHit(&trTmp,&hit,kFALSE);
  1560. ok = MpdKalmanFilter::Instance()->PropagateToHit(&trTmp,&hit,kTRUE);
  1561. if (!ok || TMath::Abs(trTmp.GetParamNew(1)) > 149.9) {
  1562. fTracks->RemoveAt(i);
  1563. cout << " *** Removing track (not going inward): " << i << " "
  1564. << track->GetTrackID() << " " << endl;
  1565. continue;
  1566. } else if (trTmp.GetPosNew() >= rad[nR-1] || rTr < rad[0]) {
  1567. // Too large DCA or inside TPC inner shell
  1568. *track = trTmp;
  1569. track->SetPos(trTmp.GetPosNew());
  1570. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  1571. track->Weight2Cov();
  1572. if (track->IsFromTpc() || track->IsMatchTpc()) track->SetLengAtHit(trTmp.GetLength()-leng);
  1573. else track->SetLengAtHit(trTmp.GetLength()-trTmp.GetLengAtHit());
  1574. track->SetLength(leng);
  1575. continue;
  1576. }
  1577. if (track->IsFromTpc() || track->IsMatchTpc()) track->SetLengAtHit(trTmp.GetLength()-leng);
  1578. else track->SetLengAtHit(trTmp.GetLength()-trTmp.GetLengAtHit());
  1579. for (Int_t j = nR-1; j >= 0; --j) {
  1580. if (rTr < rad[j]) continue;
  1581. hit.SetPos(rad[j]);
  1582. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(rad[j],0.,0.),kTRUE); // put hit in GeoScheme
  1583. //ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kTRUE);
  1584. ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,!track->IsFromTpc());
  1585. //cout << ok << " " << track->GetPos() << endl;
  1586. if (!ok) {
  1587. fTracks->RemoveAt(i);
  1588. cout << " *** Removing track (not passing TPC inner shell): " << i << " "
  1589. << track->GetTrackID() << " " << endl;
  1590. break;
  1591. }
  1592. //continue;
  1593. if (TMath::Abs(track->GetParamNew(1)) > 150.) continue; // outside TPC Z-dimensions
  1594. // Add multiple scattering
  1595. TMatrixDSym *cov = track->Weight2Cov();
  1596. Double_t th = track->GetParamNew(3);
  1597. Double_t cosTh = TMath::Cos(th);
  1598. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[j]);
  1599. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1600. (*cov)(3,3) += angle2;
  1601. //cov->Print();
  1602. Int_t iok = 0;
  1603. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1604. track->SetWeight(*cov);
  1605. }
  1606. if (!ok) continue;
  1607. /* TMatrixDSym *cov = */ track->Weight2Cov();
  1608. // Update track
  1609. track->SetParam(*track->GetParamNew());
  1610. track->SetPos(track->GetPosNew());
  1611. // Correct for energy loss
  1612. //*
  1613. Double_t theta = TMath::PiOver2() - track->GetParam(3);
  1614. Double_t pt = 1. / track->GetParam(4);
  1615. Double_t ptCor = CorrectForLoss(TMath::Abs(pt), theta, track->GetTrackID()); // ionization loss correction
  1616. track->SetParam(4,TMath::Sign(1./ptCor,pt));
  1617. track->SetParamNew(*track->GetParam());
  1618. //*/
  1619. //*
  1620. hit.SetPos(0.);
  1621. hit.SetMeas(0,track->GetParam(2));
  1622. //hit.SetRphi(track->GetParam(2)); // track Phi - check if it is necessary !!!!!!!
  1623. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(0.,0.,0.),kTRUE); // put hit in GeoScheme
  1624. //*/
  1625. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  1626. Double_t pos = track->GetPosNew();
  1627. //Double_t pos = ((MpdKalmanHitR*)track->GetHits()->Last())->GetR(); // position at last hit
  1628. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  1629. //cout << track->Momentum3().Eta() << endl;
  1630. //track->GetCovariance()->Print();
  1631. //(track->GetWeightAtHit()->Invert()).Print();
  1632. //MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  1633. MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, !track->IsFromTpc());
  1634. track->SetPos(pos); // save position after passing inner shell
  1635. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  1636. }
  1637. fTracks->Compress();
  1638. }
  1639. //__________________________________________________________________________
  1640. Double_t MpdEctTrackFinderTof::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
  1641. {
  1642. /// Apply momentum correction due to ionization loss in the beam pipe (0.5 mmof Al)
  1643. /// and TPC inner shell
  1644. //return pt; ///
  1645. const Int_t nDim = 13;
  1646. //const Double_t mom[nDim]={0.075,0.100,0.200,0.300,0.400,0.500,0.700,0.900,
  1647. //1.200,1.500,2.000,2.500,3.500};
  1648. const Double_t tkin[nDim]={0.019,0.032,0.104,0.191,0.284,0.380,0.574,0.771,
  1649. 1.069,1.367,1.865,2.364,3.363};
  1650. //const Double_t dedxp[nDim]={3.125,1.950,0.975,0.800,0.800,0.750,0.750,0.750,
  1651. //0.750,0.750,0.750,0.775,0.775}; // peak
  1652. const Double_t dedxm[nDim]={3.195,2.085,1.100,0.952,0.918,0.898,0.898,0.890,
  1653. 0.887,0.886,0.901,0.904,0.913}; // mean
  1654. const Double_t piMass = 0.13957, eMass = 0.00051;
  1655. Double_t dt, mass;
  1656. Int_t elec = GetParticleId(id);
  1657. if (elec) mass = eMass;
  1658. else mass = piMass;
  1659. Double_t mass2 = mass * mass;
  1660. Double_t p = pt / TMath::Sin(the);
  1661. Double_t t = TMath::Sqrt (p*p + mass2) - mass;
  1662. if (elec) dt = 1.45; // 1.45 MeV loss for electrons
  1663. else {
  1664. //if (t < tkin[0]) dt = dedxm[0];
  1665. if (t < tkin[0]) {
  1666. dt = dedxm[0] + (dedxm[1]-dedxm[0])/(tkin[1]-tkin[0]) * (t-tkin[0]);
  1667. dt = TMath::Min (dt,10.);
  1668. }
  1669. else if (t > tkin[nDim-1]) dt = dedxm[nDim-1];
  1670. else {
  1671. // Binary search
  1672. Int_t i1 = 0, i2 = nDim-1, i = i2;
  1673. do {
  1674. i = i1 + (i2-i1) / 2;
  1675. if (t > tkin[i]) i1 = i;
  1676. else i2 = i;
  1677. //cout << i << " " << i1 << " " << i2 << " " << t << " " << tkin[i1] << " " << tkin[i2] << endl;
  1678. } while (i2 - i1 > 1);
  1679. // Linear interpolation
  1680. dt = dedxm[i1] + (dedxm[i2]-dedxm[i1])/(tkin[i2]-tkin[i1]) * (t-tkin[i1]);
  1681. }
  1682. }
  1683. dt /= TMath::Sin(the);
  1684. t += dt * 1.e-3;
  1685. Double_t e = mass + t;
  1686. p = TMath::Sqrt (e*e - mass2);
  1687. pt = p * TMath::Sin(the);
  1688. return pt;
  1689. }
  1690. //__________________________________________________________________________
  1691. Int_t MpdEctTrackFinderTof::GetParticleId(Int_t id)
  1692. {
  1693. /// Particle ID:
  1694. /// !!! based on MC information at the moment !!!
  1695. MpdMCTrack* mcTr = (MpdMCTrack*) fMCTracks->UncheckedAt(id);
  1696. Int_t pdg = mcTr->GetPdgCode();
  1697. if (TMath::Abs(pdg) == 11) return 1;
  1698. return 0;
  1699. }
  1700. //__________________________________________________________________________
  1701. void MpdEctTrackFinderTof::MergeWithTpc(MpdEctKalmanTrack *track)
  1702. {
  1703. /// Merge track with TPC track
  1704. track->SetDirection(MpdKalmanTrack::kInward);
  1705. cout << track->GetTrackID() << " " << track->GetNofHits() << " " << track->GetChi2() << " " << track->Pt()*track->Charge() << " " ;//<< endl;
  1706. //Bool_t ok = MpdKalmanFilter::Instance()->Refit(track, 1); // from first to last hit
  1707. //*
  1708. track->GetHits()->Sort();
  1709. track->SetChi2(0.);
  1710. TMatrixDSym *w = track->GetWeight();
  1711. Int_t nHits = track->GetNofHits();
  1712. for (Int_t i = 0; i < 5; ++i) {
  1713. for (Int_t j = i; j < 5; ++j) {
  1714. if (j == i) (*w)(i,j) /= nHits;
  1715. else (*w)(i,j) = (*w)(j,i) = 0.;
  1716. }
  1717. }
  1718. TMatrixD param(5,1);
  1719. TMatrixDSym weight(5), pointWeight(5);
  1720. Int_t ibeg = 0, iend = nHits, nAcc = 0;
  1721. MpdKalmanHit *hit = 0x0;
  1722. for (Int_t i = ibeg; i != iend; ++i) {
  1723. hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
  1724. if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kFALSE)) continue;
  1725. ++nAcc;
  1726. PassWall(track, 0.001);
  1727. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hit, pointWeight, param);
  1728. track->SetChi2(track->GetChi2()+dChi2);
  1729. weight = *track->GetWeight();
  1730. weight += pointWeight;
  1731. track->SetWeight(weight);
  1732. track->SetParamNew(param);
  1733. //cout << i << " " << dChi2 << " " << 1./track->GetParamNew(4) << endl;
  1734. // Set params at hit 9-may-12
  1735. track->SetParamAtHit(*track->GetParamNew());
  1736. track->SetWeightAtHit(*track->GetWeight());
  1737. track->SetLengAtHit(track->GetLength());
  1738. track->SetPosAtHit(track->GetPosNew());
  1739. }
  1740. //*/
  1741. cout << nAcc << " " << track->GetChi2() << " " << track->Pt()*track->Charge() << endl;
  1742. // Propagate to TPC end-plate
  1743. MpdKalmanHit hitTmp;
  1744. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  1745. Double_t zEnd = 150.; // get it from geometry !!!
  1746. //hitTmp.SetPos(TMath::Sign(zEnd,track->GetPos()));
  1747. hitTmp.SetPos(TMath::Sign(zEnd,track->GetParam(3)));
  1748. MpdKalmanFilter::Instance()->PropagateToHit(track, &hitTmp, kFALSE);
  1749. Double_t thick = 0.06; // material thickness in rad. lengths
  1750. PassWall(track,thick);
  1751. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(track->GetTpcIndex());
  1752. TClonesArray *hits = tr->GetTrHits();
  1753. /*Int_t*/ nHits = hits->GetEntriesFast();
  1754. //TMatrixDSym pointWeight(5);
  1755. pointWeight = 0.;
  1756. //TMatrixD param(5,1);
  1757. TClonesArray *oldHits = track->GetTrHits();
  1758. Int_t nOld = oldHits->GetEntriesFast(), nAdd = 0, lay0 = -1, lay = -1;
  1759. nAcc = nOld;
  1760. map<Int_t,Double_t>& stepMap = track->GetSteps();
  1761. for (Int_t j = 0; j < nHits; ++j) {
  1762. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits->UncheckedAt(j);
  1763. lay0 = lay;
  1764. lay = hit1->GetLayer();
  1765. //if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit,kFALSE)) continue;
  1766. //cout << j << " " << hit->GetMeas(1) << endl;
  1767. Double_t stepBack = -1.0;
  1768. if (stepMap.find(lay) != stepMap.end() && stepMap.find(lay0) != stepMap.end())
  1769. stepBack = TMath::Abs (stepMap[lay] - stepMap[lay0]);
  1770. if (!MpdKalmanFilter::Instance()->PropagateToHit(track,hit1,kFALSE,kTRUE,stepBack)) {
  1771. track->SetNodeNew(track->GetNode()); // restore node
  1772. continue;
  1773. }
  1774. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track,hit1,pointWeight,param);
  1775. //cout << j << " " << hit->GetDist() << " " << dChi2 << endl;
  1776. if (dChi2 > fgkChi2Cut * 4) continue;
  1777. //track->SetParam(param);
  1778. track->SetParamNew(param);
  1779. (*track->GetWeight()) += pointWeight;
  1780. track->GetHits()->Add(hit1);
  1781. track->SetChi2(track->GetChi2()+dChi2);
  1782. new ((*oldHits)[nOld++]) MpdKalmanHit (*hit1);
  1783. ++nAdd;
  1784. // Set params at hit
  1785. track->SetParamAtHit(*track->GetParamNew());
  1786. track->SetWeightAtHit(*track->GetWeight());
  1787. track->SetLengAtHit(track->GetLength());
  1788. track->SetPosAtHit(track->GetPosNew());
  1789. }
  1790. if (nAdd == 0) {
  1791. cout << " Merge TPC: merging failed !!! " << endl;
  1792. track->SetFlag(track->GetFlag() ^ 1); // clear "From_TPC" flag
  1793. return;
  1794. }
  1795. /*TClonesArray *oldHits = track->GetTrHits();
  1796. Int_t nAdd = hits->GetEntriesFast(), nOld = oldHits->GetEntriesFast();
  1797. for (Int_t ih = 0; ih < nAdd; ++ih) {
  1798. new ((*oldHits)[nOld++]) MpdKalmanHit (*(MpdKalmanHit*)hits->UncheckedAt(ih));
  1799. }*/
  1800. track->SetNofHits(track->GetNofTrHits());
  1801. cout << " Merge TPC: " << nHits << " " << nAdd << " " << nAcc << " " << track->GetNofHits() << " "
  1802. << track->GetNofTrHits() << " " << tr->GetTrackID() << " " << track->GetTrackID() << " "
  1803. << tr->GetParam(3) << " " << track->GetChi2() << endl;
  1804. }
  1805. //__________________________________________________________________________
  1806. void MpdEctTrackFinderTof::Smooth()
  1807. {
  1808. /// Primary vertex constraints: smooth tracks from primary vertex (update
  1809. /// momentum and track length - covariance matrix is not updated !!!)
  1810. const Double_t cutChi2 = 15.;
  1811. MpdVertex *vtx = 0x0, *vtxOld = 0x0;
  1812. Int_t nPart = 0;
  1813. if (fPrimVtx) {
  1814. vtx = (MpdVertex*) fPrimVtx->UncheckedAt(0);
  1815. nPart = vtx->GetNTracks();
  1816. // Save original vertex
  1817. vtxOld = new MpdVertex();
  1818. *vtxOld = *vtx;
  1819. }
  1820. if (vtx == 0x0 || nPart < 20) {
  1821. // Use default vertex
  1822. TMatrixFSym cov(3);
  1823. cov(0,0) = cov(1,1) = cov(2,2) = 1.e-6;
  1824. MpdVertex vertDef("Def. vertex", "DefVertex", 0, 0, 0, 0, 0, 0, cov);
  1825. if (vtx == 0x0) vtx = new MpdVertex();
  1826. *vtx = vertDef;
  1827. if (!fPrimVtx) {
  1828. // No vertex
  1829. /*
  1830. TClonesArray tmp("MpdVertex",1);
  1831. fPrimVtx = &tmp;
  1832. MpdVertex *v = new (tmp[0]) MpdVertex();
  1833. */
  1834. TClonesArray *tmp = new TClonesArray("MpdVertex",1);
  1835. fPrimVtx = tmp;
  1836. MpdVertex *v = new ((*tmp)[0]) MpdVertex();
  1837. *v = *vtx;
  1838. }
  1839. }
  1840. MpdKfPrimaryVertexFinder vFinder;
  1841. vFinder.SetVertices(fPrimVtx);
  1842. vFinder.SetTracks(fTracks);
  1843. vFinder.Chi2Vertex();
  1844. // Select tracks close to the vertex
  1845. Int_t nTr = fTracks->GetEntriesFast(), nOk = 0;
  1846. TArrayI indsNew(nTr); //, indsOld = *(vtx->GetIndices());
  1847. Double_t *lengs = new Double_t [nTr];
  1848. for (Int_t itr = 0; itr < nTr; ++itr) {
  1849. MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(itr);
  1850. lengs[itr] = track->GetLength();
  1851. if (track->GetChi2Vertex() > cutChi2) continue;
  1852. indsNew[nOk++] = itr;
  1853. }
  1854. indsNew.Set(nOk);
  1855. vtx->SetIndices(&indsNew);
  1856. vFinder.Smooth();
  1857. // Set lower limit to the track length (distance to ETOF hit)
  1858. for (Int_t itr = 0; itr < nTr; ++itr) {
  1859. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(itr);
  1860. if (track->GetChi2Vertex() < cutChi2) track->SetLengAtHit(track->GetLengAtHit() + track->GetLength() - lengs[itr]);
  1861. //if (track->GetChi2Vertex() > cutChi2) continue;
  1862. //if (track->GetTofIndex() < 0) continue;
  1863. Double_t dist = TMath::Abs (((MpdTofHit*) fTofHits->First())->GetZ());
  1864. if (track->GetTofIndex() < 0) {
  1865. if (track->GetChi2Vertex() > cutChi2) continue;
  1866. } else {
  1867. MpdTofHit *hit = (MpdTofHit*) fTofHits->UncheckedAt(track->GetTofIndex());
  1868. Double_t dz = hit->GetZ() - vtx->GetZ();
  1869. dist = hit->GetX() * hit->GetX() + hit->GetY() * hit->GetY() + dz * dz;
  1870. dist = TMath::Sqrt (dist);
  1871. }
  1872. if (track->GetLength() < dist) track->SetLength(dist);
  1873. }
  1874. delete [] lengs;
  1875. if (vtxOld) {
  1876. // Restore old vertex
  1877. *vtx = *vtxOld;
  1878. delete vtxOld;
  1879. } else {
  1880. fPrimVtx->Delete();
  1881. fPrimVtx = 0x0;
  1882. delete vtx;
  1883. }
  1884. }
  1885. //__________________________________________________________________________
  1886. void MpdEctTrackFinderTof::GoOutward()
  1887. {
  1888. /// Propagate tracks outward (for matching with outer detectors)
  1889. const Double_t cutChi2 = 15.;
  1890. Bool_t ok = 0;
  1891. Int_t nReco = fTracks->GetEntriesFast();
  1892. FairMCEventHeader *mchead = (FairMCEventHeader*) FairRootManager::Instance()->GetObject("MCEventHeader.");
  1893. for (Int_t i = 0; i < nReco; ++i) {
  1894. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fTracks->UncheckedAt(i);
  1895. //if (track->IsFromTpc()) continue; // track from TPC
  1896. MpdEctKalmanTrack tr = *track;
  1897. tr.SetParam(*tr.GetParamAtHit());
  1898. tr.SetParamNew(*tr.GetParamAtHit());
  1899. tr.SetWeight(*tr.GetWeightAtHit());
  1900. tr.SetLength(tr.GetLengAtHit());
  1901. MpdKalmanHit *hit = (MpdKalmanHit*) tr.GetTrHits()->Last();
  1902. //2-may-12 tr.SetPos(hit->GetPos());
  1903. tr.SetPos(tr.GetPosAtHit());
  1904. tr.SetPosNew(tr.GetPos());
  1905. tr.SetDirection(MpdKalmanTrack::kOutward);
  1906. if (hit->GetType() == MpdKalmanHit::kFixedP) {
  1907. tr.SetNode(MpdKalmanFilter::Instance()->GetGeo()->Path(hit->GetDetectorID()));
  1908. tr.SetNodeNew(MpdKalmanFilter::Instance()->GetGeo()->Path(hit->GetDetectorID()));
  1909. }
  1910. //cout << "MpdEctTrackFinderTof::GoOutward(): " << i << " " << tr.GetTrackID() << " " << tr.GetPos() << " "
  1911. // << tr.GetNofHits() << " " << tr.GetNofTrHits() << " " << tr.GetChi2() << endl;
  1912. //if (tr.GetNofHits() != tr.GetNofTrHits()) exit(0);
  1913. // Just for test - from EtofMatching
  1914. /*
  1915. tr.SetPos(track->GetPosNew());
  1916. tr.SetPosNew(track->GetPosNew());
  1917. tr.SetParam(*track->GetParam());
  1918. tr.SetParamNew(*track->GetParam());
  1919. tr.ReSetWeight();
  1920. tr.SetLength(0);
  1921. tr.SetDirection(MpdKalmanTrack::kOutward);
  1922. MpdKalmanFilter::Instance()->Refit(&tr,-1,kTRUE);
  1923. tr.SetPos(tr.GetPosNew());
  1924. tr.SetParam(*tr.GetParamNew());
  1925. */
  1926. // Check if there are hits from TPC
  1927. TClonesArray *trhits = tr.GetTrHits();
  1928. trhits->Sort();
  1929. Int_t nTpc = 0, nEct = 0, nHits = trhits->GetEntriesFast();
  1930. for (Int_t jh = 0; jh < nHits; ++jh) {
  1931. MpdKalmanHit *hit1 = (MpdKalmanHit*) trhits->UncheckedAt(jh);
  1932. if (hit1->GetType() == MpdKalmanHit::kFixedR || hit1->GetType() == MpdKalmanHit::kFixedP) ++nTpc;
  1933. else ++nEct;
  1934. }
  1935. //*
  1936. TObjArray *hits = tr.GetHits();
  1937. Bool_t ok1 = 0;
  1938. if (nTpc > 0) {
  1939. // TPC
  1940. hits->Clear();
  1941. for (Int_t jh = nEct; jh < nHits; ++jh) hits->Add(trhits->UncheckedAt(jh));
  1942. //ok = MpdKalmanFilter::Instance()->Refit(&tr,-1,kFALSE);
  1943. ok1 = MpdKalmanFilter::Instance()->Refit(&tr,-1,kTRUE);
  1944. if (!ok1) cout << " Go outward: refit failed !!!" << endl;
  1945. // Go to TPC end-plate
  1946. MpdKalmanHit hitTmp;
  1947. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  1948. Double_t zEnd = 150.; // get it from geometry !!!
  1949. hitTmp.SetPos(TMath::Sign(zEnd,track->GetParam(3)));
  1950. //MpdKalmanFilter::Instance()->PropagateToHit(&tr, &hitTmp, kFALSE);
  1951. MpdKalmanFilter::Instance()->PropagateToHit(&tr, &hitTmp, kTRUE);
  1952. Double_t thick = 0.06; // material thickness in rad. lengths
  1953. PassWall(&tr,thick);
  1954. }
  1955. if (nEct > 0) {
  1956. // ECT
  1957. hits->Clear();
  1958. for (Int_t jh = 0; jh < nEct; ++jh) hits->Add(trhits->UncheckedAt(jh));
  1959. if (nTpc == 0) {
  1960. tr.SetType(MpdKalmanTrack::kEndcap);
  1961. tr.SetChi2(0.);
  1962. //ok = MpdKalmanFilter::Instance()->Refit(&tr,-1,kTRUE);
  1963. TMatrixDSym *w = tr.GetWeight();
  1964. for (Int_t k = 0; k < 5; ++k) {
  1965. for (Int_t l = k; l < 5; ++l) {
  1966. if (l == k) (*w)(k,l) /= nEct;
  1967. else (*w)(k,l) = (*w)(l,k) = 0.;
  1968. }
  1969. }
  1970. }
  1971. TMatrixDSym pointWeight(5);
  1972. TMatrixD param(5,1);
  1973. for (Int_t jh = nEct - 1; jh >= 0; --jh) {
  1974. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits->UncheckedAt(jh);
  1975. //if (!MpdKalmanFilter::Instance()->PropagateToHit(&tr,hit1,kFALSE)) continue;
  1976. if (!MpdKalmanFilter::Instance()->PropagateToHit(&tr,hit1,kTRUE)) continue;
  1977. PassWall(&tr, 0.001);
  1978. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(&tr,hit1,pointWeight,param);
  1979. //tr.SetParam(param);
  1980. tr.SetParamNew(param);
  1981. (*tr.GetWeight()) += pointWeight;
  1982. tr.SetChi2(tr.GetChi2()+dChi2);
  1983. }
  1984. }
  1985. printf("MpdEctTrackFinderTof::GoOutward(): %6d %2d %2d %10.3f %7.3f %1d \n", tr.GetTrackID(),
  1986. nTpc, nEct, tr.GetChi2(), tr.Pt(), ok);
  1987. // Set params at hit
  1988. track->SetParamAtHit(*tr.GetParamNew());
  1989. track->SetWeightAtHit(*tr.GetWeight());
  1990. track->SetLengAtHit(tr.GetLength());
  1991. track->SetPosAtHit(tr.GetPosNew());
  1992. track->GetTrHits()->Sort();
  1993. //*/
  1994. // Below is just for debugging
  1995. if (0) {
  1996. tr.SetParam(*track->GetParamAtHit());
  1997. tr.SetParamNew(*track->GetParamAtHit());
  1998. tr.SetWeight(*track->GetWeightAtHit());
  1999. tr.SetLength(track->GetLengAtHit());
  2000. tr.SetPos(((MpdKalmanHit*)track->GetTrHits()->First())->GetPos());
  2001. tr.SetPosNew(tr.GetPos());
  2002. // Extrapolate to ETOF
  2003. //Double_t zTof = ((MpdTofHit*)fTofHits->First())->GetZ();
  2004. Double_t zTof = ((FairMCPoint*)fTofPoints->First())->GetZ();
  2005. zTof = TMath::Sign (TMath::Abs(zTof),track->GetParam(3));
  2006. MpdKalmanHit hitTmp;
  2007. hitTmp.SetType(MpdKalmanHit::kFixedZ);
  2008. hitTmp.SetPos(zTof);
  2009. hitTmp.SetNofDim(2);
  2010. //MpdKalmanFilter::Instance()->PropagateToHit(&tr, &hitTmp, kFALSE);
  2011. MpdKalmanFilter::Instance()->PropagateToHit(&tr, &hitTmp, kTRUE);
  2012. // Correct min. length
  2013. if (tr.GetChi2Vertex() < cutChi2) {
  2014. Double_t dist = zTof * zTof + tr.GetParamNew(0) * tr.GetParamNew(0);
  2015. dist += (tr.GetParamNew(1) * tr.GetParamNew(1));
  2016. dist = TMath::Sqrt (dist);
  2017. //if (tr.GetLength() < dist) tr.SetLength(dist);
  2018. tr.SetLength(dist);
  2019. }
  2020. // Matching
  2021. /*
  2022. Int_t nTof = fTofHits->GetEntriesFast();
  2023. for (Int_t jh = 0; jh < nTof; ++jh) {
  2024. MpdTofHit *tofHit = (MpdTofHit*) fTofHits->UncheckedAt(jh);
  2025. Int_t nLinks = tofHit->GetNLinks();
  2026. for (Int_t lnk = 0; lnk < nLinks; ++lnk) {
  2027. FairLink link = tofHit->GetLink(lnk);
  2028. if (link.GetType() != MpdTofUtils::IsTofPointIndex) continue;
  2029. FairMCPoint *p = fTofPoints->UncheckedAt(link);
  2030. }
  2031. //hitTmp
  2032. }
  2033. */
  2034. Int_t nTof = fTofPoints->GetEntriesFast();
  2035. TVector3 pos, posTr(tr.GetParamNew(0),tr.GetParamNew(1),tr.GetPosNew());
  2036. for (Int_t jh = 0; jh < nTof; ++jh) {
  2037. FairMCPoint *p = (FairMCPoint*) fTofPoints->UncheckedAt(jh);
  2038. if (p->GetTrackID() != tr.GetTrackID()) continue;
  2039. Double_t dl = tr.GetLength() - p->GetLength();
  2040. p->Position(pos);
  2041. Double_t phi = pos.Phi();
  2042. Double_t r = pos.Pt();
  2043. Double_t phiTr = MpdKalmanFilter::Instance()->Proxim(phi,posTr.Phi());
  2044. if (lun2) fprintf(lun2, "%4d %6d %2d %10.3f %10.3f % 10.3f %10.3f %10.3f %10.3f %10.3f %1d\n", mchead->GetEventID(), tr.GetTrackID(), nTpc,
  2045. dl, tr.GetParamNew(0), tr.GetParamNew(1), tr.GetPosNew(), pos.Pt(),
  2046. posTr.Pt()-r, (phiTr-phi)*r, tr.IsFromTpc());
  2047. break;
  2048. }
  2049. } // if (0)
  2050. } //for (Int_t i = 0; i < nReco;
  2051. }
  2052. ClassImp(MpdEctTrackFinderTof);