class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; class RhoTuple; class PndPidCandidate; void psi4160at8GeV10000eventsAnalysis(int nevts=0,double pbarmom = 8.00000) { // some variables int i=0,j=0, k=0, l=0; // the output file examined TString OutFile="output_2.root"; // the files coming from the simulation TString inPidFile = "pid_complete.root"; // this file contains the PndPidCandidates and McTruth TString inParFile = "simparams.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 an output file for all histograms TFile *out = TFile::Open("psi4160_2.root","RECREATE"); // *** take constant field; needed for PocaVtx RhoCalculationTools::ForceConstantBz(20.0); // *** create some ntuples RhoTuple *nd0d0bar = new RhoTuple("nall", "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"); // // Now the analysis stuff comes... // // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nevts==0) nevts= theAnalysis->GetEntries(); // *** name of the only PidAlgo TClonesArray in fsim 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, kplus, kminus, piplus, piminus, all, mclist; // *** Mass selector for the d0 cands double m0_d0 = TDatabasePDG::Instance()->GetParticle("D0")->Mass(); // Get nominal PDG mass of the D0/anti-D0 RhoMassParticleSelector *d0MassSel=new RhoMassParticleSelector("d0",m0_d0,1.0); // *** the lorentz vector of the initial psi(4160) double m0_p = TDatabasePDG::Instance()->GetParticle("proton")->Mass(); // Get nominal PDG mass of the proton TLorentzVector ini(0, 0, pbarmom, sqrt(m0_p*m0_p + pbarmom*pbarmom) + m0_p); // *** // 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(all, "All", pidalg); PndEventShape evsh(all, ini, 0.05, 0.1); // *** Select with no PID info ('All'); type and mass are set theAnalysis->FillList(kplus, "KaonLoosePlus", pidalg); theAnalysis->FillList(kminus, "KaonLooseMinus", pidalg); theAnalysis->FillList(piplus, "PionLoosePlus", pidalg); theAnalysis->FillList(piminus, "PionLooseMinus", pidalg); /* // some mass pre selection d0.Select(d0MassSel); */ // *** combinatorics for d0 -> k- pi+ *** // d0.Combine(kminus, piplus); d0.SetType(421); int nd0mct = theAnalysis->McTruthMatch(d0); for (j=0;jColumn("ev", (Int_t) i); nd0->Column("cand", (Float_t) j); nd0->Column("ncand", (Float_t) d0.GetLength()); nd0->Column("nmct", (Float_t) nd0mct); // store info about initial 4-vector qa.qaP4("beam", ini, nd0); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("d0", d0[j], nd0); // store info about event shapes qa.qaEventShapeShort("es",&evsh, nd0); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate *truth = (d0[j]->GetMcTruth()); TLorentzVector lv; if (truth) lv = truth->P4(); qa.qaP4("trd0", lv, nd0); PndKinVtxFitter *vtxfitter=new PndKinVtxFitter(d0[j]); // instantiate a vertex fitter vtxfitter->FitAll(); RhoCandidate *d0fitvtx = d0[j]->GetFit(); //store info about vertex qa.qaVtx("vtx", d0fitvtx, nd0); Float_t chi2_vtx = vtxfitter->GetChi2(); // access chi2 of fit Float_t prob_vtx = vtxfitter->GetProb(); // access probability of fit nd0->Column("vtxprob", (Float_t) prob_vtx); nd0->Column("vtxchi2", (Float_t) chi2_vtx); PndKinFitter *mfitter=new PndKinFitter(d0[j]); mfitter->AddMassConstraint(m0_d0); // add the mass constraint 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); // // 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 mfitter; delete vtxfitter; nd0->DumpData(); } /* // some mass pre selection antid0.Select(d0MassSel); */ // *** combinatorics for antid0 -> k+ pi- *** // antid0.Combine(kplus, piminus); antid0.SetType(-421); int nantid0mct = theAnalysis->McTruthMatch(antid0); for (j=0;jColumn("ev", (Int_t) i); nantid0->Column("cand", (Float_t) j); nantid0->Column("ncand", (Float_t) antid0.GetLength()); nantid0->Column("nmct", (Float_t) nantid0mct); // store info about initial 4-vector qa.qaP4("beam", ini, nantid0); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("antid0", antid0[j], nantid0); // store info about event shapes qa.qaEventShapeShort("es",&evsh, nantid0); //store info about vertex qa.qaVtx("vtx", antid0[j], nantid0); // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate *truth = (antid0[j]->GetMcTruth()); TLorentzVector lv; if (truth) lv = truth->P4(); qa.qaP4("trantid0", lv, nantid0); PndKinVtxFitter *vtxfitter=new PndKinVtxFitter(antid0[j]); // instantiate a vertex fitter vtxfitter->FitAll(); RhoCandidate *antid0fitvtx = antid0[j]->GetFit(); //store info about vertex qa.qaVtx("vtx", antid0fitvtx, nantid0); Float_t chi2_vtx = vtxfitter->GetChi2(); // access chi2 of fit Float_t prob_vtx = vtxfitter->GetProb(); // access probability of fit nantid0->Column("vtxprob", (Float_t) prob_vtx); nantid0->Column("vtxchi2", (Float_t) chi2_vtx); PndKinFitter *mfitter=new PndKinFitter(antid0[j]); mfitter->AddMassConstraint(m0_d0); // add the mass constraint 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); // // missing particle reconstruction // TLorentzVector lvf(antid0fit->GetMomentum(),antid0fit->GetEnergy()); TLorentzVector lmissingf = ini - lvf; TLorentzVector lv(antid0[j]->GetMomentum(),antid0[j]->GetEnergy()); TLorentzVector lmissing = ini - lv; qa.qaP4("d0", lmissing, nantid0); qa.qaP4("d0fit", lmissingf, nantid0); delete mfitter; delete vtxfitter; nantid0->DumpData(); } // *** combinatorics for psi4160 -> 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()); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("psi", d0d0bar[j], nd0d0bar); 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("psifitchi2", (Float_t) d0d0barchi2); nd0d0bar->Column("psifitprob", (Float_t) d0d0barprob); nd0d0bar->Column("psifitNdf", (Float_t) d0d0barndf); // store the fitted tree of particles as well, JGM, May 2014 qa.qaComp("psifit", d0d0barfit, nd0d0bar); 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); // store the fitted tree of particles as well, JGM, May 2014 qa.qaComp("d0fit", d0fit, nd0d0bar); // store info about missing antid0 particle, before and after fit, JGM, May 2014 TLorentzVector ld0(d0d0bar[j]->Daughter(0)->P4()); TLorentzVector ld0missing = ini - ld0; TLorentzVector lfd0(d0fit->P4()); TLorentzVector lfd0missing = ini - lfd0; qa.qaP4("d0miss",ld0missing, nd0d0bar); qa.qaP4("d0missfit",lfd0missing, nd0d0bar); 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); // store the fitted tree of particles as well, JGM, May 2014 qa.qaComp("antid0fit", antid0fit, 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("antid0miss",lantid0missing, nd0d0bar); qa.qaP4("antid0missfit",lfantid0missing, nd0d0bar); // store info about initial 4-vector qa.qaP4("beam", ini, nd0d0bar); // store info about event shapes qa.qaEventShapeShort("es",&evsh, nd0d0bar); // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate *truth = d0d0bar[j]->GetMcTruth(); TLorentzVector lv; if (truth) lv = truth->P4(); qa.qaP4("trpsi", lv, nd0d0bar); nd0d0bar->DumpData(); delete kinfit; delete mfitter; delete mfitter2; } } // *** write out all the histos out->cd(); nd0->GetInternalTree()->Write(); nantid0->GetInternalTree()->Write(); nd0d0bar->GetInternalTree()->Write(); nmc->GetInternalTree()->Write(); out->Save(); }