MpdV0Finder.cxx 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441
  1. // Task for finding secondary vertices
  2. #include "MpdV0Finder.h"
  3. //__________________________________________________________________________
  4. MpdV0Finder::MpdV0Finder(const char *name, Int_t iVerbose)
  5. : FairTask(name, iVerbose),
  6. V0(0),
  7. fEvent(0),
  8. fMpdTracks(0),
  9. fKFTracks(0),
  10. fTofMatching(0),
  11. fMCTracks(0),
  12. pi(0),
  13. pr(0)
  14. {}
  15. //__________________________________________________________________________
  16. MpdV0Finder::~MpdV0Finder() {
  17. // delete V0; V0 = 0;
  18. fV0Cont->Clear("C");
  19. //fMpdTracks->Clear("C");
  20. pi->Clear("C");
  21. pr->Clear("C");
  22. }
  23. //__________________________________________________________________________
  24. InitStatus MpdV0Finder::Init() {
  25. cout << endl;
  26. cout << " MpdV0Finder::Init " << endl;
  27. npr = 0;
  28. npi = 0;
  29. PVert = 0.;
  30. vtx = 0.;
  31. fNMPDtracks = 0;
  32. fNtracks = 0;
  33. fV0Cont = 0x0;
  34. fMpdTracks = 0x0;
  35. fKFTracks = 0x0;
  36. fTofMatching = 0x0;
  37. fMCTracks = 0x0;
  38. fNMPDtracks = 0;
  39. nvert = 0;
  40. Nreco = 0;
  41. Nlam = 0;
  42. pi = new TClonesArray("MpdParticle"); //put in exec to check for mem leak
  43. pr = new TClonesArray("MpdParticle");
  44. FairRootManager *manager = FairRootManager::Instance();
  45. fV0Cont = new TClonesArray("MpdV0", nvert); //container for v0s
  46. manager->Register("V02", "MpdV0", fV0Cont, kTRUE); //Register v0s
  47. //fMpdTracks = new TClonesArray("MpdTrack"); //container for MpdTracks if used in reco chain this is for probability
  48. fPrimVertex = (TClonesArray*) manager->GetObject("Vertex");
  49. fKFTracks = (TClonesArray*) manager->GetObject("TpcKalmanTrack");
  50. fTofMatching = (TClonesArray*) manager->GetObject("TOFMatching");
  51. // fEvent = (MpdEvent*) manager->GetObject("MPDEvent.");
  52. // fEvent = (MpdEvent*) manager->ActivateBranch("MPDEvent.");
  53. fEvent = (MpdEvent*) manager->GetObjectFromInTree("MPDEvent.");
  54. //fMpdTracks = fEvent->GetGlobalTracks();
  55. //cout<<" GlobalTracks = "<<fMpdTracks->GetEntriesFast()<<endl;
  56. fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
  57. //fMpdTracks = (TClonesArray*) manager->GetObject("MpdTracks");
  58. return kSUCCESS;
  59. }
  60. //__________________________________________________________________________
  61. InitStatus MpdV0Finder::ReInit() { return kSUCCESS; }
  62. //__________________________________________________________________________
  63. void MpdV0Finder::Reset() {
  64. cout << " MpdV0Finder::Reset " << endl;
  65. //fV0Cont->Delete();
  66. //fMpdTracks->Delete();
  67. fV0Cont->Clear("C");
  68. //fMpdTracks->Clear("C");
  69. pi->Clear("C");
  70. pr->Clear("C");
  71. //fMCTracks->Clear("C");
  72. //fKalmanTrs->Clear();
  73. //fKalmanTrs->Expand(0);
  74. //fMpdTracks->Clear();
  75. Nok1 = 0;
  76. Nmc1 = 0;
  77. Nlam1 = 0;
  78. nvert = 0;
  79. npr = 0;
  80. npi = 0;
  81. // pr.resize(0);
  82. // pi.resize(0);
  83. }
  84. //__________________________________________________________________________
  85. void MpdV0Finder::Finish() {
  86. cout<<"***For all events*** "<<endl;
  87. Double_t Nbad = Nreco - Nok; //accumulates
  88. cout<<" Nreco = "<<Nreco<<" Nmc = "<<Nmc<<" Nok = "<<Nok<<" Nbad = "<<Nbad<<endl;
  89. Double_t Cont = (Double_t)Nbad/(Double_t)Nreco; //accumulates
  90. Double_t Eff = (Double_t)Nok/(Double_t)Nmc; //accumulates
  91. cout<<"Efficiency = "<<Eff*100<<"% Contamination = "<<Cont*100<<"%"<<endl;
  92. cout<<"MClambdas = "<<Nlam<<endl;
  93. }
  94. //__________________________________________________________________________
  95. void MpdV0Finder::Exec(Option_t * option) {
  96. static int eventCounter = 0;
  97. ++eventCounter;
  98. cout << endl;
  99. if (eventCounter != 796) return;
  100. time_t now = time(0); // current date/time based on current system
  101. char* dt = ctime(&now); // convert now to string form
  102. cout << "The local date and time is: " << dt << endl;
  103. cout << "Event Number " << eventCounter << endl;
  104. cout << endl;
  105. //pi = new TClonesArray("MpdParticle"); //put in exec to check for mem leak
  106. //pr = new TClonesArray("MpdParticle");
  107. Reset();
  108. //GetPrimaryVertex();
  109. MpdVertex *PrimVtx = (MpdVertex*) fPrimVertex->UncheckedAt(0);
  110. PVert.SetXYZ( PrimVtx->GetX(),PrimVtx->GetY(),PrimVtx->GetZ());
  111. FindMCVertices(); // find mc vertices with decay daughters in tpc
  112. cout<<"MC lambdas = "<<Nlam1<<endl;
  113. cout<<"Lambdas with both daughters in tpc = "<<Nmc1<<endl;
  114. //FillMpdTracks(); //like in fill dst just used for probability
  115. // MpdEvent *ev = new MpdEvent();
  116. // ev = (MpdEvent*) fEvent->UncheckedAt(0);
  117. // cout<<fEvent->GetEntries()<<endl;
  118. // MpdEvent *ev = (MpdEvent*) fEvent->UncheckedAt(1);
  119. cout<<"fEvent->GetRunInfoRunId() = "<<fEvent->GetRunInfoRunId()<<endl;
  120. fMpdTracks = (TClonesArray*) fEvent->GetGlobalTracks();
  121. fNMPDtracks = fMpdTracks->GetEntriesFast();
  122. cout<<" GlobalTracks = "<<fNMPDtracks<<endl;
  123. //cout << "MpdTracks = " << fNMPDtracks << endl;
  124. SelectTracks();
  125. //Int_t nprots = pr.size();
  126. //Int_t npions = pi.size();
  127. Int_t nprots = pr->GetEntriesFast();
  128. Int_t npions = pi->GetEntriesFast();
  129. cout<<"Protons for lambda = "<<nprots<<endl;
  130. cout<<"Pions for lambda = "<<npions<<endl;
  131. vector<MpdParticle*> cont; //container for daughters
  132. for (Int_t i = 0; i < nprots; i++) {
  133. //MpdParticle *prot = pr.at(i);
  134. MpdParticle *prot = (MpdParticle*) pr->UncheckedAt(i);
  135. Int_t ch1 = prot->GetCharge();
  136. //Int_t B = ch1; //baryon number
  137. if (ch1 > 0) prot->SetPdg(protonPDG);
  138. else prot->SetPdg(-protonPDG);
  139. prot->SetMass();
  140. for (Int_t j = 0; j < npions; j++) {
  141. //MpdParticle *pion = pi.at(j);
  142. MpdParticle *pion = (MpdParticle*) pi->UncheckedAt(j);
  143. Int_t ch2 = pion->GetCharge();
  144. if (ch1*ch2 > 0) continue;
  145. if (ch2 < 0) pion->SetPdg(-pionPDG);
  146. else pion->SetPdg(pionPDG);
  147. pion->SetMass();
  148. cont.clear();
  149. cont.push_back(prot);
  150. cont.push_back(pion);
  151. MpdParticle lambda;
  152. //cout<<i<<" "<<j<<endl;
  153. Double_t chi2 = abs(lambda.BuildMother(cont)); //why abs?
  154. if (chi2 > V0Chi2Cut) continue;
  155. Double_t mass = lambda.GetMass();
  156. if (minmassCut > mass || mass > maxmassCut)continue;
  157. Double_t DCA = lambda.Dca();
  158. if (abs(DCA) > DCAV0Cut) continue;
  159. vtx.SetXYZ(lambda.Getx()(0,0), lambda.Getx()(1,0), lambda.Getx()(2,0));
  160. Double_t xyz[3];
  161. vtx.GetXYZ(xyz);
  162. TVector3 R = vtx-PVert; //decay distance from primvert
  163. Double_t r[3];
  164. R.GetXYZ(r);
  165. if (RminCut > R.Mag() || R.Mag() > RmaxCut) continue;
  166. TVector3 P = lambda.Momentum3();
  167. if (P.Mag() < PminCut) continue;
  168. Double_t p[3];
  169. P.GetXYZ(p);
  170. Double_t cosa = (R.Dot(P)) / (P.Mag() * R.Mag()); //pointing angle
  171. if (cosa < cosaCut) continue; //negative cosa and cut or absolute?a lot of true are with positive cosa
  172. //TMatrixD covar = lambda.GetC();//MpdKfV0Fitter::Instance()->GetCovariance();
  173. MpdKalmanTrack *tr1 = (MpdKalmanTrack*) fKFTracks->UncheckedAt(prot->GetIndx()); //this is necessary because mcid and reco number should be separated
  174. MpdKalmanTrack *tr2 = (MpdKalmanTrack*) fKFTracks->UncheckedAt(pion->GetIndx());
  175. Int_t id0 = CheckMCTracks(tr1, tr2);
  176. // cout<<"nvert = "<<nvert<<endl;
  177. V0 = new ((*fV0Cont)[nvert]) MpdV0();
  178. V0->SetV0(xyz, chi2, p, DCA, cosa, mass, r , ch1, prot->GetIndx(), pion->GetIndx(), id0); // ndf = 2n-3(is there apriory? ) or 2n so its either 4 or 1??? B is baryon number
  179. ++nvert; //gets reset for every event
  180. ++Nreco; //accumulates
  181. cout<<".";
  182. }
  183. }
  184. cout<<endl;
  185. cout<<"Reconstructed "<<nvert<<" lambda particles."<<endl;
  186. Double_t Eff1 = (Double_t)Nok1/(Double_t)Nmc1;
  187. //Eff1vsZ->Fill(pvert[2],Eff1);
  188. Double_t Bad1 = nvert-Nok1;
  189. Double_t Cont1 = (Double_t)Bad1/(Double_t)nvert;
  190. cout<<"Eff1 = "<<Eff1<<endl; //pseudo-efficiency of lambda for this event
  191. cout<<"Cont1 = "<<Cont1<<endl;
  192. }
  193. //__________________________________________________________________________
  194. void MpdV0Finder::SelectTracks(){
  195. MpdTpcKalmanFilter *recoTpc = new MpdTpcKalmanFilter();//for refit where should refit?
  196. //recoTpc->SetSectorGeo(MpdTpcSectorGeo::Instance());
  197. //recoTpc->FillGeoScheme();
  198. for (Int_t i = 0; i < fNMPDtracks; i++){
  199. MpdKalmanTrack *tr = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
  200. if (tr->GetNofHits() < NofHitsCut) continue;
  201. MpdTrack *mpdtr = (MpdTrack*) fMpdTracks->UncheckedAt(i);
  202. Int_t n = GetHighestProb(mpdtr);
  203. if (n == 0 || n == 1 || n == 4) continue; //no need for electrons or kaons
  204. MpdKalmanTrack trRefit = *tr;
  205. Double_t chi2v = tr->GetChi2Vertex();
  206. Double_t Pt = tr->Pt();
  207. Double_t Mom = tr->Momentum();
  208. //protons
  209. if (n == 2){
  210. if (protonPtminCut > Pt || Pt > protonPtmaxCut) continue;
  211. if (protonPminCut > Mom || Mom > protonPmaxCut) continue;
  212. if (protonChi2vminCut > chi2v || chi2v > protonChi2vmaxCut) continue;
  213. if (PCA(tr, PVert) < protonIPCut) continue; // cut on particle or charge???
  214. recoTpc->Refit(&trRefit, 0.93827, 1); // refit // Where should I refit???
  215. MpdParticle *part = new ((*pr)[npr]) MpdParticle(*tr,i); //MC ID or reco number??? i or GetTrackID <-track id is bad because later you have to loop and check number of track
  216. // MpdParticle *part = new ((*pr)[npr]) MpdParticle(trRefit,i);
  217. npr++;
  218. //delete part;
  219. //pions
  220. }else{
  221. if (pionPtminCut > Pt || Pt > pionPtmaxCut) continue;
  222. if (pionPminCut > Mom || Mom > pionPmaxCut) continue;
  223. if (pionChi2vminCut > chi2v || chi2v > pionChi2vmaxCut) continue;
  224. if (PCA(tr, PVert) < pionIPCut) continue;
  225. //recoTpc->Refit(&trRefit, 0.93827, 1);
  226. MpdParticle *part = new ((*pi)[npi]) MpdParticle(*tr,i);
  227. //MpdParticle *part = new ((*pr)[npr]) MpdParticle(&trRefit,i);
  228. npi++;
  229. //delete part;
  230. }
  231. // delete part; //check if mem leak is fixed
  232. }//trackloop
  233. delete recoTpc; //check if mem leak is fixed
  234. }
  235. //__________________________________________________________________________
  236. Double_t MpdV0Finder::PCA(MpdKalmanTrack *tr, TVector3 point) {
  237. // Create MpdHelix for 3D DCA as done in AnalXiFull
  238. Double_t r = tr->GetPosNew();
  239. Double_t phi = tr->GetParam(0) / r;
  240. Double_t x = r * TMath::Cos(phi);
  241. Double_t y = r * TMath::Sin(phi);
  242. Double_t dip = tr->GetParam(3);
  243. Double_t cur = 0.3 * 0.01 * 5 / 10; // 5 kG
  244. cur *= TMath::Abs(tr->GetParam(4));
  245. TVector3 o(x, y, tr->GetParam(1));
  246. Int_t h = (Int_t) TMath::Sign(1.1, tr->GetParam(4));
  247. MpdHelix helix(cur, dip, tr->GetParam(2) - TMath::PiOver2() * h, o, h);
  248. TVector3 pca;
  249. Double_t s = helix.pathLength(point);
  250. pca = helix.at(s);
  251. pca -= point;
  252. Double_t dist = pca.Mag();
  253. return dist;
  254. }
  255. //__________________________________________________________________________
  256. Int_t MpdV0Finder::GetHighestProb(MpdTrack *track) {
  257. //returns the most probable particle. a cut can still be applied
  258. Double_t electron = track->GetPidProbElectron();
  259. Double_t proton = track->GetPidProbProton();
  260. Double_t pion = track->GetPidProbPion();
  261. Double_t kaon = track->GetPidProbKaon();
  262. Double_t Probs[] = {0,electron, proton, pion, kaon};
  263. Double_t maxprob = *max_element(Probs, Probs + 5);
  264. Int_t i = distance(Probs, max_element(Probs, Probs + 5));
  265. if (maxprob<ProbCut) return 0;
  266. else return i;
  267. }
  268. //__________________________________________________________________________
  269. /*
  270. void MpdV0Finder::FillMpdTracks() {
  271. //only for barrel, done for pid (needs to be a separate task in reco imho)
  272. Int_t nReco = fKFTracks ? fKFTracks->GetEntriesFast() : 0;
  273. Int_t nMatching = fTofMatching ? fTofMatching->GetEntriesFast() : 0;
  274. MpdTofMatchingData *pMatchingData;
  275. bool matchingDataExist;
  276. Int_t m = 0;
  277. for (Int_t i = 0; i < nReco; i++) {
  278. MpdKalmanTrack *kftrack = (MpdKalmanTrack*) fKFTracks->UncheckedAt(i);
  279. MpdTrack *track = new ((*fMpdTracks)[m]) MpdTrack(); //instead of AddGlobalTrack()
  280. MpdParticleIdentification *identificator = new MpdParticleIdentification();
  281. track->SetID(kftrack->GetTrackID());
  282. Float_t Ppi, Pk, Pe, Pp;
  283. if (!identificator->GetTpcProbs(kftrack->Momentum3().Mag(), kftrack->GetPartID(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 0)) {
  284. track->SetTPCpidProb(Pe, Ppi, Pk, Pp, BIT(2)); }
  285. matchingDataExist = false;
  286. for (Int_t tofIndex = 0; tofIndex < nMatching; tofIndex++) {
  287. pMatchingData = (MpdTofMatchingData*) fTofMatching->UncheckedAt(tofIndex);
  288. if (pMatchingData->GetKFTrackIndex() == i) {matchingDataExist = true; break;}
  289. } // first matching
  290. if (matchingDataExist) {
  291. track->SetTofMass2(pMatchingData->GetMass2());
  292. track->SetTofHitIndex(pMatchingData->GetTofHitIndex());
  293. if (!identificator->GetTofProbs(pMatchingData->GetMomentum().Mag(), pMatchingData->GetMass2(), kftrack->GetNofHits(), Ppi, Pk, Pp, Pe, 0)) {
  294. track->SetTOFpidProb(Pe, Ppi, Pk, Pp, BIT(1)); }
  295. }
  296. Float_t tpcProbs[4] = {track->GetTPCPidProbPion(), track->GetTPCPidProbKaon(), track->GetTPCPidProbProton(), track->GetTPCPidProbElectron()};
  297. Float_t tofProbs[4] = {track->GetTOFPidProbPion(), track->GetTOFPidProbKaon(), track->GetTOFPidProbProton(), track->GetTOFPidProbElectron()};
  298. Float_t combProbs[4]; //probabilities combined from TOF & TPC
  299. identificator->GetCombinedProbs(tofProbs, tpcProbs, combProbs, 4);
  300. Ppi = combProbs[0]; Pk = combProbs[1]; Pp = combProbs[2]; Pe = combProbs[3];
  301. track->SetCombPidProb(Pe, Ppi, Pk, Pp);
  302. m++;
  303. }
  304. }
  305. */
  306. //__________________________________________________________________________
  307. Int_t MpdV0Finder::CheckMCTracks(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2) {
  308. //checks if the found vertex of these two tracks is a MC(true) vertex primary lambda have mum -1
  309. FairMCTrack* mctr1 = (FairMCTrack*) fMCTracks->UncheckedAt(tr1->GetTrackID());
  310. FairMCTrack* mctr2 = (FairMCTrack*) fMCTracks->UncheckedAt(tr2->GetTrackID());
  311. Int_t mum1 = mctr1->GetMotherId();
  312. Int_t mum2 = mctr2->GetMotherId();
  313. if (mum1 == -1 || mum2 == -1) return -2; // we have a primary track, so vert is not true
  314. TVector3 vert1, vert2;
  315. mctr1->GetStartVertex(vert1);
  316. mctr2->GetStartVertex(vert2);
  317. if (vert1 != vert2) return -2; // vertex matches
  318. FairMCTrack* mcmumTr1 = (FairMCTrack*) fMCTracks->UncheckedAt(mum1);
  319. FairMCTrack* mcmumTr2 = (FairMCTrack*) fMCTracks->UncheckedAt(mum2);
  320. if (abs(mcmumTr1->GetPdgCode()) != lambdaPDG) return -2;
  321. ++Nok;
  322. ++Nok1;
  323. return mum1;
  324. }
  325. //__________________________________________________________________________
  326. void MpdV0Finder::FindMCVertices() {
  327. Int_t fNMC = fMCTracks->GetEntriesFast();
  328. cout << "Nmctracks = " << fNMC << endl;
  329. //find all lambdas first
  330. for (Int_t i = 0; i < fNMC; ++i) {
  331. FairMCTrack* mctr0 = (FairMCTrack*) fMCTracks->UncheckedAt(i);
  332. if (mctr0->GetPdgCode() != lambdaPDG) continue;
  333. ++Nlam;
  334. ++Nlam1;
  335. }
  336. for (Int_t i = 0; i < fNMC; ++i) { //find protons
  337. FairMCTrack* mctr1 = (FairMCTrack*) fMCTracks->UncheckedAt(i);
  338. Int_t mctrpdg1 = mctr1->GetPdgCode();
  339. if (abs(mctrpdg1) != protonPDG) continue;
  340. Int_t mom1 = mctr1->GetMotherId();
  341. if (mom1 == -1)continue;
  342. //is the mommy lambda
  343. if (abs(((FairMCTrack*) fMCTracks->UncheckedAt(mom1))->GetPdgCode()) != lambdaPDG) continue;
  344. for (Int_t j = 0; j < fNMC; ++j) {//find pions for protons
  345. FairMCTrack* mctr2 = (FairMCTrack*) fMCTracks->UncheckedAt(j);
  346. Int_t mctrpdg2 = mctr2->GetPdgCode();
  347. if (abs(mctrpdg2) != pionPDG) continue;
  348. Int_t mom2 = mctr2->GetMotherId();
  349. if (mom2 == -1)continue;
  350. if (abs(((FairMCTrack*) fMCTracks->UncheckedAt(mom2))->GetPdgCode()) != lambdaPDG) continue;
  351. //both daughters have detector points??
  352. if (mctr1->GetNPoints(kTPC) && mctr2->GetNPoints(kTPC)) {
  353. //daughters have same vertex
  354. TVector3 vert1, vert2;
  355. mctr1->GetStartVertex(vert1);
  356. mctr2->GetStartVertex(vert2);
  357. if (vert1 != vert2) continue;
  358. ++Nmc;
  359. ++Nmc1;
  360. //mcPt1->Fill(mctr1->GetPt());
  361. //mcP1 ->Fill(mctr1->GetP());
  362. //mcPt2->Fill(mctr2->GetPt());
  363. //mcP2 ->Fill(mctr2->GetP());
  364. //TVector3 PVert; //primary vertex
  365. //PVert.SetXYZ(pvert[0], pvert[1], pvert[2]);
  366. TVector3 mcR = PVert - vert1;
  367. Double_t dist = mcR.Mag();
  368. //hmcR->Fill(dist);
  369. TLorentzVector lorVec[2];
  370. FairMCTrack * trs[2] = {mctr1, mctr2};
  371. for (Int_t i = 0; i < 2; ++i) {
  372. FairMCTrack *tr = trs[i];
  373. lorVec[i].SetXYZM(tr->GetPx(), tr->GetPy(), tr->GetPz(), tr->GetMass());
  374. }
  375. TLorentzVector lorInv = lorVec[0] + lorVec[1];
  376. TVector3 mcP = lorInv.Vect();
  377. Double_t mccosa = (mcR.Dot(mcP)) / (mcP.Mag() * mcR.Mag());
  378. //mccos->Fill(mccosa);
  379. }
  380. }
  381. }
  382. }
  383. ClassImp(MpdV0Finder);