// Analyze data coming from the detector simulation // for the decay: // p+ p- (0) ⇒ aXi+ (1) + Xi-* (2) // | |→ K- (8) + Lambda0 (7) // | |→ pi- (10) + p (9) // |→ aLambda0 (3) + pi+ (4) // |→ pi+ (6) + ap (5) // class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; class RhoTuple; #include #include // define the used particle IDs for more readability enum pidNumbers { kXims = 3314, kaXip = -3312, kL0 = 3122, kaL0 = -3122, kKm = -321, kPip = 211, kPim = -211, kPp = 2212, kaPm = -2212 }; namespace andre { void qaFitter(PndRhoTupleQA * qa, TString pre, RhoCandidate * c, RhoTuple * n, RhoFitterBase * fit, bool skip = false) { if (0 == n) return; if (!skip) { qa->qaComp(pre, c, n); n->Column(pre + "Chi2", (float) fit->GetChi2(), 0.0f); n->Column(pre + "Ndf", (float) fit->GetNdf(), 0.0f); n->Column(pre + "Prob", (float) fit->GetProb(), 0.0f); } } /** * @brief Save hit information for particle candidate * @details Saves the amount of hits in the tracking detectors for the given particle candidate, if applicable (composites often don't leave any hits) * * @param pre A prefix for the leaf name * @param c Particle candidate * @param n RhoTuple the data will be stored into * @param checkDaughters Loop through all daughter particles * @param skip Skip this round */ void qaDetHits(TString pre, RhoCandidate * c, RhoTuple * n, bool checkDaughters = true, bool skip = false) { if (0 == n) return; PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate(); if ( !skip && mic ) { n->Column(pre + "hitsMVD", (Int_t) mic->GetMvdHits(), -999); n->Column(pre + "hitsSTT", (Int_t) mic->GetSttHits(), -999); n->Column(pre + "hitsGEM", (Int_t) mic->GetGemHits(), -999); // n->Column(pre + "hitsFTS", (Int_t) mic->GetFtsHits(), -999); } else { n->Column(pre + "hitsMVD", -999, -999); n->Column(pre + "hitsSTT", -999, -999); n->Column(pre + "hitsGEM", -999, -999); // n->Column(pre + "hitsFTS", -999, -999); } // loop through all daughters int nd = c->NDaughters(); for ( int d = 0; d < nd && checkDaughters; d++ ) { RhoCandidate *dau = c->Daughter(d); TString pre2 = TString::Format("%sd%d", pre.Data(), d); qaDetHits(pre2, dau, n); } } /** * @brief Save information from all particle candidates in a list. * @details Saves the production vertex and momentum information of all particles in the candidate list given. * * @param qa Pointer to the PndRhoTupleQA object * @param pre A prefix for the leaf name * @param cc Particle candidate list * @param n RhoTuple the data will be stored into * @param skip Skip this round */ void qaCandList(PndRhoTupleQA * qa, TString pre, RhoCandList &cc, RhoTuple * n, Int_t evt, bool skip = false) { if (0 == n) return; RhoCandidate * curCand, * truth; std::vector foundUids; for ( int i = 0; i < cc.GetLength(); i++ ) { curCand = cc[i]; truth = curCand->GetMcTruth(); bool mct = (0x0 != truth); bool alreadyFound = false; // check for multiplicities of tracks if (mct) { for ( int j = 0; j < foundUids.size(); j ++ ) { if ( foundUids[j] == truth->Uid() ) { alreadyFound = true; break; } } if ( !alreadyFound ) foundUids.push_back(truth->Uid()); } n->Column(pre + "evt", (Int_t) evt, -999); n->Column(pre + "mct", (Bool_t) mct, false); qa->qaCand(pre, curCand, n, skip); qa->qaCand(pre + "mc", truth, n, (!mct || skip)); n->Column(pre + "uid", (Int_t) truth->Uid(), -999); n->Column(pre + "mcmPdg", (Int_t) truth->TheMother()->PdgCode(), -999999); n->Column(pre + "mcm2Pdg", (Int_t) truth->TheMother()->TheMother()->PdgCode(), -999999); n->Column(pre + "mult", (Bool_t) alreadyFound, false); andre::qaDetHits(pre, curCand, n, false, (!mct || skip)); n->DumpData(); } } } void cascade_ana(TString prefix = "", int nEvents = 0) { bool doPreAnalysis = true; // *** the output file for FairRunAna TString OutFile = prefix + "ana_complete.root"; // *** the files coming from the simulation TString inPidFile = prefix + "pid_complete.root"; // this file contains the PndPidCandidates and McTruth TString inParFile = prefix + "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); // load Geane for back propagation FairGeane *Geane = new FairGeane(); fRun->AddTask(Geane); fRun->SetOutputFile(OutFile); fRun->SetUseFairLinks(kTRUE); fRun->Init(); // *** create an output file for all histograms TString anaOutput = prefix + TString::Format("cascade_ana-%d.root", nEvents); TFile * out = TFile::Open(anaOutput, "RECREATE"); // *** create ntuples for J/psi and psi(2S) RhoTuple * ntpMc = new RhoTuple("ntpMC", "Mctruth info"); if (doPreAnalysis) { RhoTuple * ntpPre = new RhoTuple("ntpPre", "Pre Analysis"); RhoTuple * ntpPreCand = new RhoTuple("ntpPreCand", "Pre Analysis of Particle Candidates"); } RhoTuple * ntpLambda = new RhoTuple("ntpLambda", "Lambda0 Analysis"); RhoTuple * ntpAntiLambda = new RhoTuple("ntpAntiLambda", "AntiLambda0 Analysis"); RhoTuple * ntpAntiXi = new RhoTuple("ntpAntiXi", "AntiXiplus Analysis"); RhoTuple * ntpXi = new RhoTuple("ntpXi", "Ximinus Analysis"); // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nEvents == 0 || nEvents > theAnalysis->GetEntries()) nEvents = theAnalysis->GetEntries(); // *** RhoCandLists for the analysis RhoCandList all, mclist; RhoCandList kminus, piplus, piminus, proton, antiproton; RhoCandList Lambda0List, aLambda0List, XiminusList, aXiplusList; RhoCandList Lambda0FitList, aLambda0FitList, piplus_reduced; // *** Mass selector for the Lambda cands double m0_Lambda0 = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass(); std::cout << "mass Lambda0 = " << m0_Lambda0 << std::endl; RhoMassParticleSelector * lambda0MassSelector = new RhoMassParticleSelector("lambda0", m0_Lambda0, 0.1); // 0.1 is the full width of the window! // *** Mass selector for the Xi-* cands double m0_Ximinus = 1.69; std::cout << "mass Xi-* = " << m0_Ximinus << std::endl; RhoMassParticleSelector * XiminusMassSelector = new RhoMassParticleSelector("xiplus", m0_Ximinus, 0.3); // 0.3 is the full width of the window! // *** Mass selector for the Xi+ cands double m0_aXiplus = TDatabasePDG::Instance()->GetParticle("Xi-_bar")->Mass(); std::cout << "mass Xi+ = " << m0_aXiplus << std::endl; RhoMassParticleSelector * aXiplusMassSelector = new RhoMassParticleSelector("xiplus", m0_aXiplus, 0.3); // 0.3 is the full width of the window! // *** Pid Selection Algorithms TString pidSelection = "PidAlgoIdealCharged"; // Use MC info to set P=1 for ISPARTICLE and P=0 for ISNOTPARTICLE // *** the lorentz vector of the initial reaction double pbarmom = 4.1; 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); TVector3 beamBoost = ini.BoostVector(); PndRhoTupleQA qa(theAnalysis, pbarmom); RhoCandidate * curCand; RhoCandidate * dummyCand = new RhoCandidate(); PndVtxPoca * vtxPocaFitter = new PndVtxPoca(); TVector3 vtxPoca; // my propagation object PndRhoPropDistToIP myProp(0); // *** // -------------- the event loop -------------- // *** Int_t evt = -1; while (theAnalysis->GetEvent() && ++evt < nEvents) { if ((evt%100) == 0) cout << "evt " << evt << endl; // *** get MC list and store info in ntuple theAnalysis->FillList(mclist, "McTruth"); qa.qaMcList("", mclist, ntpMc); ntpMc->DumpData(); // *** Setup event shape object theAnalysis->FillList(all, "All", pidSelection); PndEventShape evsh(all, ini, 0.05, 0.1); // *** Select with no PID info ('All'); type and mass are set // Best = take only best probablity for Kaon theAnalysis->FillList(kminus, "KaonBestMinus", pidSelection); theAnalysis->FillList(piplus, "PionBestPlus", pidSelection); theAnalysis->FillList(piminus, "PionBestMinus", pidSelection); theAnalysis->FillList(proton, "ProtonBestPlus", pidSelection); theAnalysis->FillList(antiproton, "ProtonBestMinus", pidSelection); // ***** // *** combinatorics for intermediate particles // ***** // final state RhoCandList piplus_reduced(piplus); piplus_reduced.SetType(kPip); // 1st order: Lambda -> Pi- proton Lambda0List.Combine(piminus, proton); Lambda0List.Select(lambda0MassSelector); Lambda0List.SetType(kL0); int nLambda0Mc = theAnalysis->McTruthMatch(Lambda0List); // 1st order: aLambda -> Pi+ antiproton aLambda0List.Combine(piplus, antiproton); aLambda0List.Select(lambda0MassSelector); aLambda0List.SetType(kaL0); int naLambda0Mc = theAnalysis->McTruthMatch(aLambda0List); // 2nd order: Xi- -> Lambda K- Lambda0FitList.Cleanup(); Lambda0FitList.SetType(kL0); // 2nd order: Xi+ -> aLambda Pi+ aLambda0FitList.Cleanup(); aLambda0FitList.SetType(kaL0); // *** // -------------- Lambda0 analysis -------------- // *** // get a list of fit probability, i.e. how good was the fit PndRhoFitProbPresorter Lambda_probPresorter(&Lambda0List, m0_Lambda0); Lambda_probPresorter.GenerateIndices(); for ( int l0 = 0; l0 < Lambda0List.GetLength(); l0++ ) { curCand = Lambda0List[l0]; // some general info about event (actually written for each candidate!) ntpLambda->Column("ev", (Float_t) evt); ntpLambda->Column("cand", (Float_t) l0); ntpLambda->Column("evI", (Int_t) evt); ntpLambda->Column("candI", (Int_t) l0); ntpLambda->Column("ncandL0", (Float_t) Lambda0List.GetLength()); ntpLambda->Column("ncandL0I", (Int_t) Lambda0List.GetLength()); ntpLambda->Column("nMcL0", (Float_t) nLambda0Mc); ntpLambda->Column("nMcL0I", (Int_t) nLambda0Mc); // store info about initial 4-vector qa.qaP4("beam_", ini, ntpLambda); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("L0_", curCand, ntpLambda); andre::qaDetHits("L0_", curCand, ntpLambda); // store info about event shapes qa.qaEventShapeShort("es_", &evsh, ntpLambda); // point of closest approach for the daughters qa.qaPoca("fPoca_", curCand, ntpLambda); // get the vertex fitter result PndKinVtxFitter * vertexFitter = Lambda_probPresorter.GetVtxFit(l0); RhoCandidate * lambda0Fit = curCand->GetFit(); // store info of vertex fitting ntpLambda->Column("fVtx_HowGood", (Int_t) Lambda_probPresorter.GetVtxBestFitIndex(l0)); ntpLambda->Column("fVtx_ProbIdx", (Int_t) Lambda_probPresorter.GetVtxProbIndex(l0)); andre::qaFitter(&qa, "fVtx_", lambda0Fit, ntpLambda, vertexFitter); ntpLambda->Column("fVtx_ProbRatio", (Float_t) Lambda_probPresorter.GetVtxRatioToBestProb(l0) ); ntpLambda->Column("fVtx_ProbRatio2nd", (Float_t) Lambda_probPresorter.GetVtxRatio2ndToBestProb() ); // get the mass fitter result PndKinFitter * massFitter = Lambda_probPresorter.GetMassFit(l0); RhoCandidate * lambda0Fit_mass = lambda0Fit->GetFit(); andre::qaFitter(&qa, "fMass_", lambda0Fit_mass, ntpLambda, massFitter); ntpLambda->Column("fMass_HowGood", (Int_t) Lambda_probPresorter.GetMassBestFitIndex(l0)); ntpLambda->Column("fMass_ProbIdx", (Int_t) Lambda_probPresorter.GetMassProbIndex(l0)); ntpLambda->Column("fMass_ProbRatio", (Float_t) Lambda_probPresorter.GetMassRatioToBestProb(l0) ); ntpLambda->Column("fMass_ProbRatio2nd", (Float_t) Lambda_probPresorter.GetMassRatio2ndToBestProb() ); // Differences to MC truth qa.qaMcDiff("L0_", curCand, ntpLambda); qa.qaMcDiff("fVtx_", lambda0Fit, ntpLambda); qa.qaMcDiff("fMass_", lambda0Fit_mass, ntpLambda); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate * truth = curCand->GetMcTruth(); TLorentzVector lv(-999., -999., -999., -999.); if (0x0 != truth) { lv = truth->P4(); qa.qaVtx("truth_", truth, ntpLambda); } else { qa.qaVtx("truth_", dummyCand, ntpLambda); } qa.qaP4("truth_", lv, ntpLambda); // apply a condition for combined lambdas if ( Lambda_probPresorter.GetVtxBestFitIndex(l0) == 1 && Lambda_probPresorter.GetMassBestFitIndex(l0) == 1 ) { // append the current candidate to the fitted list Lambda0FitList.Add(lambda0Fit_mass); } // // store information of boosted particles // lambda0Fit->Boost(-beamBoost); // qa.qaComp("boost_", lambda0Fit, ntpLambda); ntpLambda->DumpData(); } // loop over Lambda0 // *** // ---------- AntiLambda0 analysis ---------- // *** // get a list of fit probability, i.e. how good was the fit PndRhoFitProbPresorter aLambda_probPresorter(&aLambda0List, m0_Lambda0); aLambda_probPresorter.GenerateIndices(); for ( int al0 = 0; al0 < aLambda0List.GetLength(); al0++ ) { curCand = aLambda0List[al0]; // some general info about event (actually written for each candidate!) ntpAntiLambda->Column("ev", (Float_t) evt); ntpAntiLambda->Column("cand", (Float_t) al0); ntpAntiLambda->Column("evI", (Int_t) evt); ntpAntiLambda->Column("candI", (Int_t) al0); ntpAntiLambda->Column("ncandaL0", (Float_t) aLambda0List.GetLength()); ntpAntiLambda->Column("ncandaL0I", (Int_t) aLambda0List.GetLength()); ntpAntiLambda->Column("nMcaL0", (Float_t) naLambda0Mc); ntpAntiLambda->Column("nMcaL0I", (Int_t) naLambda0Mc); // store info about initial 4-vector qa.qaP4("beam_", ini, ntpAntiLambda); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("aL0_", curCand, ntpAntiLambda); andre::qaDetHits("aL0_", curCand, ntpAntiLambda); // store info about event shapes qa.qaEventShapeShort("es_", &evsh, ntpAntiLambda); // point of closest approach for the daughters qa.qaPoca("fPoca_", curCand, ntpAntiLambda); // get the vertex fitter result PndKinVtxFitter * vertexFitter = aLambda_probPresorter.GetVtxFit(al0); RhoCandidate * aLambda0Fit = curCand->GetFit(); // store info of vertex fitting ntpAntiLambda->Column("fVtx_HowGood", (Int_t) aLambda_probPresorter.GetVtxBestFitIndex(al0)); ntpAntiLambda->Column("fVtx_ProbIdx", (Int_t) aLambda_probPresorter.GetVtxProbIndex(al0)); andre::qaFitter(&qa, "fVtx_", aLambda0Fit, ntpAntiLambda, vertexFitter); ntpAntiLambda->Column("fVtx_ProbRatio", (Float_t) aLambda_probPresorter.GetVtxRatioToBestProb(al0) ); ntpAntiLambda->Column("fVtx_ProbRatio2nd", (Float_t) aLambda_probPresorter.GetVtxRatio2ndToBestProb() ); // get the mass fitter result PndKinFitter * massFitter = aLambda_probPresorter.GetMassFit(al0); RhoCandidate * aLambda0Fit_mass = aLambda0Fit->GetFit(); andre::qaFitter(&qa, "fMass_", aLambda0Fit_mass, ntpAntiLambda, massFitter); ntpAntiLambda->Column("fMass_HowGood", (Int_t) aLambda_probPresorter.GetMassBestFitIndex(al0)); ntpAntiLambda->Column("fMass_ProbIdx", (Int_t) aLambda_probPresorter.GetMassProbIndex(al0)); ntpAntiLambda->Column("fMass_ProbRatio", (Float_t) aLambda_probPresorter.GetMassRatioToBestProb(al0) ); ntpAntiLambda->Column("fMass_ProbRatio2nd", (Float_t) aLambda_probPresorter.GetMassRatio2ndToBestProb() ); // Differences to MC truth qa.qaMcDiff("aL0_", curCand, ntpAntiLambda); qa.qaMcDiff("fVtx_", aLambda0Fit, ntpAntiLambda); qa.qaMcDiff("fMass_", aLambda0Fit_mass, ntpAntiLambda); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate * truth = curCand->GetMcTruth(); TLorentzVector lv(-999., -999., -999., -999.); if (0x0 != truth) { lv = truth->P4(); qa.qaVtx("truth_", truth, ntpAntiLambda); } else { qa.qaVtx("truth_", dummyCand, ntpAntiLambda); } qa.qaP4("truth_", lv, ntpAntiLambda); // apply a condition for combined anti lambdas if ( aLambda_probPresorter.GetVtxBestFitIndex(al0) == 1 && aLambda_probPresorter.GetMassBestFitIndex(al0) == 1 ) { // remove pions belonging to the aLambda from the reduced list piplus_reduced.Remove(curCand->Daughter(0)); // append the current candidate to the fitted list aLambda0FitList.Add(aLambda0Fit_mass); } // // store information of boosted particles // aLambda0Fit->Boost(-beamBoost); // qa.qaComp("boost_", aLambda0Fit, ntpAntiLambda); ntpAntiLambda->DumpData(); } // loop over AntiLambda0 // *** // -------------- AntiXiplus analysis -------------- // *** // combine 2nd order: Xi+ -> aLambda(after mass fit) & Pi+ aXiplusList.Combine(aLambda0FitList, piplus_reduced); aXiplusList.Select(aXiplusMassSelector); aXiplusList.SetType(kaXip); int naXiplusMc = theAnalysis->McTruthMatch(aXiplusList); // get a list of fit probability, i.e. how good was the fit PndRhoFitProbPresorter aXiplus_probPresorter(&aXiplusList, m0_aXiplus); aXiplus_probPresorter.GenerateIndices(); for ( int i_aXi = 0; i_aXi < aXiplusList.GetLength(); i_aXi++ ) { curCand = aXiplusList[i_aXi]; // some general info about event (actually written for each candidate!) ntpAntiXi->Column("ev", (Float_t) evt); ntpAntiXi->Column("cand", (Float_t) i_aXi); ntpAntiXi->Column("evI", (Int_t) evt); ntpAntiXi->Column("candI", (Int_t) i_aXi); ntpAntiXi->Column("ncandaXip", (Float_t) aXiplusList.GetLength()); ntpAntiXi->Column("ncandaXipI", (Int_t) aXiplusList.GetLength()); ntpAntiXi->Column("nMcaXip", (Float_t) naXiplusMc); ntpAntiXi->Column("nMcaXipI", (Int_t) naXiplusMc); // store info about initial 4-vector qa.qaP4("beam_", ini, ntpAntiXi); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("aXip_", curCand, ntpAntiXi); andre::qaDetHits("aXip_", curCand, ntpAntiXi); // store info about event shapes qa.qaEventShapeShort("es_", &evsh, ntpAntiXi); // point of closest approach for the daughters qa.qaPoca("fPoca_", curCand, ntpAntiXi); // get the vertex fitter result PndKinVtxFitter * vertexFitter = aXiplus_probPresorter.GetVtxFit(i_aXi); RhoCandidate * aXiplusFit = curCand->GetFit(); // store info of vertex fitting ntpAntiXi->Column("fVtx_HowGood", (Int_t) aXiplus_probPresorter.GetVtxBestFitIndex(i_aXi)); ntpAntiXi->Column("fVtx_ProbIdx", (Int_t) aXiplus_probPresorter.GetVtxProbIndex(i_aXi)); andre::qaFitter(&qa, "fVtx_", aXiplusFit, ntpAntiXi, vertexFitter); ntpAntiXi->Column("fVtx_ProbRatio", (Float_t) aXiplus_probPresorter.GetVtxRatioToBestProb(i_aXi) ); ntpAntiXi->Column("fVtx_d0ProbRatio", (Float_t) aLambda_probPresorter.GetVtxRatioToBestProbUID(curCand->Daughter(0)->Uid()) ); ntpAntiXi->Column("fVtx_d0ProbRatio2nd", (Float_t) aLambda_probPresorter.GetMassRatio2ndToBestProb() ); ntpAntiXi->Column("fVtx_d0ProbIdx", (Int_t) aLambda_probPresorter.GetVtxProbIndexUID(curCand->Daughter(0)->Uid()) ); // propagate the candidate back to z=0 plane and get distance to IP vtxPocaFitter->GetPocaVtx(vtxPoca, curCand); if ( myProp.SetParticleTrackPar(curCand, vtxPoca) ) { myProp.PropagateWithGeane(); // .. with Geane to plane z=0 myProp.PropagateWithGeanePocaIP(); // .. with Geane until closest point to IP myProp.PropagateWithHelix(); // .. with helix to plane z=0 } // save the propagate results myProp.qaPropResult("Prop_", ntpAntiXi); // get the mass fitter result PndKinFitter * massFitter = aXiplus_probPresorter.GetMassFit(i_aXi); RhoCandidate * aXiplusFit_mass = aXiplusFit->GetFit(); andre::qaFitter(&qa, "fMass_", aXiplusFit_mass, ntpAntiXi, massFitter); ntpAntiXi->Column("fMass_HowGood", (Int_t) aXiplus_probPresorter.GetMassBestFitIndex(i_aXi)); ntpAntiXi->Column("fMass_ProbIdx", (Int_t) aXiplus_probPresorter.GetMassProbIndex(i_aXi)); ntpAntiXi->Column("fMass_ProbRatio", (Float_t) aXiplus_probPresorter.GetMassRatioToBestProb(i_aXi) ); ntpAntiXi->Column("fMass_d0ProbRatio2nd", (Float_t) aLambda_probPresorter.GetMassRatio2ndToBestProb() ); ntpAntiXi->Column("fMass_d0ProbRatio", (Float_t) aLambda_probPresorter.GetMassRatioToBestProbUID(curCand->Daughter(0)->Uid()) ); ntpAntiXi->Column("fMass_d0ProbIdx", (Int_t) aLambda_probPresorter.GetMassProbIndexUID(curCand->Daughter(0)->Uid()) ); // Differences to MC truth qa.qaMcDiff("aXip_", curCand, ntpAntiXi); qa.qaMcDiff("fVtx_", aXiplusFit, ntpAntiXi); qa.qaMcDiff("fMass_", aXiplusFit_mass, ntpAntiXi); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate * truth = curCand->GetMcTruth(); TLorentzVector lv(-999., -999., -999., -999.); if (0x0 != truth) { lv = truth->P4(); qa.qaVtx("truth_", truth, ntpAntiXi); } else { qa.qaVtx("truth_", dummyCand, ntpAntiXi); } qa.qaP4("truth_", lv, ntpAntiXi); // // store information of boosted particles // aXiplusFit->Boost(-beamBoost); // qa.qaComp("boost_", aXiplusFit, ntpAntiXi); ntpAntiXi->DumpData(); } // loop over AntiXiplus // *** // -------------- Xiplus analysis -------------- // *** // combine 2nd order: Xi- -> Lambda(after mass fit) & K- XiminusList.Combine(Lambda0FitList, kminus); XiminusList.Select(XiminusMassSelector); XiminusList.SetType(kXims); int nXiminusMc = theAnalysis->McTruthMatch(XiminusList); // get a list of fit probability, i.e. how good was the fit PndRhoFitProbPresorter Ximinus_probPresorter(&XiminusList, m0_Ximinus); Ximinus_probPresorter.GenerateIndices(); for ( int i_Xi = 0; i_Xi < XiminusList.GetLength(); i_Xi++ ) { curCand = XiminusList[i_Xi]; // some general info about event (actually written for each candidate!) ntpXi->Column("ev", (Float_t) evt); ntpXi->Column("cand", (Float_t) i_Xi); ntpXi->Column("evI", (Int_t) evt); ntpXi->Column("candI", (Int_t) i_Xi); ntpXi->Column("ncandXim", (Float_t) XiminusList.GetLength()); ntpXi->Column("ncandXimI", (Int_t) XiminusList.GetLength()); ntpXi->Column("nMcXim", (Float_t) nXiminusMc); ntpXi->Column("nMcXimI", (Int_t) nXiminusMc); // store info about initial 4-vector qa.qaP4("beam_", ini, ntpXi); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("Xim_", curCand, ntpXi); andre::qaDetHits("Xim_", curCand, ntpXi); // store info about event shapes qa.qaEventShapeShort("es_", &evsh, ntpXi); // point of closest approach for the daughters qa.qaPoca("fPoca_", curCand, ntpXi); // get the vertex fitter result PndKinVtxFitter * vertexFitter = Ximinus_probPresorter.GetVtxFit(i_Xi); RhoCandidate * XiminusFit = curCand->GetFit(); // store info of vertex fitting ntpXi->Column("fVtx_HowGood", (Int_t) Ximinus_probPresorter.GetVtxBestFitIndex(i_Xi)); ntpXi->Column("fVtx_ProbIdx", (Int_t) Ximinus_probPresorter.GetVtxProbIndex(i_Xi)); andre::qaFitter(&qa, "fVtx_", XiminusFit, ntpXi, vertexFitter); ntpXi->Column("fVtx_ProbRatio", (Float_t) Ximinus_probPresorter.GetVtxRatioToBestProb(i_Xi) ); ntpXi->Column("fVtx_d0ProbRatio", (Float_t) Lambda_probPresorter.GetVtxRatioToBestProbUID(curCand->Daughter(0)->Uid()) ); ntpXi->Column("fVtx_d0ProbRatio2nd", (Float_t) Lambda_probPresorter.GetMassRatio2ndToBestProb() ); ntpXi->Column("fVtx_d0ProbIdx", (Int_t) Lambda_probPresorter.GetVtxProbIndexUID(curCand->Daughter(0)->Uid()) ); // propagate the candidate back to z=0 plane and get distance to IP vtxPocaFitter->GetPocaVtx(vtxPoca, curCand); if ( myProp.SetParticleTrackPar(curCand, vtxPoca) ) { myProp.PropagateWithGeane(); // .. with Geane to plane z=0 myProp.PropagateWithGeanePocaIP(); // .. with Geane until closest point to IP myProp.PropagateWithHelix(); // .. with helix to plane z=0 } // save the propagate results myProp.qaPropResult("Prop_", ntpXi); // get the mass fitter result PndKinFitter * massFitter = Ximinus_probPresorter.GetMassFit(i_Xi); RhoCandidate * XiminusFit_mass = XiminusFit->GetFit(); andre::qaFitter(&qa, "fMass_", XiminusFit_mass, ntpXi, massFitter); ntpXi->Column("fMass_HowGood", (Int_t) Ximinus_probPresorter.GetMassBestFitIndex(i_Xi)); ntpXi->Column("fMass_ProbIdx", (Int_t) Ximinus_probPresorter.GetMassProbIndex(i_Xi)); ntpXi->Column("fMass_ProbRatio", (Float_t) Ximinus_probPresorter.GetMassRatioToBestProb(i_Xi) ); ntpXi->Column("fMass_d0ProbRatio2nd", (Float_t) Lambda_probPresorter.GetMassRatio2ndToBestProb() ); ntpXi->Column("fMass_d0ProbRatio", (Float_t) Lambda_probPresorter.GetMassRatioToBestProbUID(curCand->Daughter(0)->Uid()) ); ntpXi->Column("fMass_d0ProbIdx", (Int_t) Lambda_probPresorter.GetMassProbIndexUID(curCand->Daughter(0)->Uid()) ); // Differences to MC truth qa.qaMcDiff("Xim_", curCand, ntpXi); qa.qaMcDiff("fVtx_", XiminusFit, ntpXi); qa.qaMcDiff("fMass_", XiminusFit_mass, ntpXi); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate * truth = curCand->GetMcTruth(); TLorentzVector lv(-999., -999., -999., -999.); if (0x0 != truth) { lv = truth->P4(); qa.qaVtx("truth_", truth, ntpXi); } else { qa.qaVtx("truth_", dummyCand, ntpXi); } qa.qaP4("truth_", lv, ntpXi); // // store information of boosted particles // XiminusFit->Boost(-beamBoost); // qa.qaComp("boost_", XiminusFit, ntpXi); ntpXi->DumpData(); } // loop over Ximinus // store information on number of particle candidates if (doPreAnalysis) { ntpPre->Column("ev", (Int_t) evt); // final state particles ntpPre->Column("nKminus", (Int_t) kminus.GetLength()); ntpPre->Column("nKminusMCT", (Int_t) theAnalysis->McTruthMatch(kminus)); ntpPre->Column("nPiplus", (Int_t) piplus.GetLength()); ntpPre->Column("nPiplusMCT", (Int_t) theAnalysis->McTruthMatch(piplus)); ntpPre->Column("nPiminus", (Int_t) piminus.GetLength()); ntpPre->Column("nPiminusMCT", (Int_t) theAnalysis->McTruthMatch(piminus)); ntpPre->Column("nProtons", (Int_t) proton.GetLength()); ntpPre->Column("nProtonsMCT", (Int_t) theAnalysis->McTruthMatch(proton)); ntpPre->Column("nAntiProtons", (Int_t) antiproton.GetLength()); ntpPre->Column("nAntiProtonsMCT", (Int_t) theAnalysis->McTruthMatch(antiproton)); // intermediate Lambda0 ntpPre->Column("nLambda0", (Int_t) Lambda0List.GetLength()); ntpPre->Column("nLambda0MCT", (Int_t) theAnalysis->McTruthMatch(Lambda0List)); ntpPre->Column("nAntiLambda0", (Int_t) aLambda0List.GetLength()); ntpPre->Column("nAntiLambda0MCT", (Int_t) theAnalysis->McTruthMatch(aLambda0List)); // ntpPre->Column("nLambda0Fit", (Int_t) Lambda0FitList.GetLength()); // ntpPre->Column("nLambda0MCTFit", (Int_t) theAnalysis->McTruthMatch(Lambda0FitList)); ntpPre->Column("nAntiLambda0Fit", (Int_t) aLambda0FitList.GetLength()); ntpPre->Column("nAntiLambda0MCTFit", (Int_t) theAnalysis->McTruthMatch(aLambda0FitList)); // intermediate Cascades ntpPre->Column("nXiminus", (Int_t) XiminusList.GetLength()); ntpPre->Column("nXiminusMCT", (Int_t) theAnalysis->McTruthMatch(XiminusList)); ntpPre->Column("nAntiXiplus", (Int_t) aXiplusList.GetLength()); ntpPre->Column("nAntiXiplusMCT", (Int_t) theAnalysis->McTruthMatch(aXiplusList)); // final state particles' data andre::qaCandList(&qa, "", kminus, ntpPreCand, evt); andre::qaCandList(&qa, "", piplus, ntpPreCand, evt); andre::qaCandList(&qa, "", piminus, ntpPreCand, evt); andre::qaCandList(&qa, "", proton, ntpPreCand, evt); andre::qaCandList(&qa, "", antiproton, ntpPreCand, evt); ntpPre->DumpData(); } } // loop over events // *** write out all the data out->cd(); ntpMc->GetInternalTree()->Write("", TObject::kOverwrite); if (doPreAnalysis) { ntpPre->GetInternalTree()->Write("", TObject::kOverwrite); ntpPreCand->GetInternalTree()->Write("", TObject::kOverwrite); } ntpLambda->GetInternalTree()->Write("", TObject::kOverwrite); ntpAntiLambda->GetInternalTree()->Write("", TObject::kOverwrite); ntpAntiXi->GetInternalTree()->Write("", TObject::kOverwrite); ntpXi->GetInternalTree()->Write("", TObject::kOverwrite); out->Save(); }