// d0d0bar_analysis.c // Created by Alex on 14/09/16. class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; class PndPidCandidate; class RhoTuple; void NumberOfHitsInSubdetector(TString pre, RhoCandidate *c, RhoTuple *n){ PndPidCandidate *pidCand = (PndPidCandidate*)c->GetRecoCandidate(); if(pidCand){ n->Column(pre + "MvdHits", (Int_t) pidCand->GetMvdHits()); n->Column(pre + "MvdDEDX", (Float_t) pidCand->GetMvdDEDX()); n->Column(pre + "SttHits", (Int_t) pidCand->GetSttHits()); n->Column(pre + "SttDEDX", (Float_t) pidCand->GetSttMeanDEDX()); n->Column(pre + "DrcTheta", (Float_t) pidCand->GetDrcThetaC()); n->Column(pre + "DrcThetaError", (Float_t) pidCand->GetDrcThetaCErr()); n->Column(pre + "TofTime", (Float_t) pidCand->GetTofStopTime()); n->Column(pre + "TofLength", (Float_t) pidCand->GetTofTrackLength()); n->Column(pre + "TofM2", (Float_t) pidCand->GetTofM2()); n->Column(pre + "TofQuality", (Float_t) pidCand->GetTofQuality()); } } void qaFinalState(PndRhoTupleQA * qa, TString pre, RhoCandList & cl, Int_t ev, RhoTuple * n) { if (0 == n) return; for (int i = 0; i < cl.GetLength(); ++i) { n->Column("evt", (Int_t) ev); TString strIndex = TString::Format("%i", i); RhoCandidate * c = cl[i]; RhoCandidate * truth = c->GetMcTruth(); n->Column(pre + "idx", (Int_t) i); qa->qaComp(pre, c, n); qa->qaCand(pre + "tr", truth, n); n->Column(pre + "trpdg", (Int_t) truth->PdgCode()); qa->qaMcDiff(pre, c, n); n->DumpData(); } } void qaDaughters(PndRhoTupleQA * qa, TString pre, RhoCandidate * c, RhoTuple * n) { if (0 == n) return; for (int i = 0; i < c->NDaughters(); ++i) { RhoCandidate * daughter = c->Daughter(i); RhoCandidate * daughterTrue = daughter->GetMcTruth(); TString dindex = TString::Format("d%i", i); qa->qaCand(pre + dindex + "tr", daughterTrue, n); n->Column(pre + dindex + "trpdg", (Int_t) daughterTrue->PdgCode()); qa->qaMcDiff(pre + dindex, daughter, n); } } void d0d0bar_analysis_last(int nevts=1000,double pbarmom = 6.5) { // *** some variables int i=0,j=0, k=0, l=0; // *** the output file examined TString OutFile="first_output.root"; // *** the files coming from the simulation TString inPidFile =" d0d0bar_pid.root"; TString inParFile = "d0d0bar_par.root"; // *** PID table with selection thresholds; can be modified by the user TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par"; // *** initialization FairLogger::GetLogger()->SetLogToFile(kFALSE); FairRunAna* fRun = new FairRunAna(); FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); fRun->SetInputFile(inPidFile); // *** setup parameter database FairParRootFileIo* parIO = new FairParRootFileIo(); parIO->open(inParFile); FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo(); parIOPid->open(pidParFile.Data(),"in"); rtdb->setFirstInput(parIO); rtdb->setSecondInput(parIOPid); rtdb->setOutput(parIO); fRun->SetOutputFile(OutFile); fRun->Init(); // *** create output file for all histograms TFile *out = TFile::Open("d0d0bar_new.root","RECREATE"); // *** create the ntuples RhoTuple *npre = new RhoTuple("npre", "info before analysis"); RhoTuple * nprep = new RhoTuple("nprep", "Pre Pos Analysis"); RhoTuple * nprem = new RhoTuple("nprem", "Pre Neg Analysis"); // RhoTuple * nprepip = new RhoTuple("nprepip", "Pre Pi Pos Analysis"); // RhoTuple * nprepim = new RhoTuple("nprepim", "Pre Pi Neg Analysis"); RhoTuple *nd0d0bar = new RhoTuple("nd0d0bar", "exclusive d0 d0bar analysis"); RhoTuple *nd0 = new RhoTuple("nd0", "d0 analysis"); RhoTuple *nantid0 = new RhoTuple("nantid0", "antid0 analysis"); RhoTuple *nmc = new RhoTuple("nmc", "mctruth info"); // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); nevts= theAnalysis->GetEntries(); // *** name of the only PidAlgo TClonesArray in fsim // TString pidalg = "PidAlgoStt;PidAlgoMvd;PidAlgoMdtHardCuts;PidAlgoDrc;PidAlgoDisc;PidAlgoEmcBayes;PidAlgoSciT;PidAlgoRich"; // TString pidalg = "PidAlgoIdealCharged"; // *** QA tool for simple dumping of analysis results in RhoRuple PndRhoTupleQA qa(theAnalysis, pbarmom); // *** RhoCandLists for the analysis RhoCandList d0d0bar, d0, antid0, plus1, minus1, plus2, minus2, allCharged, mclist; // *** Mass selector for the d0/d0bar cands double m0_d0 = TDatabasePDG::Instance()->GetParticle("D0")->Mass(); RhoMassParticleSelector *d0MassSel=new RhoMassParticleSelector("d0",m0_d0,0.3); // *** Mass of the proton double m0_p = TDatabasePDG::Instance()->GetParticle("proton")->Mass(); // *** the lorentz vector of the initial state TLorentzVector ini(0, 0, pbarmom, sqrt(m0_p*m0_p + pbarmom*pbarmom) + m0_p); TLorentzVector mcvector(0, 0, 0, 0); // *** // the event loop // *** while (theAnalysis->GetEvent() && i++FillList(mclist, "McTruth"); nmc->Column("ev", (Int_t) i); qa.qaMcList("", mclist, nmc); nmc->DumpData(); // *** Setup event shape object theAnalysis->FillList(allCharged, "Charged"); PndEventShape evsh(allCharged, ini, 0.01, 0.01); // *** Select without PID theAnalysis->FillList(plus1, "PionAllPlus"); theAnalysis->FillList(minus1, "KaonAllMinus"); theAnalysis->FillList(plus2, "KaonAllPlus"); theAnalysis->FillList(minus2, "PionAllMinus"); npre->Column("ev", (Float_t) i); npre->Column("nd0", (Int_t) d0.GetLength()); npre->Column("nantid0", (Int_t) antid0.GetLength()); npre->Column("nplus", (Int_t) plus1.GetLength()); npre->Column("nminus", (Int_t) minus1.GetLength()); npre->Column("nplus", (Int_t) plus2.GetLength()); npre->Column("nminus", (Int_t) minus2.GetLength()); npre->DumpData(); qaFinalState(&qa, "nplus", plus1, i, nprep); qaFinalState(&qa, "nminus", minus1, i, nprem); qaFinalState(&qa, "nplus", plus2, i, nprep); qaFinalState(&qa, "nminus", minus2, i, nprem); nprep->DumpData(); nprem->DumpData(); // *** combinatorics for d0 -> k- pi+ *** // d0.Combine(minus1,plus1); d0.SetType(421); // *** combinatorics for antid0 -> k+ pi- *** // antid0.Combine(plus2,minus2); antid0.SetType(-421); for (j=0;jColumn("ev", (Int_t) i); nd0->Column("cand", (Float_t) j); nd0->Column("ncand", (Float_t) d0.GetLength()); // *** get daughters RhoCandidate *pl1 = d0[j]->Daughter(1); RhoCandidate *min1 = d0[j]->Daughter(0); NumberOfHitsInSubdetector("plus1_",pl1,nd0); NumberOfHitsInSubdetector("minus1_",min1,nd0); // *** get mctruth info // bool d0mct = theAnalysis->McTruthMatch(d0[j]); // RhoCandidate *trued0 = d0[j]->GetMcTruth(); // nd0->Column("nmct", (Float_t) d0mct); // if (trued0 !=0) // { // qa.qaP4("trd0",trued0->P4(),nd0); // } // else // { // qa.qaP4("trd0",mcvector,nd0); // } qaDaughters(&qa, "d0new", d0[j], nd0); // *** various eventshape things qa.qaEventShapeShort("evsh",&evsh, nd0); // *** instantiate a vertex fitter PndKinVtxFitter *vtxfitter=new PndKinVtxFitter(d0[j]); vtxfitter->Fit(); RhoCandidate *d0fitvtx = d0[j]->GetFit(); Float_t chi2_vtx = vtxfitter->GetChi2(); Float_t prob_vtx = vtxfitter->GetProb(); // *** store info about vertex qa.qaVtx("vtx", d0fitvtx, nd0); nd0->Column("vtxprob", (Float_t) prob_vtx); nd0->Column("vtxchi2", (Float_t) chi2_vtx); TVector3 lv2; if (d0fitvtx) lv2 = (d0fitvtx->Daughter(0))->GetPosition(); nd0->Column("trposx", (Float_t) lv2.x()); nd0->Column("trposy", (Float_t) lv2.y()); nd0->Column("trposz", (Float_t) lv2.z()); // *** poca vtx info qa.qaPoca("poca", d0[j], nd0); // *** store info about initial 4-vector qa.qaP4("beam", ini, nd0); // *** store information about composite candidate tree recursively qa.qaComp("d0", d0[j], nd0); // // *** mc differences for D0 // qa.qaMcDiff("d0mc", d0[j],nd0); // *** mass constraint fit PndKinFitter *mfitter=new PndKinFitter(d0[j]); mfitter->AddMassConstraint(m0_d0); // add the mfitter->Fit(); RhoCandidate *d0fit = d0[j]->GetFit(); Float_t chi2 = mfitter->GetChi2(); // get chi2 of fit Float_t prob = mfitter->GetProb(); // access probability of fit Int_t ndf = mfitter->GetNdf(); nd0->Column("prob", (Float_t) prob); nd0->Column("ndf", (Float_t) ndf); nd0->Column("chi2", (Float_t) chi2); qa.qaComp("d0fit", d0fit, nd0); qaDaughters(&qa, "d0fitnew", d0fit, nd0); // *** missing particle reconstruction TLorentzVector lv(d0[j]->GetMomentum(),d0[j]->GetEnergy()); TLorentzVector lmissing = ini - lv; TLorentzVector lvf(d0fit->GetMomentum(),d0fit->GetEnergy()); TLorentzVector lmissingf = ini - lvf; qa.qaP4("antid0", lmissing, nd0); qa.qaP4("antid0fit", lmissingf, nd0); // *** delete fitters and store data delete mfitter; delete vtxfitter; nd0->DumpData(); } for (j=0;jColumn("ev", (Int_t) i); nantid0->Column("cand", (Float_t) j); nantid0->Column("ncand", (Float_t) antid0.GetLength()); // *** get daughters RhoCandidate *pl1 = antid0[j]->Daughter(0); RhoCandidate *min1 = antid0[j]->Daughter(1); NumberOfHitsInSubdetector("plus1_",pl1,nantid0); NumberOfHitsInSubdetector("minus1_",min1,nantid0); // *** get mctruth info // bool antid0mct = theAnalysis->McTruthMatch(antid0[j]); // RhoCandidate *trueantid0 = antid0[j]->GetMcTruth(); // nantid0->Column("nmct", (Float_t) antid0mct); // if (trueantid0 !=0) // { // qa.qaP4("trantid0",trueantid0->P4(),nantid0); // } // else // { // qa.qaP4("trantid0",mcvector,nantid0); // } qaDaughters(&qa, "antid0new", antid0[j], nantid0); // *** various eventshape things qa.qaEventShapeShort("evsh",&evsh, nantid0); // *** instantiate a vertex fitter PndKinVtxFitter *vtxfitter=new PndKinVtxFitter(antid0[j]); vtxfitter->Fit(); RhoCandidate *antid0fitvtx = antid0[j]->GetFit(); Float_t chi2_vtx = vtxfitter->GetChi2(); Float_t prob_vtx = vtxfitter->GetProb(); // *** store info about vertex qa.qaVtx("vtx", antid0fitvtx, nantid0); nantid0->Column("vtxprob", (Float_t) prob_vtx); nantid0->Column("vtxchi2", (Float_t) chi2_vtx); TVector3 lv2; if (antid0fitvtx) lv2 = (antid0fitvtx->Daughter(0))->GetPosition(); nantid0->Column("trposx", (Float_t) lv2.x()); nantid0->Column("trposy", (Float_t) lv2.y()); nantid0->Column("trposz", (Float_t) lv2.z()); // *** poca vtx info qa.qaPoca("poca", antid0[j], nantid0); // *** store info about initial 4-vector qa.qaP4("beam", ini, nantid0); // *** store information about composite candidate tree recursively qa.qaComp("antid0", antid0[j], nantid0); // // *** mc differences for D0bar // qa.qaMcDiff("antid0mc", antid0[j],nantid0); // *** mass constraint fit PndKinFitter *mfitter=new PndKinFitter(antid0[j]); mfitter->AddMassConstraint(m0_d0); // add the mfitter->Fit(); RhoCandidate *antid0fit = antid0[j]->GetFit(); Float_t chi2 = mfitter->GetChi2(); // get chi2 of fit Float_t prob = mfitter->GetProb(); // access probability of fit Int_t ndf = mfitter->GetNdf(); nantid0->Column("prob", (Float_t) prob); nantid0->Column("ndf", (Float_t) ndf); nantid0->Column("chi2", (Float_t) chi2); qa.qaComp("antid0fit", antid0fit, nantid0); qaDaughters(&qa, "antid0fitnew", antid0fit, nantid0); // *** missing particle reconstruction TLorentzVector lv(antid0[j]->GetMomentum(),antid0[j]->GetEnergy()); TLorentzVector lmissing = ini - lv; TLorentzVector lvf(antid0fit->GetMomentum(),antid0fit->GetEnergy()); TLorentzVector lmissingf = ini - lvf; qa.qaP4("d0", lmissing, nantid0); qa.qaP4("d0fit", lmissingf, nantid0); // *** delete fitters and store data delete mfitter; delete vtxfitter; nantid0->DumpData(); } // *** combinatorics for pbarpsystem -> d0 anti-d0 *** // d0d0bar.Combine(d0, antid0); d0d0bar.SetType(88888); for (j=0;jColumn("ev", (Int_t) i); nd0d0bar->Column("cand", (Float_t) j); nd0d0bar->Column("ncand", (Float_t) d0d0bar.GetLength()); // *** get daughters RhoCandidate *pl1 = (d0d0bar[j]->Daughter(0))->Daughter(1); RhoCandidate *min1 = (d0d0bar[j]->Daughter(0))->Daughter(0); RhoCandidate *pl2 = (d0d0bar[j]->Daughter(1))->Daughter(0); RhoCandidate *min2 = (d0d0bar[j]->Daughter(1))->Daughter(1); NumberOfHitsInSubdetector("plus1_",pl1,nd0d0bar); NumberOfHitsInSubdetector("plus2_",pl2,nd0d0bar); NumberOfHitsInSubdetector("minus1_",min1,nd0d0bar); NumberOfHitsInSubdetector("minus2_",min2,nd0d0bar); // store info about initial 4-vector qa.qaP4("beam", ini, nd0d0bar); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("d0d0bar", d0d0bar[j], nd0d0bar); // // *** mc differences for D0 // qa.qaMcDiff("d0mc", (d0d0bar[j]->Daughter(0)),nd0d0bar); // // // *** mc differences for D0bar // qa.qaMcDiff("antid0mc", (d0d0bar[j]->Daughter(1)),nd0d0bar); // // // // *** get mctruth info // bool d0d0barmct = theAnalysis->McTruthMatch(d0d0bar); // RhoCandidate *true_d0d0bar = d0d0bar[j]->GetMcTruth(); // nd0d0bar->Column("nmct", (Float_t) d0d0barmct); // if (true_d0d0bar !=0) // { // qa.qaP4("trd0d0bar",true_d0d0bar->P4(),nd0d0bar); // qa.qaP4("trd0d0bard0",(true_d0d0bar->Daughter(0))->P4(),nd0d0bar); // qa.qaP4("trd0d0barantid0",(true_d0d0bar->Daughter(1))->P4(),nd0d0bar); // } // // *** various eventshape things qa.qaEventShapeShort("evsh",&evsh, nd0d0bar); // *** 4C kinematic fitter for system PndKinFitter *kinfit=new PndKinFitter(d0d0bar[j]); kinfit->Add4MomConstraint(ini); kinfit->Fit(); RhoCandidate *d0d0barfit = d0d0bar[j]->GetFit(); Float_t d0d0barchi2 = kinfit->GetChi2(); // get chi2 of fit Float_t d0d0barprob = kinfit->GetProb(); // access probability of fit Int_t d0d0barndf = kinfit->GetNdf(); nd0d0bar->Column("d0d0barfitchi2", (Float_t) d0d0barchi2); nd0d0bar->Column("d0d0barfitprob", (Float_t) d0d0barprob); nd0d0bar->Column("d0d0barfitNdf", (Float_t) d0d0barndf); qa.qaComp("d0d0barfit", d0d0barfit, nd0d0bar); // *** mass constraint fit for D0 PndKinFitter *mfitter=new PndKinFitter((d0d0bar[j]->Daughter(0))); // instantiate the PndKinFitter in psi(2S) mfitter->AddMassConstraint(m0_d0); // add the mass constraint mfitter->Fit(); RhoCandidate *d0fit = (d0d0bar[j]->Daughter(0))->GetFit(); Float_t d0chi2 = mfitter->GetChi2(); // get chi2 of fit Float_t d0prob = mfitter->GetProb(); // access probability of fit Int_t d0ndf = mfitter->GetNdf(); nd0d0bar->Column("d0fitchi2", (Float_t) d0chi2); nd0d0bar->Column("d0fitprob", (Float_t) d0prob); nd0d0bar->Column("d0fitNdf", (Float_t) d0ndf); qa.qaComp("d0fit", d0fit, nd0d0bar); // *** mass constraint fit for D0bar PndKinFitter *mfitter2=new PndKinFitter((d0d0bar[j]->Daughter(1))); // instantiate the PndKinFitter in psi(2S) mfitter2->AddMassConstraint(m0_d0); // add the mass constraint mfitter2->Fit(); RhoCandidate *antid0fit = (d0d0bar[j]->Daughter(1))->GetFit(); Float_t antid0chi2 = mfitter2->GetChi2(); // get chi2 of fit Float_t antid0prob = mfitter2->GetProb(); // access probability of fit Int_t antid0ndf = mfitter2->GetNdf(); nd0d0bar->Column("antid0fitchi2", (Float_t) antid0chi2); nd0d0bar->Column("antid0fitprob", (Float_t) antid0prob); nd0d0bar->Column("antid0fitNdf", (Float_t) antid0ndf); qa.qaComp("antid0fit", antid0fit, nd0d0bar); // *** store info about missing antid0 particle, before and after fit TLorentzVector ld0(d0d0bar[j]->Daughter(0)->P4()); TLorentzVector ld0missing = ini - ld0; TLorentzVector lfd0(d0fit->P4()); TLorentzVector lfd0missing = ini - lfd0; qa.qaP4("antid0miss",ld0missing, nd0d0bar); qa.qaP4("antid0missfit",lfd0missing, nd0d0bar); // ***store info about missing d0 particle, before and after fit, JGM, May 2014 TLorentzVector lantid0(d0d0bar[j]->Daughter(1)->P4()); TLorentzVector lantid0missing = ini - lantid0; TLorentzVector lfantid0(antid0fit->P4()); TLorentzVector lfantid0missing = ini - lfantid0; qa.qaP4("d0miss",lantid0missing, nd0d0bar); qa.qaP4("d0missfit",lfantid0missing, nd0d0bar); // *** instantiate a vertex fitter for d0 RhoCandidate *d0exc = (d0d0bar[j])->Daughter(0); PndKinVtxFitter *vtxfitter1 = new PndKinVtxFitter(d0exc); // instantiate a vertex fitter vtxfitter1->Fit(); RhoCandidate *d0excfit = d0exc->GetFit(); qa.qaVtx("d0vtx", d0excfit, nd0d0bar); Float_t chi2_vtx = vtxfitter1->GetChi2(); // access chi2 of fit Float_t prob_vtx = vtxfitter1->GetProb(); // access probability of fit nd0d0bar->Column("d0vtxprob", (Float_t) prob_vtx); nd0d0bar->Column("d0vtxchi2", (Float_t) chi2_vtx); qa.qaPoca("d0poca", d0exc, nd0d0bar); // *** instantiate a vertex fitter for d0bar RhoCandidate *d0barexc = (d0d0bar[j])->Daughter(1); PndKinVtxFitter *vtxfitter2 = new PndKinVtxFitter(d0barexc); // instantiate a vertex fitter vtxfitter2->Fit(); RhoCandidate *d0barexcfit = d0barexc->GetFit(); qa.qaVtx("d0barvtx", d0barexcfit, nd0d0bar); Float_t chi2_vtx = vtxfitter2->GetChi2(); // access chi2 of fit Float_t prob_vtx = vtxfitter2->GetProb(); // access probability of fit nd0d0bar->Column("d0barvtxprob", (Float_t) prob_vtx); nd0d0bar->Column("d0barvtxchi2", (Float_t) chi2_vtx); qa.qaPoca("d0barpoca", d0barexc, nd0d0bar); // *** get position of D0 TVector3 lv2; lv2 = (d0d0bar[j]->Daughter(0))->GetPosition(); nd0d0bar->Column("d0posx", (Float_t) lv2.x()); nd0d0bar->Column("d0posy", (Float_t) lv2.y()); nd0d0bar->Column("d0posz", (Float_t) lv2.z()); // *** get position of D0bar TVector3 lv3; lv3 = (d0d0bar[j]->Daughter(1))->GetPosition(); nd0d0bar->Column("d0barposx", (Float_t) lv3.x()); nd0d0bar->Column("d0barposy", (Float_t) lv3.y()); nd0d0bar->Column("d0barposz", (Float_t) lv3.z()); // *** delete fitters and store data delete vtxfitter1; delete vtxfitter2; delete kinfit; delete mfitter; delete mfitter2; nd0d0bar->DumpData(); } } // *** write out all the histos out->cd(); npre->GetInternalTree()->Write(); nprep->GetInternalTree()->Write(); nprem->GetInternalTree()->Write(); nd0->GetInternalTree()->Write(); nantid0->GetInternalTree()->Write(); nd0d0bar->GetInternalTree()->Write(); nmc->GetInternalTree()->Write(); out->Save(); }