/* * PndLLbarAnaTask.cxx * * Created on: July 7, 2015 * Author: Walter Ikegami Andersson */ // The header file #include "PndLLbarAnaTask.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 namespace std; using std::cout; using std::endl; // ----- Default constructor ------------------------------------------- PndLLbarAnaTask::PndLLbarAnaTask() : FairTask("PandaLLbarAnaTask"), fBeamMom(1.642), fSignalsample(true){ } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndLLbarAnaTask::~PndLLbarAnaTask() { } // ------------------------------------------------------------------------- // ----- Method to select PID candidates int PndLLbarAnaTask::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 PndLLbarAnaTask::VertexFit(RhoCandList &l){ for (int i = 0; i < l.GetLength(); 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() < 0.01){ l.Remove(l[i]); } } } } void PndLLbarAnaTask::MassFit(RhoCandList &l, Double_t mass){ for (int i = 0; i < l.GetLength(); 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() < 0.01){ l.Remove(l[i]); } } } } void PndLLbarAnaTask::Tree4CFit(RhoCandList &l, TLorentzVector ini, Double_t prob_cut){ for (int i = l.GetLength()-1; i >= 0; --i){ PndKinFitter tree4cfitter(l.Get(i)); 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 PndLLbarAnaTask::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 PndLLbarAnaTask::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 PndLLbarAnaTask::Pull(TString pre, RhoTuple *t, RhoCandidate *c, RhoCandidate *fit){ if(fit && c) { TLorentzVector difp4 = c->P4() - fit->P4(); TVector3 difpos = c->Pos() - fit->Pos(); TMatrixD cov7 = c->Cov7(); TMatrixD cov7fit = fit->Cov7(); t->Column(pre+"fitpullx", (Float_t) (difpos.X()/sqrt(cov7(0,0) - cov7fit(0,0))),0.0f ); t->Column(pre+"fitpully", (Float_t) (difpos.Y()/sqrt(cov7(1,1) - cov7fit(1,1))),0.0f ); t->Column(pre+"fitpullz", (Float_t) (difpos.Z()/sqrt(cov7(2,2) - cov7fit(2,2))),0.0f ); t->Column(pre+"fitpullpx", (Float_t) (difp4.Px()/sqrt(cov7(3,3) - cov7fit(3,3))),0.0f ); t->Column(pre+"fitpullpy", (Float_t) (difp4.Py()/sqrt(cov7(4,4) - cov7fit(4,4))),0.0f ); t->Column(pre+"fitpullpz", (Float_t) (difp4.Pz()/sqrt(cov7(5,5) - cov7fit(5,5))),0.0f ); t->Column(pre+"fitpulle", (Float_t) (difp4.E()/sqrt(cov7(6,6) - cov7fit(6,6))),0.0f ); } else{ t->Column(pre+"pullx", (Float_t) -999., 0.0f ); t->Column(pre+"pully", (Float_t) -999., 0.0f ); t->Column(pre+"pullz", (Float_t) -999., 0.0f ); t->Column(pre+"pullpx", (Float_t) -999., 0.0f ); t->Column(pre+"pullpy", (Float_t) -999., 0.0f ); t->Column(pre+"pullpz", (Float_t) -999., 0.0f ); t->Column(pre+"pulle", (Float_t) -999., 0.0f ); } } int PndLLbarAnaTask::tagHits(RhoCandidate *c){ /**@brief Tag the particle with different integers * @details Tag the particle with different integers: * 0: if there is no hit in the detector * 1: sttHits>3 or mvdHits>3 or gemHit>3 */ int tag = 0; PndPidCandidate * pidCand = (PndPidCandidate*)c->GetRecoCandidate(); int branch = trackBranch(c); if(pidCand){ int mvdHits = pidCand->GetMvdHits(); int sttHits = pidCand->GetSttHits(); int gemHits = pidCand->GetGemHits(); if(mvdHits>3 || gemHits>3 || sttHits>3) tag=1; //else if (sttHits>3 && branch==FairRootManager::Instance()->GetBranchId("SttMvdGemGenTrack")) tag=1; else tag=0; } return tag; } int PndLLbarAnaTask::trackBranch(RhoCandidate *c){ int branch=0; PndPidCandidate * pid = (PndPidCandidate*)c->GetRecoCandidate(); if(pid){ branch = pid->GetTrackBranch(); } return branch; } // ----- Public method Init -------------------------------------------- InitStatus PndLLbarAnaTask::Init() { // initialize analysis object theAnalysis = new PndAnalysis(); // reset the event counter nevts = 0; SumWeightCosThetaX = 0.; SumWeightCosThetaXBar = 0.; SumWeightCosThetaY = 0.; SumWeightCosThetaYBar = 0.; SumWeightCosThetaZ = 0.; SumWeightCosThetaZBar = 0.; SumWeightThtX_ThtX = 0.; SumWeightThtX_ThtY = 0.; SumWeightThtX_ThtZ = 0.; SumWeightThtY_ThtX = 0.; SumWeightThtY_ThtY = 0.; SumWeightThtY_ThtZ = 0.; SumWeightThtZ_ThtX = 0.; SumWeightThtZ_ThtY = 0.; SumWeightThtZ_ThtZ = 0.; SumWeight = 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(); //Lambda decay parameters 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); // 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"); fntpPbarp = new RhoTuple("ntpPbarp", "Pbarp system info"); fntpPbarp4Cfit = new RhoTuple("ntpPbarp4Cfit", "Pbarp system info 4Cfit"); // MC tuples fntpMCtruth = new RhoTuple("ntpMCtruth", "MC truth info"); cout<<"PndLLbarAnaTask - 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(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"); // *** --------------------------- *** // *** Same all Monte Carlo Tracks *** // *** --------------------------- *** RhoCandidate *mcpbarp, *mclambar, *mclam, *mcpip, *mcpim, *mcp, *mcpbar, *dummyCand; /* * Save true MC tracks * If background sample is analysed, there is no monte carlo truth * and a dummy is used instead */ dummyCand = new RhoCandidate(); if (fSignalsample) { mcpbarp = mclist[0]; // pbarp system mclambar = mclist[1]; //d0 Lambar mcpbar = mclist[1]->Daughter(0);//d0d0 Pbar mcpip = mclist[1]->Daughter(1); //d0d1 pi+ mclam = mclist[2]; //d1 Lam mcp = mclist[2]->Daughter(0); //d1d0 Proton mcpim = mclist[2]->Daughter(1); //d1d1 pi- } else { mcpbarp = dummyCand; mclambar = dummyCand; mcpbar = dummyCand; mcpip = dummyCand; mclam = dummyCand; mcp = dummyCand; mcpim = dummyCand; } /* * If a dpm background sample is analysed, llbar events should be removed * This is done here */ bool lammc = false; bool lambmc = false; if (!fSignalsample) { for (int dd = 0; dd < mclist.GetLength(); dd++ ){ if (mclist[dd]->PdgCode() == 3122) { lammc = true; } if (mclist[dd]->PdgCode() == -3122) { lambmc = true; } } if (lammc && lambmc) { return; } } /* * If background sample is analysed, there is no monte carlo truth */ // *** ----------------------------------- *** // *** 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 Lambda to CM frame TLorentzVector MC_P4lam_CM = mclam->P4(); MC_P4lam_CM.Boost(-fIni.BoostVector()); //Boost Lambdabar to CM frame TLorentzVector MC_P4lambar_CM = mclambar->P4(); MC_P4lambar_CM.Boost(-fIni.BoostVector()); //Boost proton to Lambda rest frame TLorentzVector MC_P4p_lamCM = mcp->P4(); MC_P4p_lamCM.Boost(-fIni.BoostVector()); MC_P4p_lamCM.Boost(-MC_P4lam_CM.BoostVector()); //Boost pbar to Lambdabar rest frame TLorentzVector MC_P4pbar_lambarCM = mcpbar->P4(); MC_P4pbar_lambarCM.Boost(-fIni.BoostVector()); MC_P4pbar_lambarCM.Boost(-MC_P4lambar_CM.BoostVector()); // *** -------------------------- *** // *** Construct Reference frames *** // *** -------------------------- *** //Initial frame, lab system std::vector lab_RF; lab_RF.push_back(TVector3(1,0,0)); lab_RF.push_back(TVector3(0,1,0)); lab_RF.push_back(TVector3(0,0,1)); //Lambdabar decay frame, lambdabar in rest std::vector MC_lambar_RF = GenerateUnitVectors(beamCM,MC_P4lambar_CM,MC_P4lambar_CM); //Lambda decay frame, lambda in rest std::vector MC_lam_RF = GenerateUnitVectors(beamCM,MC_P4lam_CM,MC_P4lambar_CM); TLorentzVector MC_P4pbar_rotlambarRF = TransformCoords(lab_RF,MC_lambar_RF,MC_P4pbar_lambarCM); fntpMCtruth->Column("MC_d0d0_phi", (Float_t) MC_P4pbar_rotlambarRF.Phi(), 0.0f ); fntpMCtruth->Column("MC_d0d0_CosTheta", (Float_t) TMath::Cos(MC_P4pbar_rotlambarRF.Theta()), 0.0f ); fntpMCtruth->Column("MC_d0d0_kx", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(0))), 0.0f ); fntpMCtruth->Column("MC_d0d0_ky", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(1))), 0.0f ); fntpMCtruth->Column("MC_d0d0_kz", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(2))), 0.0f ); TLorentzVector MC_P4p_rotlamRF = TransformCoords(lab_RF,MC_lam_RF,MC_P4p_lamCM); fntpMCtruth->Column("MC_d1d0_phi", (Float_t) MC_P4p_rotlamRF.Phi(), 0.0f ); fntpMCtruth->Column("MC_d1d0_CosTheta", (Float_t) TMath::Cos(MC_P4p_rotlamRF.Theta()), 0.0f ); fntpMCtruth->Column("MC_d1d0_kx", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(0))), 0.0f ); fntpMCtruth->Column("MC_d1d0_ky", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(1))), 0.0f ); fntpMCtruth->Column("MC_d1d0_kz", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(2))), 0.0f ); Double_t cmAngle = 999.; TMatrixD Cij(3,3); Double_t Py = 1.; Double_t PyBar = 1.; Double_t Py2sin = 1.; Double_t PyBar2sin = 1.; Double_t weightPy = 1.; Double_t weightPyBar = 1.; Double_t weightPy2sin = 1.; Double_t weightPyBar2sin = 1.; Double_t weightCxx = 1.; Double_t weightCyy = 1.; Double_t weightCzz = 1.; Double_t weightCxz = 1.; Double_t weightCzx = 1.; Double_t weightCxxData = 1.; Double_t weightCyyData = 1.; Double_t weightCzzData = 1.; Double_t weightCxzData = 1.; Double_t weightCzxData = 1.; Double_t eventWeight = 1.; /* * If signal sample is analysed, weights should also be calculated for each event * If background is analysed, all weights should be 1 */ if(fSignalsample){ /* * Calculate CM angle, weights, spin variables with method of moments */ cmAngle = beamCM.Vect().Angle(MC_P4lambar_CM.Vect()); Cij = GenerateSpinCorrTrigFuncs(cmAngle); //TMatrixD Cij = GenerateSpinCorrData(cmAngle); Py = TMath::Sin(cmAngle); Py2sin = TMath::Sin(2.*cmAngle); //Double_t Py = 0.; PyBar = TMath::Sin(cmAngle); PyBar2sin = TMath::Sin(2.*cmAngle); //Double_t PyBar = 0.; weightPy = 1 + alpha_lam*Py*TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); weightPyBar = 1 + alpha_lambar*PyBar*TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())); weightPy2sin = 1 + alpha_lam*Py2sin*TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); weightPyBar2sin = 1 + alpha_lambar*PyBar2sin*TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())); weightCxx = 1 + alpha_lam*alpha_lambar*Cij(0,0) * TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())) * TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())); weightCyy = 1 + alpha_lam*alpha_lambar*Cij(1,1) * TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())) * TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())); weightCzz = 1 + alpha_lam*alpha_lambar*Cij(2,2) * TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())) * TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())); weightCxz = 1 + alpha_lam*alpha_lambar*Cij(0,2) * TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())); weightCzx = 1 + alpha_lam*alpha_lambar*Cij(2,0) * TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())); weightCxxData = GenerateSpinCorrDataCxx(cmAngle); weightCyyData = GenerateSpinCorrDataCyy(cmAngle); weightCzzData = GenerateSpinCorrDataCzz(cmAngle); weightCxzData = GenerateSpinCorrDataCxz(cmAngle); weightCzxData = GenerateSpinCorrDataCxz(cmAngle); //The article assumed Cxz == Czx eventWeight = SpinWeightGeneral(MC_lam_RF,MC_lambar_RF,MC_P4p_lamCM,MC_P4pbar_lambarCM, Py, PyBar,Cij,alpha_lam, alpha_lambar); SumWeight += eventWeight; //Polarisation along x SumWeightCosThetaX += eventWeight*TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())); SumWeightCosThetaXBar += eventWeight*TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())); //Polarisation along y SumWeightCosThetaY += eventWeight*TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); SumWeightCosThetaYBar += eventWeight*TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())); //Polarisation along z SumWeightCosThetaZ += eventWeight*TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())); SumWeightCosThetaZBar += eventWeight*TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())); //Cij SumWeightThtX_ThtX += eventWeight * TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())); SumWeightThtX_ThtY += eventWeight * TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); SumWeightThtX_ThtZ += eventWeight * TMath::Cos(MC_lambar_RF.at(0).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())); SumWeightThtY_ThtX += eventWeight * TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())); SumWeightThtY_ThtY += eventWeight * TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); SumWeightThtY_ThtZ += eventWeight * TMath::Cos(MC_lambar_RF.at(1).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())); SumWeightThtZ_ThtX += eventWeight * TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(0).Angle(MC_P4p_lamCM.Vect())); SumWeightThtZ_ThtY += eventWeight * TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(1).Angle(MC_P4p_lamCM.Vect())); SumWeightThtZ_ThtZ += eventWeight * TMath::Cos(MC_lambar_RF.at(2).Angle(MC_P4pbar_lambarCM.Vect())) * TMath::Cos(MC_lam_RF.at(2).Angle(MC_P4p_lamCM.Vect())); } fntpMCtruth->Column("weightPy", (Float_t) weightPy); fntpMCtruth->Column("weightPyBar", (Float_t) weightPyBar); fntpMCtruth->Column("weightPy2sin", (Float_t) weightPy2sin); fntpMCtruth->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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); fntpMCtruth->Column("weightCxxData", (Float_t) weightCxxData); fntpMCtruth->Column("weightCyyData", (Float_t) weightCyyData); fntpMCtruth->Column("weightCzzData", (Float_t) weightCzzData); fntpMCtruth->Column("weightCxzData", (Float_t) weightCxzData); fntpMCtruth->Column("weightCzxData", (Float_t) weightCzxData); qa.qaCand("MC_d0", mclambar, fntpMCtruth); qa.qaP4Cms("MC_d0", mclambar->P4(), fntpMCtruth); qa.qaCand("MC_d0d0", mcpbar, fntpMCtruth); qa.qaP4Cms("MC_d0d0", mcpbar->P4(), fntpMCtruth); qa.qaCand("MC_d0d1", mcpip, fntpMCtruth); qa.qaP4Cms("MC_d0d1", mcpip->P4(), fntpMCtruth); qa.qaCand("MC_d1", mclam, fntpMCtruth); qa.qaP4Cms("MC_d1", mclam->P4(), fntpMCtruth); qa.qaCand("MC_d1d0", mcp, fntpMCtruth); qa.qaP4Cms("MC_d1d0", mcp->P4(), fntpMCtruth); qa.qaCand("MC_d1d1", mcpim, fntpMCtruth); qa.qaP4Cms("MC_d1d1", mcpim->P4(), fntpMCtruth); fntpMCtruth->DumpData(); // *** ---------------------------- *** // *** Now the analysis stuff comes *** // *** ---------------------------- *** // *** ---------- *** // *** Pion- loop *** // *** ---------- *** for (j=0;jP4(); //Fill Rho tuples fntpPiMinus->Column("Event", (Float_t) nevts); fntpPiMinus->Column("weightPy", (Float_t) weightPy); fntpPiMinus->Column("weightPyBar", (Float_t) weightPyBar); fntpPiMinus->Column("weightPy2sin", (Float_t) weightPy2sin); fntpPiMinus->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpPiMinus->Column("weightCyyData", (Float_t) weightCyyData); fntpPiMinus->Column("weightCzzData", (Float_t) weightCzzData); fntpPiMinus->Column("weightCxzData", (Float_t) weightCxzData); fntpPiMinus->Column("weightCzxData", (Float_t) weightCzxData); fntpPiMinus->Column("CombiId", (Float_t) j); fntpPiMinus->Column("MinHitReq", (Int_t) tagHits(piminus[j])); fntpPiMinus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piminus[j])); if(fSignalsample){ RhoCandidate* mctpbarp = piminus[j]->GetMcTruth()->TheMother()->TheMother(); RhoCandidate* mctlam = piminus[j]->GetMcTruth()->TheMother(); if (mctpbarp && mctlam && piminus[j]->GetMcTruth()->TheMother()->PdgCode() == 3122) { fntpPiMinus->Column("llbarEvent", (Int_t) 1); qa.qaCand("MC_", piminus[j]->GetMcTruth(), fntpPiMinus,false); } else { fntpPiMinus->Column("llbarEvent", (Int_t) 0); qa.qaCand("MC_", dummyCand, fntpPiMinus,false); } } //PndPidCandidate* pidpar=(PndPidCandidate*)piminus[j]->GetRecoCandidate(); FairRecoCandidate* recpar=piminus[j]->GetRecoCandidate(); fntpPiMinus->Column("TrackBranch", (Int_t) recpar->GetTrackBranch()); qa.qaCand("", piminus[j], fntpPiMinus); qa.qaMcDiff("",piminus[j],fntpPiMinus); fntpPiMinus->DumpData(); } // *** ---------- *** // *** Pion+ loop *** // *** ---------- *** for (j=0;jP4(); //Fill Rho tuples fntpPiPlus->Column("Event", (Float_t) nevts); fntpPiPlus->Column("weightPy", (Float_t) weightPy); fntpPiPlus->Column("weightPyBar", (Float_t) weightPyBar); fntpPiPlus->Column("weightPy2sin", (Float_t) weightPy2sin); fntpPiPlus->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpPiPlus->Column("weightCyyData", (Float_t) weightCyyData); fntpPiPlus->Column("weightCzzData", (Float_t) weightCzzData); fntpPiPlus->Column("weightCxzData", (Float_t) weightCxzData); fntpPiPlus->Column("weightCzxData", (Float_t) weightCzxData); fntpPiPlus->Column("CombiId", (Float_t) j); fntpPiPlus->Column("MinHitReq", (Int_t) tagHits(piplus[j])); fntpPiPlus->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(piplus[j])); if(fSignalsample){ RhoCandidate* mctpbarp = piplus[j]->GetMcTruth()->TheMother()->TheMother(); RhoCandidate* mctlam = piplus[j]->GetMcTruth()->TheMother(); if (mctpbarp && mctlam && piplus[j]->GetMcTruth()->TheMother()->PdgCode() == -3122) { fntpPiPlus->Column("llbarEvent", (Int_t) 1); qa.qaCand("MC_", piplus[j]->GetMcTruth(), fntpPiPlus,false); } else { fntpPiPlus->Column("llbarEvent", (Int_t) 0); qa.qaCand("MC_", dummyCand, fntpPiPlus,false); } } //PndPidCandidate* pidpar=(PndPidCandidate*)piminus[j]->GetRecoCandidate(); FairRecoCandidate* recpar=piplus[j]->GetRecoCandidate(); fntpPiPlus->Column("TrackBranch", (Int_t) recpar->GetTrackBranch()); qa.qaCand("", piplus[j], fntpPiPlus); qa.qaMcDiff("",piplus[j],fntpPiPlus); fntpPiPlus->DumpData(); } // *** ----------- *** // *** Proton loop *** // *** ----------- *** for (j=0;jP4(); //Fill Rho tuples fntpProton->Column("Event", (Float_t) nevts); fntpProton->Column("weightPy", (Float_t) weightPy); fntpProton->Column("weightPyBar", (Float_t) weightPyBar); fntpProton->Column("weightPy2sin", (Float_t) weightPy2sin); fntpProton->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpProton->Column("weightCyyData", (Float_t) weightCyyData); fntpProton->Column("weightCzzData", (Float_t) weightCzzData); fntpProton->Column("weightCxzData", (Float_t) weightCxzData); fntpProton->Column("weightCzxData", (Float_t) weightCzxData); fntpProton->Column("CombiId", (Float_t) j); fntpProton->Column("MinHitReq", (Int_t) tagHits(p[j])); fntpProton->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(p[j])); if(fSignalsample){ RhoCandidate* mctpbarp = p[j]->GetMcTruth()->TheMother()->TheMother(); RhoCandidate* mctlam = p[j]->GetMcTruth()->TheMother(); if (mctpbarp && mctlam && p[j]->GetMcTruth()->TheMother()->PdgCode() == 3122) { fntpProton->Column("llbarEvent", (Int_t) 1); qa.qaCand("MC_", p[j]->GetMcTruth(), fntpProton,false); } else { fntpProton->Column("llbarEvent", (Int_t) 0); qa.qaCand("MC_", dummyCand, fntpProton,false); } } //PndPidCandidate* pidpar=(PndPidCandidate*)piminus[j]->GetRecoCandidate(); FairRecoCandidate* recpar=p[j]->GetRecoCandidate(); fntpProton->Column("TrackBranch", (Int_t) recpar->GetTrackBranch()); qa.qaCand("", p[j], fntpProton); qa.qaMcDiff("",p[j],fntpProton); fntpProton->DumpData(); } // *** --------------- *** // *** AntiProton loop *** // *** --------------- *** for (j=0;jP4(); //Fill Rho tuples fntpAntiProton->Column("Event", (Float_t) nevts); fntpAntiProton->Column("weightPy", (Float_t) weightPy); fntpAntiProton->Column("weightPyBar", (Float_t) weightPyBar); fntpAntiProton->Column("weightPy2sin", (Float_t) weightPy2sin); fntpAntiProton->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpAntiProton->Column("weightCyyData", (Float_t) weightCyyData); fntpAntiProton->Column("weightCzzData", (Float_t) weightCzzData); fntpAntiProton->Column("weightCxzData", (Float_t) weightCxzData); fntpAntiProton->Column("weightCzxData", (Float_t) weightCzxData); fntpAntiProton->Column("CombiId", (Float_t) j); fntpAntiProton->Column("MinHitReq", (Int_t) tagHits(pbar[j])); fntpAntiProton->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(pbar[j])); if(fSignalsample){ RhoCandidate* mctpbarp = pbar[j]->GetMcTruth()->TheMother()->TheMother(); RhoCandidate* mctlam = pbar[j]->GetMcTruth()->TheMother(); if (mctpbarp && mctlam && pbar[j]->GetMcTruth()->TheMother()->PdgCode() == -3122) { fntpAntiProton->Column("llbarEvent", (Int_t) 1); qa.qaCand("MC_", pbar[j]->GetMcTruth(), fntpAntiProton,false); } else { fntpAntiProton->Column("llbarEvent", (Int_t) 0); qa.qaCand("MC_", dummyCand, fntpAntiProton,false); } } //PndPidCandidate* pidpar=(PndPidCandidate*)piminus[j]->GetRecoCandidate(); FairRecoCandidate* recpar=pbar[j]->GetRecoCandidate(); fntpAntiProton->Column("TrackBranch", (Int_t) recpar->GetTrackBranch()); qa.qaCand("", pbar[j], fntpAntiProton); qa.qaMcDiff("",pbar[j],fntpAntiProton); fntpAntiProton->DumpData(); } // *** ----------------------------- *** // *** Lambda (-> Proton Pion-) loop *** // *** ----------------------------- *** lam0.Combine(p,piminus); lam0.Select(lambdaMassSelector); std::map bestVtxFitLam0, bestMassFitLam0; bestVtxFitLam0 = VertexQaIndex(&lam0, 0.01); bestMassFitLam0 = MassFitQaIndex(&lam0, m0_lam, 0.01); //VertexFit(lam0, 0.01); //MassFit(lam0,m0_lam,0.01); lam0.SetType(3122); piminus.Cleanup(); p.Cleanup(); for (j=0;jP4(); //Fill Rho tuples fntpLambda->Column("Event", (Float_t) nevts); fntpLambda->Column("weightPy", (Float_t) weightPy); fntpLambda->Column("weightPyBar", (Float_t) weightPyBar); fntpLambda->Column("weightPy2sin", (Float_t) weightPy2sin); fntpLambda->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpLambda->Column("weightCyyData", (Float_t) weightCyyData); fntpLambda->Column("weightCzzData", (Float_t) weightCzzData); fntpLambda->Column("weightCxzData", (Float_t) weightCxzData); fntpLambda->Column("weightCzxData", (Float_t) weightCzxData); fntpLambda->Column("CombiId", (Float_t) j); fntpLambda->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(lam0[j])); 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("VtxFit_",lam0[j]->GetFit(),fntpLambda); qa.qaPull("VtxFit_pull_",lam0[j]->GetFit(),fntpLambda,false); qa.qaComp("VtxFit_Comp_",lam0[j]->GetFit(),fntpLambda,false,false); qa.qaMcDiff("d0_",lam0[j]->Daughter(0),fntpLambda); qa.qaMcDiff("d1_",lam0[j]->Daughter(1),fntpLambda); qa.qaMcDiff("d0_fit",lam0[j]->Daughter(0)->GetFit(),fntpLambda); qa.qaMcDiff("d1_fit",lam0[j]->Daughter(1)->GetFit(),fntpLambda); /* RhoCandidate *d = lam0[j]->GetFit(); PndKinFitter massfitter(d); //PndKinFitter massfitter(lam0[j]->GetFit()); massfitter.AddMassConstraint(m0_lam); massfitter.Fit(); fntpLambda->Column("MassFit_HowGood", (Int_t) bestMassFitLam0[j]); qa.qaFitter("MassFit_",&massfitter,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){ lam0_best.Append(lam0[j]->GetFit()); //lam0_best.Append(d->GetFit()); lam0_raw.Append(lam0[j]); } fntpLambda->DumpData(); } lam0.Cleanup(); // *** ------------------------------------ *** // *** LambdaBar (-> AntiProton Pion+) loop *** // *** ------------------------------------ *** lam0bar.Combine(pbar,piplus); lam0bar.Select(lambdaMassSelector); std::map bestVtxFitLam0bar, bestMassFitLam0bar; bestVtxFitLam0bar = VertexQaIndex(&lam0bar, 0.01); bestMassFitLam0bar = MassFitQaIndex(&lam0bar, m0_lam, 0.01); //VertexFit(lam0bar, 0.01); //MassFit(lam0bar,m0_lam,0.01); lam0bar.SetType(-3122); pbar.Cleanup(); piplus.Cleanup(); for (j=0;jP4(); //Fill Rho tuples fntpLambdaBar->Column("Event", (Float_t) nevts); fntpLambdaBar->Column("weightPy", (Float_t) weightPy); fntpLambdaBar->Column("weightPyBar", (Float_t) weightPyBar); fntpLambdaBar->Column("weightPy2sin", (Float_t) weightPy2sin); fntpLambdaBar->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpLambdaBar->Column("weightCyyData", (Float_t) weightCyyData); fntpLambdaBar->Column("weightCzzData", (Float_t) weightCzzData); fntpLambdaBar->Column("weightCxzData", (Float_t) weightCxzData); fntpLambdaBar->Column("weightCzxData", (Float_t) weightCzxData); fntpLambdaBar->Column("CombiId", (Float_t) j); fntpLambdaBar->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(lam0bar[j])); 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("VtxFit_",lam0bar[j]->GetFit(),fntpLambdaBar); qa.qaPull("VtxFit_pull_",lam0bar[j]->GetFit(),fntpLambdaBar,false); qa.qaComp("VtxFit_Comp_",lam0bar[j]->GetFit(),fntpLambdaBar,false,false); qa.qaMcDiff("d0_",lam0bar[j]->Daughter(0),fntpLambdaBar); qa.qaMcDiff("d1_",lam0bar[j]->Daughter(1),fntpLambdaBar); qa.qaMcDiff("d0_fit",lam0bar[j]->Daughter(0)->GetFit(),fntpLambdaBar); qa.qaMcDiff("d1_fit",lam0bar[j]->Daughter(1)->GetFit(),fntpLambdaBar); /* RhoCandidate *d = lam0bar[j]->GetFit(); PndKinFitter massfitter(d); //PndKinFitter massfitter(lam0bar[j]->GetFit()); massfitter.AddMassConstraint(m0_lam); massfitter.Fit(); fntpLambdaBar->Column("MassFit_HowGood", (Int_t) bestMassFitLam0bar[j]); qa.qaFitter("MassFit_",&massfitter,fntpLambdaBar,false); */ //if (bestVtxFitLam0bar[j]==1 && bestMassFitLam0bar[j]>0){ if (bestVtxFitLam0bar[j]==1){ lam0bar_best.Append(lam0bar[j]->GetFit()); //lam0bar_best.Append(d->GetFit()); lam0bar_raw.Append(lam0bar[j]); } fntpLambdaBar->DumpData(); } lam0bar.Cleanup(); // *** --------------------------------------- *** // *** PbarP system (-> Lambda LambdaBar) loop *** // *** --------------------------------------- *** pbarpsystem.Combine(lam0bar_best,lam0_best); //pbarpsystem.Combine(lam0bar_raw,lam0_raw); pbarpsystem.SetType(88888); pbarp_raw.Combine(lam0bar_raw,lam0_raw); pbarp_raw.SetType(88888); lam0_best.Cleanup(); lam0bar_best.Cleanup(); if (pbarpsystem.GetLength() > 1) cout<<"WARNING: More than one pbarp combi! Total: " <Column("Event", (Float_t) nevts); fntpPbarp->Column("weightPy", (Float_t) weightPy); fntpPbarp->Column("weightPyBar", (Float_t) weightPyBar); fntpPbarp->Column("weightPy2sin", (Float_t) weightPy2sin); fntpPbarp->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); 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("weightCxxData", (Float_t) weightCxxData); fntpPbarp->Column("weightCyyData", (Float_t) weightCyyData); fntpPbarp->Column("weightCzzData", (Float_t) weightCzzData); fntpPbarp->Column("weightCxzData", (Float_t) weightCxzData); fntpPbarp->Column("weightCzxData", (Float_t) weightCzxData); fntpPbarp->Column("CombiId", (Float_t) j); fntpPbarp->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(pbarpsystem[j])); //Fill Rho tuples fntpPbarp4Cfit->Column("Event", (Float_t) nevts); fntpPbarp4Cfit->Column("weightPy", (Float_t) weightPy); fntpPbarp4Cfit->Column("weightPyBar", (Float_t) weightPyBar); fntpPbarp4Cfit->Column("weightPy2sin", (Float_t) weightPy2sin); fntpPbarp4Cfit->Column("weightPyBar2sin", (Float_t) weightPyBar2sin); fntpPbarp4Cfit->Column("weightCxx", (Float_t) weightCxx); fntpPbarp4Cfit->Column("weightCyy", (Float_t) weightCyy); fntpPbarp4Cfit->Column("weightCzz", (Float_t) weightCzz); fntpPbarp4Cfit->Column("weightCxz", (Float_t) weightCxz); fntpPbarp4Cfit->Column("weightCzx", (Float_t) weightCzx); fntpPbarp4Cfit->Column("weightCxxData", (Float_t) weightCxxData); fntpPbarp4Cfit->Column("weightCyyData", (Float_t) weightCyyData); fntpPbarp4Cfit->Column("weightCzzData", (Float_t) weightCzzData); fntpPbarp4Cfit->Column("weightCxzData", (Float_t) weightCxzData); fntpPbarp4Cfit->Column("weightCzxData", (Float_t) weightCzxData); fntpPbarp4Cfit->Column("CombiId", (Float_t) j); fntpPbarp4Cfit->Column("McTruthMatch", (Bool_t) theAnalysis->McTruthMatch(pbarpsystem[j])); TLorentzVector lam0bar_p4 = pbarpsystem[j]->Daughter(0)->P4(); TLorentzVector lam0_p4 = pbarpsystem[j]->Daughter(1)->P4(); TLorentzVector P4pbar_lambarCM = pbarpsystem[j]->Daughter(0)->Daughter(0)->P4(); TLorentzVector piP4p_lamCM = pbarpsystem[j]->Daughter(0)->Daughter(1)->P4(); TLorentzVector P4p_lamCM = pbarpsystem[j]->Daughter(1)->Daughter(0)->P4(); TLorentzVector pim_p4 = pbarpsystem[j]->Daughter(1)->Daughter(1)->P4(); //d0 Lambar //d0d0 Pbar //d0d1 pi+ //d1 Lam //d1d0 Proton //d1d1 pi- //Boost Lambda to CM frame lam0_p4.Boost(-fIni.BoostVector()); //Boost Lambdabar to CM frame lam0bar_p4.Boost(-fIni.BoostVector()); //Boost proton to Lambda rest frame P4p_lamCM.Boost(-fIni.BoostVector()); P4p_lamCM.Boost(-lam0_p4.BoostVector()); //Boost pbar to Lambdabar rest frame P4pbar_lambarCM.Boost(-fIni.BoostVector()); P4pbar_lambarCM.Boost(-lam0bar_p4.BoostVector()); //Lambda decay frame, lambda in rest std::vector lam_RF = GenerateUnitVectors(beamCM,lam0_p4,lam0bar_p4); //Lambdabar decay frame, lambdabar in rest std::vector lambar_RF = GenerateUnitVectors(beamCM,lam0bar_p4,lam0bar_p4); fntpPbarp->Column("d0d0_kx", (Float_t) TMath::Cos(P4pbar_lambarCM.Vect().Angle(lambar_RF.at(0))), 0.0f ); fntpPbarp->Column("d0d0_ky", (Float_t) TMath::Cos(P4pbar_lambarCM.Vect().Angle(lambar_RF.at(1))), 0.0f ); fntpPbarp->Column("d0d0_kz", (Float_t) TMath::Cos(P4pbar_lambarCM.Vect().Angle(lambar_RF.at(2))), 0.0f ); fntpPbarp->Column("d1d0_kx", (Float_t) TMath::Cos(P4p_lamCM.Vect().Angle(lam_RF.at(0))), 0.0f ); fntpPbarp->Column("d1d0_ky", (Float_t) TMath::Cos(P4p_lamCM.Vect().Angle(lam_RF.at(1))), 0.0f ); fntpPbarp->Column("d1d0_kz", (Float_t) TMath::Cos(P4p_lamCM.Vect().Angle(lam_RF.at(2))), 0.0f ); /* * Fill Monte Carlo cosine angles */ fntpPbarp->Column("MC_d0d0_kx", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(0))), 0.0f ); fntpPbarp->Column("MC_d0d0_ky", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(1))), 0.0f ); fntpPbarp->Column("MC_d0d0_kz", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(2))), 0.0f ); fntpPbarp->Column("MC_d1d0_kx", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(0))), 0.0f ); fntpPbarp->Column("MC_d1d0_ky", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(1))), 0.0f ); fntpPbarp->Column("MC_d1d0_kz", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(2))), 0.0f ); qa.qaComp("",pbarpsystem[j],fntpPbarp,false,true); qa.qaMcDiff("",pbarpsystem[j],fntpPbarp); qa.qaMcDiff("d0_",pbarpsystem[j]->Daughter(0),fntpPbarp); qa.qaMcDiff("d1_",pbarpsystem[j]->Daughter(1),fntpPbarp); qa.qaMcDiff("d0d0_",pbarpsystem[j]->Daughter(0)->Daughter(0),fntpPbarp); qa.qaMcDiff("d0d1_",pbarpsystem[j]->Daughter(0)->Daughter(1),fntpPbarp); qa.qaMcDiff("d1d0_",pbarpsystem[j]->Daughter(1)->Daughter(0),fntpPbarp); qa.qaMcDiff("d1d1_",pbarpsystem[j]->Daughter(1)->Daughter(1),fntpPbarp); qa.qaCand("MC_d0", mclambar, fntpPbarp); qa.qaP4Cms("MC_d0", mclambar->P4(), fntpPbarp); qa.qaCand("MC_d1", mclam, fntpPbarp); qa.qaP4Cms("MC_d1", mclam->P4(), fntpPbarp); /* * Adding new code to test fix to RhoFitter by Xsong */ PndKinVtxFitter vtxfitter(pbarpsystem[j]); // instantiate a vertex fitter vtxfitconv = false; //vtxfitconv = vtxfitter.Fit(); //TMatrixD Tal0_t=vtxfitter.GetAl0(); //TMatrixD TV_al0_t=vtxfitter.GetVal0(); //RhoCandidate *pbarpvtxfit = pbarpsystem[j]->GetFit(); fntpPbarp->Column("VtxFit_Convergence", (Bool_t) vtxfitconv); fntpPbarp4Cfit->Column("VtxFit_Convergence", (Bool_t) vtxfitconv); qa.qaFitter("VtxFit_",&vtxfitter,fntpPbarp,false); qa.qaFitter("VtxFit_",&vtxfitter,fntpPbarp4Cfit,false); /* * 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]); //with vtx fit FourMomfitter.Add4MomConstraint(fIni); FourMomfitter.SetVerbose(true); //FourMomfitter.SetInputMatrix(Tal0_t, TV_al0_t); //Use only with Xsong fix FourMomfitter.SetInputMatrix(v0, comp); //Use only with Xsong fix Bool_t is4fit = FourMomfitter.Fit(); RhoCandidate * pbarp4cfit; pbarp4cfit = pbarpsystem[j]->GetFit(); //with vtx fit /* * Fill ntuples with angles from 4C fit */ fntpPbarp->Column("FourMomFit_Convergence", (Bool_t) is4fit); fntpPbarp4Cfit->Column("FourMomFit_Convergence", (Bool_t) is4fit); qa.qaFitter("FourMomFit_",&FourMomfitter,fntpPbarp,false); qa.qaFitter("FourMomFit_",&FourMomfitter,fntpPbarp4Cfit,false); TLorentzVector fit_lam0bar_p4 = pbarp4cfit->Daughter(0)->P4(); TLorentzVector fit_lam0_p4 = pbarp4cfit->Daughter(1)->P4(); TLorentzVector fit_P4pbar_lambarCM = pbarp4cfit->Daughter(0)->Daughter(0)->P4(); TLorentzVector fit_piP4p_lamCM = pbarp4cfit->Daughter(0)->Daughter(1)->P4(); TLorentzVector fit_P4p_lamCM = pbarp4cfit->Daughter(1)->Daughter(0)->P4(); TLorentzVector fit_pim_p4 = pbarp4cfit->Daughter(1)->Daughter(1)->P4(); //Boost Lambda to CM frame fit_lam0_p4.Boost(-fIni.BoostVector()); //Boost Lambdabar to CM frame fit_lam0bar_p4.Boost(-fIni.BoostVector()); //Boost proton to Lambda rest frame fit_P4p_lamCM.Boost(-fIni.BoostVector()); fit_P4p_lamCM.Boost(-lam0_p4.BoostVector()); //Boost pbar to Lambdabar rest frame fit_P4pbar_lambarCM.Boost(-fIni.BoostVector()); fit_P4pbar_lambarCM.Boost(-lam0bar_p4.BoostVector()); //Lambda decay frame, lambda in rest std::vector fit_lam_RF = GenerateUnitVectors(beamCM,fit_lam0_p4,fit_lam0bar_p4); //Lambdabar decay frame, lambdabar in rest std::vector fit_lambar_RF = GenerateUnitVectors(beamCM,fit_lam0bar_p4,fit_lam0bar_p4); fntpPbarp4Cfit->Column("d0d0_kx", (Float_t) TMath::Cos(fit_P4pbar_lambarCM.Vect().Angle(fit_lambar_RF.at(0))), 0.0f ); fntpPbarp4Cfit->Column("d0d0_ky", (Float_t) TMath::Cos(fit_P4pbar_lambarCM.Vect().Angle(fit_lambar_RF.at(1))), 0.0f ); fntpPbarp4Cfit->Column("d0d0_kz", (Float_t) TMath::Cos(fit_P4pbar_lambarCM.Vect().Angle(fit_lambar_RF.at(2))), 0.0f ); fntpPbarp4Cfit->Column("d1d0_kx", (Float_t) TMath::Cos(fit_P4p_lamCM.Vect().Angle(fit_lam_RF.at(0))), 0.0f ); fntpPbarp4Cfit->Column("d1d0_ky", (Float_t) TMath::Cos(fit_P4p_lamCM.Vect().Angle(fit_lam_RF.at(1))), 0.0f ); fntpPbarp4Cfit->Column("d1d0_kz", (Float_t) TMath::Cos(fit_P4p_lamCM.Vect().Angle(fit_lam_RF.at(2))), 0.0f ); /* * Fill Monte Carlo cosine angles */ fntpPbarp4Cfit->Column("MC_d0d0_kx", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(0))), 0.0f ); fntpPbarp4Cfit->Column("MC_d0d0_ky", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(1))), 0.0f ); fntpPbarp4Cfit->Column("MC_d0d0_kz", (Float_t) TMath::Cos(MC_P4pbar_lambarCM.Vect().Angle(MC_lambar_RF.at(2))), 0.0f ); fntpPbarp4Cfit->Column("MC_d1d0_kx", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(0))), 0.0f ); fntpPbarp4Cfit->Column("MC_d1d0_ky", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(1))), 0.0f ); fntpPbarp4Cfit->Column("MC_d1d0_kz", (Float_t) TMath::Cos(MC_P4p_lamCM.Vect().Angle(MC_lam_RF.at(2))), 0.0f ); qa.qaComp("",pbarp4cfit,fntpPbarp4Cfit,false,false); qa.qaMcDiff("",pbarp4cfit,fntpPbarp4Cfit); qa.qaMcDiff("d0_",pbarp4cfit->Daughter(0),fntpPbarp4Cfit); qa.qaMcDiff("d1_",pbarp4cfit->Daughter(1),fntpPbarp4Cfit); qa.qaMcDiff("d0d0_",pbarp4cfit->Daughter(0)->Daughter(0),fntpPbarp4Cfit); qa.qaMcDiff("d0d1_",pbarp4cfit->Daughter(0)->Daughter(1),fntpPbarp4Cfit); qa.qaMcDiff("d1d0_",pbarp4cfit->Daughter(1)->Daughter(0),fntpPbarp4Cfit); qa.qaMcDiff("d1d1_",pbarp4cfit->Daughter(1)->Daughter(1),fntpPbarp4Cfit); qa.qaCand("MC_d0", mclambar, fntpPbarp4Cfit); qa.qaP4Cms("MC_d0", mclambar->P4(), fntpPbarp4Cfit); qa.qaCand("MC_d1", mclam, fntpPbarp4Cfit); qa.qaP4Cms("MC_d1", mclam->P4(), fntpPbarp4Cfit); Pull("d0_",fntpPbarp4Cfit,pbarpsystem[j]->Daughter(0),pbarp4cfit->Daughter(0)); Pull("d1_",fntpPbarp4Cfit,pbarpsystem[j]->Daughter(1),pbarp4cfit->Daughter(1)); TMatrixD d0d0cov = pbarp4cfit->Daughter(0)->Daughter(0)->Cov7(); Float_t fds = pbarp4cfit->Daughter(0)->Daughter(0)->GetPosition().Z(); fntpPbarp4Cfit->Column("FDS", (Float_t) fds, 0.0f ); /* * Write data into root files */ fntpPbarp->DumpData(); fntpPbarp4Cfit->DumpData(); } pbarpsystem.Cleanup(); mclist.Cleanup(); } //End PndLLbarAnaTask::Exec() void PndLLbarAnaTask::Finish() { Double_t Px_rec = (3/alpha_lam)*(SumWeightCosThetaX/SumWeight); Double_t Pxbar_rec = (3/alpha_lambar)*(SumWeightCosThetaXBar/SumWeight); Double_t Py_rec = (3/alpha_lam)*(SumWeightCosThetaY/SumWeight); Double_t PyBar_rec = (3/alpha_lambar)*(SumWeightCosThetaYBar/SumWeight); Double_t Pz_rec = (3/alpha_lam)*(SumWeightCosThetaZ/SumWeight); Double_t Pzbar_rec = (3/alpha_lambar)*(SumWeightCosThetaZBar/SumWeight); Double_t Cxx = (9/(alpha_lam*alpha_lambar))*(SumWeightThtX_ThtX/SumWeight); Double_t Cxy = (9/(alpha_lam*alpha_lambar))*(SumWeightThtX_ThtY/SumWeight); Double_t Cxz = (9/(alpha_lam*alpha_lambar))*(SumWeightThtX_ThtZ/SumWeight); Double_t Cyx = (9/(alpha_lam*alpha_lambar))*(SumWeightThtY_ThtX/SumWeight); Double_t Cyy = (9/(alpha_lam*alpha_lambar))*(SumWeightThtY_ThtY/SumWeight); Double_t Cyz = (9/(alpha_lam*alpha_lambar))*(SumWeightThtY_ThtZ/SumWeight); Double_t Czx = (9/(alpha_lam*alpha_lambar))*(SumWeightThtZ_ThtX/SumWeight); Double_t Czy = (9/(alpha_lam*alpha_lambar))*(SumWeightThtZ_ThtY/SumWeight); Double_t Czz = (9/(alpha_lam*alpha_lambar))*(SumWeightThtZ_ThtZ/SumWeight); cout<<"MC Polarization for p along x: "<GetInternalTree()->Write(); fntpPiPlus->GetInternalTree()->Write(); fntpProton->GetInternalTree()->Write(); fntpAntiProton->GetInternalTree()->Write(); fntpLambda->GetInternalTree()->Write(); fntpLambdaBar->GetInternalTree()->Write(); fntpPbarp->GetInternalTree()->Write(); fntpPbarp4Cfit->GetInternalTree()->Write(); fntpMCtruth->GetInternalTree()->Write(); } //End PndLLbarAnaTask::Finish() ClassImp(PndLLbarAnaTask)