MpdFillDstTask.cxx 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518
  1. // Author: Oleg Rogachevsky
  2. // Update: 2009-10-07 17:56:17+0400
  3. // Copyright: 2009 (C) MPD coll.
  4. //
  5. // fill dst task
  6. #include "MpdFillDstTask.h"
  7. #include "MpdKalmanTrack.h"
  8. #include "MpdTpcKalmanTrack.h"
  9. #include "MpdEctKalmanTrack.h"
  10. #include "MpdTofMatchingData.h"
  11. #include "MpdVertex.h"
  12. #include "MpdParticleIdentification.h"
  13. #include "MpdHelix.h"
  14. #include "MpdFieldCreator.h"
  15. #include "MpdMapPar.h"
  16. #include "MpdMCTrack.h"
  17. #include "FairRootManager.h"
  18. #include "FairRunAna.h"
  19. #include "FairRuntimeDb.h"
  20. #include <TMath.h>
  21. #include <TMatrixD.h>
  22. #include "TGeoManager.h"
  23. #include "TGeoBBox.h"
  24. #include "TGeoTube.h"
  25. #include <assert.h>
  26. #include <set>
  27. #include <iostream>
  28. #include <algorithm>
  29. #include <vector>
  30. #include <bitset>
  31. using namespace std;
  32. //using std::cout;
  33. //using std::flush;
  34. //using std::endl;
  35. // ----- constructor with names ------------------------------------
  36. MpdFillDstTask::MpdFillDstTask(const char *name, const char *title)
  37. : FairTask(name),fEvent(nullptr),
  38. fKFTracks(nullptr),fKFEctTracks(nullptr),fMCTracks(nullptr), fGenTracks(nullptr), fTpcHits(nullptr),
  39. fMCEventHeader(nullptr),fTofMatching(nullptr),fEtofMatching(nullptr),
  40. fVertex(nullptr),fZdcSkeletonesSaved(kFALSE),
  41. fHistZdc1En(nullptr),fHistZdc2En(nullptr),
  42. fELossZdc1Histo(nullptr),fELossZdc2Histo(nullptr),
  43. fELossZdc1Value(nullptr),fELossZdc2Value(nullptr),
  44. fhTrackMotherId(nullptr),fhTrackPrimaryPDG(nullptr),
  45. fhTrackVertex(nullptr),fhTruthVertex(nullptr),
  46. fPID(nullptr),fSharedHitArraySize(0),fSharedHitArray(nullptr)
  47. {
  48. }
  49. // ----- Destructor -----------------------------------------------
  50. MpdFillDstTask::~MpdFillDstTask() {
  51. if (fEvent) delete fEvent;
  52. if(fSharedHitArraySize>0) delete []fSharedHitArray;
  53. }
  54. // -------------------------------------------------------------------
  55. InitStatus MpdFillDstTask::Init() {
  56. fEvent = new MpdEvent();
  57. FairRootManager *manager = FairRootManager::Instance();
  58. fKFTracks = (TClonesArray *) manager->GetObject("TpcKalmanTrack");
  59. fKFEctTracks = (TClonesArray *) manager->GetObject("EctTrack");
  60. fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
  61. fGenTracks = (TClonesArray *) manager->GetObject("GenTracks");
  62. fTpcHits = (TClonesArray*) manager->GetObject("TpcRecPoint");
  63. fMCEventHeader = (FairMCEventHeader*) manager->GetObject("MCEventHeader.");
  64. fTofMatching = (TClonesArray *) manager->GetObject("TOFMatching");
  65. fEtofMatching = (TClonesArray *) manager->GetObject("ETOFMatching");
  66. fVertex = (TClonesArray *) manager->GetObject("Vertex"); //AZ
  67. fELossZdc1Value = (TClonesArray *) manager->GetObject("ELossZdc1Value"); //EL
  68. fELossZdc2Value = (TClonesArray *) manager->GetObject("ELossZdc2Value"); //EL
  69. fELossZdc1Histo = (TClonesArray *) manager->GetObject("ELossZdc1Histo"); //EL
  70. fELossZdc2Histo = (TClonesArray *) manager->GetObject("ELossZdc2Histo"); //EL
  71. fHistZdc1En = (TH2F *) manager->GetObject("HistZdc1En"); //EL
  72. fHistZdc2En = (TH2F *) manager->GetObject("HistZdc2En"); //EL
  73. fZdcSkeletonesSaved = 0;
  74. if(fPID==nullptr){
  75. fPID = new MpdPid(0, 0, 0, 0, "DEFAULT", "CF", "");
  76. //PID with cluster finder
  77. // fPID->Init("DEFAULT","CF","");
  78. }
  79. if(fTpcHits==nullptr){
  80. cout<<"WARNING: no TpcRec array, can't fill hit maps!"<<endl;
  81. }
  82. fSharedHitArraySize = 1000;
  83. fSharedHitArray = new Short_t[fSharedHitArraySize];
  84. FairRootManager::Instance()->Register("MCEventHeader.", "MC", fMCEventHeader, kTRUE);
  85. FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
  86. if (fGenTracks)
  87. FairRootManager::Instance()->Register("GenTracks", "MC", fGenTracks, kTRUE);
  88. //FairRootManager::Instance()->Register("MpdEvent","MpdEvents", fEvents, kTRUE);
  89. FairRootManager::Instance()->Register("MPDEvent.", "MpdEvent", fEvent, kTRUE);
  90. // FairRootManager::Instance()->Register("EZdc1","EZdc1", fELossZdc1Histo, kTRUE);
  91. // FairRootManager::Instance()->Register("EZdc2","EZdc2", fELossZdc2Histo, kTRUE);
  92. //fhTrackMotherId = new TH1F("TrackMotherId","",1000,-5,995);
  93. //fhTrackPrimaryPDG = new TH1F("TrackPrimaryPDG","",1000,-500,500);
  94. //fhTrackVertex = new TH2F("TrackVertex","",1000,-500,500,1000,-500,500);
  95. //fhTruthVertex = new TH2F("TruthVertex","",1000,-500,500,1000,-500,500);
  96. return kSUCCESS;
  97. }
  98. // -------------------------------------------------------------------
  99. void MpdFillDstTask::Exec(Option_t * option) {
  100. Reset(); // Clear previos event information
  101. if (!fZdcSkeletonesSaved) { // empty skeletones saved only once
  102. FairRootManager* ioman = FairRootManager::Instance();
  103. if (fHistZdc1En) {
  104. fHistZdc1En->SetDirectory((TFile*) ioman->GetOutFile());
  105. fHistZdc2En->SetDirectory((TFile*) ioman->GetOutFile());
  106. fHistZdc1En->Write();
  107. fHistZdc2En->Write();
  108. fZdcSkeletonesSaved = 1;
  109. }
  110. }
  111. Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
  112. Int_t nEctReco = fKFEctTracks ? fKFEctTracks->GetEntriesFast() : 0;
  113. cout << "\n-I- [MpdFillDstTask::Exec] " << nReco + nEctReco << " reconstruced tracks to write" << endl;
  114. if(fTpcHits)
  115. CalculateSharedArrayMap();
  116. FairRunAna *fRun = FairRunAna::Instance();
  117. fEvent->SetRunInfoRunId(fRun->GetRunId());
  118. FairField *fPar = fRun->GetField();
  119. fEvent->SetRunInfoMagneticFieldZ(fPar->GetType());
  120. if (fVertex) {
  121. MpdVertex *vertex = (MpdVertex*) fVertex->UncheckedAt(0);
  122. fEvent->SetPrimaryVerticesX(vertex->GetX());
  123. fEvent->SetPrimaryVerticesY(vertex->GetY());
  124. fEvent->SetPrimaryVerticesZ(vertex->GetZ());
  125. } else {
  126. // Vertex info is not available
  127. fEvent->SetPrimaryVerticesX(0.);
  128. fEvent->SetPrimaryVerticesY(0.);
  129. fEvent->SetPrimaryVerticesZ(0.);
  130. }
  131. // helix works with meters, primary vertices are in cm
  132. TVector3 recoVertex(fEvent->GetPrimaryVerticesX(),
  133. fEvent->GetPrimaryVerticesY(),
  134. fEvent->GetPrimaryVerticesZ());
  135. TVector3 mcVertex;
  136. fMCEventHeader->GetVertex(mcVertex);
  137. // check clone track into ECT
  138. typedef std::set<Int_t> trackSet;
  139. trackSet EctTrackSet;
  140. for (Int_t index = 0; index < nEctReco; index++) // cycle by ECT KF tracks
  141. {
  142. MpdEctKalmanTrack *track = (MpdEctKalmanTrack*) fKFEctTracks->UncheckedAt(index);
  143. if (track->IsFromTpc()) EctTrackSet.insert(trackSet::value_type(track->GetTpcIndex()));
  144. }
  145. Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
  146. MpdTofMatchingData *pMatchingData;
  147. bool matchingDataExist;
  148. //AZ MpdParticleIdentification *identificator = new MpdParticleIdentification();
  149. MpdParticleIdentification identificator;
  150. for (Int_t i = 0; i < nReco; i++) {
  151. // check clone track into ECT
  152. if (nEctReco > 0) {
  153. trackSet::iterator it = EctTrackSet.find(i);
  154. if (it != EctTrackSet.end()) continue;
  155. }
  156. MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
  157. MpdTpcKalmanTrack *kfTPCtrack = (MpdTpcKalmanTrack*) fKFTracks->UncheckedAt(i);
  158. MpdTrack *track = fEvent->AddGlobalTrack();
  159. if (kfTPCtrack->GetRecoQuality()) track->SetEdgeCut(kTRUE);
  160. track->SetID(kftrack->GetTrackID());
  161. track->SetNofHits(kftrack->GetNofHits());
  162. track->SetdEdXTPC(kftrack->GetPartID());
  163. Float_t Ppi, Pk, Pe, Pp;
  164. if (!identificator.GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 0)) { //0 - equal bayesian coefficients
  165. track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2));
  166. }
  167. matchingDataExist = false;
  168. for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++) {
  169. pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
  170. if (pMatchingData->GetKFTrackIndex() == i) {
  171. matchingDataExist = true;
  172. break;
  173. } // first matching
  174. }
  175. if (matchingDataExist) {
  176. track->SetTofBeta(pMatchingData->GetBeta());
  177. track->SetTofMass2(pMatchingData->GetMass2());
  178. track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
  179. if (!identificator.GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetBeta(), Ppi, Pk, Pp, Pe, 0)) {
  180. track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1));
  181. }
  182. }
  183. Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
  184. Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
  185. Float_t combProbs[4]; //probabilities combined from TOF & TPC
  186. identificator.GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
  187. Ppi = combProbs[0];
  188. Pk = combProbs[1];
  189. Pp = combProbs[2];
  190. Pe = combProbs[3];
  191. track->SetCombPidProb(Pe, Ppi, Pk, Pp);
  192. if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
  193. else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
  194. track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
  195. track->SetPhi(kftrack->GetParam(2)); // Phi
  196. TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
  197. track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
  198. track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
  199. track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
  200. track->SetChi2(kftrack->GetChi2());
  201. Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
  202. track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
  203. track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
  204. track->SetFirstPointZ(kftrack->GetParam(1));
  205. track->SetLastPointX(0.); // AZ - currently not available
  206. track->SetLastPointY(0.); // AZ - currently not available
  207. track->SetLastPointZ(0.); // AZ - currently not available
  208. FillTrackDCA(track, &recoVertex, &mcVertex);
  209. FillTrackPID(track);
  210. FillTrackTpcHits(i, track);
  211. }
  212. Int_t nEMatching = fEtofMatching ? fEtofMatching->GetEntries() : 0;
  213. MpdTofMatchingData *pEMatchingData;
  214. for (Int_t i = 0; i < nEctReco; i++) {
  215. MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFEctTracks->UncheckedAt(i);
  216. MpdTrack *track = fEvent->AddGlobalTrack();
  217. track->SetID(kftrack->GetTrackID());
  218. track->SetNofHits(kftrack->GetNofHits());
  219. matchingDataExist = false;
  220. for (Int_t etofIndex = 0; etofIndex < nEMatching; etofIndex++) {
  221. pEMatchingData = (MpdTofMatchingData*) fEtofMatching->UncheckedAt(etofIndex);
  222. if (pEMatchingData->GetKFTrackIndex() == i) {
  223. matchingDataExist = true;
  224. break;
  225. } // first matching
  226. }
  227. if (matchingDataExist) {
  228. track->SetTofMass2(pEMatchingData->GetMass2());
  229. track->SetTofHitIndex(pEMatchingData->GetTofHitIndex());
  230. }
  231. if (kftrack->GetParam(4) == 0.) track->SetPt(TMath::Sqrt(-1)); /*NaN*/
  232. else track->SetPt(1. / kftrack->GetParam(4)); /*signed Pt*/
  233. track->SetTheta(TMath::PiOver2() - kftrack->GetParam(3)); // Theta: angle from beam line
  234. track->SetPhi(kftrack->GetParam(2)); // Phi
  235. TMatrixD Cov = *kftrack->GetCovariance(); // Error matrix
  236. track->SetPtError(TMath::Sqrt(Cov(4, 4))); // Pt error
  237. track->SetThetaError(TMath::Sqrt(Cov(3, 3))); // Theta error
  238. track->SetPhiError(TMath::Sqrt(Cov(2, 2))); // Phi error
  239. track->SetChi2(kftrack->GetChi2());
  240. Double_t phi = kftrack->GetParam(0) / kftrack->GetPosNew();
  241. track->SetFirstPointX(kftrack->GetPosNew() * TMath::Cos(phi)); // closest to beam line
  242. track->SetFirstPointY(kftrack->GetPosNew() * TMath::Sin(phi));
  243. track->SetFirstPointZ(kftrack->GetParam(1));
  244. track->SetLastPointX(0.); // AZ - currently not available
  245. track->SetLastPointY(0.); // AZ - currently not available
  246. track->SetLastPointZ(0.); // AZ - currently not available
  247. FillTrackDCA(track, &recoVertex, &mcVertex);
  248. FillTrackPID(track);
  249. }
  250. }
  251. // -------------------------------------------------------------------
  252. //Delete MC tracks being outside the MPD
  253. void MpdFillDstTask::CleanMC() {
  254. cout << "-I- [MpdFillDstTask::Exec] Cleaning from outer decays..." << flush;
  255. Int_t nMC = fMCTracks->GetEntriesFast(), motherId;
  256. MpdMCTrack *track;
  257. //reading outer radius of TOF
  258. FairRunAna::Instance()->GetRuntimeDb()->getContainer("FairBaseParSet");
  259. assert(gGeoManager);
  260. TGeoTube* geoTube = (TGeoTube*) gGeoManager->GetVolume("tof1")->GetShape();
  261. Double_t rMax = geoTube->GetRmax();
  262. //cout <<"rMax="<<rMax<<endl;
  263. vector<int> removedMother;
  264. Int_t i;//, j;
  265. for (i = 0; i < nMC; i++) {
  266. track = (MpdMCTrack*) fMCTracks->UncheckedAt(i);
  267. //filling hostograms
  268. motherId = track->GetMotherId();
  269. //fhTrackMotherId->Fill(motherId);
  270. if (motherId == -1) {
  271. //fhTrackPrimaryPDG->Fill(track->GetPdgCode());
  272. } else {
  273. //fhTrackVertex->Fill(track->GetStartX(),track->GetStartY());
  274. //if decay was out of MPD
  275. if (sqrt(track->GetStartX() * track->GetStartX() + track->GetStartY() * track->GetStartY()) > rMax) {
  276. fMCTracks->Remove(track);
  277. removedMother.push_back(i);
  278. } else {
  279. //motherId < childId (i) REQUIRED
  280. /*for (j = 0; j < cntRemovedMother; j++){
  281. if (motherId == removedMother[j]){
  282. fMCTracks->Remove(track);
  283. removedMother[cntRemovedMother] = i;
  284. cntRemovedMother++;
  285. break;
  286. }
  287. }
  288. if (j == cntRemovedMother) fhTruthVertex->Fill(track->GetStartX(),track->GetStartY());
  289. else cnt++;*/
  290. /*Int_t left = 0, right = cntRemovedMother-1, middle;
  291. while (left <= right){
  292. middle = (left + right)/2;
  293. if (motherId > removedMother[middle])
  294. left = middle + 1;
  295. else{
  296. if (motherId < removedMother[middle])
  297. right = middle - 1;
  298. else
  299. break;
  300. }
  301. }
  302. if (left > right) fhTruthVertex->Fill(track->GetStartX(),track->GetStartY());
  303. else{
  304. fMCTracks->Remove(track);
  305. removedMother[cntRemovedMother] = i;
  306. cntRemovedMother++;
  307. }*/
  308. //binary search
  309. if (binary_search(removedMother.begin(), removedMother.end(), motherId)) {
  310. fMCTracks->Remove(track);
  311. removedMother.push_back(i);
  312. }
  313. //else
  314. // fhTruthVertex->Fill(track->GetStartX(),track->GetStartY());
  315. }
  316. }
  317. }
  318. fMCTracks->Compress();
  319. cout << endl;
  320. }
  321. // -------------------------------------------------------------------
  322. void MpdFillDstTask::Reset() {
  323. fEvent->Reset();
  324. }
  325. // -------------------------------------------------------------------
  326. void MpdFillDstTask::Finish() {
  327. //cout<<"\n-I- [MpdFillDstTask::Finish] "<< endl;
  328. //fEvents->Dump();
  329. //cout << "\n";
  330. /* !this code gives some problems! -- !add extra event and so segfault!
  331. FairRootManager *fManager = FairRootManager::Instance();
  332. fManager->Fill();
  333. */
  334. //fhTrackMotherId->Write();
  335. //fhTrackPrimaryPDG->Write();
  336. //fhTrackVertex->Write();
  337. //fhTruthVertex->Write();
  338. delete fEvent;
  339. fEvent = nullptr;
  340. }
  341. void MpdFillDstTask::CalculateSharedArrayMap() {
  342. if(fTpcHits->GetEntriesFast()>fSharedHitArraySize){
  343. delete []fSharedHitArray;
  344. fSharedHitArraySize = fTpcHits->GetEntriesFast();
  345. fSharedHitArray = new Short_t [fSharedHitArraySize];
  346. }
  347. for(int i=0;i<fTpcHits->GetEntriesFast();i++){
  348. fSharedHitArray[i]=0;
  349. }
  350. for(int i=0;i<fKFTracks->GetEntriesFast();i++){
  351. MpdTpcKalmanTrack *kalman = (MpdTpcKalmanTrack*) fKFTracks->UncheckedAt(i);
  352. TObjArray *khits = kalman->GetTrHits();
  353. if(khits)
  354. for(int j=0;j<khits->GetEntriesFast();j++){
  355. MpdKalmanHit *hit = (MpdKalmanHit*)khits->UncheckedAt(j);
  356. Int_t id = hit->GetIndex();// index of MpdTpcHit
  357. if(fSharedHitArray[id]<3)
  358. fSharedHitArray[id]++;
  359. }
  360. }
  361. }
  362. void MpdFillDstTask::FillTrackTpcHits(Int_t particle_index, MpdTrack *track) {
  363. MpdTpcKalmanTrack *kalman = (MpdTpcKalmanTrack*) fKFTracks->UncheckedAt(particle_index);
  364. ULong64_t layerHit = 0;
  365. ULong64_t sharedHit = 0;
  366. if(fTpcHits){// calculate both maps
  367. TObjArray *khits = kalman->GetTrHits();
  368. if(khits)
  369. for(int j=0;j<khits->GetEntriesFast();j++){
  370. MpdKalmanHit *hit = (MpdKalmanHit*)khits->UncheckedAt(j);
  371. Int_t id = hit->GetIndex();// index of MpdTpcHit
  372. Int_t layer = hit->GetLayer();
  373. if(layer>=0){
  374. if(fSharedHitArray[id]>1)
  375. SETBIT(sharedHit,layer);
  376. SETBIT(layerHit,layer);
  377. }
  378. }
  379. }else{//calculate only hit layer map
  380. TObjArray *khits = kalman->GetTrHits();
  381. if(khits)
  382. for(int j=0;j<khits->GetEntriesFast();j++){
  383. MpdKalmanHit *hit = (MpdKalmanHit*)khits->UncheckedAt(j);
  384. Int_t layer = hit->GetLayer();
  385. if(layer>=0){
  386. SETBIT(layerHit,layer);
  387. }
  388. }
  389. }
  390. track->SetLayerHitMap(layerHit);
  391. track->SetSharedHitMap(sharedHit);
  392. }
  393. // -------------------------------------------------------------------
  394. /*MpdEvent *MpdFillDstTask::AddEvent(Option_t * option)
  395. {
  396. TClonesArray &events = *fEvents;
  397. Int_t size = events.GetEntriesFast(); // It is really always only one record
  398. MpdEvent* event = new(events[size]) MpdEvent();
  399. return event;
  400. }*/
  401. MpdTrack *MpdFillDstTask::AddPrimaryTrack() {
  402. return nullptr;
  403. }
  404. void MpdFillDstTask::FillTrackDCA(MpdTrack* track, TVector3 *recoVertex, TVector3 *mcVertex) {
  405. MpdHelix helix = track->GetHelix();
  406. Double_t path_at_mcVertex;
  407. Double_t path_at_recoVertex;
  408. path_at_mcVertex = helix.pathLength(*mcVertex);
  409. path_at_recoVertex = helix.pathLength(*recoVertex);
  410. TVector3 DCA_MC = helix.at(path_at_mcVertex);
  411. TVector3 DCA_RECO = helix.at(path_at_recoVertex);
  412. // set dca global as dca to MC vertex DW
  413. track->SetDCAGlobalX(DCA_MC.X()-mcVertex->X());
  414. track->SetDCAGlobalY(DCA_MC.Y()-mcVertex->Y());
  415. track->SetDCAGlobalZ(DCA_MC.Z()-mcVertex->Z());
  416. // set dca as dca to reconstructed vertex DW
  417. track->SetDCAX(DCA_RECO.X()-recoVertex->X());
  418. track->SetDCAY(DCA_RECO.Y()-recoVertex->Y());
  419. track->SetDCAZ(DCA_RECO.Z()-recoVertex->Z());
  420. }
  421. void MpdFillDstTask::FillTrackPID(MpdTrack* track) {
  422. TVector3 mom(track->GetPx(), track->GetPy(), track->GetPz());
  423. Double_t p = mom.Mag();
  424. Double_t dedx = track->GetdEdXTPC();
  425. Double_t dedx_el = fPID->GetDedxElParam(p);
  426. Double_t dedx_pi = fPID->GetDedxPiParam(p);
  427. Double_t dedx_ka = fPID->GetDedxKaParam(p);
  428. Double_t dedx_pr = fPID->GetDedxPrParam(p);
  429. Double_t sigma_el = fPID->GetDedxWidthValue(p, MpdPidUtils::kElectron) * dedx_el;
  430. Double_t sigma_pi = fPID->GetDedxWidthValue(p, MpdPidUtils::kPion) * dedx_pi;
  431. Double_t sigma_ka = fPID->GetDedxWidthValue(p, MpdPidUtils::kKaon) * dedx_ka;
  432. Double_t sigma_pr = fPID->GetDedxWidthValue(p, MpdPidUtils::kProton) * dedx_pr;
  433. if (sigma_el > FLT_EPSILON)
  434. track->SetNSigmaElectron((dedx - dedx_el) / sigma_el);
  435. if (sigma_pi > FLT_EPSILON)
  436. track->SetNSigmaPion((dedx - dedx_pi) / sigma_pi);
  437. if (sigma_ka > FLT_EPSILON)
  438. track->SetNSigmaKaon((dedx - dedx_ka) / sigma_ka);
  439. if (sigma_pr > FLT_EPSILON)
  440. track->SetNSigmaProton((dedx - dedx_pr) / sigma_pr);
  441. }
  442. // -------------------------------------------------------------------
  443. ClassImp(MpdFillDstTask);