MpdEctTrackFinderCpc.cxx 101 KB

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