MpdTpcKalmanFilter.cxx 108 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717271827192720
  1. /// \class MpdTpcKalmanFilter
  2. ///
  3. /// Kalman filter track reconstructor in the MPD central detector
  4. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  5. #include "MpdTpcKalmanFilter.h"
  6. #include "MpdCodeTimer.h"
  7. #include "MpdTpcHit.h"
  8. #include "MpdTpcSectorGeo.h"
  9. #include "MpdTpcKalmanTrack.h"
  10. #include "MpdKalmanFilter.h"
  11. #include "MpdKalmanTrack.h"
  12. #include "MpdKalmanGeoScheme.h"
  13. #include "MpdKalmanHit.h"
  14. #include "MpdTpcDedxTask.h"
  15. #include "MpdVertexZfinder.h"
  16. #include "MpdTpcFoundHit.h"
  17. #include "MpdTpcSectorGeo.h"
  18. #include "TpcCluster.h"
  19. #include "TpcGeoPar.h"
  20. #include "TpcPoint.h"
  21. #include "TpcLheTrack.h"
  22. #include "MpdMCTrack.h"
  23. #include "FairEventHeader.h"
  24. #include "FairField.h"
  25. #include "FairMCEventHeader.h"
  26. #include "FairRootManager.h"
  27. #include "FairRunAna.h"
  28. #include "FairRuntimeDb.h"
  29. #include "FairTask.h"
  30. //#include "MpdTofGeoPar.h" //AZ
  31. #include <TClonesArray.h>
  32. #include <TGeoManager.h>
  33. #include <TGeoMatrix.h>
  34. #include <TGeoPgon.h>
  35. #include <TGeoTrd1.h>
  36. #include <TGeoTube.h>
  37. #include <TRandom.h>
  38. #include <TROOT.h>
  39. #ifdef XROOTD
  40. #include "TXNetFile.h"
  41. #endif
  42. #ifdef _OPENMP
  43. #include "omp.h"
  44. #include "sys/time.h"
  45. #endif
  46. const Double_t MpdTpcKalmanFilter::fgkChi2Cut = 20; //50; //20; //100;
  47. FILE *lunTpc = 0x0; //fopen("dl.dat","w");
  48. //__________________________________________________________________________
  49. MpdTpcKalmanFilter::MpdTpcKalmanFilter()
  50. : FairTask("TPC Kalman filter"),
  51. fNofEvents(0),
  52. fNTracks(0),
  53. fNPass(1),
  54. fHits(0x0),
  55. fKHits(new TClonesArray("MpdKalmanHit")),
  56. fTracks(new TClonesArray("MpdTpcKalmanTrack")),
  57. fTrackCand(new TClonesArray("MpdTpcKalmanTrack")),
  58. fLHEtracks(0x0),
  59. fMCtracks(0x0),
  60. fTpcPoints(0x0),
  61. fLayPointers(0x0),
  62. fhLays(0x0),
  63. fVertZ(0.0),
  64. fZflag(0),
  65. fModular(0),
  66. //fPadPlane(TpcPadPlane::Instance())
  67. //fSecGeo(MpdTpcSectorGeo::Instance())
  68. //fCache(new std::map<Double_t,tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> >)
  69. fCache(new std::map<Double_t,matrix4>)
  70. {
  71. /// Default constructor
  72. fUseMCHit = kTRUE;
  73. }
  74. //__________________________________________________________________________
  75. MpdTpcKalmanFilter::MpdTpcKalmanFilter(const char *name,
  76. const char *title)
  77. //: FairTask(name),
  78. : FairTask("TPC Kalman filter"),
  79. fNofEvents(0),
  80. fNTracks(0),
  81. fNPass(1),
  82. fHits(0x0),
  83. fKHits(new TClonesArray("MpdKalmanHit")),
  84. fTracks(new TClonesArray("MpdTpcKalmanTrack")),
  85. fTrackCand(new TClonesArray("MpdTpcKalmanTrack")),
  86. fLHEtracks(0x0),
  87. fMCtracks(0x0),
  88. fTpcPoints(0x0),
  89. fLayPointers(0x0),
  90. fhLays(0x0),
  91. fVertZ(0.0),
  92. fZflag(0),
  93. fModular(0),
  94. //fPadPlane(TpcPadPlane::Instance())
  95. //fSecGeo(MpdTpcSectorGeo::Instance())
  96. //fCache(new std::map<Double_t,tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> >)
  97. fCache(new std::map<Double_t,matrix4>)
  98. {
  99. /// Constructor
  100. FairTask *dedx = new MpdTpcDedxTask();
  101. Add(dedx);
  102. fUseMCHit = kTRUE;
  103. }
  104. //__________________________________________________________________________
  105. MpdTpcKalmanFilter::~MpdTpcKalmanFilter()
  106. {
  107. /// Destructor
  108. //new delete fHits;
  109. delete fKHits;
  110. delete fTracks;
  111. delete fTrackCand;
  112. delete [] fLayPointers;
  113. delete fhLays;
  114. }
  115. //__________________________________________________________________________
  116. InitStatus MpdTpcKalmanFilter::Init() {
  117. cout << "InitStatus MpdTpcKalmanFilter::Init\n\n";
  118. //fVerbose = 10;
  119. FairRootManager *manager = FairRootManager::Instance();
  120. fSecGeo = MpdTpcSectorGeo::Instance();
  121. fhLays = new TH1F("hLays","TPC layers",150,0,150);
  122. //fHits = (TClonesArray*) manager->GetObject("TpcCluster");
  123. if (fUseMCHit) fHits = (TClonesArray*) manager->GetObject("TpcHit");
  124. //else fHits = (TClonesArray*) manager->GetObject("MpdTpcHit");
  125. else fHits = (TClonesArray*) manager->GetObject("TpcRecPoint");
  126. if ( ! fHits ) {
  127. cout << "-I- "<< GetName() << "::Init: No MpdTpcHit array!" << endl;
  128. return kERROR;
  129. }
  130. //fKHits->SetOwner(kTRUE);
  131. fLHEtracks = (TClonesArray *)manager->GetObject("LheGeantTrack");
  132. fMCtracks = (TClonesArray *)manager->GetObject("MCTrack");
  133. fTpcPoints = (TClonesArray *)manager->GetObject("TpcPoint");
  134. #ifdef GEANT
  135. //fListMCtracks = (TClonesArray *)fManager->GetObject("MCTrack");
  136. #endif
  137. fNTracks = fNofEvents = 0;
  138. fNPass = 2; //2;
  139. Register();
  140. // Get event number offset
  141. /* TFile* f = new TFile(FairRootManager::Instance()->GetInFile()->GetName(),"READ"); // EL
  142. #ifdef XROOTD
  143. if (!f->IsOpen()) f = new TXNetFile(FairRootManager::Instance()->GetInFile()->GetName(),"READ");
  144. #endif
  145. cout << f->GetName() << " " << manager->GetObject("MCEventHeader.") << endl;
  146. TTree *cbmsim = (TTree*) f->FindObjectAny("cbmsim");
  147. FairMCEventHeader *head = 0x0; //new FairMCEventHeader();
  148. cbmsim->SetBranchAddress("MCEventHeader.",&head);
  149. cbmsim->GetEntry(0);
  150. Int_t event = head->GetEventID() - 1;
  151. f->Close();
  152. delete f;
  153. fNofEvents += event;
  154. */
  155. /*TH1F *hZ = */new TH1F ("hZ","Z of vertex",150,-75,75);
  156. /*TH1F *hPt = */new TH1F ("hPt","Pt of cand.",100,0,50);
  157. /*TH1F *hChi2 = */new TH1F ("hChi2","Chi2 of track",100,0,50);
  158. /*TH1F *hNhits = */new TH1F ("hNhits","Points per track",100,0,100);
  159. return kSUCCESS;
  160. }
  161. // ----- Private method SetParContainers -------------------------------
  162. void MpdTpcKalmanFilter::SetParContainers() {
  163. //return;
  164. // Get run and runtime database
  165. FairRunAna* run = FairRunAna::Instance();
  166. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  167. FairRuntimeDb* db = run->GetRuntimeDb();
  168. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  169. // Get TPC geometry parameter container
  170. db->getContainer("TpcGeoPar");
  171. db->getContainer("MpdTofGeoPar");
  172. }
  173. //__________________________________________________________________________
  174. InitStatus MpdTpcKalmanFilter::ReInit()
  175. {
  176. return kSUCCESS;
  177. }
  178. //__________________________________________________________________________
  179. void MpdTpcKalmanFilter::Reset()
  180. {
  181. ///
  182. cout << " MpdTpcKalmanFilter::Reset " << endl;
  183. ProcInfo_t proc;
  184. gSystem->GetProcInfo(&proc);
  185. cout << " User CPU time: " << proc.fCpuUser << " Memory: resident " << proc.fMemResident << ", virtual " << proc.fMemVirtual << endl;
  186. fhLays->Reset();
  187. fKHits->Delete();
  188. fTracks->Delete();
  189. fTrackCand->Delete();
  190. delete [] fLayPointers;
  191. fLayPointers = 0x0;
  192. for (Int_t i = 0; i < fgkLays; ++i)
  193. for (Int_t j = 0; j < fgkSecs; ++j) fLaySecBegPointers[i][j] = fLaySecEndPointers[i][j] = -1;
  194. fNTracks = 0;
  195. }
  196. //__________________________________________________________________________
  197. void MpdTpcKalmanFilter::Register()
  198. {
  199. ///
  200. FairTask *itsHF = (FairTask*) FairRun::Instance()->GetTask("Ideal STS hit Producer");
  201. TClonesArray *itsP = (TClonesArray*) FairRootManager::Instance()->GetObject("StsPoint");
  202. if (itsHF && itsP) FairRootManager::Instance()->
  203. Register("TpcKalmanTrack","MpdKalmanTrack", fTracks, kFALSE);
  204. else FairRootManager::Instance()->
  205. Register("TpcKalmanTrack","MpdKalmanTrack", fTracks, kTRUE);
  206. FairRootManager::Instance()->Register("TpcKalmanHit","KfHit", fKHits, kFALSE);
  207. }
  208. //__________________________________________________________________________
  209. void MpdTpcKalmanFilter::Finish()
  210. {
  211. ///
  212. fKHits->Delete();
  213. fTracks->Delete();
  214. fTrackCand->Delete();
  215. //((TH1F*)gROOT->FindObjectAny("hChi2"))->Write();
  216. //((TH1F*)gROOT->FindObjectAny("hNhits"))->Write();
  217. //((TH1F*)gROOT->FindObjectAny("hZ"))->Write();
  218. }
  219. //__________________________________________________________________________
  220. void MpdTpcKalmanFilter::Exec(Option_t * option)
  221. {
  222. ///
  223. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  224. static Int_t first = 1;
  225. if (first && fHits->GetEntriesFast()) {
  226. first = 0;
  227. // Check geometry
  228. if (fHits->First()->GetUniqueID()) fModular = 1;
  229. cout << "\n !!!!! ***** ***** ***** !!!!! " << endl;
  230. if (fModular) {
  231. FillGeoScheme();
  232. cout << " ***** TPC with modular design of readout chambers ***** " << endl;
  233. } else cout << " ***** TPC with cylindrical design ***** " << endl;
  234. cout << " !!!!! ***** ***** ***** !!!!! " << endl;
  235. }
  236. cout << "\n\n *** Event # " << ++fNofEvents << endl;
  237. cout << " ===== MpdTpcKalmanFilter =====\n";
  238. Reset();
  239. // Create Kalman hits
  240. if (fHits->GetEntriesFast() == 0) return;
  241. TString name = fHits->UncheckedAt(0)->GetName();
  242. cout << " Name: " << name << endl;
  243. if (name.Contains("TpcCluster")) Cluster2KalmanHits();
  244. else if (fModular == 0) MakeKalmanHits();
  245. else MakeKalmanHitsModul();
  246. // Evaluate vertex Z
  247. MpdVertexZfinder *vertexZ = (MpdVertexZfinder*) FairRun::Instance()->GetTask("MpdVertexZfinder");
  248. if (vertexZ) {
  249. vertexZ->SetHits(fKHits,fhLays);
  250. fVertZ = vertexZ->FindZ(fLayPointers, fZflag);
  251. }
  252. for (Int_t i = 0; i < fNPass; ++i) {
  253. fTrackCand->Delete();
  254. GetTrackSeeds(i);
  255. if (!fUseMCHit && i == 0) GetTrackSeedsEndCaps();
  256. cout << " Total number of hits for tracking: " << fHits->GetEntriesFast() << endl;
  257. cout << " Total number of track seeds: " << fTrackCand->GetEntriesFast() << endl;
  258. DoTracking(i);
  259. StoreTracks();
  260. cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
  261. if (i != fNPass - 1) ExcludeHits(); // exclude used hits
  262. }
  263. RemoveShorts(); // remove short tracks (Nhits < 4)
  264. cout << " Total number of found long tracks: " << fTracks->GetEntriesFast() << endl;
  265. if (0) {
  266. // To improve dL-resolution
  267. Int_t ntr = fTracks->GetEntriesFast();
  268. for (Int_t i = 0; i < ntr; ++i) {
  269. MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
  270. //track->GetCovariance()->Print();
  271. //track->GetParamNew()->Print();
  272. MpdKalmanFilter::Instance()->Refit(track, -1);
  273. track->SetLength(0);
  274. MpdKalmanFilter::Instance()->Refit(track, 1, kTRUE);
  275. //track->GetParamNew()->Print();
  276. //track->GetCovariance()->Print();
  277. }
  278. }
  279. AddHits(); // add hit objects to tracks
  280. // Apply track merging procedure (for clusters) - 11.10.2015
  281. if (!fUseMCHit) {
  282. for (Int_t ipass = 0; ipass < 3; ++ipass) {
  283. MergeTracks(ipass);
  284. }
  285. }
  286. GoToBeamLine(); // propagate tracks to the beam line
  287. GoOut(); // go outward for matching with outer detectors
  288. #ifdef GEANT
  289. //CheckTracks();
  290. #ifdef TRACK
  291. //PrintTracks(0);
  292. #endif
  293. #endif
  294. #ifdef MASSPROD
  295. //fTpcPoints->Clear("C");
  296. // fTstPoints->Clear("C");
  297. #endif
  298. /*
  299. FairMCEventHeader *mchead = (FairMCEventHeader*) FairRootManager::Instance()->GetObject("MCEventHeader.");
  300. Int_t event = mchead->GetEventID();
  301. //FairEventHeader *head = (FairEventHeader*) FairRootManager::Instance()->GetObject("EventHeader.");
  302. FairEventHeader *head = FairRunAna::Instance()->GetEvtHeader();
  303. head->SetEventId(event);
  304. */
  305. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  306. }
  307. //__________________________________________________________________________
  308. void MpdTpcKalmanFilter::FillGeoScheme()
  309. {
  310. /// Fill Kalman filter geometry manager info
  311. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  312. TVector2 p2;
  313. Int_t nRows[2] = {fSecGeo->NofRowsReg(0), fSecGeo->NofRowsReg(1)};
  314. Int_t nSec = fSecGeo->NofSectors();
  315. Double_t dPhi = fSecGeo->Dphi();
  316. // Get size along Z
  317. TString volName = "tpc01sv";
  318. TGeoVolume *sv = gGeoManager->GetVolume(volName);
  319. TGeoShape *shape = sv->GetShape();
  320. Double_t dZ = ((TGeoTube*)shape)->GetDZ();
  321. if (TString(shape->ClassName()) == "TGeoTube") {
  322. // Old sensitive volume
  323. TGeoVolume *endPlateG10 = gGeoManager->GetVolume("tpc01ppS");
  324. dZ -= ((TGeoTube*)endPlateG10->GetShape())->GetDZ() * 2;
  325. TGeoVolume *endPlateAl = gGeoManager->GetVolume("tpc01bpS");
  326. dZ -= ((TGeoTube*)endPlateAl->GetShape())->GetDZ() * 2;
  327. } else if (TString(shape->ClassName()) == "TGeoPgon") {
  328. // New sensitive volume
  329. Double_t phi1 = ((TGeoPgon*)shape)->Phi1();
  330. Double_t loc[3] = {100, 0, 10}, glob[3] = {0};
  331. gGeoManager->FindNode(loc[0],loc[1],loc[2]);
  332. if (TString(gGeoManager->GetPath()).Contains("row")) gGeoManager->CdUp(); // layered geometry
  333. if (TString(gGeoManager->GetPath()).Contains("sec")) gGeoManager->CdUp();
  334. gGeoManager->LocalToMaster(loc,glob);
  335. phi1 += -TMath::ATan2 (glob[1],glob[0]) * TMath::RadToDeg(); // due to rotation
  336. if (TMath::Abs(phi1-dPhi*TMath::RadToDeg()/2) > 0.001)
  337. Fatal("MpdTpcKalmanFilter::FillGeoScheme()"," !!! Inconsistent sens. volume parameters !!! ");
  338. dZ = ((TGeoPgon*)shape)->GetDZ() * 2;
  339. //FairGeoNode* membrane = (FairGeoNode*) (passNodes->FindObject("tpc01mb"));
  340. //if (!membrane) membrane = (FairGeoNode*) (sensNodes->FindObject("tpc01mb")); // it was a test version !
  341. TGeoVolume *membrane = gGeoManager->GetVolume("tpc01mb");
  342. dZ += ((TGeoTube*)membrane->GetShape())->GetDZ(); // membrane
  343. } else Fatal("MpdTpcKalmanFilter::FillGeoScheme()"," !!! Unknown sensitive volume shape !!! ");
  344. // Use TGeoManager to describe readout chamber geometry
  345. //TGeoMedium *medium = gGeoManager->GetMedium("TPCmixture");
  346. //TGeoVolume *tpcVolOrig = gGeoManager->GetVolume("tpc01sv");
  347. TGeoMedium *medium = gGeoManager->GetMedium("air");
  348. TGeoVolume *tpcVolOrig = gGeoManager->GetVolume("tpcChamber1");
  349. TGeoVolume *tpcVol = new TGeoVolume("tpc01sv_new", tpcVolOrig->GetShape(), medium);
  350. //TGeoVolume* tpcsec[2];
  351. //Double_t rads[2];
  352. Double_t dx1 = fSecGeo->GetMinY() * TMath::Tan(dPhi/2);
  353. Double_t dx2 = fSecGeo->GetMaxY() * TMath::Tan(dPhi/2);
  354. Double_t dy = dZ, dz = (fSecGeo->GetMaxY() - fSecGeo->GetMinY()) / 2.;
  355. Double_t rad = (fSecGeo->GetMaxY() + fSecGeo->GetMinY()) / 2.;
  356. TGeoVolume *tpcsec = gGeoManager->MakeTrd1("tpcsec", medium, dx1, dx2, dy, dz), *vol = NULL;
  357. //TGeoVolume * Division(const char* name, const char* mother, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step, Int_t numed = 0, Option_t* option = "");
  358. //gGeoManager->Division("tpcrow", "tpcsec", 3, nRows, 0., 0., 0, "N");
  359. // The code below is based on TGeoTrd1::Divide()
  360. Int_t ndiv = nRows[0] + nRows[1];
  361. Double_t start = -dz, dx1n, dx2n, zmax, step = fSecGeo->PadHeight(), zmin = start - step; // end = dz,
  362. for (Int_t id = 0; id < ndiv; ++id) {
  363. zmin += step;
  364. if (id == nRows[0]) step = fSecGeo->PadHeight(1);
  365. zmax = zmin + step;
  366. dx1n = 0.5 * (dx1 * (dz - zmin) + dx2 * (dz + zmin)) / dz;
  367. dx2n = 0.5 * (dx1 * (dz - zmax) + dx2 * (dz + zmax)) / dz;
  368. shape = new TGeoTrd1(dx1n, dx2n, dy, step/2.);
  369. vol = new TGeoVolume("tpcrow", shape, tpcsec->GetMedium());
  370. TGeoHMatrix *matr = new TGeoHMatrix();
  371. Double_t transl[3] = {0.0, 0.0, zmin+step/2};
  372. matr->SetTranslation(transl);
  373. tpcsec->AddNode(vol, id+1, matr);
  374. }
  375. //
  376. for (Int_t isec = 0; isec < nSec; ++isec) {
  377. Double_t phi = dPhi * isec;
  378. TGeoTranslation t(rad*TMath::Cos(phi), rad*TMath::Sin(phi), 0.);
  379. TGeoRotation r;
  380. Double_t phiDeg = TMath::RadToDeg() * phi;
  381. //r.SetAngles(90,phiDeg+90,0,0,90,phiDeg);
  382. r.SetAngles(90,phiDeg-90,0,0,90,phiDeg);
  383. TGeoCombiTrans c(t,r);
  384. TGeoHMatrix *matr = new TGeoHMatrix(c);
  385. //gGeoManager->RegisterMatrix(matr);
  386. // AddNode(const TGeoVolume* vol, Int_t copy_no, TGeoMatrix* mat = 0, Option_t* option = "")
  387. tpcVol->AddNode(tpcsec, isec+1, matr);
  388. }
  389. //gGeoManager->ClearPhysicalNodes();
  390. //gGeoManager->ClearThreadData();
  391. //gGeoManager->ClearThreadsMap();
  392. gGeoManager->ReplaceVolume(tpcVolOrig, tpcVol);
  393. //gGeoManager->RefreshPhysicalNodes();
  394. // Fill GeoScheme
  395. TVector3 p3norm, p3pos;
  396. Double_t padHeight = fSecGeo->PadHeight();
  397. Double_t rStart = fSecGeo->GetMinY() - padHeight / 2;
  398. Int_t nRows2 = fSecGeo->NofRows();
  399. for (Int_t isec = 0; isec < nSec; ++isec) {
  400. Double_t phi = dPhi * isec;
  401. p3norm.SetXYZ(TMath::Cos(phi), TMath::Sin(phi), 0.);
  402. Double_t x0 = rStart * p3norm.X(), y0 = rStart * p3norm.Y();
  403. padHeight = fSecGeo->PadHeight();
  404. for (Int_t lay = 0; lay < nRows2; ++lay) {
  405. Int_t detId = lay * 1000000 + isec;
  406. geo->SetNormal(detId, p3norm); // normal vector
  407. if (lay >= nRows[0]) padHeight = fSecGeo->PadHeight(1); // outer ROC region
  408. x0 += padHeight * p3norm.X();
  409. y0 += padHeight * p3norm.Y();
  410. TGeoNode *node = gGeoManager->FindNode(x0, y0, 0.);
  411. geo->SetPath(detId, gGeoManager->GetPath()); // TGeo path
  412. geo->SetDetId(gGeoManager->GetPath(), detId); // path to detId
  413. const Double_t *transl = gGeoManager->GetCurrentMatrix()->GetTranslation();
  414. p3pos.SetXYZ(transl[0], transl[1], 0.);
  415. geo->SetGlobalPos(detId, p3pos); // padrow position
  416. TGeoTrd1 *shape1 = (TGeoTrd1*) node->GetVolume()->GetShape();
  417. if (TString(shape1->Class()->GetName()) != "TGeoTrd1") { cout <<
  418. " !!! MpdTpcKalmanFilter::FillGeoScheme(): Wrong shape " << endl; exit(0); }
  419. p2.Set((shape1->GetDx1()+shape1->GetDx2())/2, shape1->GetDy());
  420. geo->SetSize(detId, p2); // sector size
  421. }
  422. //exit(0);
  423. }
  424. }
  425. //__________________________________________________________________________
  426. void MpdTpcKalmanFilter::RemoveShorts()
  427. {
  428. /// Remove short tracks (Nhits < 4)
  429. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  430. Int_t nReco = fTracks->GetEntriesFast();
  431. for (Int_t i = 0; i < nReco; ++i) {
  432. MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
  433. if (track->GetNofHits() < 4) fTracks->RemoveAt(i);
  434. }
  435. fTracks->Compress();
  436. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  437. }
  438. //__________________________________________________________________________
  439. void MpdTpcKalmanFilter::GoToBeamLine()
  440. {
  441. /// Propagate tracks to the beam line (mult. scat. and energy loss corrections)
  442. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  443. const Int_t nR = 7;
  444. //const Double_t rad[nR] = {27.035, 27.57, 28.105, 30.64, 33.15, 33.665, 34.178};
  445. //const Double_t dx[nR] = {0.0031, 0.00085, 0.0031, 0.0003, 0.001, 0.00085, 0.001};
  446. const Double_t rad[nR] = {27.00, 27.01, 27.31, 27.32, 33.82, 33.83, 34.13};
  447. 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};
  448. Bool_t ok = 0;
  449. Int_t nReco = fTracks->GetEntriesFast();
  450. for (Int_t i = 0; i < nReco; ++i) {
  451. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  452. if (track->GetNodeNew() == "") continue;
  453. MpdKalmanHit hit;
  454. hit.SetDetectorID(-999);
  455. hit.SetType(MpdKalmanHit::kFixedR);
  456. for (Int_t j = nR-1; j >= 0; --j) {
  457. hit.SetPos(rad[j]);
  458. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(rad[j],0.,0.),kTRUE); // put hit in GeoScheme
  459. ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kTRUE);
  460. if (!ok) {
  461. fTracks->RemoveAt(i);
  462. break;
  463. }
  464. //continue;
  465. // Add multiple scattering
  466. TMatrixDSym *cov = track->Weight2Cov();
  467. Double_t th = track->GetParamNew(3);
  468. Double_t cosTh = TMath::Cos(th);
  469. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[j]);
  470. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  471. (*cov)(3,3) += angle2;
  472. //cov->Print();
  473. Int_t iok = 0;
  474. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  475. track->SetWeight(*cov);
  476. }
  477. if (!ok) continue;
  478. TMatrixDSym *cov = track->Weight2Cov();
  479. // Update track
  480. track->SetParam(*track->GetParamNew());
  481. track->SetPos(track->GetPosNew());
  482. // Correct for energy loss
  483. //*
  484. Double_t theta = TMath::PiOver2() - track->GetParam(3);
  485. Double_t pt = 1. / track->GetParam(4);
  486. Double_t ptCor = CorrectForLoss(TMath::Abs(pt), theta, track->GetTrackID()); // ionization loss correction
  487. track->SetParam(4,TMath::Sign(1./ptCor,pt));
  488. track->SetParamNew(*track->GetParam());
  489. Double_t sigmaPt = CorrectForLossFluct(TMath::Abs(ptCor), theta, 0.1396, 1); // ionization loss correction fluct.
  490. (*cov)(4,4) += sigmaPt*sigmaPt;
  491. //*/
  492. //*
  493. hit.SetPos(0.);
  494. hit.SetMeas(0,track->GetParam(2));
  495. //hit.SetRphi(track->GetParam(2)); // track Phi - check if it is necessary !!!!!!!
  496. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(0.,0.,0.),kTRUE); // put hit in GeoScheme
  497. //*/
  498. /*
  499. MpdKalmanHitZ hit;
  500. hit.SetZ(0.);
  501. */
  502. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  503. Double_t pos = track->GetPosNew();
  504. //Double_t pos = ((MpdKalmanHitR*)track->GetHits()->Last())->GetR(); // position at last hit
  505. //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
  506. //track->GetCovariance()->Print();
  507. //(track->GetWeightAtHit()->Invert()).Print();
  508. MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
  509. //track->Weight2Cov();
  510. //track->GetCovariance()->Print();
  511. track->SetPos(pos); // save position after passing inner shell
  512. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  513. //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
  514. track->SetLengAtHit(track->GetLength()-track->GetLengAtHit()); // track length from PCA to the nearest hit
  515. }
  516. fTracks->Compress();
  517. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  518. }
  519. //__________________________________________________________________________
  520. void MpdTpcKalmanFilter::StoreTracks()
  521. {
  522. /// Transfer tracks from fTrackCand to fTracks
  523. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  524. Int_t nFound = fTracks->GetEntriesFast();
  525. for (Int_t i = 0; i < fNTracks; ++i) {
  526. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTrackCand->UncheckedAt(i);
  527. track->Weight2Cov();
  528. new ((*fTracks)[nFound++]) MpdTpcKalmanTrack(*track);
  529. //fTrackCand->RemoveAt(i);
  530. }
  531. fNTracks = 0;
  532. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  533. }
  534. //__________________________________________________________________________
  535. void MpdTpcKalmanFilter::AddHits(Int_t indx0)
  536. {
  537. /// Add hit objects to tracks and compute number of wrongly assigned hits
  538. /// (hits with ID different from ID of the track)
  539. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  540. Int_t ib = 0, ie = fTracks->GetEntriesFast() - 1;
  541. if (indx0 >= 0) ib = ie = indx0; // process one track
  542. for (Int_t i = ib; i <= ie; ++i) {
  543. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  544. Int_t nHits = track->GetNofHits();
  545. track->SetNofHits(nHits);
  546. TClonesArray &trHits = *track->GetTrHits();
  547. //Int_t idOld = track->GetTrackID(); // 24.12.2018
  548. SetTrackID(track); // set track ID as ID of majority of its hits
  549. Int_t nWrong = 0, motherID = track->GetTrackID();
  550. // Get track mother ID
  551. MpdMCTrack *mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(motherID);
  552. while (mctrack->GetMotherId() >= 0) {
  553. motherID = mctrack->GetMotherId();
  554. mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(mctrack->GetMotherId());
  555. }
  556. TObjArray *hits = track->GetHits();
  557. for (Int_t j = 0; j < nHits; ++j) {
  558. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  559. MpdKalmanHit *newHit = new (trHits[j]) MpdKalmanHit(*hit);
  560. MpdTpcHit *p = GetTpcHit(hit);
  561. Int_t motherID1 = p->GetTrackID();
  562. //AZ 19-07-15 if (!fUseMCHit) motherID1 = p->GetRefIndex();
  563. // Get point mother ID
  564. mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(motherID1);
  565. while (mctrack->GetMotherId() >= 0) {
  566. motherID1 = mctrack->GetMotherId();
  567. mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(mctrack->GetMotherId());
  568. }
  569. if (motherID1 != motherID) nWrong++;
  570. hits->AddAt(newHit, j); // replace content - 13.10.2015
  571. }
  572. track->SetNofWrong(nWrong);
  573. track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
  574. }
  575. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  576. }
  577. //__________________________________________________________________________
  578. void MpdTpcKalmanFilter::SetTrackID(MpdTpcKalmanTrack* track)
  579. {
  580. /// Set track ID as ID of majority of its hits
  581. const Int_t idMax = 9999;
  582. vector<Int_t> ids(idMax);
  583. Int_t nHits = track->GetNofHits(), locMax = 0, size = idMax;
  584. TObjArray *hits = track->GetHits();
  585. for (Int_t i = 0; i < nHits; ++i) {
  586. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
  587. Int_t id = GetHitID(hit);
  588. if (id >= size) { size = id + 1; ids.resize(size); }
  589. ++ids[id];
  590. if (ids[id] > ids[locMax]) locMax = id;
  591. }
  592. if (ids[locMax] > 1) track->SetTrackID(locMax);
  593. }
  594. //__________________________________________________________________________
  595. void MpdTpcKalmanFilter::GetTrackSeeds(Int_t iPass)
  596. {
  597. /// Build track seeds from central detector hits
  598. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  599. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  600. Int_t layMax0 = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  601. TH1F *hPt = (TH1F*) gROOT->FindObjectAny("hPt");
  602. TH1F *hZ = (TH1F*) gROOT->FindObjectAny("hZ");
  603. // Estimate Z-position of vertex
  604. // Loop over layers
  605. Int_t layBeg = layMax0, layEnd = layBeg - 2, iDir = -1, dLays = 6; //10;
  606. //Int_t layMax = layMax0, dLays = 6; //10;
  607. Double_t zToler = 10;
  608. if (iPass > 0) {
  609. //25.12.2018 layBeg = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
  610. layBeg = 6;
  611. //25.12.2018 layEnd = layBeg + 2;
  612. layEnd = layBeg + 4;
  613. //25.12.2018 iDir = 1;
  614. iDir = 2;
  615. //25.12.2018 dLays = 1;
  616. dLays = 2;
  617. zToler = 30;
  618. }
  619. if (fZflag < 0) zToler = 50; // bad quality of vertex estimate
  620. MpdTpcKalmanTrack *track = 0x0;
  621. MpdTpcHit *h1 = 0x0;//, *h2 = 0x0;
  622. TVector3 pos1, pos2, posLoc;
  623. Double_t rad1, phi1, rad2, phi2, params[3];
  624. for (Int_t lay = layBeg; lay != layEnd; lay+=iDir) {
  625. Int_t nHits1 = (Int_t) fhLays->GetBinContent(lay+1,0), isec = -1;
  626. Int_t nHits2 = (Int_t) fhLays->GetBinContent(lay+dLays*TMath::Sign(1,iDir)+1,0);
  627. //cout << "Hits: " << nHits1 << " " << nHits2 << endl;
  628. // Loop over hits in first layer
  629. for (Int_t ihit1 = 0; ihit1 < nHits1; ++ihit1) {
  630. MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay]+ihit1);
  631. if (hit1->GetIndex() >= 0) h1 = (MpdTpcHit*) fHits->UncheckedAt(hit1->GetIndex());
  632. if (hit1->GetFlag() & MpdKalmanHit::kUsed) continue;
  633. if (fModular) {
  634. //posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  635. posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  636. isec = hit1->GetDetectorID() % 1000000;
  637. //pos1 = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
  638. fSecGeo->Local2Global(isec, posLoc, pos1);
  639. rad1 = pos1.Pt();
  640. phi1 = pos1.Phi();
  641. //cout << posLoc.X() << " " << posLoc.Y() << " " << posLoc.Z() << endl;
  642. //cout << pos1.X() << " " << pos1.Y() << " " << pos1.Z() << endl;
  643. //cout << rad1 << " " << phi1 << endl;
  644. } else {
  645. rad1 = hit1->GetPos();
  646. phi1 = hit1->GetMeas(0) / rad1;
  647. }
  648. // Loop over hits in second (closer to the beam for the first pass) layer
  649. for (Int_t ihit2 = 0; ihit2 < nHits2; ++ihit2) {
  650. MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay+dLays*TMath::Sign(1,iDir)]+ihit2);
  651. //if (hit2->GetIndex() >= 0) h2 = (MpdTpcHit*) fHits->UncheckedAt(hit2->GetIndex());
  652. //if (h2->GetTrackID() != h1->GetTrackID()) continue; //!!! for test - exact ID matching
  653. if (hit2->GetFlag() & MpdKalmanHit::kUsed) continue;
  654. if (fModular) {
  655. //posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
  656. posLoc.SetXYZ(hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
  657. Int_t isec1 = hit2->GetDetectorID() % 1000000;
  658. //if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fPadPlane->nSectors()-1) continue;
  659. //pos2 = fPadPlane->fromSectorReferenceFrame(posLoc, isec1);
  660. if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fSecGeo->NofSectors()-1) continue;
  661. fSecGeo->Local2Global(isec1, posLoc, pos2);
  662. rad2 = pos2.Pt();
  663. if ((rad2 - rad1) * TMath::Sign(1,iDir) < 0.5) continue;
  664. phi2 = pos2.Phi();
  665. } else {
  666. rad2 = hit2->GetPos();
  667. phi2 = hit2->GetMeas(0) / rad2;
  668. }
  669. Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi1,phi2) - phi1;
  670. if (TMath::Abs(dPhi) > TMath::PiOver4()) continue;
  671. Double_t pt = EvalPt(hit1, hit2, params);
  672. //Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1)-hit2->GetMeas(1)) / (rad1-rad2) * rad2;
  673. Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1) - hit2->GetMeas(1)) / params[0] * params[1];
  674. if (hZ && iPass) hZ->Fill(zvert);
  675. if (TMath::Abs(zvert-fVertZ) > zToler) continue;
  676. //printf("%4d %4d %4d %2d %2d %5.3f %5.3f %7.3f %7.3f %4d %4d %7.3f %6.3f\n", fNTracks, ihit1, ihit2,
  677. // hit1->GetLayer(), hit2->GetLayer(), rad1, rad2, hit1->GetMeas(1), hit2->GetMeas(1),
  678. // h1->GetTrackID(), h2->GetTrackID(), zvert, pt);
  679. if (hPt) hPt->Fill(TMath::Abs(pt));
  680. //if (TMath::Abs(pt) > 5) cout << " ***!!! " << pt << endl;
  681. //if (TMath::Abs(pt) < 0.05) continue; // skip low-Pt tracks
  682. if (TMath::Abs(pt) < 0.02) continue; // skip low-Pt tracks
  683. if (TMath::Abs(pt) > 10.) continue; // skip high-Pt tracks
  684. TVector3 vertex(0,0,0);
  685. if (fModular) {
  686. if (iPass == 0) track = new ((*fTrackCand)[fNTracks++])
  687. MpdTpcKalmanTrack(hit1, hit2, vertex, pos1, pos2, pt, params);
  688. else track = new ((*fTrackCand)[fNTracks++])
  689. MpdTpcKalmanTrack(hit2, hit1, vertex, pos2, pos1, pt, params);
  690. TString s1 = geo->Path (hit1->GetDetectorID());
  691. track->SetNodeNew(s1);
  692. //cout << hit1->GetDetectorID() << " " << track->GetNodeNew() << endl;
  693. } else {
  694. if (iPass == 0) track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2,
  695. vertex, pt, params);
  696. else track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit2, hit1, vertex,
  697. pt, params);
  698. }
  699. if (lay <= layMax0 - 1) track->SetDirection(MpdKalmanTrack::kOutward);
  700. if (iPass > 0) {
  701. track->GetHits()->Clear();
  702. track->GetHits()->Add(hit1);
  703. //track->GetHits()->Add(hit2);
  704. //if (fModular) track->SetNodeNew(geo->Path(hit2->GetDetectorID())); // 29-apr-12
  705. // 7-may-12
  706. track->SetParam(1,hit1->GetMeas(1)); // hit closest to the beam line
  707. track->SetPos(TMath::Min(rad1,rad2));
  708. //31.12.2018 track->SetParam(0,track->GetParam(0)/track->GetPosNew()*track->GetPos());
  709. //31.12.2018
  710. if (rad1 < rad2) track->SetParam(0,MpdKalmanFilter::Instance()->Proxim(track->GetParam(0)/rad2,pos1.Phi())*rad1);
  711. else track->SetParam(0,MpdKalmanFilter::Instance()->Proxim(track->GetParam(0)/rad1,pos2.Phi())*rad2);
  712. //
  713. track->SetPosNew(track->GetPos());
  714. track->SetParamNew(*track->GetParam());
  715. }
  716. //*
  717. if (fModular) {
  718. /*TString s1, s2;
  719. s1 += hit1->GetDetectorID() % 1000000; // sector No.
  720. s2 += hit2->GetDetectorID() % 1000000; // sector No.
  721. if (iPass > 0) { track->SetNode(s1); track->SetNodeNew(s2); }
  722. else { track->SetNode(s2); track->SetNodeNew(s1); }
  723. */
  724. /*
  725. TString s1 = geo->Path(hit1->GetDetectorID());
  726. TString s2 = geo->Path(hit2->GetDetectorID());
  727. if (iPass > 0) track->SetNodeNew(s2); // trick
  728. else track->SetNodeNew(s1);
  729. */
  730. }
  731. //*/
  732. if (h1) track->SetTrackID(h1->GetTrackID());
  733. //((MpdKalmanTrack*)fTrackCand->UncheckedAt(fNTracks-1))->GetParam()->Print();
  734. //cout << fNTracks-1 << endl;
  735. //((TpcLheKalmanTrack*)fTrackCand->UncheckedAt(fNTracks-1))->Weight2Cov()->Print();
  736. }
  737. }
  738. }
  739. //hLays->Write();
  740. //hZ->Write();
  741. //hPt->Write();
  742. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  743. }
  744. //__________________________________________________________________________
  745. void MpdTpcKalmanFilter::GetTrackSeedsEndCaps()
  746. {
  747. /// Build track seeds from end-cap regions (inward tracks from hits near endcaps)
  748. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  749. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  750. Int_t layBeg = fSecGeo->NofRows() - 3, layEnd = 4, iDir = -1, dLays = 5;
  751. Double_t zToler = 10, zMax = 160;
  752. if (fZflag < 0) zToler = 50; // bad quality of vertex estimate
  753. MpdTpcKalmanTrack *track = 0x0;
  754. MpdTpcHit *h1 = 0x0;//, *h2 = 0x0;
  755. TVector3 pos1, pos2, posLoc;
  756. Double_t rad1, phi1, rad2, phi2, params[3];
  757. // Loop over layers
  758. for (Int_t lay = layBeg; lay != layEnd; lay+=iDir) {
  759. Int_t nHits1 = (Int_t) fhLays->GetBinContent(lay+1,0), isec = -1;
  760. Int_t nHits2 = (Int_t) fhLays->GetBinContent(lay+dLays*iDir+1,0);
  761. if (nHits1 <= 0 || nHits2 <= 0) continue;
  762. // Loop over hits in first layer
  763. for (Int_t ihit1 = 0; ihit1 < nHits1; ++ihit1) {
  764. MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay]+ihit1);
  765. if (hit1->GetIndex() >= 0) h1 = (MpdTpcHit*) fHits->UncheckedAt(hit1->GetIndex());
  766. if (hit1->GetFlag() & MpdKalmanHit::kUsed) continue;
  767. if (TMath::Abs(hit1->GetMeas(1)) > zMax) continue; // too close to end-cap
  768. if (fModular) {
  769. //posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  770. posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  771. isec = hit1->GetDetectorID() % 1000000;
  772. fSecGeo->Local2Global(isec, posLoc, pos1);
  773. rad1 = pos1.Pt();
  774. phi1 = pos1.Phi();
  775. Double_t tanThe = rad1 / TMath::Abs(pos1.Z()-fVertZ);
  776. Int_t row = fSecGeo->PadRow(hit1->GetDetectorID());
  777. Double_t padHeight = fSecGeo->PadHeight();
  778. if (row >= fSecGeo->NofRowsReg(0)) padHeight = fSecGeo->PadHeight(1);
  779. if (TMath::Abs(pos1.Z()) < zMax - padHeight / tanThe) continue; // too far from end-cap
  780. //cout << posLoc.X() << " " << posLoc.Y() << " " << posLoc.Z() << endl;
  781. //cout << pos1.X() << " " << pos1.Y() << " " << pos1.Z() << endl;
  782. //cout << rad1 << " " << phi1 << endl;
  783. } else {
  784. rad1 = hit1->GetPos();
  785. phi1 = hit1->GetMeas(0) / rad1;
  786. }
  787. // Loop over hits in second (closer to the beam for the first pass) layer
  788. for (Int_t ihit2 = 0; ihit2 < nHits2; ++ihit2) {
  789. MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay+dLays*iDir]+ihit2);
  790. //if (hit2->GetIndex() >= 0) h2 = (MpdTpcHit*) fHits->UncheckedAt(hit2->GetIndex());
  791. //if (h2->GetTrackID() != h1->GetTrackID()) continue; //!!! for test - exact ID matching
  792. if (hit2->GetFlag() & MpdKalmanHit::kUsed) continue;
  793. if (fModular) {
  794. //posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
  795. posLoc.SetXYZ(hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
  796. Int_t isec1 = hit2->GetDetectorID() % 1000000;
  797. //if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fPadPlane->nSectors()-1) continue;
  798. if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fSecGeo->NofSectors()-1) continue;
  799. fSecGeo->Local2Global(isec1, posLoc, pos2);
  800. rad2 = pos2.Pt();
  801. if ((rad2 - rad1) * iDir < 0.5) continue;
  802. phi2 = pos2.Phi();
  803. } else {
  804. rad2 = hit2->GetPos();
  805. phi2 = hit2->GetMeas(0) / rad2;
  806. }
  807. Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi1,phi2) - phi1;
  808. if (TMath::Abs(dPhi) > TMath::PiOver4()) continue;
  809. Double_t pt = EvalPt(hit1, hit2, params);
  810. // Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1)-hit2->GetMeas(1)) / (rad1-rad2) * rad2;
  811. Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1) - hit2->GetMeas(1)) / params[0] * params[1];
  812. if (TMath::Abs(zvert-fVertZ) > zToler) continue;
  813. //cout << fNTracks << " " << ihit1 << " " << ihit2 << " " << hit1->GetLayer() << " " << hit2->GetLayer() << " " << rad1 << " " << rad2 << " " << hit1->GetMeas(1) << " " << hit2->GetMeas(1) << " " << h1->GetTrackID() << " " << h2->GetTrackID() << " " << zvert << " " /*<< hit1->GetUsage() << " " << hit2->GetUsage() << " "*/ << pt << endl;
  814. if (TMath::Abs(pt) < 0.02) continue; // skip low-Pt tracks
  815. if (TMath::Abs(pt) > 10.) continue; // skip high-Pt tracks
  816. TVector3 vertex(0,0,0);
  817. if (fModular) {
  818. track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2, vertex,
  819. pos1, pos2, pt, params);
  820. TString s1 = geo->Path (hit1->GetDetectorID());
  821. track->SetNodeNew(s1);
  822. //cout << hit1->GetDetectorID() << " " << track->GetNodeNew() << endl;
  823. } else {
  824. track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2, vertex,
  825. pt, params);
  826. }
  827. track->SetDirection(MpdKalmanTrack::kOutward);
  828. //*
  829. if (fModular) {
  830. /*TString s1, s2;
  831. s1 += hit1->GetDetectorID() % 1000000; // sector No.
  832. s2 += hit2->GetDetectorID() % 1000000; // sector No.
  833. if (iPass > 0) { track->SetNode(s1); track->SetNodeNew(s2); }
  834. else { track->SetNode(s2); track->SetNodeNew(s1); }
  835. */
  836. /*
  837. TString s1 = geo->Path(hit1->GetDetectorID());
  838. TString s2 = geo->Path(hit2->GetDetectorID());
  839. if (iPass > 0) track->SetNodeNew(s2); // trick
  840. else track->SetNodeNew(s1);
  841. */
  842. }
  843. //*/
  844. if (h1) track->SetTrackID(h1->GetTrackID());
  845. }
  846. }
  847. }
  848. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  849. }
  850. //__________________________________________________________________________
  851. void MpdTpcKalmanFilter::Cluster2KalmanHits()
  852. {
  853. /// Create Kalman hits from Rec. Points found from clusters
  854. Int_t lay, layMax = 0, nClus = fHits->GetEntriesFast();
  855. // !!!!!!!!! Geometry dependent values
  856. //Double_t rMin = 20.5, rMax = 109.5, dR = (rMax-rMin)/50; // dR = 1.78; // old version
  857. //Double_t rMin = 29.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.4123; // new version (with dead material)
  858. Double_t rMin = 34.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.1962; // new version (with dead material)
  859. //Double_t rMin = 35.00, rMax = 100.00, dR = (rMax-rMin)/50; // 1.1962; // new version (with dead material)
  860. // !!!!!!!!!
  861. Double_t xyz[3], rPhiErr = 0.05, zErr = 0.1;
  862. //TClonesArray hits("TpcLheHit");
  863. for (int j = 0; j < nClus; ++j ) {
  864. TpcCluster* clus = (TpcCluster*) fHits->UncheckedAt(j);
  865. xyz[0] = clus->ret_x();
  866. xyz[1] = clus->ret_y();
  867. xyz[2] = clus->ret_z();
  868. //TpcLheHit *hit = new (hits[j]) TpcLheHit(xyz, 0);
  869. //hit->SetR (TMath::Sqrt (xyz[0]*xyz[0]+xyz[1]*xyz[1]));
  870. Double_t rad = TMath::Sqrt (xyz[0] * xyz[0] + xyz[1] * xyz[1]);
  871. Double_t rPhi = TMath::ATan2(xyz[1],xyz[0]) * rad;
  872. lay = (Int_t) ((rad - rMin) / dR);
  873. lay = TMath::Max (lay, 0);
  874. layMax = TMath::Max (lay, layMax);
  875. //hit->SetLayer(lay);
  876. fhLays->Fill(lay+0.1);
  877. // Create Kalman hits
  878. //MpdKalmanHit *hitK = new ((*fKHits)[j]) MpdKalmanHit(rad, rPhi, xyz[2],
  879. // rPhiErr, zErr, lay, -1);
  880. //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
  881. Double_t meas[2] = {rPhi, xyz[2]};
  882. Double_t err[2] = {rPhiErr, zErr};
  883. Double_t cossin[2] = {1., 0.};
  884. MpdKalmanHit *hitK = new ((*fKHits)[j]) MpdKalmanHit(0, 2, MpdKalmanHit::kFixedR, meas, err, cossin, 0., rad, -1);
  885. }
  886. cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << endl;
  887. fKHits->Sort();
  888. fLayPointers = new Int_t [layMax+1];
  889. Int_t ipos = 0;
  890. for (Int_t i = layMax; i >= 0; --i) {
  891. //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
  892. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  893. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  894. fLayPointers[i] = ipos;
  895. ipos += (Int_t) fhLays->GetBinContent(i+1,0);
  896. }
  897. }
  898. //__________________________________________________________________________
  899. void MpdTpcKalmanFilter::MakeKalmanHits()
  900. {
  901. /// Create Kalman hits.
  902. Int_t nHits = fHits->GetEntriesFast();
  903. cout << " Number of merged hits: " << nHits << endl;
  904. Int_t nKh = 0, lay = 0, layMax = 0;
  905. for (int j = 0; j < nHits; ++j ) {
  906. MpdTpcHit* hit = (MpdTpcHit*) fHits->UncheckedAt(j);
  907. //if (hit->GetTrackID() != 1783) continue; // keep only one track
  908. lay = hit->GetLayer();
  909. layMax = TMath::Max (lay, layMax);
  910. fhLays->Fill(lay+0.1);
  911. // Add errors
  912. Double_t dRPhi = 0, dZ = 0;
  913. gRandom->Rannor(dRPhi,dZ);
  914. Double_t z = hit->GetZ() + dZ * hit->GetDz(); //add error
  915. Double_t rPhi = hit->GetRphi() + dRPhi * hit->GetDx(); //add error
  916. //Double_t phi = rPhi / hit->GetR();
  917. //hit->SetX(hit->GetR()*TMath::Cos(phi));
  918. //hit->SetY(hit->GetR()*TMath::Sin(phi));
  919. // Create Kalman hits
  920. //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
  921. Double_t meas[2] = {rPhi, z};
  922. Double_t err[2] = {hit->GetDx(), hit->GetDz()};
  923. Double_t cossin[2] = {1., 0.};
  924. MpdKalmanHit *hitK = new ((*fKHits)[nKh++])
  925. MpdKalmanHit(lay*1000000+nKh-1, 2, MpdKalmanHit::kFixedR, meas, err, cossin, hit->GetEnergyLoss()/hit->GetStep(), hit->GetR(), j);
  926. hitK->SetLength(hit->GetLength());
  927. //hitK->SetDedx (hit->GetEdep()/hit->GetStep());
  928. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(hitK,TVector3(hit->GetX(),hit->GetY(),hit->GetZ()));
  929. }
  930. fLayPointers = new Int_t [layMax+1];
  931. Int_t ipos = 0;
  932. for (Int_t i = layMax; i >= 0; --i) {
  933. //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
  934. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  935. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  936. fLayPointers[i] = ipos;
  937. fLaySecBegPointers[i][0] = ipos;
  938. ipos += (Int_t) fhLays->GetBinContent(i+1,0);
  939. fLaySecEndPointers[i][0] = ipos;
  940. }
  941. }
  942. //__________________________________________________________________________
  943. void MpdTpcKalmanFilter::MakeKalmanHitsModul()
  944. {
  945. /// Create Kalman hits for the modular geometry of readout chambers.
  946. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  947. Int_t nHits = fHits->GetEntriesFast();
  948. cout << " Number of merged hits: " << nHits << endl;
  949. Int_t nKh = 0, lay = 0, layMax = 0;
  950. for (int j = 0; j < nHits; ++j ) {
  951. MpdTpcHit* hit = (MpdTpcHit*) fHits->UncheckedAt(j);
  952. lay = hit->GetLayer();
  953. layMax = TMath::Max (lay, layMax);
  954. fhLays->Fill(lay+0.1);
  955. // Add errors
  956. Double_t dX = 0, dZ = 0;
  957. gRandom->Rannor(dX, dZ);
  958. // Create Kalman hits
  959. //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
  960. Double_t err[2] = {hit->GetDx(), hit->GetDz()};
  961. //Double_t meas[2] = {-hit->GetLocalX()-dX*err[0], hit->GetLocalZ()+dZ*err[1]};
  962. Double_t meas[2] = {hit->GetLocalX(), hit->GetLocalZ()};
  963. if (fUseMCHit) { meas[0] += dX * err[0]; meas[1] += dZ * err[1]; }
  964. Double_t cossin[2] = {1., 0.};
  965. Int_t padID = hit->GetDetectorID();
  966. //TpcPadID padIDobj = TpcPadID::numberToPadID(padID);
  967. //padID = padIDobj.sector();
  968. padID = fSecGeo->Sector(padID);
  969. Double_t xsec = (hit->GetLocalY() + fSecGeo->GetMinY()) * TMath::Tan(fSecGeo->Dphi()/2);
  970. Double_t xloc = hit->GetLocalX(); Double_t edge = 0.0; // distance to sector boundary
  971. if (xloc < 0) edge = xloc + xsec;
  972. else edge = xloc - xsec;
  973. MpdKalmanHit *hitK = new ((*fKHits)[nKh++])
  974. MpdKalmanHit(lay*1000000+padID, 2, MpdKalmanHit::kFixedP, meas, err, cossin, hit->GetEnergyLoss()/hit->GetStep(), hit->GetLocalY(), j, edge);
  975. hitK->SetLength(hit->GetLength());
  976. hitK->SetFlag (hitK->GetFlag() | hit->GetFlag());
  977. //hitK->SetDedx (hit->GetEdep()/hit->GetStep());
  978. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(hitK,TVector3(hit->GetX(),hit->GetY(),hit->GetZ()));
  979. }
  980. fKHits->Sort();
  981. fLayPointers = new Int_t [layMax+1];
  982. Int_t ipos = 0;
  983. for (Int_t i = layMax; i >= 0; --i) {
  984. //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
  985. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  986. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  987. fLayPointers[i] = ipos;
  988. ipos += (Int_t) fhLays->GetBinContent(i+1,0);
  989. }
  990. // Fill indices for different sectors
  991. /*
  992. Int_t iseco = -1, isec = 0;
  993. for (Int_t i = layMax; i >= 0; --i) {
  994. Int_t nLay = GetNofHitsInLayer(lay);
  995. Int_t indx0 = GetHitsInLayer(lay);
  996. ipos = indx0;
  997. for (Int_t indx = 0; indx != nLay; ++indx) {
  998. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  999. isec = hit->GetDetectorID() % 1000000;
  1000. if (isec != iseco) {
  1001. fLaySecBegPointers[i].insert(pair<Int_t,Int_t>(isec,ipos+indx));
  1002. if (iseco > 0) fLaySecEndPointers[i].insert(pair<Int_t,Int_t>(iseco,ipos+indx));
  1003. iseco = isec;
  1004. }
  1005. }
  1006. }
  1007. fLaySecEndPointers[0].insert(pair<Int_t,Int_t>(isec,nHits));
  1008. */
  1009. nHits = fKHits->GetEntriesFast();
  1010. for (Int_t j = 0; j < nHits; ++j ) {
  1011. MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(j);
  1012. lay = hit->GetLayer();
  1013. Int_t isec = hit->GetDetectorID() % 1000000;
  1014. if (fLaySecBegPointers[lay][isec] < 0) {
  1015. fLaySecBegPointers[lay][isec] = j;
  1016. fLaySecEndPointers[lay][isec] = j + 1;
  1017. } else {
  1018. if (fLaySecBegPointers[lay][isec] > j) fLaySecBegPointers[lay][isec] = j;
  1019. if (fLaySecEndPointers[lay][isec] < j+1) fLaySecEndPointers[lay][isec] = j + 1;
  1020. }
  1021. }
  1022. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1023. }
  1024. //__________________________________________________________________________
  1025. Double_t MpdTpcKalmanFilter::EvalPt(const MpdKalmanHit *hit1, const MpdKalmanHit *hit2, Double_t *params)
  1026. {
  1027. /// Evaluate signed track Pt (curvature) assuming the track coming from the
  1028. /// primary vertex
  1029. TVector3 pos, posLoc;
  1030. Double_t rad1, phi1, rad2, phi2;
  1031. if (fModular) {
  1032. posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), 0.);
  1033. Int_t isec = hit1->GetDetectorID() % 1000000;
  1034. //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
  1035. fSecGeo->Local2Global(isec, posLoc, pos);
  1036. rad1 = pos.Pt();
  1037. phi1 = pos.Phi();
  1038. posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), 0.);
  1039. isec = hit2->GetDetectorID() % 1000000;
  1040. //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
  1041. fSecGeo->Local2Global(isec, posLoc, pos);
  1042. rad2 = pos.Pt();
  1043. phi2 = pos.Phi();
  1044. } else {
  1045. rad1 = hit1->GetPos();
  1046. rad2 = hit2->GetPos();
  1047. phi1 = hit1->GetMeas(0) / rad1;
  1048. phi2 = hit2->GetMeas(0) / rad2;
  1049. }
  1050. TVector2 vec1(rad1*TMath::Cos(phi1)-0.,rad1*TMath::Sin(phi1)-0.);
  1051. TVector2 vec2(rad2*TMath::Cos(phi2)-0.,rad2*TMath::Sin(phi2)-0.);
  1052. TVector2 vec21 = vec1 - vec2;
  1053. Double_t dist2 = vec2.Mod();
  1054. Double_t dist12 = vec21.Mod();
  1055. Double_t cosAlpha = vec2 * vec21 / dist2 / dist12;
  1056. Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
  1057. params[0] = TMath::ASin (dist12 / 2 / rad) * rad * TMath::Sign(2.0,rad1-rad2); // arc12 length
  1058. params[1] = TMath::ASin (dist2 / 2 / rad) * 2 * rad; // arc02 length
  1059. params[2] = params[0] / rad; // arc12 angle
  1060. Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
  1061. Double_t factor = 0.0003 * bz; // 0.3 * 0.01 * 5kG / 10
  1062. Double_t charge = phi1 - MpdKalmanFilter::Instance()->Proxim(phi1,phi2);
  1063. if (hit1->GetLayer() < hit2->GetLayer()) charge = -charge;
  1064. return factor * TMath::Abs(rad) * TMath::Sign(1., -charge);
  1065. }
  1066. //__________________________________________________________________________
  1067. void MpdTpcKalmanFilter::DoTracking(Int_t iPass)
  1068. {
  1069. /// Run Kalman tracking
  1070. //struct timeval tvStart, tvEnd;
  1071. //struct timezone tz;
  1072. //gettimeofday(&tvStart, &tz);
  1073. //fTracks->Sort(); // in descending order in Pt
  1074. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1075. Int_t nTracks = fTrackCand->GetEntriesFast();
  1076. //TH1F *hChi2 = (TH1F*) gROOT->FindObjectAny("hChi2");
  1077. //TH1F *hNhits = (TH1F*) gROOT->FindObjectAny("hNhits");
  1078. //#pragma omp parallel for
  1079. for (Int_t i = 0; i < nTracks; ++i) {
  1080. Int_t iok = 0;
  1081. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTrackCand->UncheckedAt(i);
  1082. MpdKalmanTrack::TrackDir dir0 = MpdKalmanTrack::kInward;
  1083. if (((MpdKalmanHit*)track->GetHits()->First())->GetLayer() < 10) dir0 = MpdKalmanTrack::kOutward;
  1084. //cout << " Track seed No. " << i << " " << track->GetTrackID() << " " << ((MpdKalmanHit*)track->GetHits()->First())->GetLayer() << endl;
  1085. // << track->GetPosNew() << endl;
  1086. //track->GetParamNew()->Print();
  1087. // Check if the seed hit was used already
  1088. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetHits()->First();
  1089. if (hit->GetFlag() & MpdKalmanHit::kUsed) iok = -1;
  1090. else {
  1091. if (track->GetDirection() == MpdKalmanTrack::kOutward) {
  1092. //cout << " !!! " << track->GetNode() << " !!! " << track->GetNodeNew() <<endl;
  1093. GoOutward(track);
  1094. /*cout << i << " " << track->GetChi2() << endl;*/
  1095. if (track->GetChi2()>1000) {
  1096. //cout << " !!! " << " " << track->GetNofHits() << " " << track->GetNode() << " " << track->GetNodeNew() << endl;
  1097. }
  1098. }
  1099. iok = RunKalmanFilter(track);
  1100. //if (fNofEvents > 54)
  1101. //cout << iok << " " << track->GetChi2() << " " << track->GetNofHits() << endl;
  1102. }
  1103. //if (iOK != 0) track->SetStatus(-1);
  1104. if (iok != 0) {
  1105. //fTracks->RemoveAt(i);
  1106. //cout << track->GetNofHits() << endl;
  1107. //#pragma omp critical
  1108. {
  1109. fTrackCand->RemoveAt(i);
  1110. }
  1111. //delete track; // !!! Object has been deleted already
  1112. } else {
  1113. track->SetDirection(dir0); // restore initial direction flag
  1114. /*
  1115. Int_t nHits = track->GetNofHits();
  1116. //#pragma omp critical
  1117. {
  1118. hChi2->Fill(track->GetChi2());
  1119. }
  1120. //#pragma omp critical
  1121. {
  1122. hNhits->Fill(nHits+0.1);
  1123. }
  1124. */
  1125. //if (track->GetChi2() < 1.) cout << " Strange: " << track->GetChi2() << " " << nHits << " " << track->GetDirection() << endl;
  1126. // Mark hits as being used
  1127. /*
  1128. TObjArray *hits = track->GetHits();
  1129. for (Int_t j = 0; j < nHits; ++j) {
  1130. hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1131. hit->SetFlag(-1);
  1132. }
  1133. */
  1134. }
  1135. if (iPass == 0 && track->GetNofHits() < 10) fTrackCand->RemoveAt(i); //AZ - 190820
  1136. }
  1137. //fTracks->Compress();
  1138. fTrackCand->Compress();
  1139. /*
  1140. nTracks = fTrackCand->GetEntriesFast();
  1141. for (Int_t i = 0; i < nTracks; ++i) {
  1142. TpcLheKalmanTrack *track = (TpcLheKalmanTrack*) fTrackCand->UncheckedAt(i);
  1143. Int_t nHits = track->GetNofHits(), iCharge = 0;
  1144. TObjArray *hits = track->GetHits();
  1145. TVector3 pmom, pmom1;
  1146. for (Int_t j = 0; j < nHits; ++j) {
  1147. TpcLheHit *hit = (TpcLheHit*) hits->UncheckedAt(j);
  1148. if (j == nHits-1) {
  1149. // Compare with MC track
  1150. Int_t nMC = fLHEtracks->GetEntriesFast();
  1151. for (Int_t itr = 0; itr < nMC; ++itr) {
  1152. TpcLheTrack *trMC = (TpcLheTrack*) fLHEtracks->UncheckedAt(itr);
  1153. if (trMC->GetTrackNumber() != hit->GetTrackID()) continue;
  1154. pmom = trMC->GetMomentum();
  1155. break;
  1156. }
  1157. TpcPoint *point = (TpcPoint*) fTpcPoints->UncheckedAt(hit->GetUniqueID());
  1158. point->Momentum(pmom1);
  1159. iCharge = -hit->GetRefIndex();
  1160. }
  1161. printf(" %4d",hit->GetTrackID());
  1162. }
  1163. printf(" hits\n");
  1164. cout << " Momentum: " << i << " " << 1./track->GetParamNew(4) << " " << pmom1.Pt()*iCharge << " " << pmom.Pt() << " " << track->GetChi2() << endl;
  1165. } // for (Int_t i = 0; i < nTracks;
  1166. */
  1167. // Remove doubles
  1168. RemoveDoubles(fTrackCand);
  1169. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1170. //gettimeofday(&tvEnd, &tz);
  1171. //int dif = 1000000*((int)tvEnd.tv_sec - (int)tvStart.tv_sec) + (tvEnd.tv_usec - tvStart.tv_usec);
  1172. //cout<<"TPC::DoTracking "<<dif<<" ns\n";
  1173. }
  1174. //__________________________________________________________________________
  1175. Int_t MpdTpcKalmanFilter::RunKalmanFilter(MpdTpcKalmanTrack *track)
  1176. {
  1177. /// Run Kalman filter
  1178. //cout << fHits->GetEntriesFast() << endl;
  1179. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1180. const Double_t dStepMax = 9; // max track length distance between consecutive hits
  1181. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  1182. Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  1183. MpdKalmanHit *hitOK = (MpdKalmanHit*) track->GetHits()->Last();
  1184. Int_t layOK = hitOK->GetLayer();
  1185. MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
  1186. //Int_t layBeg = layOK - 1, layEnd = -1, dLay = -1;
  1187. Int_t layBeg = layOK - 1, layEnd = ((MpdKalmanHit*)fKHits->Last())->GetLayer() - 1, dLay = -1;
  1188. if (trackDir == MpdKalmanTrack::kOutward) {
  1189. layBeg = layOK + 1;
  1190. layEnd = layMax + 1;
  1191. dLay = 1;
  1192. }
  1193. //Int_t indxOK = hits->IndexOf(hitOK);
  1194. //Int_t nHits = hits->GetEntriesFast();
  1195. Int_t miss = 0, iOK = 1;
  1196. if (fModular) iOK = 0;
  1197. TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5), weightOK = *track->GetWeight();
  1198. TMatrixD param(5,1), paramTmp(5,1), paramOK = *track->GetParamNew();
  1199. Double_t saveR = 0., rOK = track->GetPosNew(), saveLeng = 0., lengOK = 0.;
  1200. TString saveNode, nodeOK = track->GetNodeNew();
  1201. //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
  1202. // Loop over layers
  1203. Int_t nLay;
  1204. for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  1205. //#pragma omp critical
  1206. {
  1207. //nLay = GetNofHitsInLayer(lay);
  1208. }
  1209. //Int_t indx0 = GetHitsInLayer(lay);
  1210. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Time1");
  1211. Double_t dChi2Min = 1.e+6, padH = lay < MpdTpcSectorGeo::Instance()->NofRowsReg(0) ?
  1212. MpdTpcSectorGeo::Instance()->PadHeight() : MpdTpcSectorGeo::Instance()->PadHeight(1);
  1213. MpdKalmanHit *hitMin = 0x0;
  1214. //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
  1215. Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1, secFirst = -1, isecHit = -1, isecDif = 0;
  1216. if (trackDir == MpdKalmanTrack::kOutward) {
  1217. indxBeg = nLay - 1;
  1218. indxEnd = -1;
  1219. dIndx = -1;
  1220. }
  1221. Int_t isec = 0;
  1222. if (fModular) isec = SectorNo(track->GetNodeNew()); // sector number
  1223. Int_t sector3[3] = {-1, -1, -1};
  1224. MpdTpcKalmanTrack* trackBr[3] = {NULL, NULL, NULL};
  1225. // Loop over hits of 3 sectors (central (where track is located) and +- 1)
  1226. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Hits");
  1227. Int_t i3sec[3] = { 0, -1, 1};
  1228. for (Int_t iadj = 0; iadj < 3; ++iadj) {
  1229. Int_t isecLoop = isec + i3sec[iadj];
  1230. //if (isecLoop < 0) isecLoop += fPadPlane->nSectors();
  1231. //else if (isecLoop > fPadPlane->nSectors() - 1) isecLoop -= fPadPlane->nSectors();
  1232. if (isecLoop < 0) isecLoop += fSecGeo->NofSectors();
  1233. else if (isecLoop > fSecGeo->NofSectors() - 1) isecLoop -= fSecGeo->NofSectors();
  1234. if (fLaySecBegPointers[lay][isecLoop] < 0) continue;
  1235. indxBeg = fLaySecBegPointers[lay][isecLoop];
  1236. indxEnd = fLaySecEndPointers[lay][isecLoop];
  1237. if (trackDir == MpdKalmanTrack::kOutward) {
  1238. Int_t tmp = indxBeg;
  1239. indxBeg = indxEnd - 1;
  1240. indxEnd = tmp - 1;
  1241. }
  1242. MpdKalmanHit *hit = NULL;
  1243. for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
  1244. //MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
  1245. hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx);
  1246. // Exclude used hits
  1247. if (hit->GetFlag() & MpdKalmanHit::kUsed) continue;
  1248. // !!! Exact ID match
  1249. //if (GetPoint(hit)->GetTrackID() != track->GetTrackID()) continue;
  1250. //if (GetTpcHit(hit)->GetTrackID() != track->GetTrackID()) continue;
  1251. // Check if the hit within some window (15x9cm for the moment - check!!!)
  1252. if (!fModular) {
  1253. if (trackBr[0] == NULL) trackBr[0] = new MpdTpcKalmanTrack(*track);
  1254. if (TMath::Abs(hit->GetMeas(1)-trackBr[0]->GetParamNew(1)) > 9) continue;
  1255. if (TMath::Abs(MpdKalmanFilter::Instance()->Proxim(trackBr[0],hit)
  1256. -trackBr[0]->GetParamNew(0)) > 15) continue;
  1257. } else {
  1258. isecHit = hit->GetDetectorID() % 1000000;
  1259. isecDif = isecHit - isec;
  1260. // Track and hit not in adjacent sectors
  1261. //if (TMath::Abs(isecDif) > 1 && TMath::Abs(isecDif) < fPadPlane->nSectors()-1) break;
  1262. //if (TMath::Abs(isecDif) > 1) isecDif -= TMath::Sign(fPadPlane->nSectors(),isecDif);
  1263. if (TMath::Abs(isecDif) > 1 && TMath::Abs(isecDif) < fSecGeo->NofSectors()-1) break;
  1264. if (TMath::Abs(isecDif) > 1) isecDif -= TMath::Sign(fSecGeo->NofSectors(),isecDif);
  1265. ++isecDif;
  1266. if (isecHit == secFirst && sector3[isecDif] < 0) break; // excluded sector
  1267. if (isecHit != secFirst && trackBr[isecDif] == NULL) trackBr[isecDif] = new MpdTpcKalmanTrack(*track);
  1268. if (TMath::Abs(hit->GetMeas(1)-trackBr[isecDif]->GetParamNew(1)) > 9) continue;
  1269. }
  1270. Double_t leng = trackBr[isecDif]->GetLength();
  1271. if (fModular && isecHit != secFirst) {
  1272. // Modular geometry
  1273. secFirst = isecHit;
  1274. sector3[isecDif] = 1;
  1275. if (trackBr[isecDif]->GetNode() == "") trackBr[isecDif]->SetNodeNew("");
  1276. //if (!MpdKalmanFilter::Instance()->PropagateToHit(&tr1,hit,kTRUE,kTRUE) ||
  1277. // TMath::Abs(leng - tr1.GetLength()) < MpdKalmanFilter::Instance()->fgkEpsilon) {
  1278. if (!MpdKalmanFilter::Instance()->PropagateToHit(trackBr[isecDif],hit,kTRUE,kTRUE)) {
  1279. // Not reaching padrow
  1280. sector3[isecDif] = -1;
  1281. //iOK = 0;
  1282. break; //continue;
  1283. }
  1284. //tr1.GetParamNew()->Print();
  1285. //cout << " New node: " << tr1.GetNodeNew() << endl;
  1286. iOK = 1;
  1287. // Check for edge effect
  1288. if (TMath::Abs(trackBr[isecDif]->GetParamNew(0)) > geo->Size(hit).X() + 0.5) {
  1289. // Outside target sector
  1290. sector3[isecDif] = -1;
  1291. break; //continue;
  1292. }
  1293. trackBr[isecDif]->SetNode(trackBr[isecDif]->GetNodeNew());
  1294. } else if (fModular == 0) {
  1295. if (!MpdKalmanFilter::Instance()->PropagateToHit(trackBr[0],hit)) { iOK = 0; break; }
  1296. sector3[0] = 1;
  1297. } // if (fModular && isecHit != secFirst)
  1298. // Add multiple scattering in TPC gas (TPCmixture)
  1299. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Scat");
  1300. Double_t step = trackBr[isecDif]->GetLength() - leng;
  1301. Double_t th = trackBr[isecDif]->GetParamNew(3);
  1302. Double_t cosTh = TMath::Cos(th);
  1303. if (1 && step > 1.e-4) {
  1304. Double_t x0 = 13363.6; // rad. length - TPCMixture
  1305. TMatrixDSym *cov = trackBr[isecDif]->Weight2Cov();
  1306. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(trackBr[isecDif], x0, step);
  1307. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  1308. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1309. (*cov)(3,3) += angle2;
  1310. Int_t iok = 0;
  1311. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1312. trackBr[isecDif]->SetWeight(*cov);
  1313. }
  1314. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Scat");
  1315. Double_t dChi2 = 999999.0;
  1316. if (trackBr[isecDif]->GetLength() - track->GetLength() < padH / TMath::Abs(cosTh) * 1.8 * (miss + 1))
  1317. dChi2 = MpdKalmanFilter::Instance()->FilterHit(trackBr[isecDif],hit,pointWeight,param);
  1318. //cout << lay << " " << isecHit << " " << isec << " " << dChi2 << " " << track->GetTrackID()
  1319. // << " " << GetTpcHit(hit)->GetTrackID() << endl;
  1320. if (dChi2 < dChi2Min) {
  1321. dChi2Min = dChi2;
  1322. hitMin = hit;
  1323. saveWeight = *trackBr[isecDif]->GetWeight();
  1324. saveR = trackBr[isecDif]->GetPosNew();
  1325. saveLeng = trackBr[isecDif]->GetLength();
  1326. saveNode = trackBr[isecDif]->GetNodeNew();
  1327. // temporary storage for the current track
  1328. paramTmp = param;
  1329. pointWeightTmp = pointWeight;
  1330. /*
  1331. cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetZ() << endl;
  1332. cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
  1333. cout << hit->GetRphi() << " " << hit->GetZ() << endl;
  1334. cout << param(0,0) << " " << param(1,0) << endl;
  1335. */
  1336. //paramTmp.Print();
  1337. }
  1338. } // for (Int_t indx = indxBeg; indx != indxEnd;
  1339. // No hits in the sector where the track is located (or its padrow can't be reached)- try both neighbours
  1340. //26-dec-2012 if (fModular && (trackBr[1] == NULL || sector3[1] < 0)) continue;
  1341. if (fModular && trackBr[1] == NULL) continue;
  1342. if (fModular && sector3[1] < 0) {
  1343. if (trackBr[isecDif]->GetParam(0) > 0) ++iadj; // skip one step in the loop
  1344. continue;
  1345. }
  1346. if (!fModular || isecHit < 0 || TMath::Abs(trackBr[isecDif]->GetParamNew(0)) < geo->Size(hit).X() - 1.0) break;
  1347. // Close to sector edge - try adjacent sector
  1348. //i3sec = (trackBr[isecDif]->GetParamNew(0) < 0) ? -1 : 1;
  1349. //if (-trackBr[isecDif]->GetParamNew(0) > 0) ++iadj; // skip one step in the loop
  1350. if (trackBr[isecDif]->GetParamNew(0) > 0) ++iadj; // skip one step in the loop
  1351. } //for (Int_t iadj = 0; iadj < 3; ++iadj)
  1352. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Time1");
  1353. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Hits");
  1354. if (!iOK) {
  1355. // Curling track
  1356. //*
  1357. track->SetPos(rOK); //
  1358. track->SetPosNew(rOK); //
  1359. track->SetParam(paramOK); //
  1360. track->SetParamNew(paramOK); //
  1361. track->SetWeight(weightOK); //
  1362. track->SetLength(lengOK); //
  1363. track->SetNode(nodeOK);
  1364. track->SetNodeNew(nodeOK);
  1365. //*/
  1366. for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
  1367. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1368. return -1;
  1369. } else if (dChi2Min < fgkChi2Cut) {
  1370. //if (fNofEvents > 54)
  1371. //cout << lay << " " << hitMin->GetDetectorID() % 1000000 << " " << isec << " " << dChi2Min << " "
  1372. // << track->GetTrackID() << " " << GetTpcHit(hitMin)->GetTrackID() << " " << saveLeng-track->GetLength() << endl;
  1373. //layOK = lay;
  1374. hitOK = hitMin;
  1375. track->GetHits()->Add(hitOK);
  1376. miss = 0;
  1377. // Restore track params at the best hit
  1378. track->SetChi2(track->GetChi2()+dChi2Min);
  1379. //if (track->GetTrackID() >= 0) cout << " *** Adding hit: " << track->GetTrackID() << " " << GetHitID(hitOK) << " " << hitOK->GetLayer() << " " << hitOK->GetDetectorID() << " "/*<< hitOK->GetTrackID() << " " << hitOK->GetUsage() << " "*/ << dChi2Min << " " << track->GetChi2() << " " << track->GetLength() << " " << saveLeng-track->GetLength() << endl;
  1380. saveWeight += pointWeightTmp;
  1381. track->SetWeight(saveWeight);
  1382. track->SetPos(saveR);
  1383. track->SetPosNew(saveR);
  1384. track->SetLength(saveLeng);
  1385. track->SetParam(paramTmp);
  1386. track->SetParamNew(paramTmp);
  1387. track->SetNode(saveNode);
  1388. track->SetNodeNew(saveNode);
  1389. if (track->GetDirection() == MpdKalmanTrack::kInward) {
  1390. // Save track params at last hit
  1391. track->SetLengAtHit(saveLeng);
  1392. track->SetParamAtHit(paramTmp);
  1393. track->SetWeightAtHit(saveWeight);
  1394. } else {
  1395. // 14-may-12
  1396. track->GetSteps().insert(pair<Int_t,Double_t>(lay,saveLeng));
  1397. }
  1398. // Check if the accepted hit is the same as the seeded hit
  1399. //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
  1400. paramOK = paramTmp;
  1401. rOK = saveR;
  1402. lengOK = saveLeng;
  1403. weightOK = saveWeight;
  1404. nodeOK = saveNode;
  1405. } else {
  1406. ++miss;
  1407. Int_t ibr = -1;
  1408. for (Int_t jj = 0; jj < 3; ++jj) {
  1409. if (sector3[jj] > 0) { ibr = jj; break; }
  1410. }
  1411. if (ibr < 0 || TMath::Abs(trackBr[ibr]->GetLength() - track->GetLength()) > dStepMax) {
  1412. //cout << " !!! Too few branches " << sector3.size() << endl;
  1413. //exit(0);
  1414. track->SetPos(rOK); //
  1415. track->SetPosNew(rOK); //
  1416. track->SetParam(paramOK); //
  1417. track->SetParamNew(paramOK); //
  1418. track->SetWeight(weightOK); //
  1419. track->SetLength(lengOK); //
  1420. //27-04-2012 track->SetNode(nodeOK);
  1421. track->SetNodeNew(nodeOK);
  1422. } else {
  1423. Int_t ip = ibr;
  1424. Double_t leng = 999999.;
  1425. for (Int_t jj = ip ; jj < 3; ++jj) {
  1426. if (sector3[jj] < 0) continue;
  1427. if (trackBr[jj]->GetLength() < leng) { ip = jj; leng = trackBr[ip]->GetLength(); }
  1428. }
  1429. // 12-may-12
  1430. Int_t secPos1 = trackBr[ip]->GetNodeNew().Index("tpcsec");
  1431. Int_t pos1 = trackBr[ip]->GetNodeNew().Index("/",secPos1);
  1432. Int_t secPos2 = track->GetNodeNew().Index("tpcsec");
  1433. Int_t pos2 = track->GetNodeNew().Index("/",secPos2);
  1434. if (trackBr[ip]->GetNodeNew()(secPos1,pos1-secPos1) != track->GetNodeNew()(secPos2,pos2-secPos2)) {
  1435. track->SetPos(rOK); //
  1436. track->SetPosNew(rOK); //
  1437. track->SetParam(paramOK); //
  1438. track->SetParamNew(paramOK); //
  1439. track->SetWeight(weightOK); //
  1440. track->SetLength(lengOK); //
  1441. track->SetNodeNew(nodeOK);
  1442. } else {
  1443. track->SetPos(trackBr[ip]->GetPosNew()); //
  1444. track->SetPosNew(trackBr[ip]->GetPosNew()); //
  1445. track->SetParam(*trackBr[ip]->GetParamNew()); //
  1446. track->SetParamNew(*trackBr[ip]->GetParamNew()); //
  1447. track->SetWeight(*trackBr[ip]->GetWeight()); //
  1448. track->SetLength(trackBr[ip]->GetLength()); //
  1449. track->SetNode(trackBr[ip]->GetNodeNew());
  1450. track->SetNodeNew(trackBr[ip]->GetNodeNew());
  1451. }
  1452. }
  1453. //if (miss > 1) {
  1454. if (miss > 5) {
  1455. /*
  1456. track->SetPos(rOK); //
  1457. track->SetPosNew(rOK); //
  1458. track->SetParam(paramOK); //
  1459. track->SetParamNew(paramOK); //
  1460. */
  1461. //return -1;
  1462. track->SetPos(rOK); //
  1463. track->SetPosNew(rOK); //
  1464. track->SetParam(paramOK); //
  1465. track->SetParamNew(paramOK); //
  1466. track->SetWeight(weightOK); //
  1467. track->SetLength(lengOK); //
  1468. track->SetNodeNew(nodeOK);
  1469. for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
  1470. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1471. return 0; // keep short tracks
  1472. }
  1473. }
  1474. //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
  1475. for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
  1476. } // for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
  1477. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1478. return 0;
  1479. }
  1480. //__________________________________________________________________________
  1481. Int_t MpdTpcKalmanFilter::SectorNo(const char* cpath)
  1482. {
  1483. /// Extract sector number from TGeo path
  1484. TString path(cpath);
  1485. Int_t ipos = path.Index("tpcrow_");
  1486. if (ipos < 0) return -1;
  1487. TString sub = TString(path(ipos-3, 2)); // check for 2 digits
  1488. if (!sub.IsDigit()) sub = TString(path(ipos-2, 1)); // 1 digit
  1489. return sub.Atoi() - 1;
  1490. }
  1491. //__________________________________________________________________________
  1492. MpdTpcHit* MpdTpcKalmanFilter::GetTpcHit(MpdKalmanHit *hit)
  1493. {
  1494. /// Get MpdTpcHit pointer for the Kalman hit
  1495. return (MpdTpcHit*) fHits->UncheckedAt(hit->GetIndex());
  1496. }
  1497. //__________________________________________________________________________
  1498. Int_t MpdTpcKalmanFilter::GetHitID(MpdKalmanHit *hit)
  1499. {
  1500. /// Get hit ID from MpdTpcHit ID
  1501. MpdTpcHit *h = GetTpcHit(hit);
  1502. //AZ 19-07-15 return (fUseMCHit) ? h->GetTrackID() : h->GetRefIndex();
  1503. return h->GetTrackID();
  1504. }
  1505. //__________________________________________________________________________
  1506. void MpdTpcKalmanFilter::GoOutward(MpdTpcKalmanTrack *track)
  1507. {
  1508. /// Propagate track in the outward direction
  1509. TMatrixDSym weight = *track->GetWeight(); // save starting weight matrix
  1510. Int_t ok = RunKalmanFilter(track);
  1511. Int_t lastLay = ((MpdKalmanHit*)track->GetHits()->Last())->GetLayer();
  1512. //cout << ok << " " << track->GetChi2() << " " << track->GetNofHits() << " " << lastLay << " " << ((MpdKalmanHit*)track->GetHits()->Last())->GetMeas(1) << endl;
  1513. track->SetLastLay(lastLay);
  1514. track->SetLength(0.);
  1515. BackTrace(track, weight);
  1516. //cout << track->GetChi2() << " " << track->GetTrackID() << " " << track->GetNofHits() << endl;
  1517. }
  1518. //__________________________________________________________________________
  1519. Bool_t MpdTpcKalmanFilter::BackTrace(MpdTpcKalmanTrack *track, TMatrixDSym &weight0, Int_t iDir, Bool_t corDedx)
  1520. {
  1521. /// Propagate track using already found hits
  1522. track->StartBack();
  1523. if (iDir > 0) track->SetDirection(MpdKalmanTrack::kInward);
  1524. else track->SetDirection(MpdKalmanTrack::kOutward);
  1525. Int_t nHits = track->GetNofHits();
  1526. //28-apr-2012 if (nHits == 1) track->SetWeight(weight0); // restore original weight
  1527. //if (fModular && track->GetNode() == "") track->SetNodeNew(""); // 27-04-2012
  1528. TMatrixD param(5,1);
  1529. TMatrixDSym weight(5), pointWeight(5);
  1530. Int_t ibeg = iDir > 0 ? 0 : nHits-1;
  1531. Int_t iend = iDir > 0 ? nHits-1 : 0;
  1532. Int_t lay0 = -1, lay = -1, ifirst = 1;
  1533. iend += iDir;
  1534. Bool_t ok = kTRUE;
  1535. MpdKalmanHit *hit = 0x0;
  1536. map<Int_t,Double_t>& stepMap = track->GetSteps();
  1537. for (Int_t i = ibeg; i != iend; i+=iDir) {
  1538. //for (Int_t i = 1; i < nHits; ++i) {
  1539. if (i == ibeg && iDir > 0) track->SetLength(0.); // for correct track length
  1540. hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
  1541. lay0 = lay;
  1542. lay = hit->GetLayer();
  1543. /*AZ-100620 - W.t.f. is this???
  1544. if (!fUseMCHit && corDedx && i == ibeg) {
  1545. Double_t cosTh = TMath::Cos(track->GetParam(3));
  1546. hit->SetSignal(hit->GetSignal() * cosTh);
  1547. }
  1548. */
  1549. if (i == ibeg) continue; // skip first hit
  1550. Double_t leng = track->GetLength(), stepBack = -1;
  1551. if (stepMap.find(lay) != stepMap.end() && stepMap.find(lay0) != stepMap.end())
  1552. stepBack = TMath::Abs (stepMap[lay] - stepMap[lay0]);
  1553. ok = MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE, stepBack);
  1554. if (!ok) continue;
  1555. //cout << track->GetLength() << endl;
  1556. // Add multiple scattering in TPC gas (TPCmixture)
  1557. Double_t step = track->GetLength() - leng;
  1558. if (1 && step > 1.e-4) {
  1559. Double_t x0 = 13363.6; // rad. length
  1560. TMatrixDSym *cov = track->Weight2Cov();
  1561. Double_t th = track->GetParamNew(3);
  1562. Double_t cosTh = TMath::Cos(th);
  1563. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, step);
  1564. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  1565. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  1566. (*cov)(3,3) += angle2;
  1567. Int_t iok = 0;
  1568. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  1569. track->SetWeight(*cov);
  1570. }
  1571. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track,hit,pointWeight,param);
  1572. if (ifirst && dChi2 > 1000) {
  1573. ifirst = 0;
  1574. //cout << " !!! Back: " << dChi2 << " " << nHits << " " << hit->GetLayer() << " " << ok << " " << step << " " << track->GetNode() << " " << track->GetNodeNew() << endl;
  1575. }
  1576. track->SetChi2(track->GetChi2()+dChi2);
  1577. weight = *track->GetWeight();
  1578. weight += pointWeight;
  1579. track->SetWeight(weight);
  1580. track->SetParamNew(param);
  1581. //if (i == 0 && iDir > 0) track->SetLength(0.); // for correct track length
  1582. // Correct dE/dx for step length (for clusters) - due to track angles
  1583. /* Angle corrections moved to Dedx task
  1584. if (!fUseMCHit && corDedx) {
  1585. Double_t cosTh = TMath::Cos(track->GetParamNew(3));
  1586. Double_t phi = track->GetParamNew(2);
  1587. //Double_t secPhi = fSecGeo->SectorAngle(fSecGeo->Sector(hit->GetDetectorID()));
  1588. Double_t secPhi = phi; // !!! it seems to work better !!!
  1589. Double_t cosPhi = TMath::Cos(phi-secPhi);
  1590. if (cosPhi < 0.1) { cout << " oops " << secPhi << " " << phi << " " << track->Pt() << endl; cosPhi = 1; }
  1591. hit->SetSignal(hit->GetSignal() * cosTh * cosPhi);
  1592. }
  1593. */
  1594. }
  1595. // Save track params at last hit
  1596. track->SetLengAtHit(track->GetLength());
  1597. track->SetParamAtHit(*track->GetParamNew());
  1598. track->SetWeightAtHit(*track->GetWeight());
  1599. track->GetHits()->UnSort(); // 20-02-2012
  1600. track->GetHits()->Sort(); // 20-02-2012
  1601. return ok;
  1602. }
  1603. //__________________________________________________________________________
  1604. void MpdTpcKalmanFilter::GoOut()
  1605. {
  1606. /// Backpropagate tracks in the outward direction for matching with outer detectors
  1607. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
  1608. Int_t nReco = fTracks->GetEntriesFast(), isec = 0;
  1609. TMatrixDSym tmp;
  1610. TVector3 posLoc, pos;
  1611. for (Int_t i = 0; i < nReco; ++i) {
  1612. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  1613. //track->GetSteps().clear();
  1614. MpdTpcKalmanTrack tr = *track;
  1615. //Int_t id = tr.GetTrackID();
  1616. tr.SetParam(*tr.GetParamAtHit());
  1617. tr.SetParamNew(*tr.GetParamAtHit());
  1618. tr.SetWeight(*tr.GetWeightAtHit());
  1619. tr.SetLength(tr.GetLengAtHit());
  1620. MpdKalmanHit *hit = (MpdKalmanHit*) tr.GetTrHits()->Last();
  1621. tr.SetPos(hit->GetPos());
  1622. if (fModular) {
  1623. tr.SetNode(MpdKalmanFilter::Instance()->GetGeo()->Path(hit->GetDetectorID()));
  1624. tr.SetNodeNew(tr.GetNode());
  1625. //posLoc.SetXYZ(-hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
  1626. posLoc.SetXYZ(hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
  1627. isec = hit->GetDetectorID() % 1000000;
  1628. //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
  1629. fSecGeo->Local2Global(isec, posLoc, pos);
  1630. gGeoManager->cd(tr.GetNodeNew());
  1631. gGeoManager->CdUp();
  1632. Double_t v7[3], v77[3];
  1633. pos.GetXYZ(v7);
  1634. gGeoManager->MasterToLocal(v7,v77);
  1635. tr.SetPos(v77[2]);
  1636. }
  1637. tr.SetPosNew(tr.GetPos());
  1638. //cout << tr.GetTrackID() << " " << tr.GetLengAtHit() << " " << tr.GetPos() << endl;
  1639. Bool_t ok = BackTrace(&tr,tmp,-1,kTRUE);
  1640. //cout << " aaaaa " << tr.Momentum3().Eta() << " " << tr.GetLength() << endl;
  1641. if (ok) {
  1642. MpdKalmanHit hitTmp;
  1643. hitTmp.SetType(MpdKalmanHit::kFixedR);
  1644. if (fModular) {
  1645. // Transform to global frame
  1646. /*
  1647. hit = (MpdKalmanHit*) tr.GetTrHits()->First();
  1648. posLoc.SetXYZ(-hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
  1649. isec = hit->GetDetectorID() % 1000000;
  1650. pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
  1651. */
  1652. //Double_t v7[3] = {tr.GetParamNew(0), tr.GetParamNew(1), tr.GetPosNew()}, v77[3];
  1653. Double_t v7[3] = {-tr.GetParamNew(0), tr.GetParamNew(1), tr.GetPosNew()}, v77[3];
  1654. gGeoManager->cd(tr.GetNodeNew());
  1655. gGeoManager->CdUp();
  1656. gGeoManager->LocalToMaster(v7,v77);
  1657. //MpdKalmanHit hitTmp;
  1658. hitTmp.SetDist(TMath::Sqrt(v77[0]*v77[0]+v77[1]*v77[1]) + 0.2);
  1659. //hitTmp.SetType(MpdKalmanHit::kFixedR);
  1660. } else hitTmp.SetDist(tr.GetPosNew() + 0.2);
  1661. ok = MpdKalmanFilter::Instance()->PropagateToHit(&tr,&hitTmp);
  1662. }
  1663. if (!ok) track->SetFlag(MpdKalmanTrack::kNoOut); //AZ
  1664. track->SetLengAtHit(tr.GetLength());
  1665. track->SetParamAtHit(*tr.GetParamNew());
  1666. track->SetWeightAtHit(*tr.GetWeight());
  1667. track->SetPosAtHit(tr.GetPosNew());
  1668. }
  1669. // Below is just for debugging
  1670. if (0) {
  1671. TClonesArray *tofPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("TOFPoint");
  1672. Int_t nTof = tofPoints->GetEntriesFast();
  1673. //cout << nTof << endl;
  1674. for (Int_t i = 0; i < nReco; ++i) {
  1675. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  1676. MpdTpcKalmanTrack tr = *track;
  1677. //cout << i << " " << tr.GetTrackID() << endl;
  1678. tr.SetParam(*tr.GetParamAtHit());
  1679. tr.SetParamNew(*tr.GetParamAtHit());
  1680. tr.SetWeight(*tr.GetWeightAtHit());
  1681. tr.SetLength(tr.GetLengAtHit());
  1682. tr.SetPos(((MpdKalmanHit*)tr.GetTrHits()->First())->GetPos());
  1683. tr.SetPosNew(tr.GetPos());
  1684. for (Int_t j = 0; j < nTof; ++j) {
  1685. FairMCPoint *p = (FairMCPoint*) tofPoints->UncheckedAt(j);
  1686. if (p->GetTrackID() != tr.GetTrackID()) continue;
  1687. p->Position(pos);
  1688. MpdKalmanHit hitTmp;
  1689. hitTmp.SetDist(pos.Pt());
  1690. hitTmp.SetType(MpdKalmanHit::kFixedR);
  1691. MpdKalmanFilter::Instance()->PropagateToHit(&tr,&hitTmp);
  1692. Double_t dl = tr.GetLength() - p->GetLength();
  1693. //cout << dl << " " << tr.GetPosNew()-pos.Pt() << " " << tr.GetParamNew(1)-pos.Z() << " " << pos.Z() << endl;
  1694. if (lunTpc) fprintf(lunTpc, "%10.3f %10.3f %10.3f %10.3f \n",
  1695. dl, tr.GetPosNew()-pos.Pt(), tr.GetParamNew(1)-pos.Z(), pos.Z());
  1696. break;
  1697. }
  1698. }
  1699. }
  1700. if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
  1701. }
  1702. //__________________________________________________________________________
  1703. void MpdTpcKalmanFilter::RemoveDoubles(TClonesArray *tracks)
  1704. {
  1705. /// Remove double tracks (keep the ones with better quality)
  1706. Int_t ntracks = tracks->GetEntriesFast();
  1707. cout << " Total tracks: " << ntracks << endl;
  1708. // Fill map with track quality
  1709. multimap<Double_t,MpdTpcKalmanTrack*> qualMap;
  1710. multimap<Double_t,MpdTpcKalmanTrack*>::iterator mit1, mit2;
  1711. for (Int_t i = 0; i < ntracks; i++) {
  1712. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) tracks->UncheckedAt(i);
  1713. //AZ-190820 Double_t quality = tr->GetNofHits() + TMath::Min(tr->GetChi2(),999.0) / 1000.0;
  1714. Double_t quality = tr->GetNofHits() + (1 - TMath::Min(tr->GetChi2(),999.0) / 1000.0);
  1715. qualMap.insert(pair<Double_t,MpdTpcKalmanTrack*>(-quality,tr));
  1716. }
  1717. MpdTpcKalmanTrack *tr1, *tr2;
  1718. for (mit1 = qualMap.begin(); mit1 != qualMap.end(); ++mit1) {
  1719. tr1 = mit1->second;
  1720. if (tr1 == NULL) continue; // killed track
  1721. mit2 = mit1;
  1722. ++mit2;
  1723. for ( ; mit2 != qualMap.end(); ++mit2) {
  1724. tr2 = mit2->second;
  1725. if (tr2 == NULL) continue; // killed track
  1726. //Int_t nHitsCommon = GetNofCommonHits(tr1, tr2);
  1727. //if ((float)nHitsCommon / TMath::Min(tr1->GetNofHits(),tr2->GetNofHits()) < 0.5) continue;
  1728. if (!AreTracksDoubles(tr1, tr2)) continue;
  1729. if (mit1->first < mit2->first) { tracks->Remove(tr2); mit2->second = NULL; }
  1730. else {
  1731. if (tr1->GetChi2() <= tr2->GetChi2()) { tracks->Remove(tr2); mit2->second = NULL; }
  1732. else { tracks->Remove(tr1); mit1->second = NULL; break; }
  1733. }
  1734. }
  1735. }
  1736. tracks->Compress();
  1737. fNTracks = tracks->GetEntriesFast();
  1738. }
  1739. //__________________________________________________________________________
  1740. Int_t MpdTpcKalmanFilter::GetNofCommonHits(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2)
  1741. {
  1742. /// Compute number of common hits in 2 tracks
  1743. TObjArray *hits1 = tr1->GetHits(), *hits2 = tr2->GetHits();
  1744. Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast();
  1745. Int_t nCom = 0;//, id = -1;
  1746. for (Int_t i = 0; i < nh1; ++i) {
  1747. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
  1748. for (Int_t j = 0; j < nh2; ++j) {
  1749. MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
  1750. if (hit1 == hit2) {
  1751. ++nCom;
  1752. //id = hit1->GetTrackID();
  1753. break;
  1754. }
  1755. if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam
  1756. }
  1757. }
  1758. /*
  1759. if (nCom > 0) cout << " nCom: " << nCom << " " << nh1 << " " << nh2 << " " << id
  1760. << " " << 1./tr1->GetParamNew(4) << " " << 1./tr2->GetParamNew(4)
  1761. << " " << tr1->GetChi2() << " " << tr2->GetChi2() << endl;
  1762. */
  1763. return nCom;
  1764. }
  1765. //__________________________________________________________________________
  1766. Bool_t MpdTpcKalmanFilter::AreTracksDoubles(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2)
  1767. {
  1768. /// Searching common hits in 2 tracks to determine doubles
  1769. //track1 consists of less count of hits than track2
  1770. MpdKalmanTrack *track1 = tr1, *track2 = tr2;
  1771. if (tr1->GetNofHits() > tr2->GetNofHits())
  1772. track1 = tr2, track2 = tr1;
  1773. Int_t limCommonPoint = (track1->GetNofHits()+1) / 2; // at least so many common hits should be found
  1774. TObjArray *hits1 = track1->GetHits(), *hits2 = track2->GetHits();
  1775. Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast(), nHitsCommon = 0, j = 0;
  1776. for (Int_t i = 0; i < nh1; i++){
  1777. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
  1778. for ( ; j < nh2; j++){
  1779. MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
  1780. // Is the hit common for two tracks?
  1781. if (hit1 == hit2) {
  1782. nHitsCommon++;
  1783. if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
  1784. break;
  1785. }
  1786. if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam
  1787. }
  1788. if (i+limCommonPoint-nHitsCommon > nh1) return kFALSE; // there'll be not enough common hits already
  1789. }
  1790. if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
  1791. return kTRUE;
  1792. }
  1793. //__________________________________________________________________________
  1794. void MpdTpcKalmanFilter::ExcludeHits()
  1795. {
  1796. /// Exclude hits, already used for tracking, from consideration during the next passes
  1797. Int_t nReco = fTracks->GetEntriesFast();
  1798. cout << " nReco: " << nReco << endl;
  1799. for (Int_t i = 0; i < nReco; ++i) {
  1800. MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
  1801. Int_t nhitsKF = track->GetNofHits();
  1802. TObjArray *hits = track->GetHits();
  1803. for (Int_t j = 0; j < nhitsKF; ++j) {
  1804. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
  1805. hit->SetFlag (hit->GetFlag() | MpdKalmanHit::kUsed);
  1806. }
  1807. }
  1808. }
  1809. //__________________________________________________________________________
  1810. Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
  1811. {
  1812. /// Apply momentum correction due to ionization loss in the beam pipe (1 mm of Be)
  1813. /// and TPC inner shell (tpc_v6) - for pions and electrons
  1814. /// Attention!!! Neither pion nor muon with p = 0.050 GeV/c survive - they decay
  1815. /// before entering TPC
  1816. const Double_t piMass = 0.13957, eMass = 0.00051;
  1817. Double_t mass;
  1818. Int_t elec = GetParticleId(id);
  1819. if (elec) mass = eMass;
  1820. else mass = piMass;
  1821. return CorrectForLoss(TMath::Abs(pt), the, mass, 1); // ionization loss correction
  1822. }
  1823. //__________________________________________________________________________
  1824. /*
  1825. Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
  1826. {
  1827. /// Apply momentum correction due to ionization loss in the beam pipe (1 mm of Be)
  1828. /// and TPC inner shell (tpc_v6) - for pions and electrons
  1829. /// Attention!!! Neither pion nor muon with p = 0.050 GeV/c survive - they decay
  1830. /// before entering TPC
  1831. const Int_t nDim = 13;
  1832. const Double_t mom[nDim]= {0.075, 0.100, 0.200, 0.300, 0.400, 0.500, 0.700, 0.900,
  1833. 1.200, 1.500, 2.000, 2.500, 3.500};
  1834. const Double_t tkin[nDim]= {0.019, 0.032, 0.104, 0.191, 0.284, 0.380, 0.574, 0.771,
  1835. 1.069, 1.367, 1.865, 2.364, 3.363};
  1836. const Double_t dedxp[nDim]={8.700, 4.900, 2.350, 1.950, 1.850, 1.850, 1.750, 1.750,
  1837. 1.750, 1.750, 1.850, 1.850, 1.850}; // peak
  1838. const Double_t dedxm[nDim]={8.773, 5.002, 2.597, 2.259, 2.175, 2.143, 2.129, 2.131,
  1839. 2.152, 2.154, 2.182, 2.193, 2.239}; // mean of fit Gauss*Landau
  1840. const Double_t piMass = 0.13957, eMass = 0.00051;
  1841. Double_t dt, mass;
  1842. Int_t elec = GetParticleId(id);
  1843. if (elec) mass = eMass;
  1844. else mass = piMass;
  1845. Double_t mass2 = mass * mass;
  1846. Double_t p = pt / TMath::Sin(the);
  1847. Double_t t = TMath::Sqrt (p*p + mass2) - mass;
  1848. if (elec) dt = 2.9; // 2.9 MeV loss for electrons
  1849. else {
  1850. if (t < tkin[0]) dt = dedxm[0];
  1851. //if (t < tkin[0]) {
  1852. //dt = dedxm[0] + (dedxm[1]-dedxm[0])/(tkin[1]-tkin[0]) * (t-tkin[0]);
  1853. //dt = TMath::Min (dt,10.);
  1854. //}
  1855. else if (t > tkin[nDim-1]) dt = dedxm[nDim-1];
  1856. else {
  1857. // Binary search
  1858. Int_t i1 = 0, i2 = nDim-1, i = i2;
  1859. do {
  1860. i = i1 + (i2-i1) / 2;
  1861. if (t > tkin[i]) i1 = i;
  1862. else i2 = i;
  1863. //cout << i << " " << i1 << " " << i2 << " " << t << " " << tkin[i1] << " " << tkin[i2] << endl;
  1864. } while (i2 - i1 > 1);
  1865. // Linear interpolation
  1866. dt = dedxm[i1] + (dedxm[i2]-dedxm[i1])/(tkin[i2]-tkin[i1]) * (t-tkin[i1]);
  1867. }
  1868. }
  1869. dt /= TMath::Sin(the);
  1870. t += dt * 1.e-3;
  1871. Double_t e = mass + t;
  1872. p = TMath::Sqrt (e*e - mass2);
  1873. pt = p * TMath::Sin(the);
  1874. return pt;
  1875. }
  1876. */
  1877. /*
  1878. //__________________________________________________________________________
  1879. Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
  1880. {
  1881. /// Apply momentum correction due to ionization loss in the beam pipe (0.5 mm of Al)
  1882. /// and TPC inner shell
  1883. //return pt; ///
  1884. const Int_t nDim = 13;
  1885. const Double_t mom[nDim]={0.075,0.100,0.200,0.300,0.400,0.500,0.700,0.900,
  1886. 1.200,1.500,2.000,2.500,3.500};
  1887. const Double_t tkin[nDim]={0.019,0.032,0.104,0.191,0.284,0.380,0.574,0.771,
  1888. 1.069,1.367,1.865,2.364,3.363};
  1889. const Double_t dedxp[nDim]={3.125,1.950,0.975,0.800,0.800,0.750,0.750,0.750,
  1890. 0.750,0.750,0.750,0.775,0.775}; // peak
  1891. const Double_t dedxm[nDim]={3.195,2.085,1.100,0.952,0.918,0.898,0.898,0.890,
  1892. 0.887,0.886,0.901,0.904,0.913}; // mean
  1893. const Double_t piMass = 0.13957, eMass = 0.00051;
  1894. Double_t dt, mass;
  1895. Int_t elec = GetParticleId(id);
  1896. if (elec) mass = eMass;
  1897. else mass = piMass;
  1898. Double_t mass2 = mass * mass;
  1899. Double_t p = pt / TMath::Sin(the);
  1900. Double_t t = TMath::Sqrt (p*p + mass2) - mass;
  1901. if (elec) dt = 1.45; // 1.45 MeV loss for electrons
  1902. else {
  1903. //if (t < tkin[0]) dt = dedxm[0];
  1904. if (t < tkin[0]) {
  1905. dt = dedxm[0] + (dedxm[1]-dedxm[0])/(tkin[1]-tkin[0]) * (t-tkin[0]);
  1906. dt = TMath::Min (dt,10.);
  1907. }
  1908. else if (t > tkin[nDim-1]) dt = dedxm[nDim-1];
  1909. else {
  1910. // Binary search
  1911. Int_t i1 = 0, i2 = nDim-1, i = i2;
  1912. do {
  1913. i = i1 + (i2-i1) / 2;
  1914. if (t > tkin[i]) i1 = i;
  1915. else i2 = i;
  1916. //cout << i << " " << i1 << " " << i2 << " " << t << " " << tkin[i1] << " " << tkin[i2] << endl;
  1917. } while (i2 - i1 > 1);
  1918. // Linear interpolation
  1919. dt = dedxm[i1] + (dedxm[i2]-dedxm[i1])/(tkin[i2]-tkin[i1]) * (t-tkin[i1]);
  1920. }
  1921. }
  1922. dt /= TMath::Sin(the);
  1923. t += dt * 1.e-3;
  1924. Double_t e = mass + t;
  1925. p = TMath::Sqrt (e*e - mass2);
  1926. pt = p * TMath::Sin(the);
  1927. return pt;
  1928. }
  1929. */
  1930. //__________________________________________________________________________
  1931. Int_t MpdTpcKalmanFilter::GetParticleId(Int_t id)
  1932. {
  1933. /// Particle ID:
  1934. /// !!! based on MC information at the moment !!!
  1935. MpdMCTrack* mcTr = (MpdMCTrack*) fMCtracks->UncheckedAt(id);
  1936. Int_t pdg = mcTr->GetPdgCode();
  1937. if (TMath::Abs(pdg) == 11) return 1;
  1938. return 0;
  1939. }
  1940. //__________________________________________________________________________
  1941. //Bool_t MpdTpcKalmanFilter::SameOrigin(TpcLheHit *hit, Int_t idKF, Int_t *mcTracks)
  1942. Bool_t MpdTpcKalmanFilter::SameOrigin(TpcPoint *hit, Int_t idKF, Int_t *mcTracks)
  1943. {
  1944. /// Check if two hits belong to tracks from the same mother
  1945. //cout << " Wrongs: " << hit->GetTrackID() << " " << idKF << " " << mcTracks[hit->GetTrackID()] << " " << mcTracks[idKF] << endl;
  1946. if (mcTracks[hit->GetTrackID()] <= 0 && mcTracks[idKF] <= 0) return kFALSE; // both hits are from primaries
  1947. else if (mcTracks[hit->GetTrackID()] > 0) {
  1948. // Check mother
  1949. MpdMCTrack *track = (MpdMCTrack*) fMCtracks->UncheckedAt(hit->GetTrackID());
  1950. while (track->GetMotherId() > 0) {
  1951. //cout << " Wrongs1: " << track->GetMotherId() << endl;
  1952. if (track->GetMotherId() == idKF) return kTRUE;
  1953. track = (MpdMCTrack*) fMCtracks->UncheckedAt(track->GetMotherId());
  1954. }
  1955. } else {
  1956. MpdMCTrack *track = (MpdMCTrack*) fMCtracks->UncheckedAt(idKF);
  1957. while (track->GetMotherId() > 0) {
  1958. //cout << " Wrongs2: " << track->GetMotherId() << endl;
  1959. if (track->GetMotherId() == hit->GetTrackID()) return kTRUE;
  1960. track = (MpdMCTrack*) fMCtracks->UncheckedAt(track->GetMotherId());
  1961. }
  1962. }
  1963. return kFALSE;
  1964. }
  1965. //__________________________________________________________________________
  1966. Bool_t MpdTpcKalmanFilter::Refit(MpdKalmanTrack *track, Double_t mass, Int_t charge, Bool_t skip, Int_t iDir,
  1967. Bool_t exclude, std::map<Double_t,matrix4> *cache)
  1968. // Bool_t exclude, std::map<Double_t,std::tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> > *cache)
  1969. {
  1970. /// Refit track in TPC using track hits (toward beam line) for some
  1971. /// particle mass and charge hypothesis
  1972. //#pragma omp critical
  1973. {
  1974. track->GetHits()->UnSort();
  1975. track->GetHits()->Sort();
  1976. }
  1977. track->SetChi2(0.);
  1978. if (iDir == 1) track->SetDirection(MpdKalmanTrack::kInward);
  1979. //if (GetNofHits() == 1) return;
  1980. track->SetWeight(*track->GetWeightAtHit());
  1981. TMatrixDSym *w = track->GetWeight();
  1982. Int_t nHits = track->GetNofHits(), nHits2 = nHits * nHits;;
  1983. for (Int_t i = 0; i < 5; ++i) {
  1984. for (Int_t j = i; j < 5; ++j) {
  1985. //if (j == i) (*w)(i,j) /= 1000.;
  1986. if (j == i) (*w)(i,j) /= nHits2;
  1987. //if (j == i) (*w)(i,j) /= 1;
  1988. else (*w)(i,j) = (*w)(j,i) = 0.;
  1989. }
  1990. }
  1991. TMatrixD param(5,1);
  1992. TMatrixDSym weight(5), pointWeight(5);
  1993. //Int_t iDir = 1;
  1994. Int_t ibeg = (iDir > 0) ? 0 : nHits-1;
  1995. Int_t iend = (iDir > 0) ? nHits-1 : 0;
  1996. iend += iDir;
  1997. MpdKalmanHit *hit = 0x0;
  1998. track->SetPosNew(track->GetPosAtHit()+0.2);
  1999. track->SetParamNew(*track->GetParamAtHit());
  2000. track->SetNodeNew("");
  2001. Double_t leng = 0.0;
  2002. track->SetLength(leng);
  2003. TString mass2 = "";
  2004. mass2 += mass * mass;
  2005. for (Int_t i = ibeg; i != iend; i+=iDir) {
  2006. hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
  2007. if (hit->GetUniqueID()) continue; // 5-jan-2018 - skip ITS hits
  2008. //if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE)) return kFALSE;
  2009. if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE)) {
  2010. if (fVerbose > 1) cout << " III " << i << " " << hit->GetDist() << " " << track->Pt() << " " << track->GetDirection() << endl;
  2011. return kFALSE;
  2012. }
  2013. Double_t step = track->GetLength() - leng;
  2014. Double_t th = track->GetParamNew(3);
  2015. Double_t cosTh = TMath::Cos(th);
  2016. if (step > 1.e-4) {
  2017. Double_t x0 = 13363.6; // rad. length - TPCMixture
  2018. TMatrixDSym *cov = track->Weight2Cov();
  2019. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, step, mass2, charge);
  2020. //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
  2021. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  2022. (*cov)(3,3) += angle2;
  2023. Int_t iok = 0;
  2024. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  2025. track->SetWeight(*cov);
  2026. }
  2027. Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hit, pointWeight, param);
  2028. leng = track->GetLength();
  2029. if (exclude && dChi2 > fgkChi2Cut) {
  2030. track->GetHits()->Remove(hit);
  2031. continue; // too high Chi2
  2032. }
  2033. track->SetChi2(track->GetChi2()+dChi2);
  2034. weight = *track->GetWeight();
  2035. weight += pointWeight;
  2036. if (cache) {
  2037. // For smoother
  2038. TMatrixD fc(MpdKalmanFilter::Instance()->GetJacob(),TMatrixD::kMult,*track->GetWeight()); // DW
  2039. (*cache)[track->GetPosNew()] = matrix4(*track->GetParamNew(),param,fc,weight);
  2040. }
  2041. track->SetWeight(weight);
  2042. track->SetParamNew(param);
  2043. //cout << i << " " << dChi2 << " " << 1./track->GetParamNew(4) << endl;
  2044. }
  2045. if (exclude) track->GetHits()->Compress();
  2046. if (iDir < 0 || skip) return kTRUE; // going outward or do not go to the beam line
  2047. // Go to the beam line
  2048. const Int_t nR = 7;
  2049. const Double_t rad[nR] = {27.00, 27.01, 27.31, 27.32, 33.82, 33.83, 34.13};
  2050. 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};
  2051. Bool_t ok = 0;
  2052. MpdKalmanHit hitK;
  2053. hitK.SetType(MpdKalmanHit::kFixedR);
  2054. for (Int_t j = nR-1; j >= 0; --j) {
  2055. hitK.SetPos(rad[j]);
  2056. ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hitK,kTRUE);
  2057. // Add multiple scattering
  2058. TMatrixDSym *cov = track->Weight2Cov();
  2059. Double_t th = track->GetParamNew(3);
  2060. Double_t cosTh = TMath::Cos(th);
  2061. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[j], mass2, charge);
  2062. (*cov)(2,2) += (angle2 / cosTh / cosTh);
  2063. (*cov)(3,3) += angle2;
  2064. //cov->Print();
  2065. Int_t iok = 0;
  2066. MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
  2067. track->SetWeight(*cov);
  2068. }
  2069. //if (!ok) continue;
  2070. /* TMatrixDSym *cov = */track->Weight2Cov();
  2071. // Update track
  2072. track->SetParam(*track->GetParamNew());
  2073. track->SetPos(track->GetPosNew());
  2074. // Correct for energy loss
  2075. //*
  2076. Double_t theta = TMath::PiOver2() - track->GetParam(3);
  2077. Double_t pt = 1. / track->GetParam(4);
  2078. Double_t ptCor = CorrectForLoss(TMath::Abs(pt), theta, mass, charge); // ionization loss correction
  2079. track->SetParam(4,TMath::Sign(1./ptCor,pt));
  2080. track->SetParamNew(*track->GetParam());
  2081. //*/
  2082. //*
  2083. hitK.SetPos(0.);
  2084. hitK.SetMeas(0,track->GetParam(2));
  2085. //hit.SetRphi(track->GetParam(2)); // track Phi - check if it is necessary !!!!!!!
  2086. Double_t pos = track->GetPosNew();
  2087. MpdKalmanFilter::Instance()->PropagateToHit(track, &hitK, kTRUE);
  2088. track->SetPos(pos); // save position after passing inner shell
  2089. track->SetParam(*track->GetParamNew()); // !!! track params at PCA
  2090. //track->SetLengAtHit(track->GetLength()-track->GetLengAtHit()); // track length from PCA to the nearest hit
  2091. if (0) {
  2092. TVector3 pos3 = MpdKalmanFilter::Instance()->GetGeo()->GlobalPos(hit);
  2093. cout << " Hit: " << pos3.Pt()*pos3.Phi() << " " << pos3.Pt() << " " << pos3.Z() << endl;
  2094. cout << " Track: " << nHits << " " << track->GetParamNew(0) << " " << track->GetPosNew() << " " << track->GetParamNew(1) << " " << track->GetChi2() << endl;
  2095. }
  2096. return kTRUE;
  2097. }
  2098. //__________________________________________________________________________
  2099. Double_t MpdTpcKalmanFilter::CorrectForLossFluct(Double_t pt, Double_t the, Double_t mass, Int_t charge)
  2100. {
  2101. // Ionization loss correction fluctuations - only for low momentum pions
  2102. // Pion
  2103. const Int_t nP = 4, nThe = 5;
  2104. const Double_t moms[nP] = {0.1, 0.2, 0.3, 0.4}; // p
  2105. const Double_t thes[nThe] = {0.015, 0.315, 0.615, 0.915, 1.215}; // dip angle
  2106. static vector<vector<Double_t> > sigmas(nP,vector<Double_t>(nThe));
  2107. static Int_t first = 1;
  2108. if (first) {
  2109. first = 0;
  2110. sigmas[0][0] = 3.2; sigmas[0][1] = 3.3; sigmas[0][2] = 3.6; sigmas[0][3] = 5.3; sigmas[0][4] = 10.7;
  2111. sigmas[1][0] = 2.8; sigmas[1][1] = 2.4; sigmas[1][2] = 2.1; sigmas[1][3] = 2.3; sigmas[1][4] = 2.8;
  2112. sigmas[2][0] = 1.6; sigmas[2][1] = 1.7; sigmas[2][2] = 1.7; sigmas[2][3] = 1.5; sigmas[2][4] = 2.0;
  2113. sigmas[3][0] = 2.5; sigmas[3][1] = 2.4; sigmas[3][2] = 2.2; sigmas[3][3] = 1.8; sigmas[3][4] = 1.9;
  2114. for (Int_t i = 0; i < nP; ++i) {
  2115. for (Int_t j = 0; j < nThe; ++j) sigmas[i][j] *= 1.e-4;
  2116. }
  2117. }
  2118. if (mass < 0.1 || mass > 0.2) return 0;
  2119. Double_t p = pt / TMath::Sin(the) * charge;
  2120. if (p > 0.4) return 0;
  2121. Double_t e = TMath::Sqrt (p*p + mass*mass), sigt = 0.0;
  2122. sigt = Interp2d(moms, thes, sigmas, p, TMath::Abs(TMath::PiOver2()-the));
  2123. //fprintf(lunTpc,"%f %f %f\n",p,TMath::Abs(TMath::PiOver2()-the),sigt);
  2124. sigt /= TMath::Sin(the);
  2125. return e * sigt / p / p / p / TMath::Sin(the); // Delta (1/pt)
  2126. }
  2127. //__________________________________________________________________________
  2128. Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Double_t mass, Int_t charge)
  2129. {
  2130. // Ionization loss correction
  2131. // Pion
  2132. const Int_t nPi = 13; // pi-
  2133. //const Double_t momPi[nPi] = {0.075, 0.100, 0.200, 0.300, 0.400, 0.500, 0.700, 0.900,
  2134. // 1.200, 1.500, 2.000, 2.500, 3.500};
  2135. const Double_t momPi[nPi] = {0.055, 0.091, 0.197, 0.297, 0.398, 0.498, 0.698, 0.898,
  2136. 1.198, 1.498, 1.998, 2.498, 3.498};
  2137. //const Double_t dedxPi[nPi] = {8.773, 5.002, 2.597, 2.259, 2.175, 2.143, 2.129, 2.131,
  2138. const Double_t dedxPi[nPi] = {8.476, 5.050, 2.625, 2.304, 2.175, 2.143, 2.129, 2.131,
  2139. 2.152, 2.154, 2.182, 2.193, 2.239}; // mean of fit Gauss*Landau
  2140. //const Double_t sigmaPi[nPi] = {4.2e-4, 3.4e-4, 2.2e-4, 2.0e-4}; // gaussian fit
  2141. // Kaon
  2142. const Int_t nK = 12; // K-
  2143. //const Double_t momK[nK] = {0.200, 0.250, 0.300, 0.400, 0.500, 0.700, 0.900, 1.200,
  2144. // 1.500, 2.000, 2.500, 3.500};
  2145. const Double_t momK[nK] = {0.167, 0.233, 0.289, 0.393, 0.495, 0.697, 0.897, 1.198,
  2146. 1.498, 1.998, 2.498, 3.498};
  2147. const Double_t dedxK[nK] = {11.40, 7.584, 5.858, 4.128, 3.307, 2.635, 2.394, 2.253,
  2148. 2.314, 2.256, 2.253, 2.249}; // mean of fit Gauss*Landau
  2149. // Proton
  2150. const Int_t nP = 11;
  2151. //const Double_t momP[nP] = {0.300, 0.350, 0.400, 0.500, 0.700, 0.900, 1.200, 1.500, 2.000,
  2152. // 2.500, 3.500};
  2153. const Double_t momP[nP] = {0.239, 0.313, 0.375, 0.485, 0.693, 0.895, 1.196, 1.498, 1.997,
  2154. 2.498, 3.498};
  2155. const Double_t dedxP[nP] = {16.87, 12.31, 9.751, 6.948, 4.460, 3.483, 2.808, 2.508, 2.295,
  2156. 2.221, 2.164}; // mean of fit Gauss*Landau
  2157. // Deuteron
  2158. const Int_t nD = 10;
  2159. //const Double_t momD[nD] = {0.450, 0.500, 0.600, 0.700, 0.900, 1.200, 1.500, 2.000,
  2160. // 2.500, 3.500};
  2161. const Double_t momD[nD] = {0.303, 0.406, 0.547, 0.666, 0.882, 1.190, 1.493, 1.996,
  2162. 2.497, 3.497};
  2163. const Double_t dedxD[nD] = {28.85, 22.02, 15.40, 11.73, 7.962, 5.369, 4.159, 3.178,
  2164. 2.705, 2.298}; // mean of fit Gauss*Landau
  2165. // Triton
  2166. const Int_t nT = 10;
  2167. //const Double_t momT[nT] = {0.600, 0.650, 0.700, 0.800, 1.000, 1.200, 1.500, 2.000,
  2168. // 2.500, 3.500};
  2169. const Double_t momT[nT] = {0.400, 0.512, 0.594, 0.730, 0.962, 1.176, 1.485, 1.992,
  2170. 2.494, 3.496};
  2171. const Double_t dedxT[nT] = {35.02, 27.97, 23.83, 18.43, 12.58, 9.341, 6.835, 4.730,
  2172. 3.662, 2.814}; // mean of fit Gauss*Landau
  2173. // He3
  2174. const Int_t nHe3 = 11;
  2175. //const Double_t momHe3[nHe3] = {0.850, 0.900, 1.000, 1.100, 1.200, 1.400, 1.600, 2.000, 2.500,
  2176. // 3.000, 3.500};
  2177. const Double_t momHe3[nHe3] = {0.455, 0.623, 0.819, 0.967, 1.096, 1.330, 1.548, 1.967, 2.478,
  2178. 2.983, 3.485};
  2179. const Double_t dedxHe3[nHe3] = {89.13, 72.42, 55.85, 45.98, 39.47, 30.53, 25.15, 19.07, 14.73,
  2180. 12.76, 11.42}; // mean of fit Gauss*Landau
  2181. // He4
  2182. const Int_t nHe4 = 9;
  2183. //const Double_t momHe4[nHe4] = {1.050, 1.100, 1.200, 1.400, 1.600, 2.000, 2.500,
  2184. // 3.000, 3.500};
  2185. const Double_t momHe4[nHe4] = {0.606, 0.759, 0.962, 1.256, 1.500, 1.941, 2.464,
  2186. 2.974, 3.479};
  2187. const Double_t dedxHe4[nHe4] = {96.18, 82.70, 66.24, 48.44, 38.43, 27.48, 20.09,
  2188. 16.48, 14.14}; // mean of fit Gauss*Landau
  2189. charge = TMath::Abs(charge);
  2190. Double_t p = pt / TMath::Sin(the) * charge, mass2 = mass * mass;
  2191. Double_t t = TMath::Sqrt (p*p + mass2) - mass, dt = 0.0;
  2192. if (mass < 0.1) dt = 2.9; // 2.9 MeV loss for electrons
  2193. else if (mass < 0.2) {
  2194. if (p < 0.45) {
  2195. dt = 0.0152867 - 0.0905572*p - 1.98303*p*p + 35.3591*p*p*p - 245.426*p*p*p*p + 915.386*p*p*p*p*p
  2196. - 1926.86*p*p*p*p*p*p + 2157.13*p*p*p*p*p*p*p - 999.496*p*p*p*p*p*p*p*p;
  2197. dt *= 1000;
  2198. } else dt = MpdKalmanFilter::Instance()->Interp(nPi, momPi, dedxPi, p); // Pion
  2199. }
  2200. else if (mass < 0.6) dt = MpdKalmanFilter::Instance()->Interp(nK, momK, dedxK, p); // Kaon
  2201. else if (mass < 1.0) dt = MpdKalmanFilter::Instance()->Interp(nP, momP, dedxP, p); // Proton
  2202. else if (mass < 2.0) dt = MpdKalmanFilter::Instance()->Interp(nD, momD, dedxD, p); // Deuteron
  2203. else if (mass < 2.9 && charge == 1) dt = MpdKalmanFilter::Instance()->Interp(nT, momT, dedxT, p); // Triton
  2204. else if (mass < 2.9) dt = MpdKalmanFilter::Instance()->Interp(nHe3, momHe3, dedxHe3, p); // He3
  2205. else if (mass < 3.9) dt = MpdKalmanFilter::Instance()->Interp(nHe4, momHe4, dedxHe4, p); // He4
  2206. dt /= TMath::Sin(the);
  2207. t += dt * 1.e-3;
  2208. Double_t e = mass + t;
  2209. p = TMath::Sqrt (e*e - mass2);
  2210. pt = p * TMath::Sin(the);
  2211. if (charge) pt /= charge;
  2212. return pt;
  2213. }
  2214. //__________________________________________________________________________
  2215. void MpdTpcKalmanFilter::MergeTracks(Int_t ipass)
  2216. {
  2217. // Merge tracks
  2218. multimap<Int_t,Int_t> secMaps[24]; // first hit row# - track index
  2219. Int_t nReco = fTracks->GetEntriesFast();
  2220. for (Int_t i = 0; i < nReco; ++i) {
  2221. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  2222. MpdKalmanHit *hit = (MpdKalmanHit*) track->GetTrHits()->Last();
  2223. Int_t isec = fSecGeo->Sector(hit->GetDetectorID());
  2224. if (hit->GetMeas(1) < 0) isec += 12;
  2225. Int_t lay = hit->GetLayer();
  2226. secMaps[isec].insert(pair<Int_t,Int_t>(lay,i));
  2227. }
  2228. // Loop over tracks and try to merge them: at first within one sector,
  2229. // then in adjacent sectors
  2230. TVector3 mom;
  2231. multimap<Int_t,Int_t>::iterator it, it1;
  2232. // Loop over sectors
  2233. for (Int_t isec = 0; isec < 24; ++isec) {
  2234. //cout << " Sector No: " << isec << ", tracks: " << secMaps[isec].size() << endl;
  2235. multimap<Int_t,Int_t> ids;
  2236. Int_t isec1 = isec;
  2237. if (ipass == 2) {
  2238. // Merge sectors
  2239. isec1 = isec + 1;
  2240. if (isec1 == 12) isec1 = 0;
  2241. else if (isec1 == 24) isec1 = 12;
  2242. }
  2243. // Loop over tracks in the sector
  2244. for (it = secMaps[isec].begin(); it != secMaps[isec].end(); ++it) {
  2245. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(it->second);
  2246. if (tr->GetTrackID() < 0) continue; // already processed
  2247. Int_t hitLays[70] = {0}, nh = tr->GetNofTrHits();
  2248. if (nh >= 50) continue; // exclude long tracks
  2249. if (ipass == 0 && nh < 10) continue; // at first exclude short tracks
  2250. TClonesArray *hits = tr->GetTrHits();
  2251. for (Int_t ih = 0; ih < nh; ++ih) {
  2252. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(ih);
  2253. hitLays[hit->GetLayer()] = ih + 1;
  2254. }
  2255. // Loop over tracks in the sector
  2256. it1 = it;
  2257. if (isec1 == isec) ++it1;
  2258. else it1 = secMaps[isec1].begin();
  2259. for ( ; it1 != secMaps[isec1].end(); ++it1) {
  2260. //if (it1 == it) continue;
  2261. MpdTpcKalmanTrack *tr1 = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(it1->second);
  2262. if (tr1->GetTrackID() < 0) continue; // already processed
  2263. Int_t hitLays1[70] = {0}, nh1 = tr1->GetNofTrHits(), nSameLay = 0;
  2264. if (nh1 >= 50) continue; // exclude long tracks
  2265. if (ipass == 0 && nh1 < 10) continue; // at first exclude short tracks
  2266. //if (TMath::Abs(tr->Momentum3().Eta()) < 1.5) {
  2267. if (TMath::Abs(tr->Momentum3().Eta()) < 1.3) {
  2268. if (nh < 10 && nh1 < 10) continue; // both tracks too short
  2269. } else if (nh < 5 && nh1 < 5) continue; // both tracks too short
  2270. TClonesArray *hits1 = tr1->GetTrHits();
  2271. // Do not merge tracks with too different angles or going in the same direction or having different charges
  2272. Double_t dphi = tr->GetParam(2) - tr1->GetParam(2);
  2273. if (TMath::Abs(dphi) > TMath::Pi()) dphi -= TMath::Sign (TMath::TwoPi(),1.*dphi);
  2274. if (ipass == 0 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.05 ||
  2275. tr->GetDirection() + tr1->GetDirection() != 1 ||
  2276. TMath::Abs(dphi) > 0.5 || tr1->Charge() - tr->Charge() != 0)) continue;
  2277. else if (ipass == 1 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.1 ||
  2278. tr->GetDirection() + tr1->GetDirection() != 1 ||
  2279. TMath::Abs(dphi) > 0.75)) continue;
  2280. else if (ipass == 2 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.1 ||
  2281. tr->GetDirection() + tr1->GetDirection() != 1 ||
  2282. TMath::Abs(dphi) > 0.75)) continue;
  2283. for (Int_t ih = 0; ih < nh1; ++ih) {
  2284. MpdKalmanHit *hit = (MpdKalmanHit*) hits1->UncheckedAt(ih);
  2285. hitLays1[hit->GetLayer()] = ih + 1;
  2286. if (hitLays[hit->GetLayer()]) nSameLay += 1; // number of hits on the same rows
  2287. }
  2288. // If there are hits on the same rows, compute average distance between them
  2289. if (nSameLay) {
  2290. // Compute average distance
  2291. Double_t dist = 0;
  2292. TVector3 posLoc, pos1, pos2;
  2293. for (Int_t ih = 0; ih < 70; ++ih) {
  2294. if (hitLays[ih] == 0 || hitLays1[ih] == 0) continue;
  2295. MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(hitLays[ih]-1);
  2296. MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(hitLays1[ih]-1);
  2297. if (isec1 == isec) {
  2298. Double_t distx = hit->GetMeas(0) - hit1->GetMeas(0);
  2299. distx *= distx;
  2300. Double_t errx2 = 1; //hit->GetErr(0) * hit->GetErr(0) + hit1->GetErr(0) * hit1->GetErr(0);
  2301. Double_t distz = hit->GetMeas(1) - hit1->GetMeas(1);
  2302. distz *= distz;
  2303. Double_t errz2 = 1; //hit->GetErr(1) * hit->GetErr(1) + hit1->GetErr(1) * hit1->GetErr(1);
  2304. //chi2 += (distx / errx2 + distz / errz2);
  2305. dist += TMath::Sqrt (distx / errx2 + distz / errz2);
  2306. } else {
  2307. posLoc.SetXYZ(hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
  2308. fSecGeo->Local2Global(isec, posLoc, pos1);
  2309. posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  2310. fSecGeo->Local2Global(isec1, posLoc, pos2);
  2311. pos2 -= pos1;
  2312. dist += pos2.Mag();
  2313. }
  2314. }
  2315. dist /= nSameLay;
  2316. if (dist > 5) continue; // do not merge tracks
  2317. }
  2318. // Try to merge tracks
  2319. // Take track with more hits and try to add the second one
  2320. MpdKalmanHit *hit, *hits2[70] = {NULL};
  2321. MpdTpcKalmanTrack *trL = tr, *trS = tr1;
  2322. Int_t *hitLaysL = hitLays, *hitLaysS = hitLays1;
  2323. if (nh < nh1) { trL = tr1; trS = tr; hitLaysL = hitLays1; hitLaysS = hitLays; }
  2324. for (Int_t ih = 0; ih < 70; ++ih) {
  2325. if (hitLaysL[ih]) hit = (MpdKalmanHit*) trL->GetTrHits()->UncheckedAt(hitLaysL[ih]-1);
  2326. else if (hitLaysS[ih]) hit = (MpdKalmanHit*) trS->GetTrHits()->UncheckedAt(hitLaysS[ih]-1);
  2327. else hit = NULL;
  2328. hits2[ih] = hit;
  2329. }
  2330. //cout << " ---------------- " << trL->GetNofTrHits() << " " << trS->GetNofTrHits() << endl;
  2331. MpdTpcKalmanTrack track = *trL;
  2332. MpdKalmanHit hitTmp;
  2333. hitTmp.SetType(MpdKalmanHit::kFixedR);
  2334. //if (track.GetDirection() == MpdKalmanTrack::kInward) MpdKalmanFilter::Instance()->Refit(&track, 1);
  2335. if (track.GetDirection() == MpdKalmanTrack::kInward) {
  2336. track.SetDirection(MpdKalmanTrack::kOutward);
  2337. MpdKalmanFilter::Instance()->Refit(&track, -1);
  2338. }
  2339. //else if (track.GetDirection() == MpdKalmanTrack::kOutward) MpdKalmanFilter::Instance()->Refit(&track, -1);
  2340. //cout << " Refit chi2: " << track.GetChi2() << " " << track.GetNofHits() << " " << trL->GetDirection() << endl;
  2341. if (track.GetChi2() > 1000) continue; // do not merge tracks
  2342. Int_t ok = 0;
  2343. if (trL->GetDirection() == MpdKalmanTrack::kInward) {
  2344. //hitTmp.SetDist(fSecGeo->GetMaxY()/TMath::Cos(fSecGeo->Dphi()/2));
  2345. Double_t radMax = TMath::Max (((MpdKalmanHit*)trL->GetHits()->First())->GetDist(),
  2346. ((MpdKalmanHit*)trS->GetHits()->First())->GetDist());
  2347. radMax += fSecGeo->GetMinY() + 5.0;
  2348. hitTmp.SetDist(radMax);
  2349. ok = MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
  2350. } else {
  2351. hitTmp.SetDist(fSecGeo->GetMinY());
  2352. ok = MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
  2353. }
  2354. if (!ok) continue; // do not merge tracks
  2355. // Save info for further use
  2356. track.SetWeightAtHit(*track.GetWeight());
  2357. track.SetParamAtHit(*track.GetParamNew());
  2358. track.SetPosAtHit(track.GetPosNew());
  2359. TObjArray *hitss = track.GetHits();
  2360. hitss->Clear();
  2361. for (Int_t ih = 0; ih < 70; ++ih) {
  2362. if (hits2[ih]) hitss->Add(hits2[ih]);
  2363. //if (hits2[ih]) { hitss->Add(hits2[ih]); cout << " dist: " << hits2[ih]->GetLayer() << " " << hits2[ih]->GetDist() << " "
  2364. // << hits2[ih]->GetMeas(0) << " " << hits2[ih]->GetMeas(1) << endl; }
  2365. }
  2366. // Refit with multiple scattering
  2367. if (trL->GetDirection() == MpdKalmanTrack::kInward) {
  2368. if (!Refit(&track, 0.13957, 1, kTRUE, 1, kTRUE)) continue; // do not merge tracks
  2369. //cout << " *** RefitIn ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
  2370. // << " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
  2371. if (1.0*(track.GetNofHits() - trL->GetNofTrHits()) / trS->GetNofTrHits() < 0.3) continue; // do not merge tracks
  2372. } else if (trL->GetDirection() == MpdKalmanTrack::kOutward) {
  2373. // Refit outward
  2374. if (!Refit(&track, 0.13957, 1, kTRUE, -1, kTRUE) || track.GetNofHits() == 0) continue; // do not merge tracks
  2375. //cout << " *** RefitOut ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
  2376. //<< " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
  2377. //hitTmp.SetDist(fSecGeo->GetMaxY()/TMath::Cos(fSecGeo->Dphi()/2));
  2378. hitTmp.SetDist(((MpdKalmanHit*)track.GetHits()->First())->GetDist() + fSecGeo->GetMinY() + 5.0);
  2379. MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
  2380. track.SetPosAtHit(track.GetPosNew());
  2381. track.SetWeightAtHit(*track.GetWeight());
  2382. track.SetParamAtHit(*track.GetParamNew());
  2383. track.SetDirection(MpdKalmanTrack::kInward);
  2384. Refit(&track, 0.13957, 1, kTRUE);
  2385. //cout << " *** RefitOut1 ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
  2386. // << " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
  2387. if (1.0*(track.GetNofHits() - trL->GetNofTrHits()) / trS->GetNofTrHits() < 0.3) continue; // do not merge tracks
  2388. }
  2389. // Merge tracks
  2390. //cout << track.GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetNofWrong() << " " << endl;
  2391. MpdTpcKalmanTrack *newTrack = new ((*fTracks)[fTracks->GetEntriesFast()]) MpdTpcKalmanTrack(track);
  2392. //newTrack->SetLengAtHit(saveLeng);
  2393. newTrack->SetParamAtHit(*track.GetParamNew());
  2394. newTrack->SetWeightAtHit(*track.GetWeight());
  2395. newTrack->SetLengAtHit(track.GetLength());
  2396. AddHits(fTracks->GetEntriesFast()-1);
  2397. //cout << newTrack->GetNofTrHits() << " " << newTrack->GetNofWrong() << " " << newTrack->GetNofHits() << endl;
  2398. trS->SetTrackID(-trS->GetTrackID()-1);
  2399. trL->SetTrackID(-trL->GetTrackID()-1);
  2400. break;
  2401. } // for (it1 = it; it1 != secMaps[isec1].end();
  2402. } // for (it = secMaps[isec].begin(); it != secMaps[isec].end();
  2403. } // for (Int_t isec = 0; isec < 24;
  2404. for (Int_t i = 0; i < nReco; ++i) {
  2405. MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
  2406. if (track->GetTrackID() < 0) fTracks->Remove(track);
  2407. }
  2408. fTracks->Compress();
  2409. }
  2410. //__________________________________________________________________________
  2411. Double_t MpdTpcKalmanFilter::Interp2d(const Double_t *moms, const Double_t *thes, vector<vector<Double_t> > &sigmas,
  2412. Double_t p, Double_t dip)
  2413. {
  2414. // 2-d linear interpolation
  2415. Int_t nx = sigmas.size(), ny = sigmas[0].size(), i0 = 0, j0 = 0, i1 = 0, j1 = 0;
  2416. for ( ; i1 < nx; ++i1) if (moms[i1] >= p) break;
  2417. if (i1 == nx) --i1;
  2418. if (i1 == 0) i0 = 1;
  2419. else i0 = i1 - 1;
  2420. Double_t u = (p - moms[i0]) / (moms[i1] - moms[i0]);
  2421. for ( ; j1 < ny; ++j1) if (thes[j1] >= dip) break;
  2422. if (j1 == ny) --j1;
  2423. if (j1 == 0) j0 = 1;
  2424. else j0 = j1 - 1;
  2425. Double_t v = (dip - thes[j0]) / (thes[j1] - thes[j0]);
  2426. Double_t uv = u * v;
  2427. Double_t res = sigmas[i0][j0] * (1 - u - v + uv) + sigmas[i1][j0] * (u - uv) +
  2428. sigmas[i0][j1] * (v - uv) + sigmas[i1][j1] * uv;
  2429. return res;
  2430. }
  2431. //__________________________________________________________________________
  2432. void MpdTpcKalmanFilter::Smooth(MpdTpcKalmanTrack *track, std::vector<std::pair<Double_t,TMatrixD> >& vecSmooth)
  2433. {
  2434. // Smooth track
  2435. fCache->clear();
  2436. Refit(track, 0.13957, 1, kTRUE, 1, kFALSE, fCache); // refit toward beam line
  2437. map<Double_t,matrix4>::iterator mit = fCache->begin(), mit1 = mit;
  2438. for ( ; mit != fCache->end(); ++mit) {
  2439. if (mit == fCache->begin()) {
  2440. vecSmooth.push_back(pair<Double_t,TMatrixD>(mit->first,mit->second.xfilt));
  2441. mit1 = mit;
  2442. continue;
  2443. }
  2444. //*
  2445. TMatrixD wfilt(mit->second.wfilt); // filtered weight
  2446. wfilt.Invert(); // covar
  2447. TMatrixD ak(wfilt,TMatrixD::kMult,mit1->second.jacob); // smoother gain matrix
  2448. TMatrixD dx(vecSmooth.back().second); // smoothed previous
  2449. dx -= (mit1->second.xextr); // -extrap. previous
  2450. TMatrixD xsmooth(ak,TMatrixD::kMult,dx);
  2451. xsmooth += (mit->second.xfilt);
  2452. vecSmooth.push_back(pair<Double_t,TMatrixD>(mit->first,xsmooth));
  2453. //*/
  2454. mit1 = mit;
  2455. }
  2456. }
  2457. //__________________________________________________________________________
  2458. ClassImp(MpdTpcKalmanFilter)