MpdKfParticleFinder.cxx 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649
  1. #include "MpdKfParticleFinder.h"
  2. MpdKfParticleFinder::MpdKfParticleFinder(TString dstName) :
  3. fTpcTracks(nullptr),
  4. fMcTracks(nullptr),
  5. fMiniTracks(nullptr),
  6. fMiniCovMatrices(nullptr),
  7. fPrimaryVertices(nullptr),
  8. fTopoReconstructor(nullptr),
  9. fPerformance(nullptr),
  10. fSecondaries(nullptr),
  11. fKFVertices(nullptr),
  12. fMCHeader(nullptr),
  13. fMCEvent(nullptr),
  14. fHeader(nullptr),
  15. isUseCuts(kFALSE),
  16. isMini(kFALSE),
  17. isPerformance(kFALSE),
  18. fCuts(nullptr), fPdgDB(nullptr), fInput("") {
  19. fVerbose = 0;
  20. // Setting default PID hypo to be used ...
  21. fPidHypo.first = 2212;
  22. fPidHypo.second = -211;
  23. // Initialize FairRoot manager
  24. ioman = FairRootManager::Instance();
  25. FairFileSource* fFileSource = nullptr;
  26. FairSource* fSource = nullptr;
  27. MpdMiniDstFileSource* fMiniSource = nullptr;
  28. // It can be a single file or a list of files (*.lis or *.list) to process ...
  29. if (dstName.Contains(".lis") || dstName.Contains(".list")) {
  30. // Reading passed list with inputs ...
  31. ReadList(dstName);
  32. if (isMini) {
  33. fMiniSource = new MpdMiniDstFileSource(*fInputs.begin());
  34. for (auto it = next(fInputs.begin(), 1); it != fInputs.end(); it++)
  35. fMiniSource->AddFile(*it);
  36. fSource = (FairSource*) fMiniSource;
  37. } else {
  38. fFileSource = new FairFileSource(*fInputs.begin());
  39. for (auto it = next(fInputs.begin(), 1); it != fInputs.end(); it++)
  40. fFileSource->AddFile(*it);
  41. }
  42. } else if (dstName.Contains(".root")) {
  43. fInput = dstName;
  44. if (dstName.Contains("MiniDst")) {
  45. isMini = kTRUE;
  46. if (checkBranchStatus(fInput))
  47. fSource = new MpdMiniDstFileSource(fInput);
  48. }
  49. else
  50. if (checkBranchStatus(fInput))
  51. fFileSource = new FairFileSource(fInput);
  52. } else
  53. Fatal("MpdKfParticleFinder::MpdKfParticleFinder(TString dstName)", "Provided input seems to be incorrect!");
  54. // Setting already specified file source ...
  55. if (fFileSource)
  56. FairRunAna::Instance()->SetSource(fFileSource);
  57. else if (fSource)
  58. FairRunAna::Instance()->SetSource(fSource);
  59. }
  60. void MpdKfParticleFinder::ReadList(TString dstName) {
  61. string line;
  62. ifstream f(dstName.Data(), ios::in);
  63. TString dst = "";
  64. Int_t id = 0;
  65. const Int_t nMaxFiles = 5000;
  66. Int_t iCurrentFile = 0;
  67. while (!f.eof()) {
  68. f >> id >> dst;
  69. fInputs.push_back(dst);
  70. fIds[dst] = id;
  71. getline(f, line);
  72. iCurrentFile++;
  73. if (iCurrentFile > nMaxFiles)
  74. Fatal("MpdKfParticleFinder::ReadList()", "Exceeded number of files or something wrong with passed list!");
  75. }
  76. for (auto it = fInputs.begin(); it != fInputs.end(); it++) {
  77. dst = *it;
  78. if (!dst.Contains(".root"))
  79. fInputs.erase(it--);
  80. }
  81. auto it = fInputs.begin();
  82. isMini = ((*it).Contains("MiniDst") ? kTRUE : kFALSE);
  83. // Checking branch status ...
  84. for (auto it0 = fInputs.begin(); it0 != fInputs.end(); it0++)
  85. if (!checkBranchStatus(*it0))
  86. fInputs.erase(it0--);
  87. }
  88. InitStatus MpdKfParticleFinder::Init() {
  89. if (isUseCuts)
  90. fCuts = new TrackCuts();
  91. // Getting necessary branches to be used ...
  92. // It depends on a dst filename passed as argument
  93. if (isMini) {
  94. fMiniTracks = (TClonesArray*) ioman->GetObject("Track");
  95. fMiniCovMatrices = (TClonesArray*) ioman->GetObject("TrackCovMatrix");
  96. fMCEvent = (TClonesArray*) ioman->GetObject("McEvent");
  97. } else {
  98. // Tpc tracks ...
  99. fTpcTracks = (TClonesArray*) ioman->GetObject("TpcKalmanTrack");
  100. // Monte Carlo tracks ...
  101. fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  102. // Reconstructed PV's ...
  103. fPrimaryVertices = (TClonesArray*) ioman->GetObject("Vertex");
  104. // MC event header ...
  105. fMCHeader = (FairMCEventHeader*) ioman->GetObject("MCEventHeader.");
  106. }
  107. fHeader = (FairEventHeader*) ioman->GetObject("EventHeader.");
  108. // Register output branch ...
  109. fSecondaries = new TClonesArray("KFParticle");
  110. ioman->Register("Secondaries", "Secondaries_", fSecondaries, kTRUE);
  111. fKFVertices = new TClonesArray("KFVertex");
  112. ioman->Register("Vertex", "Vertex_", fKFVertices, kTRUE);
  113. // Calling topoReconstructor ...
  114. fTopoReconstructor = new KFParticleTopoReconstructor();
  115. fTopoReconstructor->SetField(5.); // FIXME !!! (to be read automatically)
  116. // Setting a decay to be reconstructed ...
  117. map <Int_t, Bool_t> decayList;
  118. Int_t decayPdg = -1;
  119. if (fPidHypo.first == 2212 && fPidHypo.second == -211)
  120. decayPdg = 3122; // \Lambda^{0} -> p + \pion^{-}
  121. else if (fPidHypo.first == 211 && fPidHypo.second == -211)
  122. decayPdg = 310; // \Kaon^{0}_{s} -> \pion^{+} + \pion^{-}
  123. if (decayPdg != -1)
  124. decayList[decayPdg] = kTRUE;
  125. fTopoReconstructor->GetKFParticleFinder()->SetReconstructionList(decayList);
  126. map <Int_t, Bool_t> decayListToCheck = fTopoReconstructor->GetKFParticleFinder()->GetReconstructionList();
  127. if (decayListToCheck.size() != 1)
  128. Fatal("MpdKfParticleFinder::Init()", "Only one decay in the list is supported !!!"); // FIXME !!!
  129. for (auto pdg : decayListToCheck)
  130. fDecPdgs.push_back(pdg.first);
  131. // Let's see a list of decays to be reconstructed (if necessary) ...
  132. if (fVerbose)
  133. for (auto pdgDecay : decayListToCheck)
  134. cout << "PDG# " << pdgDecay.first << " isActive# " << pdgDecay.second << endl;
  135. // Right now KFPerformance does not support miniDst!
  136. if (isMini)
  137. isPerformance = kFALSE;
  138. if (isPerformance) {
  139. // Involving performance tools ...
  140. fPerformance = new KFTopoPerformance();
  141. // fPerformance->DoNotStoreMCHistograms();
  142. TString outPerf = (TString) FairRunAna::Instance()->GetOutputFile()->GetName();
  143. outPerf = "KFPerformance_" + outPerf;
  144. fOutFile = new TFile(outPerf.Data(), "recreate");
  145. fPerformance->CreateHistos("KFHistos", (TDirectory*) fOutFile, decayListToCheck);
  146. fPdgDB = new TDatabasePDG();
  147. }
  148. return kSUCCESS;
  149. }
  150. void MpdKfParticleFinder::Exec(Option_t* option) {
  151. if (evCounter % 1000 == 0)
  152. cout << "Processing Event# " << evCounter << endl;
  153. if (fInputs.size()) {
  154. auto it = fIds.find((TString) ioman->GetRootFile()->GetName());
  155. MpdMiniMcEvent* mcEv = (MpdMiniMcEvent*) fMCEvent->UncheckedAt(0);
  156. if (isMini)
  157. fHeader->SetMCEntryNumber(mcEv->eventId());
  158. if (it != fIds.end())
  159. fHeader->SetInputFileId(it->second);
  160. else
  161. Fatal("MpdKfParticleFinder::Exec", "Input file is not associated with its id!");
  162. }
  163. fTopoReconstructor->Clear();
  164. fSecondaries->Delete();
  165. fKFVertices->Delete();
  166. ProcessDst();
  167. evCounter++;
  168. }
  169. void MpdKfParticleFinder::Finish() {
  170. if (isPerformance) {
  171. // Performance tools ...
  172. fPerformance->GetHistosDirectory()->Write();
  173. fOutFile->Close();
  174. }
  175. // PDG DB ...
  176. if (fPdgDB)
  177. delete fPdgDB;
  178. }
  179. void MpdKfParticleFinder::Mini2Kalman(TClonesArray* arr) {
  180. for (Int_t iTrack = 0; iTrack < fMiniTracks->GetEntriesFast(); iTrack++) {
  181. MpdMiniTrack* mini = (MpdMiniTrack*) fMiniTracks->UncheckedAt(iTrack);
  182. MpdMiniTrackCovMatrix* cMatrix = (MpdMiniTrackCovMatrix*) fMiniCovMatrices->UncheckedAt(iTrack);
  183. MpdTpcKalmanTrack* kalTrack = new ((*arr)[arr->GetEntriesFast()]) MpdTpcKalmanTrack();
  184. kalTrack->SetPos(cMatrix->R0());
  185. kalTrack->SetPosNew(cMatrix->R());
  186. kalTrack->SetNofHits(mini->nHits());
  187. TMatrixD stateVector = cMatrix->stateVector();
  188. kalTrack->SetParam(stateVector);
  189. TMatrixDSym covMatrix = cMatrix->covarianceMatrix();
  190. kalTrack->SetCovariance(covMatrix);
  191. }
  192. }
  193. void MpdKfParticleFinder::ProcessDst() {
  194. vector <ExtendedMpdTpcKalmanTrack> tpcTracks;
  195. TClonesArray* tpcKalmanTracks = nullptr;
  196. if (isMini) {
  197. tpcKalmanTracks = new TClonesArray("MpdTpcKalmanTrack");
  198. Mini2Kalman(tpcKalmanTracks);
  199. }
  200. // 1. Loop over found TPC tracks ...
  201. TClonesArray* inTrackArray = isMini ? tpcKalmanTracks : fTpcTracks;
  202. for (Int_t iTpcTrack = 0; iTpcTrack < inTrackArray->GetEntriesFast(); iTpcTrack++) {
  203. MpdTpcKalmanTrack* tr = (MpdTpcKalmanTrack*) inTrackArray->UncheckedAt(iTpcTrack);
  204. ExtendedMpdTpcKalmanTrack track = *(ExtendedMpdTpcKalmanTrack*) tr;
  205. if (isUseCuts) {
  206. // Skipping tracks not satisfying to the established cuts ...
  207. Double_t eta = Log(Tan(track.Theta() / 2.));
  208. if (track.Pt() < fCuts->GetPtMin() || track.Pt() > fCuts->GetPtMax() || Abs(eta) > fCuts->GetAbsEta())
  209. continue;
  210. }
  211. track.fIdx = iTpcTrack;
  212. tpcTracks.push_back(track);
  213. }
  214. if (fVerbose)
  215. cout << "Number of found TPC tracks# " << tpcTracks.size() << endl;
  216. if (tpcTracks.size() == 0)
  217. return;
  218. // 2. Getting and adding information on primary vertex ...
  219. vector <Int_t> pvTrackIds;
  220. if (isMini)
  221. for (Int_t iMiniTrack = 0; iMiniTrack < fMiniTracks->GetEntriesFast(); iMiniTrack++) {
  222. MpdMiniTrack* mini = (MpdMiniTrack*) fMiniTracks->UncheckedAt(iMiniTrack);
  223. if (mini->isPrimary())
  224. pvTrackIds.push_back(iMiniTrack);
  225. } else {
  226. MpdVertex* recoVp = (MpdVertex*) fPrimaryVertices->UncheckedAt(0);
  227. Int_t* indcs = recoVp->GetIndices()->GetArray();
  228. Int_t size = recoVp->GetIndices()->GetSize();
  229. pvTrackIds.assign(indcs, indcs + size);
  230. }
  231. // 3. Filling KFPTrackVector (xyz, PxPyPz, Q, PDG) with data ...
  232. KFPTrackVector trVector;
  233. FillKFPTrackVector(tpcTracks, &trVector, pvTrackIds);
  234. // 4. Getting KFPTracks from KFPTrackVector ...
  235. vector <KFParticle> kfParticles;
  236. vector <Int_t>* pdgHypos = new vector <Int_t>;
  237. for (Int_t iTrack = 0; iTrack < trVector.Size(); iTrack++) {
  238. KFPTrack track;
  239. trVector.GetTrack(track, iTrack);
  240. Int_t q = track.Charge();
  241. Int_t hypo = (q > 0) ? fPidHypo.first : fPidHypo.second;
  242. KFParticle particle(track, hypo);
  243. particle.SetId(track.Id());
  244. kfParticles.push_back(particle);
  245. pdgHypos->push_back(hypo);
  246. }
  247. // 5. Initializing topoReconstructor calling one of possible user constructors ...
  248. fTopoReconstructor->Init(kfParticles, pdgHypos);
  249. // 6. Reconstruct primary vertex ...
  250. fTopoReconstructor->ReconstructPrimVertex();
  251. // 7. Reconstruct secondary particles and select best candidates ...
  252. fTopoReconstructor->SortTracks();
  253. fTopoReconstructor->ReconstructParticles();
  254. fTopoReconstructor->SelectParticleCandidates();
  255. // 8. Writing found candidates to output file ...
  256. for (KFParticle recoSecondary : fTopoReconstructor->GetParticles())
  257. new ((*fSecondaries) [fSecondaries->GetEntriesFast()]) KFParticle(recoSecondary);
  258. delete pdgHypos;
  259. // 9. Getting reconstructed primary vertex to write it to output ...
  260. KFVertex vertex = fTopoReconstructor->GetPrimKFVertex();
  261. new ((*fKFVertices) [fKFVertices->GetEntriesFast()]) KFVertex(vertex);
  262. if (isPerformance && !isMini) {
  263. // Involving KFParticlePerformance ...
  264. fPerformance->SetTopoReconstructor(fTopoReconstructor);
  265. // Filling KFMCTrack for current event ....
  266. vector <KFMCTrack> mcTracksInEvent = FillKFMCTrack();
  267. // Setting MC tracks ...
  268. fPerformance->SetMCTracks(mcTracksInEvent);
  269. // Checking MC tracks ...
  270. fPerformance->CheckMCTracks();
  271. // Doing KFParticle --> KFMCParticle (nDaughters = 1)
  272. vector <Int_t> kf2kfmc;
  273. kf2kfmc.resize(fTpcTracks->GetEntriesFast());
  274. for (Int_t iRP = 0; iRP < fTopoReconstructor->GetParticles().size(); iRP++) {
  275. const KFParticle& rPart = fTopoReconstructor->GetParticles()[iRP];
  276. if (rPart.NDaughters() != 1)
  277. continue;
  278. MpdKalmanTrack* tr = (MpdKalmanTrack*) fTpcTracks->UncheckedAt(rPart.DaughterIds().at(0));
  279. kf2kfmc.at(rPart.DaughterIds().at(0)) = tr->GetTrackID();
  280. }
  281. fPerformance->SetTrackMatch(kf2kfmc);
  282. // Matching mc and reco tracks ...
  283. fPerformance->MatchTracks();
  284. // Filling a set of QA histos ...
  285. fPerformance->FillHistos();
  286. }
  287. if (tpcKalmanTracks)
  288. delete tpcKalmanTracks;
  289. }
  290. vector <KFMCTrack> MpdKfParticleFinder::FillKFMCTrack() {
  291. vector <KFMCTrack> tracks;
  292. tracks.clear();
  293. // Getting Monte Carlo ids of reconstructed tracks ...
  294. vector <Int_t> recoMcIdx;
  295. recoMcIdx.clear();
  296. for (Int_t iMc = 0; iMc < fMcTracks->GetEntriesFast(); iMc++) {
  297. MpdMCTrack* mcTrack = (MpdMCTrack*) fMcTracks->UncheckedAt(iMc);
  298. KFMCTrack track;
  299. track.SetX(mcTrack->GetStartX());
  300. track.SetY(mcTrack->GetStartY());
  301. track.SetZ(mcTrack->GetStartZ());
  302. track.SetPx(mcTrack->GetPx());
  303. track.SetPy(mcTrack->GetPy());
  304. track.SetPz(mcTrack->GetPz());
  305. track.SetPDG(mcTrack->GetPdgCode());
  306. track.SetMotherId(mcTrack->GetMotherId());
  307. TParticlePDG* pdgParticle = fPdgDB->GetParticle(mcTrack->GetPdgCode());
  308. Int_t q = (pdgParticle) ? pdgParticle->Charge() / 3 : 0;
  309. track.SetQP(q / mcTrack->GetP());
  310. track.SetNMCPoints(mcTrack->GetNPoints(kTPC));
  311. // Trying to know whether it was reconstructed ...
  312. auto it = find(recoMcIdx.begin(), recoMcIdx.end(), iMc);
  313. if (it != recoMcIdx.end())
  314. track.SetReconstructed();
  315. else
  316. track.SetNotReconstructed();
  317. // Is it in the TPC acceptance? ...
  318. if (!isInAcceptance(mcTrack))
  319. track.SetOutOfDetector();
  320. tracks.push_back(track);
  321. }
  322. return tracks;
  323. }
  324. Bool_t MpdKfParticleFinder::isInAcceptance(MpdMCTrack* track) {
  325. Double_t eta = TMath::ATanH(track->GetPz() / track->GetP());
  326. Double_t pt = track->GetPt();
  327. if (TMath::Abs(eta) < 1.2 && pt > .15 && pt < 2.) // FIXME !!!
  328. return kTRUE;
  329. else
  330. return kFALSE;
  331. }
  332. void MpdKfParticleFinder::tpcInnerShell(MpdTpcKalmanTrack in, MpdTpcKalmanTrack& out) {
  333. // Let's get covariance matrix and track params. approximately corresponding to the inner TPC shell ...
  334. // We need to have a copy of considering track to be able to fill fParamNew
  335. // to be instantiated when doing the copy ...
  336. MpdTpcKalmanTrack tmp(in);
  337. tmp.SetDirection(MpdKalmanTrack::kOutward);
  338. /// Propagate TPC track to 27 cm inward (fPos == 27)
  339. MpdKalmanHit hitTmp;
  340. hitTmp.SetType(MpdKalmanHit::kFixedR);
  341. hitTmp.SetPos(tmp.GetPos()); // fPos ~= 27 cm
  342. tmp.SetParamNew(*tmp.GetParam());
  343. tmp.SetPos(tmp.GetPosNew());
  344. tmp.ReSetWeight();
  345. TMatrixDSym w = *tmp.GetWeight(); // save current weight matrix
  346. Bool_t ok = MpdKalmanFilter::Instance()->PropagateToHit(&tmp, &hitTmp, kFALSE, kFALSE);
  347. if (!ok) {
  348. out = tmp;
  349. return;
  350. }
  351. // tmp.SetWeight(w); // restore original weight matrix (near TPC inner shell)
  352. // После них GetCovariance() и GetParamNew() дадут ков. матр. и параметры в точке GetPosNew(), которая около R=27 см (на внутренней оболочке TPC).
  353. TMatrixDSym cov(*tmp.Weight2Cov());
  354. TMatrixD param(*tmp.GetParamNew());
  355. // Resetting params and cov. matrix by the propagated ones ...
  356. tmp.SetParam(param);
  357. tmp.SetCovariance(cov);
  358. out = tmp;
  359. }
  360. void MpdKfParticleFinder::lastHit(MpdTpcKalmanTrack in, MpdTpcKalmanTrack& out) {
  361. MpdTpcKalmanTrack tmp(in);
  362. // Resetting params and cov. matrix to those ones at the last hit ...
  363. tmp.SetParam(*in.GetParamAtHit());
  364. TMatrixDSym* weightAtLastHit = in.GetWeightAtHit();
  365. TMatrixDSym& covarAtLastHit = weightAtLastHit->InvertFast();
  366. tmp.SetCovariance(covarAtLastHit);
  367. out = tmp;
  368. }
  369. void MpdKfParticleFinder::dcaToBeamline(MpdTpcKalmanTrack in, MpdTpcKalmanTrack& out) {
  370. MpdTpcKalmanTrack tmp(in);
  371. // Resetting params and cov. matrix to those ones at dca2beamline ...
  372. tmp.SetParam(*in.GetParam());
  373. tmp.SetCovariance(*in.GetCovariance());
  374. out = tmp;
  375. }
  376. void MpdKfParticleFinder::FillKFPTrackVector(vector <ExtendedMpdTpcKalmanTrack> tracksFromTpc, KFPTrackVector* tracks, vector <Int_t> trIdxPv) {
  377. tracks->Resize(tracksFromTpc.size());
  378. for (Int_t iTrack = 0; iTrack < tracksFromTpc.size(); iTrack++) {
  379. ExtendedMpdTpcKalmanTrack in = tracksFromTpc[iTrack];
  380. MpdTpcKalmanTrack out;
  381. if (fVerbose) {
  382. cout << "Track params. at 2D-DCA to beamline: " << endl;
  383. in.GetParam()->Print();
  384. in.GetCovariance()->Print();
  385. }
  386. // Parameters are propagated to the inner TPC shell at R of order of 27 cm
  387. tpcInnerShell(in, out);
  388. // lastHit(in, out);
  389. // dcaToBeamline(in, out);
  390. if (fVerbose) {
  391. cout << "Propagated parameters: " << endl;
  392. out.GetParam()->Print();
  393. out.GetCovariance()->Print();
  394. }
  395. // Getting propagated track coordinates / angles / momenta ...
  396. Double_t phi = out.GetParam(0) / out.GetPosNew();
  397. Double_t X = out.GetPosNew() * Cos(phi);
  398. Double_t Y = out.GetPosNew() * Sin(phi);
  399. Double_t Z = out.GetZ();
  400. Double_t q = out.Charge();
  401. Double_t Px = out.Momentum3().X();
  402. Double_t Py = out.Momentum3().Y();
  403. Double_t Pz = out.Momentum3().Z();
  404. // Setting state vector ....
  405. vector <Double_t> parVectorToSet{X, Y, Z, Px, Py, Pz};
  406. for (Int_t iParam = 0; iParam < parVectorToSet.size(); iParam++)
  407. tracks->SetParameter(parVectorToSet[iParam], iParam, iTrack);
  408. tracks->SetQ(q, iTrack);
  409. tracks->SetId(in.fIdx, iTrack);
  410. // Setting PID hypo ...
  411. if (q > 0)
  412. tracks->SetPDG(fPidHypo.first, iTrack);
  413. else
  414. tracks->SetPDG(fPidHypo.second, iTrack);
  415. // Setting PVIndex = 0 to tracks considered as primaries ...
  416. auto it = find(trIdxPv.begin(), trIdxPv.end(), iTrack);
  417. if (it != trIdxPv.end())
  418. tracks->SetPVIndex(0, iTrack);
  419. // Probably, looks as a secondary ...
  420. else
  421. tracks->SetPVIndex(-1, iTrack);
  422. vector <Double_t> cov21 = DoErrorPropagationToXYZPxPyPz(out.GetCovariance(), out.GetParam(), &out);
  423. for (Int_t iCov = 0; iCov < cov21.size(); iCov++)
  424. tracks->SetCovariance(cov21.at(iCov), iCov, iTrack);
  425. }
  426. }
  427. vector <Double_t> MpdKfParticleFinder::DoErrorPropagationToXYZPxPyPz(TMatrixDSym* cov, TMatrixD * params, MpdTpcKalmanTrack* track) {
  428. // State track vector in MPD:
  429. // 0) r\Phi_{T}
  430. // 1) z
  431. // 2) \Phi
  432. // 3) \Lambda (dip angle)
  433. // 4) -q / Pt
  434. // J : A --> B
  435. // A : x, y, z, px, py, pz
  436. // B : x, y, z, \Phi, \Lambda, -q / Pt
  437. // Getting Jacobian ...
  438. Float_t J[6][6];
  439. for (Int_t i = 0; i < 6; i++)
  440. for (Int_t j = 0; j < 6; j++)
  441. J[i][j] = 0;
  442. // Getting corresponding track parameters ...
  443. Int_t q = track->Charge();
  444. Double_t pt = track->Pt();
  445. Double_t Phi = (*params)(2, 0);
  446. Double_t Lambda = (*params)(3, 0);
  447. // Coordinate transformations ...
  448. J[0][0] = 1.; // dx / dx
  449. J[1][1] = 1.; // dy / dy
  450. J[2][2] = 1.; // dz / dz
  451. // Momentum transformations ...
  452. J[3][3] = -pt * Sin(Phi); // dPx / d\Phi
  453. J[3][5] = pt * pt * Cos(Phi) / q; // dPx / d(-q / Pt)
  454. J[4][3] = pt * Cos(Phi); // dPy / d\Phi
  455. J[4][5] = pt * pt * Sin(Phi) / q; // dPy / d(-q / Pt)
  456. J[5][4] = pt / (Cos(Lambda) * Cos(Lambda)); // dPz / dLambda
  457. J[5][5] = Tan(Lambda) * pt * pt / q; // dPz / d(-q / pt)
  458. // Extending track covariance matrix by one row and one column ...
  459. Float_t CovIn[6][6]; // triangular -> symmetric matrix
  460. CovIn[0][0] = 1e-4; // dx. start from nowhere
  461. for (Int_t i = 1; i < 6; i++) {
  462. CovIn[i][0] = 0;
  463. CovIn[0][i] = 0;
  464. }
  465. for (Int_t i = 1; i < 6; i++) {
  466. for (Int_t j = 1; j <= i; j++) {
  467. CovIn[i][j] = (*cov)(i - 1, j - 1);
  468. CovIn[j][i] = (*cov)(i - 1, j - 1);
  469. }
  470. }
  471. Float_t CovInJt[6][6]; // CovInJt = CovIn * J^t
  472. for (Int_t i = 0; i < 6; i++)
  473. for (Int_t j = 0; j < 6; j++) {
  474. CovInJt[i][j] = 0;
  475. for (Int_t k = 0; k < 6; k++)
  476. CovInJt[i][j] += CovIn[i][k] * J[j][k];
  477. }
  478. Float_t CovOut[6][6]; // CovOut = J * CovInJt
  479. for (Int_t i = 0; i < 6; i++)
  480. for (Int_t j = 0; j < 6; j++) {
  481. CovOut[i][j] = 0;
  482. for (Int_t k = 0; k < 6; k++)
  483. CovOut[i][j] += J[i][k] * CovInJt[k][j];
  484. }
  485. vector <Double_t> KFPCov; // symmetric matrix -> triangular
  486. for (Int_t i = 0; i < 6; i++)
  487. for (Int_t j = 0; j <= i; j++)
  488. KFPCov.push_back(CovOut[i][j]);
  489. return KFPCov;
  490. }