class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; class RhoTuple; class PndPidCandidate; void psi4160at8GeV10000eventsAnalysis(int nevts=0) { // some variables int i=0,j=0, k=0, l=0; // the output file examined TString InFile="psi4160_fast.root"; TString OutFile="psi4160_out.root"; // initialization FairLogger::GetLogger()->SetLogToFile(kFALSE); FairRunAna* fRun = new FairRunAna(); fRun->SetWriteRunInfoFile(kFALSE); fRun->SetInputFile(InFile); fRun->SetOutputFile(OutFile); fRun->Init(); RhoCalculationTools::ForceConstantBz(20.0); // create an output file for all histograms TFile *out = TFile::Open("psi4160.root","RECREATE"); // create ntuples for psi(4160), D0 and anti-D0 RhoTuple *npsi4160 = new RhoTuple("npsi4160","npsi4160 Analysis"); RhoTuple *nd0 = new RhoTuple("nd0","nD0 Analysis"); RhoTuple *nantid0 = new RhoTuple("nantid0","nanti-D0 Analysis"); // *** Now the analysis stuff *** // // the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if(nevts==0) nevts= theAnalysis->GetEntries(); // RhoCandLists for the analysis RhoCandList psi4160, d0, antid0, kplus, kminus, piplus, piminus, all; // Mass selector for the psi4160, do/anti-d0, K+/K-, pi+/pi- 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); // Pid Selection Algorithms TString pidSelection = "PidChargedProbability"; // the lorentz vector of the initial psi(4160) TLorentzVector ini(0, 0, 8.000, 9.016); // *** the event loop *** // while (theAnalysis->GetEvent() && i++FillList(all, "All", pidSelection); PndEventShape evsh(all, ini, 0.05, 0.1); theAnalysis->FillList(kplus, "KaonLoosePlus", pidSelection); theAnalysis->FillList(kminus, "KaonLooseMinus", pidSelection); theAnalysis->FillList(piplus, "PionLoosePlus", pidSelection); theAnalysis->FillList(piminus, "PionLooseMinus", pidSelection); /* // some mass pre selection d0.Select(d0MassSel); antid0.Select(d0MassSel); */ // *** combinatorics for d0 -> k- pi+ *** // d0.Combine(kminus, piplus); d0.SetType(421); // *** do all kind of analysis and store in N-tuple *** // for (j=0;jDaughter(0); RhoCandidate *pipl = d0[j]->Daughter(1); PndPidCandidate *kmin_rec = (PndPidCandidate*)kmin->GetRecoCandidate(); PndPidCandidate *pipl_rec = (PndPidCandidate*)pipl->GetRecoCandidate(); // get truth information bool mct = theAnalysis->McTruthMatch(d0[j]); RhoCandidate *true_d0 = d0[j]->GetMcTruth(); // perform vertex fitter PndKinVtxFitter vtxfitter(d0[j]); // instantiate a vertex fitter vtxfitter.Fit(); RhoCandidate *fitvtx_d0 = d0[j]->GetFit(); double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit double prob_vtx = vtxfitter.GetProb(); // access probability of fit TVector3 vtxpos(-999.,-999.,-999.); if (fitvtx_d0) vtxpos = fitvtx_d0->Daughter(0)->Pos(); // perform mass fit PndKinFitter mfitter(d0[j]); // instantiate the PndKinFitter in d0 mfitter.AddMassConstraint(m0_d0); // add the mass constraint mfitter.Fit(); // do fit RhoCandidate *fitmass_d0 = d0[j]->GetFit(); double chi2_mass = mfitter.GetChi2(); // get chi2 of fit double prob_mass = mfitter.GetProb(); // access probability of fit // now write ntuple information // general event info nd0->Column("ev", (Float_t) i, -999.9f); nd0->Column("cand", (Float_t) j, -999.9f); // basic d0 info nd0->Column("d0pdg", (Float_t) d0[j]->PdgCode(), -999.9f); nd0->Column("d0e", (Float_t) d0[j]->E(), -999.9f); nd0->Column("d0m", (Float_t) d0[j]->M(), -999.9f); nd0->Column("d0p", (Float_t) d0[j]->P(), -999.9f); nd0->Column("d0pt", (Float_t) d0[j]->P3().Pt(), -999.9f); nd0->Column("d0tht", (Float_t) d0[j]->P3().Theta()*57.30, -999.9f); nd0->Column("d0missm", (Float_t) (ini-(d0[j]->P4())).M(), -999.9f); nd0->Column("d0track", (Float_t) d0[j]->GetTrackNumber(), -999.9f); // MC truth info nd0->Column("mct", (Float_t) mct, -999.9f); if (true_d0) { nd0->Column("td0m", (Float_t) true_d0->M(), -999.9f); nd0->Column("td0p", (Float_t) true_d0->P(), -999.9f); nd0->Column("td0tht", (Float_t) true_d0->P3().Theta()*57.30, -999.9f); } else { nd0->Column("td0m", (Float_t) d0[j]->M(), -999.9f); nd0->Column("td0p", (Float_t) d0[j]->P(), -999.9f); nd0->Column("td0tht", (Float_t) d0[j]->P3().Theta()*57.30, -999.9f); } // fitting info nd0->Column("d0mvtx", (Float_t) fitvtx_d0->M(), -999.9f); nd0->Column("chi2vtx", (Float_t) chi2_vtx, -999.9f); nd0->Column("probvtx", (Float_t) prob_vtx, -999.9f); nd0->Column("vtxx", vtxpos.X(), -999.9f); nd0->Column("vtxy", vtxpos.Y(), -999.9f); nd0->Column("vtxz", vtxpos.Z(), -999.9f); nd0->Column("d0mmass", (Float_t) fitmass_d0->M(), -999.9f); nd0->Column("chi2mass", (Float_t) chi2_mass, -999.9f); nd0->Column("probmass", (Float_t) prob_mass, -999.9f); // kinematic info of daughters nd0->Column("kmip", (Float_t) kmin->P(), -999.9f); nd0->Column("kmipt", (Float_t) kmin->P3().Pt(), -999.9f); nd0->Column("kmitht", (Float_t) kmin->P3().Theta()*57.30, -999.9f); nd0->Column("pipp", (Float_t) pipl->P(), -999.9f); nd0->Column("pippt", (Float_t) pipl->P3().Pt(), -999.9f); nd0->Column("piptht", (Float_t) pipl->P3().Theta()*57.30, -999.9f); // PID info of daughters nd0->Column("kmipid", (Float_t) kmin->GetPidInfo(3), -999.9f); nd0->Column("pippid", (Float_t) pipl->GetPidInfo(2), -999.9f); // and finally FILL Ntuple nd0->DumpData(); } // *** combinatorics for anti-d0 -> k+ pi- *** // antid0.Combine(kplus, piminus); antid0.SetType(-421); // *** do all kind of analysis and store in N-tuple *** // for (j=0;jDaughter(0); RhoCandidate *pimi = antid0[j]->Daughter(1); PndPidCandidate *kplu_rec = (PndPidCandidate*)kplu->GetRecoCandidate(); PndPidCandidate *pimi_rec = (PndPidCandidate*)pimi->GetRecoCandidate(); // get truth information bool mct = theAnalysis->McTruthMatch(antid0[j]); RhoCandidate *true_antid0 = antid0[j]->GetMcTruth(); // perform vertex fitter PndKinVtxFitter vtxfitter(antid0[j]); // instantiate a vertex fitter vtxfitter.Fit(); RhoCandidate *fitvtx_antid0 = antid0[j]->GetFit(); double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit double prob_vtx = vtxfitter.GetProb(); // access probability of fit TVector3 vtxpos(-999.,-999.,-999.); if (fitvtx_antid0) vtxpos = fitvtx_antid0->Daughter(0)->Pos(); // perform mass fit PndKinFitter mfitter(antid0[j]); // instantiate the PndKinFitter in antid0 mfitter.AddMassConstraint(m0_d0); // add the mass constraint mfitter.Fit(); // do fit RhoCandidate *fitmass_antid0 = antid0[j]->GetFit(); double chi2_mass = mfitter.GetChi2(); // get chi2 of fit double prob_mass = mfitter.GetProb(); // access probability of fit // now write ntuple information // general event info nantid0->Column("ev", (Float_t) i, -999.9f); nantid0->Column("cand", (Float_t) j, -999.9f); // basic antid0 info nantid0->Column("antid0pdg", (Float_t) antid0[j]->PdgCode(), -999.9f); nantid0->Column("antid0e", (Float_t) antid0[j]->E(), -999.9f); nantid0->Column("antid0m", (Float_t) antid0[j]->M(), -999.9f); nantid0->Column("antid0p", (Float_t) antid0[j]->P(), -999.9f); nantid0->Column("antid0tht", (Float_t) antid0[j]->P3().Theta()*57.30, -999.9f); nantid0->Column("antid0pt", (Float_t) antid0[j]->P3().Pt(), -999.9f); nantid0->Column("antid0missm", (Float_t) (ini-(antid0[j]->P4())).M(), -999.9f); nantid0->Column("antid0track", (Float_t) antid0[j]->GetTrackNumber(), -999.9f); // MC truth info nantid0->Column("mct", (Float_t) mct, -999.9f); if (true_antid0) { nantid0->Column("tantid0m", (Float_t) true_antid0->M(), -999.9f); nantid0->Column("tantid0p", (Float_t) true_antid0->P(), -999.9f); nantid0->Column("tantid0tht", (Float_t) true_antid0->P3().Theta()*57.30, -999.9f); } else { nantid0->Column("tantid0m", (Float_t) antid0[j]->M(), -999.9f); nantid0->Column("tantid0p", (Float_t) antid0[j]->P(), -999.9f); nantid0->Column("tantid0tht", (Float_t) antid0[j]->P3().Theta()*57.30, -999.9f); } // fitting info nantid0->Column("antid0mvtx", (Float_t) fitvtx_antid0->M(), -999.9f); nantid0->Column("chi2vtx", (Float_t) chi2_vtx, -999.9f); nantid0->Column("probvtx", (Float_t) prob_vtx, -999.9f); nantid0->Column("vtxx", vtxpos.X(), -999.9f); nantid0->Column("vtxy", vtxpos.Y(), -999.9f); nantid0->Column("vtxz", vtxpos.Z(), -999.9f); nantid0->Column("antid0mmass", (Float_t) fitmass_antid0->M(), -999.9f); nantid0->Column("chi2mass", (Float_t) chi2_mass, -999.9f); nantid0->Column("probmass", (Float_t) prob_mass, -999.9f); // kinematic info of daughters nantid0->Column("kplp", (Float_t) kplu->P(), -999.9f); nantid0->Column("kplpt", (Float_t) kplu->P3().Pt(), -999.9f); nantid0->Column("kpltht", (Float_t) kplu->P3().Theta()*57.30, -999.9f); nantid0->Column("pimp", (Float_t) pimi->P(), -999.9f); nantid0->Column("pimpt", (Float_t) pimi->P3().Pt(), -999.9f); nantid0->Column("pimtht", (Float_t) pimi->P3().Theta()*57.30, -999.9f); // PID info of daughters nantid0->Column("kplpid", (Float_t) kplu->GetPidInfo(3), -999.9f); nantid0->Column("pimpid", (Float_t) pimi->GetPidInfo(2), -999.9f); // and finally FILL Ntuple nantid0->DumpData(); } // *** combinatorics for psi4160 -> d0 anti-d0 *** // psi4160.Combine(d0, antid0); psi4160.SetType(60443); // *** do all kind of analysis and store in N-tuple *** // for (j=0;jDaughter(0); RhoCandidate *dantid0 = psi4160[j]->Daughter(1); //PndPidCandidate *dd0_rec = (PndPidCandidate*)dd0->GetRecoCandidate(); //PndPidCandidate *dantid0_rec = (PndPidCandidate*)dantid0->GetRecoCandidate(); // get truth information bool mct = theAnalysis->McTruthMatch(psi4160[j]); RhoCandidate *true_psi4160 = psi4160[j]->GetMcTruth(); // do 4C fit PndKinFitter fitter(psi4160[j]); // instantiate the kin fitter in psi(2S) fitter.Add4MomConstraint(ini); // set 4 constraint fitter.Fit(); // do fit RhoCandidate *fit4c_psi4160 = psi4160[j]->GetFit(); // get fitted psi4160 double chi2_4c = fitter.GetChi2(); // get chi2 of fit double prob_4c = fitter.GetProb(); // access probability of fit // general event info npsi4160->Column("ev", (Float_t) i, -999.9f); npsi4160->Column("cand", (Float_t) j, -999.9f); // basic psi4160 info npsi4160->Column("psi4160m", (Float_t) psi4160[j]->M(), -999.9f); npsi4160->Column("psi4160p", (Float_t) psi4160[j]->P(), -999.9f); npsi4160->Column("psi4160pt", (Float_t) psi4160[j]->P3().Pt(), -999.9f); npsi4160->Column("psi4160tht", (Float_t) psi4160[j]->P3().Theta()*57.30, -999.9f); npsi4160->Column("psi41604c", (Float_t) fit4c_psi4160->M(), -999.9f); // MC truth info npsi4160->Column("mct", (Float_t) mct, -999.9f); if (true_psi4160) { npsi4160->Column("tpsi4160m", (Float_t) true_psi4160->M(), -999.9f); npsi4160->Column("tpsi4160p", (Float_t) true_psi4160->P(), -999.9f); npsi4160->Column("tpsi4160tht", (Float_t) true_psi4160->P3().Theta()*57.30, -999.9f); } else { npsi4160->Column("tpsi4160m", (Float_t) psi4160[j]->M(), -999.9f); npsi4160->Column("tpsi4160p", (Float_t) psi4160[j]->P(), -999.9f); npsi4160->Column("tpsi4160tht", (Float_t) psi4160[j]->P3().Theta()*57.30, -999.9f); } // kinematic info of daughters npsi4160->Column("dd0p", (Float_t) dd0->P(), -999.9f); npsi4160->Column("dd0pt", (Float_t) dd0->P3().Pt(), -999.9f); npsi4160->Column("dd0tht", (Float_t) dd0->P3().Theta()*57.30, -999.9f); npsi4160->Column("dantid0p", (Float_t) dantid0->P(), -999.9f); npsi4160->Column("dantid0pt", (Float_t) dantid0->P3().Pt(), -999.9f); npsi4160->Column("dantid0tht", (Float_t) dantid0->P3().Theta()*57.30, -999.9f); // and finally FILL Ntuple npsi4160->DumpData(); } } // write out all the histos out->cd(); nd0->GetInternalTree()->Write(); nantid0->GetInternalTree()->Write(); npsi4160->GetInternalTree()->Write(); out->Save(); }