/* * PndXiXibarAnaTask.cxx * * Created on: Sep 23, 2016 * Author: walan603 */ //The header file #include "PndXiXibarAnaTask.h" // C++ headers #include #include // FAIR headers #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairRun.h" #include "FairRuntimeDb.h" // ROOT headers #include "TClonesArray.h" #include "TLorentzVector.h" #include "TVector3.h" #include "TH1F.h" #include "TH2F.h" #include "TParticlePDG.h" // RHO headers #include "RhoCandidate.h" #include "RhoHistogram/RhoTuple.h" #include "RhoFactory.h" #include "RhoMassParticleSelector.h" #include "RhoDecayTreeFitter.h" #include "PndRhoTupleQA.h" // PandaRoot headers #include "PndAnalysis.h" #include "Pnd4CFitter.h" #include "PndKinVtxFitter.h" #include "PndKinFitter.h" #include "PndVtxPoca.h" #include "PndPidCandidate.h" // Custom headers #include "PndSpinObsTools.h" using std::cout; using std::endl; // ----- Default constructor ------------------------------------------- PndXiXibarAnaTask::PndXiXibarAnaTask() : FairTask("PandaXiXibarAnaTask"), fBeamMom(4.6){ } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndXiXibarAnaTask::~PndXiXibarAnaTask() { } // ------------------------------------------------------------------------- // ----- Method to select PID candidates int PndXiXibarAnaTask::SelectTruePid(PndAnalysis *ana, RhoCandList &l) { int removed = 0; for (int ii=l.GetLength()-1;ii>=0;--ii) { if ( !(ana->McTruthMatch(l[ii])) ) { l.Remove(l[ii]); removed++; } } return removed; } void PndXiXibarAnaTask::VertexFit(RhoCandList &l, Double_t prob_cut){ std::map IndexToRank; for (int i = l.GetLength()-1; i >= 0; --i){ PndKinVtxFitter vtxfitter(l.Get(i)); vtxfitter.Fit(); bool failedchi2 = TMath::IsNaN(vtxfitter.GetChi2()); bool failedprob = TMath::IsNaN(vtxfitter.GetProb()); if(!failedchi2 && !failedprob){ if (vtxfitter.GetProb() < prob_cut){ l.Remove(l[i]); } } } } void PndXiXibarAnaTask::MassFit(RhoCandList &l, Double_t mass, Double_t prob_cut){ for (int i = l.GetLength()-1; i >= 0; --i){ PndKinFitter massfitter(l.Get(i)); massfitter.AddMassConstraint(mass); massfitter.Fit(); bool failedchi2 = TMath::IsNaN(massfitter.GetChi2()); bool failedprob = TMath::IsNaN(massfitter.GetProb()); if(!failedchi2 && !failedprob){ if (massfitter.GetProb() < prob_cut){ l.Remove(l[i]); } } } } void PndXiXibarAnaTask::Tree4CFit(RhoCandList &l, TLorentzVector ini, Double_t prob_cut){ for (int i = l.GetLength()-1; i >= 0; --i){ PndKinFitter tree4cfitter(l.Get(i)); //massfitter.AddMassConstraint(mass); tree4cfitter.Add4MomConstraint(ini); tree4cfitter.Fit(); bool failedchi2 = TMath::IsNaN(tree4cfitter.GetChi2()); bool failedprob = TMath::IsNaN(tree4cfitter.GetProb()); if(!failedchi2 && !failedprob){ if (tree4cfitter.GetProb() < prob_cut){ l.Remove(l[i]); } } } } std::map PndXiXibarAnaTask::VertexQaIndex(RhoCandList* candList, Double_t prob_cut=0.01){ /** @brief give back the order of the best chi2 * @details give back the order of the best chi2! 1 means best, 2: second best (same with negative values for bad chi2 ) * @details: Credit to Jenny Puetz for the code */ std::map chi2_good, chi2_bad; for (int j=0; jGetLength(); ++j){ PndKinVtxFitter vtxfitter(candList->Get(j)); vtxfitter.Fit(); bool failedchi2 = TMath::IsNaN(vtxfitter.GetChi2()); bool failedprob = TMath::IsNaN(vtxfitter.GetProb()); if(!failedchi2 && !failedprob){ if (vtxfitter.GetProb() > prob_cut){ //Prob > 0.01 chi2_good[vtxfitter.GetChi2()]=j; } else{ //Prob <= 0.01 chi2_bad[vtxfitter.GetChi2()]=j; } } } std::map::iterator is_good, is_bad; std::map indexBestFit; int running = 0; for (is_good = chi2_good.begin(); is_good != chi2_good.end(); is_good++, running++){ indexBestFit[is_good->second] = running + 1; } running =0; for (is_bad = chi2_bad.begin(); is_bad != chi2_bad.end(); is_bad++, running++){ indexBestFit[is_bad->second] = - (running + 1); } return indexBestFit; } std::map PndXiXibarAnaTask::MassFitQaIndex(RhoCandList* candList, Double_t m0, Double_t prob_cut=0.01){ /** @brief give back the order of the best chi2 for MassFit * @details give back the order of the best chi2 for the MassFit! 1 means best, 2: second best (analoge for bad chi2 with negative values) * @details: Credit to Jenny Puetz for the code */ if(m0==0) std::cout << "Mass is missing for mass fit" << std::endl; std::map chi2_good, chi2_bad; for (int i=0; iGetLength(); i++){ PndKinFitter massfitter(candList->Get(i)); massfitter.AddMassConstraint(m0); massfitter.Fit(); bool failedchi2 = TMath::IsNaN(massfitter.GetChi2()); bool failedprob = TMath::IsNaN(massfitter.GetProb()); if(!failedchi2 && !failedprob){ if (massfitter.GetProb() > prob_cut){ chi2_good[massfitter.GetChi2()]=i; } else{ chi2_bad[massfitter.GetChi2()]=i; } } } std::map::iterator is_good, is_bad; std::map bestMassFit; int run =0; for (is_good = chi2_good.begin(); is_good != chi2_good.end(); is_good++, run++){ bestMassFit[is_good->second] = run + 1; } run = 0; for (is_bad = chi2_bad.begin(); is_bad != chi2_bad.end(); is_bad++, run++){ bestMassFit[is_bad->second] = - (run + 1); } return bestMassFit; } void PndXiXibarAnaTask::CombinedList(RhoCandidate* cand, RhoCandList* combinedList, int pdg_id){ /** * @brief: gives back a list of already combined particles * @details: The function creates a list of already combined particles for the analysis * @details: Credit to Jenny Puetz for the code * * @CHECK IF THIS CODE IS COMPATIBLE WITH AN ANALYSIS WITHOUD PID!! */ for (int daughter=0; daughterNDaughters(); daughter++){ RhoCandidate * daughterCand = cand->Daughter(daughter); if (daughterCand->PdgCode()==pdg_id){ combinedList->Append(daughterCand); } } combinedList->RemoveClones(); } void PndXiXibarAnaTask::NotCombinedList(RhoCandList combinedList, RhoCandList *candList){ for (int j=0; jRemove(combinedCand); } } // ----- Public method Init -------------------------------------------- InitStatus PndXiXibarAnaTask::Init() { // initialize analysis object theAnalysis = new PndAnalysis(); // reset the event counter nevts = 0; //PDG database pdg = new TDatabasePDG(); //Particle masses m0_pi = ((TParticlePDG*)pdg->GetParticle(-211))->Mass(); m0_p = ((TParticlePDG*)pdg->GetParticle(-2212))->Mass(); m0_lam = ((TParticlePDG*)pdg->GetParticle(-3122))->Mass(); m0_xi = ((TParticlePDG*)pdg->GetParticle(-3312))->Mass(); //Decay parameters alpha_xi = -0.458; alpha_xibar = 0.458; alpha_lam = 0.642; alpha_lambar = -0.642; //4-momentum vectors fTarg.SetXYZT(0,0,0,m0_p); fBeam.SetXYZT(0,0,fBeamMom,sqrt(m0_p*m0_p + fBeamMom*fBeamMom)); fIni = fTarg+fBeam; //Mass selector lambdaMassSelector = new RhoMassParticleSelector("lambda0", m0_lam, 0.3); xiMassSelector = new RhoMassParticleSelector("xi", m0_xi, 0.3); // tuples fntpPiMinus = new RhoTuple("ntpPiMinus", "PiMinus info"); fntpPiPlus = new RhoTuple("ntpPiPlus", "PiPlus info"); fntpProton = new RhoTuple("ntpProton", "Proton info"); fntpAntiProton = new RhoTuple("ntpAntiProton", "AntiProton info"); fntpLambda = new RhoTuple("ntpLambda", "Lambda info"); fntpLambdaBar = new RhoTuple("ntpLambdaBar", "LambdaBar info"); fntpXi_PiMinus = new RhoTuple("ntpXiPiMinus", "Xi PiMinus info"); fntpXiBar_PiPlus = new RhoTuple("ntpXiBarPiPlus", "XiBar PiPlus info"); fntpXi = new RhoTuple("ntpXi", "Xi info"); fntpXiBar = new RhoTuple("ntpXiBar", "Xi info"); fntpPbarp = new RhoTuple("ntpPbarp", "Pbarp system info"); // MC tuples fntpMCtruth = new RhoTuple("ntpMCtruth", "MC truth info"); // *** ------------------------- *** // *** Create some 1D histograms *** // *** ------------------------- *** // *** --------------------- *** // *** True MC 1D histograms *** // *** --------------------- *** // *** ------------------------- *** // *** Create some 2D histograms *** // *** ------------------------- *** return kSUCCESS; cout<<"PndXiXibarAnaTask - Initalization successful!"<GetEvent(); if (!(++nevts%100)) cout << "evt "<FillList(p,"ProtonBestPlus","PidAlgoIdealCharged"); theAnalysis->FillList(pbar,"ProtonBestMinus","PidAlgoIdealCharged"); theAnalysis->FillList(piplus,"PionBestPlus","PidAlgoIdealCharged"); theAnalysis->FillList(piminus,"PionBestMinus","PidAlgoIdealCharged"); theAnalysis->FillList(piplus_noncombined,"PionBestPlus","PidAlgoIdealCharged"); theAnalysis->FillList(piminus_noncombined,"PionBestMinus","PidAlgoIdealCharged"); theAnalysis->FillList(mclist,"McTruth"); // Select with no PID info ('All'); type and mass are set //theAnalysis->FillList(p,"ProtonAllPlus"); //theAnalysis->FillList(pbar,"ProtonAllMinus"); //theAnalysis->FillList(piplus,"PionAllPlus"); //theAnalysis->FillList(piminus,"PionAllMinus"); //theAnalysis->FillList(mclist,"McTruth"); //theAnalysis->FillList(piplus_noncombined,"PionAllPlus"); //theAnalysis->FillList(piminus_noncombined,"PionAllMinus"); // *** --------------------------- *** // *** Same all Monte Carlo Tracks *** // *** --------------------------- *** RhoCandidate *mcxi, *mcxibar, *mcxipim, *mcxibarpip, *mclam, *mclambar, *mclambarpip, *mclampim, *mcp, *mcpbar; mcxi = mclist[2]; mcxibar = mclist[1]; mcxipim = mclist[2]->Daughter(1); mcxibarpip = mclist[1]->Daughter(1); mclambar = mclist[1]->Daughter(0); mclam = mclist[2]->Daughter(0); mclambarpip = mclist[1]->Daughter(0)->Daughter(1); mclampim = mclist[2]->Daughter(0)->Daughter(1); mcp = mclist[2]->Daughter(0)->Daughter(0); mcpbar = mclist[1]->Daughter(0)->Daughter(0); // *** ----------------------------------- *** // *** Boost particles to relevant systems *** // *** ----------------------------------- *** //Boost beam and target to CM frame TLorentzVector beamCM = fBeam; beamCM.Boost(-fIni.BoostVector()); TLorentzVector targCM = fTarg; targCM.Boost(-fIni.BoostVector()); //Boost Xi to CM frame TLorentzVector XiP4CM = mcxi->P4(); XiP4CM.Boost(-fIni.BoostVector()); //Boost XiBar to CM frame TLorentzVector XiBarP4CM = mcxibar->P4(); XiBarP4CM.Boost(-fIni.BoostVector()); //Boost Lambda to Xi rest frame TLorentzVector lambdaindecayframe = mclam->P4(); lambdaindecayframe.Boost(-fIni.BoostVector()); lambdaindecayframe.Boost(-XiP4CM.BoostVector()); //Boost LambdaBar to XiBar rest frame TLorentzVector lambdabarindecayframe = mclambar->P4(); lambdabarindecayframe.Boost(-fIni.BoostVector()); lambdabarindecayframe.Boost(-XiBarP4CM.BoostVector()); // *** ------------------------ *** // *** Construct Xi decay frame *** // *** ------------------------ *** std::vector final = GenerateUnitVectors(beamCM,XiP4CM,XiBarP4CM); std::vector init_fr; init_fr.push_back(TVector3(1,0,0)); init_fr.push_back(TVector3(0,1,0)); init_fr.push_back(TVector3(0,0,1)); // *** --------------------------- *** // *** Construct XiBar decay frame *** // *** --------------------------- *** std::vector LambdaBarfinal = GenerateUnitVectors(beamCM,XiBarP4CM,XiBarP4CM); TLorentzVector rot_lam_decayframe = TransformCoords(init_fr,final,lambdaindecayframe); fntpMCtruth->Column("lam0_DecFrame_theta", (Float_t) rot_lam_decayframe.Theta(), 0.0f ); fntpMCtruth->Column("lam0_DecFrame_phi", (Float_t) rot_lam_decayframe.Phi(), 0.0f ); fntpMCtruth->Column("lam0_DecFrame_kx", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(0))), 0.0f ); fntpMCtruth->Column("lam0_DecFrame_ky", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(1))), 0.0f ); fntpMCtruth->Column("lam0_DecFrame_kz", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(2))), 0.0f ); TLorentzVector rot_lambar_decayframe = TransformCoords(init_fr,LambdaBarfinal,lambdabarindecayframe); fntpMCtruth->Column("lam0bar_DecFrame_theta", (Float_t) rot_lambar_decayframe.Theta(), 0.0f ); fntpMCtruth->Column("lam0bar_DecFrame_phi", (Float_t) rot_lambar_decayframe.Phi(), 0.0f ); fntpMCtruth->Column("lam0bar_DecFrame_kx", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(0))), 0.0f ); fntpMCtruth->Column("lam0bar_DecFrame_ky", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(1))), 0.0f ); fntpMCtruth->Column("lam0bar_DecFrame_kz", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(2))), 0.0f ); //Calculate CM angle, weights, spin variables with method of moments Double_t cmAngle = beamCM.Vect().Angle(XiBarP4CM.Vect()); //Use only one weight TMatrixD Cij = GenerateSpinCorrTrigFuncs(cmAngle); Double_t Py = TMath::Sin(cmAngle); Double_t Pybar = TMath::Sin(cmAngle); Double_t weightPy = 1 + alpha_xi*Py*TMath::Cos(final.at(1).Angle(lambdaindecayframe.Vect())); Double_t weightPyBar = 1 + alpha_xibar*Pybar*TMath::Cos(LambdaBarfinal.at(1).Angle(lambdabarindecayframe.Vect())); Double_t weightCxx = 1 + alpha_xi*alpha_xibar*Cij(0,0) * TMath::Cos(final.at(0).Angle(lambdaindecayframe.Vect())) * TMath::Cos(LambdaBarfinal.at(0).Angle(lambdabarindecayframe.Vect())); Double_t weightCyy = 1 + alpha_xi*alpha_xibar*Cij(1,1) * TMath::Cos(final.at(1).Angle(lambdaindecayframe.Vect())) * TMath::Cos(LambdaBarfinal.at(1).Angle(lambdabarindecayframe.Vect())); Double_t weightCzz = 1 + alpha_xi*alpha_xibar*Cij(2,2) * TMath::Cos(final.at(2).Angle(lambdaindecayframe.Vect())) * TMath::Cos(LambdaBarfinal.at(2).Angle(lambdabarindecayframe.Vect())); Double_t weightCxz = 1 + alpha_xi*alpha_xibar*Cij(0,2) * TMath::Cos(final.at(2).Angle(lambdaindecayframe.Vect())) * TMath::Cos(LambdaBarfinal.at(0).Angle(lambdabarindecayframe.Vect())); Double_t weightCzx = 1 + alpha_xi*alpha_xibar*Cij(2,0) * TMath::Cos(final.at(0).Angle(lambdaindecayframe.Vect())) * TMath::Cos(LambdaBarfinal.at(2).Angle(lambdabarindecayframe.Vect())); Double_t eventWeight = SpinWeightGeneral(final,LambdaBarfinal,lambdaindecayframe,lambdabarindecayframe, Py, Pybar,Cij,alpha_xi, alpha_xibar); fntpMCtruth->Column("weight", (Float_t) eventWeight); fntpMCtruth->Column("weightPy", (Float_t) weightPy); fntpMCtruth->Column("weightPyBar", (Float_t) weightPyBar); fntpMCtruth->Column("weightCxx", (Float_t) weightCxx); fntpMCtruth->Column("weightCyy", (Float_t) weightCyy); fntpMCtruth->Column("weightCzz", (Float_t) weightCzz); fntpMCtruth->Column("weightCxz", (Float_t) weightCxz); fntpMCtruth->Column("weightCzx", (Float_t) weightCzx); qa.qaP4Cms("xi_", mcxi->P4(), fntpMCtruth); qa.qaCand("xi_", mcxi, fntpMCtruth); qa.qaP4Cms("xibar_", mcxibar->P4(), fntpMCtruth); qa.qaCand("xibar_", mcxibar, fntpMCtruth); qa.qaP4Cms("xipim_", mcxipim->P4(), fntpMCtruth); qa.qaCand("xipim_", mcxipim, fntpMCtruth); qa.qaP4Cms("xibarpip_", mcxipim->P4(), fntpMCtruth); qa.qaCand("xibarpip_", mcxipim, fntpMCtruth); qa.qaP4Cms("lam0_", mclam->P4(), fntpMCtruth); qa.qaCand("lam0_", mclam, fntpMCtruth); qa.qaP4Cms("lam0bar_", mclambar->P4(), fntpMCtruth); qa.qaCand("lam0bar_", mclambar, fntpMCtruth); qa.qaP4Cms("lam0pim_", mclampim->P4(), fntpMCtruth); qa.qaCand("lam0pim_", mclampim, fntpMCtruth); qa.qaP4Cms("lam0barpip_", mclambarpip->P4(), fntpMCtruth); qa.qaCand("lam0barpip_", mclambarpip, fntpMCtruth); qa.qaP4Cms("p_", mcp->P4(), fntpMCtruth); qa.qaCand("p_", mcp, fntpMCtruth); qa.qaP4Cms("pbar_", mcpbar->P4(), fntpMCtruth); qa.qaCand("pbar_", mcpbar, fntpMCtruth); fntpMCtruth->DumpData(); // *** ---------------------------- *** // *** Now the analysis stuff comes *** // *** ---------------------------- *** // *** ---------- *** // *** Pion- loop *** // *** ---------- *** for (int j=0;jP4(); //Fill Rho tuples fntpPiMinus->Column("Event", (Float_t) nevts); fntpPiMinus->Column("weight", (Float_t) eventWeight); fntpPiMinus->Column("weightCxx", (Float_t) weightCxx); fntpPiMinus->Column("weightCyy", (Float_t) weightCyy); fntpPiMinus->Column("weightCzz", (Float_t) weightCzz); fntpPiMinus->Column("weightCxz", (Float_t) weightCxz); fntpPiMinus->Column("weightCzx", (Float_t) weightCzx); fntpPiMinus->Column("CombiId", (Float_t) j); fntpPiMinus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piminus[j])); qa.qaP4Cms("P4Cms_", piminus4, fntpPiMinus); qa.qaCand("Cand_", piminus[j], fntpPiMinus); fntpPiMinus->DumpData(); } // *** ---------- *** // *** Pion+ loop *** // *** ---------- *** for (int j=0;jP4(); //Construct ROOT TLorentzVectors //Fill Rho tuples fntpPiPlus->Column("Event", (Float_t) nevts); fntpPiPlus->Column("weight", (Float_t) eventWeight); fntpPiPlus->Column("weightCxx", (Float_t) weightCxx); fntpPiPlus->Column("weightCyy", (Float_t) weightCyy); fntpPiPlus->Column("weightCzz", (Float_t) weightCzz); fntpPiPlus->Column("weightCxz", (Float_t) weightCxz); fntpPiPlus->Column("weightCzx", (Float_t) weightCzx); fntpPiPlus->Column("CombiId", (Float_t) j); fntpPiPlus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piplus[j])); qa.qaP4Cms("P4Cms_", piplus4, fntpPiPlus); qa.qaCand("Cand_", piplus[j], fntpPiPlus); fntpPiPlus->DumpData(); } // *** ----------- *** // *** Proton loop *** // *** ----------- *** for (int j=0;jP4(); //Fill Rho tuples fntpProton->Column("Event", (Float_t) nevts); fntpProton->Column("weight", (Float_t) eventWeight); fntpProton->Column("weightCxx", (Float_t) weightCxx); fntpProton->Column("weightCyy", (Float_t) weightCyy); fntpProton->Column("weightCzz", (Float_t) weightCzz); fntpProton->Column("weightCxz", (Float_t) weightCxz); fntpProton->Column("weightCzx", (Float_t) weightCzx); fntpProton->Column("CombiId", (Float_t) j); fntpProton->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(p[j])); qa.qaP4Cms("P4Cms_", p4, fntpProton); qa.qaCand("Cand_", p[j], fntpProton); fntpProton->DumpData(); } // *** --------------- *** // *** AntiProton loop *** // *** --------------- *** for (int j=0;jP4(); //Fill Rho tuples fntpAntiProton->Column("Event", (Float_t) nevts); fntpAntiProton->Column("weight", (Float_t) eventWeight); fntpAntiProton->Column("weightCxx", (Float_t) weightCxx); fntpAntiProton->Column("weightCyy", (Float_t) weightCyy); fntpAntiProton->Column("weightCzz", (Float_t) weightCzz); fntpAntiProton->Column("weightCxz", (Float_t) weightCxz); fntpAntiProton->Column("weightCzx", (Float_t) weightCzx); fntpAntiProton->Column("CombiId", (Float_t) j); fntpAntiProton->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(pbar[j])); qa.qaP4Cms("P4Cms_", pbar4, fntpAntiProton); qa.qaCand("Cand_", pbar[j], fntpAntiProton); fntpAntiProton->DumpData(); } // *** ----------------------------- *** // *** Lambda (-> Proton Pion-) loop *** // *** ----------------------------- *** lam0.Combine(p,piminus); std::map bestVtxFitLam0, bestMassFitLam0; bestVtxFitLam0 = VertexQaIndex(&lam0, 0.001); bestMassFitLam0 = MassFitQaIndex(&lam0, m0_lam, 0.001); //lam0.Select(lambdaMassSelector); //VertexFit(lam0, 0.001); //MassFit(lam0,m0_lam,0.001); lam0.SetType(3122); p.Cleanup(); piminus.Cleanup(); for (int j=0;jP4(); //Fill Rho tuples fntpLambda->Column("Event", (Float_t) nevts); fntpLambda->Column("weight", (Float_t) eventWeight); fntpLambda->Column("weightCxx", (Float_t) weightCxx); fntpLambda->Column("weightCyy", (Float_t) weightCyy); fntpLambda->Column("weightCzz", (Float_t) weightCzz); fntpLambda->Column("weightCxz", (Float_t) weightCxz); fntpLambda->Column("weightCzx", (Float_t) weightCzx); fntpLambda->Column("CombiId", (Float_t) j); fntpLambda->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(lam0[j])); qa.qaP4Cms("P4Cms_", lam04, fntpLambda); qa.qaCand("Cand_", lam0[j], fntpLambda); PndKinVtxFitter vtxfitter(lam0[j]); vtxfitter.Fit(); fntpLambda->Column("VtxFit_HowGood", (Int_t) bestVtxFitLam0[j]); qa.qaFitter("VtxFit_",&vtxfitter,fntpLambda,false); qa.qaVtx("Vtx_",lam0[j]->GetFit(),fntpLambda); qa.qaPull("pull_",lam0[j]->GetFit(),fntpLambda,false); //Select candidates with best vertex fit and mass fit which passed prob cut //if (bestVtxFitLam0[j]==1 && bestMassFitLam0[j]>0){ if (bestVtxFitLam0[j]==1){ CombinedList(lam0[j], &piminus_combined, -211); lam0_best.Append(lam0[j]->GetFit()); //lam0_best.Append(lam0[j]); } fntpLambda->DumpData(); } lam0.Cleanup(); // *** ------------------------------------ *** // *** LambdaBar (-> AntiProton Pion+) loop *** // *** ------------------------------------ *** lam0bar.Combine(pbar,piplus); std::map bestVtxFitLam0bar, bestMassFitLam0bar; bestVtxFitLam0bar = VertexQaIndex(&lam0bar, 0.001); bestMassFitLam0bar = MassFitQaIndex(&lam0bar, m0_lam, 0.001); //lam0bar.Select(lambdaMassSelector); //VertexFit(lam0bar, 0.001); //MassFit(lam0bar,m0_lam,0.001); lam0bar.SetType(-3122); pbar.Cleanup(); piplus.Cleanup(); for (int j=0;jP4(); //Fill Rho tuples fntpLambdaBar->Column("Event", (Float_t) nevts); fntpLambdaBar->Column("weight", (Float_t) eventWeight); fntpLambdaBar->Column("weightCxx", (Float_t) weightCxx); fntpLambdaBar->Column("weightCyy", (Float_t) weightCyy); fntpLambdaBar->Column("weightCzz", (Float_t) weightCzz); fntpLambdaBar->Column("weightCxz", (Float_t) weightCxz); fntpLambdaBar->Column("weightCzx", (Float_t) weightCzx); fntpLambdaBar->Column("CombiId", (Float_t) j); fntpLambdaBar->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(lam0bar[j])); qa.qaP4Cms("P4Cms_", lam0bar4, fntpLambdaBar); qa.qaCand("Cand_", lam0bar[j], fntpLambdaBar); PndKinVtxFitter vtxfitter(lam0bar[j]); vtxfitter.Fit(); fntpLambdaBar->Column("VtxFit_HowGood", (Int_t) bestVtxFitLam0bar[j]); qa.qaFitter("VtxFit_",&vtxfitter,fntpLambdaBar,false); qa.qaVtx("Vtx_",lam0bar[j]->GetFit(),fntpLambdaBar); qa.qaPull("pull_",lam0bar[j]->GetFit(),fntpLambdaBar,false); //Select candidates with best vertex fit and mass fit which passed prob cut //if (bestVtxFitLam0bar[j]==1 && bestMassFitLam0bar[j]>0){ if (bestVtxFitLam0bar[j]==1){ CombinedList(lam0bar[j], &piplus_combined, 211); lam0bar_best.Append(lam0bar[j]->GetFit()); //lam0bar_best.Append(lam0bar[j]); } fntpLambdaBar->DumpData(); } lam0bar.Cleanup(); // *** ------------------------- *** // *** Xi (-> Lambda Pion-) loop *** // *** ------------------------- *** NotCombinedList(piminus_combined, &piminus_noncombined); //xi.Combine(lam0,piminus); xi.Combine(lam0_best,piminus_noncombined); std::map bestVtxFitXi, bestMassFitXi; bestVtxFitXi = VertexQaIndex(&xi, 0.001); bestMassFitXi = MassFitQaIndex(&xi, m0_xi, 0.001); //xi.Select(xiMassSelector); //VertexFit(xi,0.001); //MassFit(xi,m0_xi,0.001); xi.SetType(3312); for (int j=0;jP4(); //Fill Rho tuples fntpXi->Column("Event", (Float_t) nevts); fntpXi->Column("weight", (Float_t) eventWeight); fntpXi->Column("weightCxx", (Float_t) weightCxx); fntpXi->Column("weightCyy", (Float_t) weightCyy); fntpXi->Column("weightCzz", (Float_t) weightCzz); fntpXi->Column("weightCxz", (Float_t) weightCxz); fntpXi->Column("weightCzx", (Float_t) weightCzx); fntpXi->Column("CombiId", (Float_t) j); fntpXi->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(xi[j])); qa.qaP4Cms("P4Cms_", xi04, fntpXi); qa.qaCand("Cand_", xi[j], fntpXi); PndKinVtxFitter vtxfitter(xi[j]); vtxfitter.Fit(); fntpXi->Column("VtxFit_HowGood", (Int_t) bestVtxFitXi[j]); qa.qaFitter("VtxFit_",&vtxfitter,fntpXi,false); qa.qaVtx("Vtx_",xi[j]->GetFit(),fntpXi); qa.qaPull("pull_",xi[j]->GetFit(),fntpXi,false); //Select candidates with best vertex fit and mass fit which passed prob cut //if (bestVtxFitXi[j]==1 && bestMassFitXi[j]>0){ if (bestVtxFitXi[j]==1){ xi_best.Append(xi[j]->GetFit()); //xi_best.Append(xi[j]); } fntpXi->DumpData(); } // *** ------------------------------- *** // *** XiBar (-> LambdaBar Pion+) loop *** // *** ------------------------------- *** NotCombinedList(piplus_combined, &piplus_noncombined); //xibar.Combine(lam0bar,piplus); xibar.Combine(lam0bar_best,piplus_noncombined); std::map bestVtxFitXibar, bestMassFitXibar; bestVtxFitXibar = VertexQaIndex(&xibar, 0.001); bestMassFitXibar = MassFitQaIndex(&xibar, m0_xi, 0.001); //xibar.Select(xiMassSelector); //VertexFit(xibar,0.001); //MassFit(xibar,m0_xi,0.001); xibar.SetType(-3312); for (int j=0;jP4(); //Fill Rho tuples fntpXiBar->Column("Event", (Float_t) nevts); fntpXiBar->Column("weight", (Float_t) eventWeight); fntpXiBar->Column("weightCxx", (Float_t) weightCxx); fntpXiBar->Column("weightCyy", (Float_t) weightCyy); fntpXiBar->Column("weightCzz", (Float_t) weightCzz); fntpXiBar->Column("weightCxz", (Float_t) weightCxz); fntpXiBar->Column("weightCzx", (Float_t) weightCzx); fntpXiBar->Column("CombiId", (Float_t) j); fntpXiBar->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(xibar[j])); qa.qaP4Cms("P4Cms_", xibar04, fntpXiBar); qa.qaCand("Cand_", xibar[j], fntpXiBar); PndKinVtxFitter vtxfitter(xibar[j]); vtxfitter.Fit(); fntpXiBar->Column("VtxFit_HowGood", (Int_t) bestVtxFitXibar[j]); qa.qaFitter("VtxFit_",&vtxfitter,fntpXiBar,false); qa.qaVtx("Vtx_",xibar[j]->GetFit(),fntpXiBar); qa.qaPull("pull_",xibar[j]->GetFit(),fntpXiBar,false); //Select candidates with best vertex fit and mass fit which passed prob cut //if (bestVtxFitXibar[j]==1 && bestMassFitXibar[j]>0){ if (bestVtxFitXi[j]==1){ xibar_best.Append(xibar[j]->GetFit()); //xibar_best.Append(xibar[j]); } fntpXiBar->DumpData(); } // *** ---------------------------- *** // *** Pion- (from Xi->Pi- Lam) loop *** // *** ---------------------------- *** for (int j=0;jP4(); //Fill Rho tuples fntpXi_PiMinus->Column("Event", (Float_t) nevts); fntpXi_PiMinus->Column("weight", (Float_t) eventWeight); fntpXi_PiMinus->Column("weightCxx", (Float_t) weightCxx); fntpXi_PiMinus->Column("weightCyy", (Float_t) weightCyy); fntpXi_PiMinus->Column("weightCzz", (Float_t) weightCzz); fntpXi_PiMinus->Column("weightCxz", (Float_t) weightCxz); fntpXi_PiMinus->Column("weightCzx", (Float_t) weightCzx); fntpXi_PiMinus->Column("CombiId", (Float_t) j); fntpXi_PiMinus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piminus_noncombined[j])); qa.qaP4Cms("P4Cms_", piminus4, fntpXi_PiMinus); qa.qaCand("Cand_", piminus_noncombined[j], fntpXi_PiMinus); fntpXi_PiMinus->DumpData(); } // *** ---------------------------- *** // *** Pion+ (from Xibar->Pi+ Lam0bar) loop *** // *** ---------------------------- *** for (int j=0;jP4(); //Fill Rho tuples fntpXiBar_PiPlus->Column("Event", (Float_t) nevts); fntpXiBar_PiPlus->Column("weight", (Float_t) eventWeight); fntpXiBar_PiPlus->Column("weightCxx", (Float_t) weightCxx); fntpXiBar_PiPlus->Column("weightCyy", (Float_t) weightCyy); fntpXiBar_PiPlus->Column("weightCzz", (Float_t) weightCzz); fntpXiBar_PiPlus->Column("weightCxz", (Float_t) weightCxz); fntpXiBar_PiPlus->Column("weightCzx", (Float_t) weightCzx); fntpXiBar_PiPlus->Column("CombiId", (Float_t) j); fntpXiBar_PiPlus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piplus_noncombined[j])); qa.qaP4Cms("P4Cms_", piplus4, fntpXiBar_PiPlus); qa.qaCand("Cand_fromXi_", piplus_noncombined[j], fntpXiBar_PiPlus); fntpXiBar_PiPlus->DumpData(); } // *** --------------------------------------- *** // *** PbarP system (-> Xi XiBar) loop *** // *** --------------------------------------- *** pbarpsystem.Combine(xi_best,xibar_best); pbarpsystem.SetType(88888); if (pbarpsystem.GetLength() > 1) cout<<"More than one pbarp combi! Total: " <P4(); //Fill Rho tuples fntpPbarp->Column("Event", (Float_t) nevts); fntpPbarp->Column("weight", (Float_t) eventWeight); fntpPbarp->Column("weightPy", (Float_t) weightPy); fntpPbarp->Column("weightPyBar", (Float_t) weightPyBar); fntpPbarp->Column("weightCxx", (Float_t) weightCxx); fntpPbarp->Column("weightCyy", (Float_t) weightCyy); fntpPbarp->Column("weightCzz", (Float_t) weightCzz); fntpPbarp->Column("weightCxz", (Float_t) weightCxz); fntpPbarp->Column("weightCzx", (Float_t) weightCzx); fntpPbarp->Column("CombiId", (Float_t) j); fntpPbarp->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(pbarpsystem[j])); /* * Create covariance matrix */ TMatrixD lamb = pbarpsystem[j]->Daughter(0)->Cov7(); TMatrixD lam = pbarpsystem[j]->Daughter(1)->Cov7(); TMatrixD cov7lam(7,7); TMatrixD cov7lamb(7,7); TMatrixD comp(14,14); TMatrixD v0(14,1); /* * The covariance matrices in RhoCandidates (x,y,z,px,py,pz,E) * are not in the same format as in PndKinFitter (px,py,pz,E,x,y,z) * The following code rearranges the covariance matrix to fitter format. */ for(int l=0; l<7; l++) { for(int m=0; m<7; m++) { if(l>=3) { if(m>=3) { cov7lam[l-3][m-3] = lam[l][m]; } else { cov7lam[l-3][m+4] = lam[l][m]; } } else { if(m>=3) { cov7lam[l+4][m-3] = lam[l][m]; } else { cov7lam[l+4][m+4] = lam[l][m]; } } } } for(int l=0; l<7; l++) { for(int m=0; m<7; m++) { if(l>=3) { if(m>=3) { cov7lamb[l-3][m-3] = lamb[l][m]; } else { cov7lamb[l-3][m+4] = lamb[l][m]; } } else { if(m>=3) { cov7lamb[l+4][m-3] = lamb[l][m]; } else { cov7lamb[l+4][m+4] = lamb[l][m]; } } } } /* * Get starting values for the 4C fit */ v0[0][0] = pbarpsystem[j]->Daughter(0)->Px(); v0[1][0] = pbarpsystem[j]->Daughter(0)->Py(); v0[2][0] = pbarpsystem[j]->Daughter(0)->Pz(); v0[3][0] = pbarpsystem[j]->Daughter(0)->E(); v0[4][0] = pbarpsystem[j]->Daughter(0)->GetPosition().X(); v0[5][0] = pbarpsystem[j]->Daughter(0)->GetPosition().Y(); v0[6][0] = pbarpsystem[j]->Daughter(0)->GetPosition().Z(); v0[7][0] = pbarpsystem[j]->Daughter(1)->Px(); v0[8][0] = pbarpsystem[j]->Daughter(1)->Py(); v0[9][0] = pbarpsystem[j]->Daughter(1)->Pz(); v0[10][0] = pbarpsystem[j]->Daughter(1)->E(); v0[11][0] = pbarpsystem[j]->Daughter(1)->GetPosition().X(); v0[12][0] = pbarpsystem[j]->Daughter(1)->GetPosition().Y(); v0[13][0] = pbarpsystem[j]->Daughter(1)->GetPosition().Z(); /* * Combine the covariance matrices from both lambda and lambdabar into one 14x14 */ comp.SetSub(0,0,cov7lamb); comp.SetSub(7,7,cov7lam); /* * End new code for Xinying part */ PndKinFitter FourMomfitter(pbarpsystem[j]); FourMomfitter.Add4MomConstraint(fIni); FourMomfitter.SetVerbose(true); FourMomfitter.SetInputMatrix(v0, comp); //Use only with Xsong fix FourMomfitter.Fit(); qa.qaFitter("Treefit_",&FourMomfitter,fntpPbarp,false); qa.qaComp("comp_",pbarpsystem[j],fntpPbarp,false,true); TLorentzVector xi_p4 = pbarpsystem[j]->Daughter(0)->P4(); TLorentzVector xibar_p4 = pbarpsystem[j]->Daughter(1)->P4(); TLorentzVector lam0_p4 = pbarpsystem[j]->Daughter(0)->Daughter(0)->P4(); TLorentzVector lam0bar_p4 = pbarpsystem[j]->Daughter(1)->Daughter(0)->P4(); //Boost Lambda to CM frame xi_p4.Boost(-fIni.BoostVector()); //Boost Lambdabar to CM frame xibar_p4.Boost(-fIni.BoostVector()); //Boost proton to Lambda rest frame lam0_p4.Boost(-fIni.BoostVector()); lam0_p4.Boost(-xi_p4.BoostVector()); //Boost pbar to Lambdabar rest frame lam0bar_p4.Boost(-fIni.BoostVector()); lam0bar_p4.Boost(-xibar_p4.BoostVector()); //Lambda decay frame, lambda in rest std::vector xi_decfr = GenerateUnitVectors(beamCM,xi_p4,xibar_p4); //Lambdabar decay frame, lambdabar in rest std::vector xibar_decfr = GenerateUnitVectors(beamCM,xibar_p4,xibar_p4); fntpPbarp->Column("lam0_DecFrame_kx", (Float_t) TMath::Cos(lam0_p4.Vect().Angle(xi_decfr.at(0))), 0.0f ); fntpPbarp->Column("lam0_DecFrame_ky", (Float_t) TMath::Cos(lam0_p4.Vect().Angle(xi_decfr.at(1))), 0.0f ); fntpPbarp->Column("lam0_DecFrame_kz", (Float_t) TMath::Cos(lam0_p4.Vect().Angle(xi_decfr.at(2))), 0.0f ); fntpPbarp->Column("lam0bar_DecFrame_kx", (Float_t) TMath::Cos(lam0bar_p4.Vect().Angle(xibar_decfr.at(0))), 0.0f ); fntpPbarp->Column("lam0bar_DecFrame_ky", (Float_t) TMath::Cos(lam0bar_p4.Vect().Angle(xibar_decfr.at(1))), 0.0f ); fntpPbarp->Column("lam0bar_DecFrame_kz", (Float_t) TMath::Cos(lam0bar_p4.Vect().Angle(xibar_decfr.at(2))), 0.0f ); fntpPbarp->Column("MClam0_DecFrame_kx", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(0))), 0.0f ); fntpPbarp->Column("MClam0_DecFrame_ky", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(1))), 0.0f ); fntpPbarp->Column("MClam0_DecFrame_kz", (Float_t) TMath::Cos(lambdaindecayframe.Vect().Angle(final.at(2))), 0.0f ); fntpPbarp->Column("MClam0bar_DecFrame_kx", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(0))), 0.0f ); fntpPbarp->Column("MClam0bar_DecFrame_ky", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(1))), 0.0f ); fntpPbarp->Column("MClam0bar_DecFrame_kz", (Float_t) TMath::Cos(lambdabarindecayframe.Vect().Angle(LambdaBarfinal.at(2))), 0.0f ); qa.qaP4Cms("xi_", mcxi->P4(), fntpPbarp); qa.qaCand("xi_", mcxi, fntpPbarp); qa.qaP4Cms("xibar_", mcxibar->P4(), fntpPbarp); qa.qaCand("xibar_", mcxibar, fntpPbarp); fntpPbarp->DumpData(); } pbarpsystem.Cleanup(); mclist.Cleanup(); } //End PndXiXibarAnaTask::Exec() void PndXiXibarAnaTask::Finish() { //Write RhoTuple data fntpPiMinus->GetInternalTree()->Write(); fntpPiPlus->GetInternalTree()->Write(); fntpProton->GetInternalTree()->Write(); fntpAntiProton->GetInternalTree()->Write(); fntpLambda->GetInternalTree()->Write(); fntpLambdaBar->GetInternalTree()->Write(); fntpXi_PiMinus->GetInternalTree()->Write(); fntpXiBar_PiPlus->GetInternalTree()->Write(); fntpXi->GetInternalTree()->Write(); fntpXiBar->GetInternalTree()->Write(); fntpPbarp->GetInternalTree()->Write(); fntpMCtruth->GetInternalTree()->Write(); } //End PndXiXibarAnaTask::Finish() ClassImp(PndXiXibarAnaTask)