class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; class RhoTuple; #include #include namespace andi { /** * @brief qa extension for fitters (vtx fit, mass fit, *) * @details A useful routine to save information of the currently employed fit * * Should be part of the official QA tools by now * * @param qa The QA object, needs to be given as this method is not part of the QA tools yet; this parameter is going to be removed once the method is assimilated into the QA class * @param pre Prefixing string for the TTree, something like "fitVtx" * @param c The RhoCandidate * @param n The TTree where information should be saved * @param fit The object of the fitter * @param skip I actually don't know what this parameter is, but it's there for every call in the QA tools */ 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 qa extension for saving MC truth information of all daughters of a composite particle * @details Goes through all daughters of a RhoCandidate and saves MC truth information. This information is missing in the qaComp method */ void qaDaughters(PndRhoTupleQA * qa, TString pre, RhoCandidate * c, RhoTuple * n) { if (0 == n) return; for (int i = 0; i < c->NDaughters(); ++i) { // std::cout << "i: " <Daughter(i); // theAnalysis->McTruthMatch(daughter); RhoCandidate * daughterTrue = daughter->GetMcTruth(); TString dindex = TString::Format("d%i", i); // std::cout << dindex << std::endl; qa->qaCand(pre + dindex + "tr", daughterTrue, n); n->Column(pre + dindex + "trpdg", (Int_t) daughterTrue->PdgCode()); qa->qaMcDiff(pre + dindex, daughter, n); } } /** * @brief qa extension for saving a lot of information for a list of particles * @details Shorthand method to save a basic amount of information (using other qa methods) for every particle in a RhoCandList of particles */ void qaFinalState(PndRhoTupleQA * qa, TString pre, RhoCandList & cl, Int_t ev, RhoTuple * n) { if (0 == n) return; for (int i = 0; i < cl.GetLength(); ++i) { n->Column("evt", (Int_t) ev); TString strIndex = TString::Format("%i", i); RhoCandidate * c = cl[i]; RhoCandidate * truth = c->GetMcTruth(); n->Column(pre + "idx", (Int_t) i); qa->qaComp(pre, c, n); qa->qaCand(pre + "tr", truth, n); n->Column(pre + "trpdg", (Int_t) truth->PdgCode()); qa->qaMcDiff(pre, c, n); n->DumpData(); } } /** * @brief Is a certain composite particle covered by the particles in a list of particles? * @details Sorts a given list of particles into two lists. The first one contains all the particles, which are covered by a given composite particle. The second one contains all particles, which are not. Have a look at the code, it's pretty self-explanatory. * * Purpose of the method is to compare all particles of the isNotPartList particle list with the composite candidate for a semin-inclusive analysis. */ void isPartOfComposite(const RhoCandidate & compositeParticle, const RhoCandList & particleList, RhoCandList & isPartList, RhoCandList & isNotPartList) { for (int i_entry = 0; i_entry < particleList.GetLength(); i_entry++) { bool isPartOfCompositeParticle = false; for (int i_composite = 0; i_composite < compositeParticle.NDaughters(); i_composite++) { if (compositeParticle.Daughter(i_composite)->IsCloneOf(* particleList[i_entry])) { isPartOfCompositeParticle = true; } } if (isPartOfCompositeParticle) { // std::cout << "It Is Part!" << " PDG ID: " << particleList[i_entry]->PdgCode() << std::endl; isPartList.Add(particleList[i_entry]); } else { isNotPartList.Add(particleList[i_entry]); } } } /** * @name Tools for distances for a semi-inclusive analysis * @{ */ double getSimpleDistZ(const RhoCandidate & vertex, const RhoCandidate & testParticle) { TVector3 vtx = vertex.GetPosition(); TVector3 momTest = testParticle.GetMomentum(); // vtx.Print(); // momTest.Print(); return vtx.Z() - TMath::Sqrt(vtx.X() * vtx.X() + vtx.Y() * vtx.Y()) * momTest.Pz() / momTest.Pt(); } std::vector getSimpleDistZ(const RhoCandidate & vertex, const RhoCandList & testList) { std::vector result; for (int i_test = 0; i_test < testList.GetLength(); i_test++) { result.push_back(getSimpleDistZ(vertex, * testList[i_test])); } return result; } double getMinDistZ(const RhoCandidate & vertex, const RhoCandidate & testParticle) { TVector3 vtx = vertex.GetPosition(); TVector3 posTest = testParticle.GetPosition(); TVector3 momTest = testParticle.GetMomentum(); Double_t charge = testParticle.GetCharge(); PndHelixPropagator prop(2.0, posTest, momTest, charge); FairTrackPar closestPosition = prop.PropagateToXYPos(vtx); // vtx.Print(); // momTest.Print(); // std::cout << "closestPosition.GetZ(): " << closestPosition.GetZ() << std::endl; // std::cout << "Vtx: "; // vtx.Print(); // std::cout << std::endl; return closestPosition.GetZ(); } std::vector getMinDistZ(const RhoCandidate & vertex, const RhoCandList & testList) { std::vector result; for (int i_test = 0; i_test < testList.GetLength(); i_test++) { result.push_back(getMinDistZ(vertex, * testList[i_test])); } return result; } double getMean(std::vector vec) { double all = std::accumulate(vec.begin(), vec.end(), 0.0); return all / (double) vec.size(); } int op_abs(double val) {return fabs(val);}; double getAbsMean(std::vector vec) { std::vector resultVec; // resultVec.resize(vec.size()); // std::transform(vec.begin(), vec.end(), resultVec.begin(), andi::op_abs); for (std::vector::iterator it = vec.begin(); it != vec.end(); ++it) { resultVec.push_back(op_abs(*it)); } double all = std::accumulate(resultVec.begin(), resultVec.end(), 0.0); return all / (double) resultVec.size(); } /** * @} */ } void ana(TString prefix = "", int nevts = 0) { // *** some variables int i = 0, j = 0, k = 0, l = 0; gStyle->SetOptFit(1011); // Some flags I set for running specific parts of the code - could also be done via parameters to the method call bool doPreAnalysis = true; bool doPreDaughterAnalysis = true; bool doRecoAnalysis = false; bool isDpmAnalysis = false; // *** the output file for FairRunAna TString OutFile = prefix + "output.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); fRun->SetOutputFile(OutFile); fRun->Init(); // *** create an output file for all histograms TString anaOutput = prefix + "output_ana"; anaOutput += "-"; anaOutput += nevts; anaOutput += ".root"; TFile * out = TFile::Open(anaOutput, "RECREATE"); // *** NTuples, one for every composite particle (also the composite of composites); additionally and optionally, a tuple for every final state partcile // The MC truth tuple is only filled when not running with DPM, because this makes the files veeeery large RhoTuple * ntpDp = new RhoTuple("ntpDp", "Dp Analysis"); if (doPreAnalysis) RhoTuple * ntpPreParticles = new RhoTuple("ntpPre", "Pre Analysis"); if (doPreAnalysis && doPreDaughterAnalysis) { RhoTuple * ntpPreKp = new RhoTuple("ntpPreKp", "Pre Kaon Pos Analysis"); RhoTuple * ntpPreKm = new RhoTuple("ntpPreKm", "Pre Kaon Neg Analysis"); RhoTuple * ntpPrePip = new RhoTuple("ntpPrePip", "Pre Pi Pos Analysis"); RhoTuple * ntpPrePim = new RhoTuple("ntpPrePim", "Pre Pi Neg Analysis"); } RhoTuple * ntpDm = new RhoTuple("ntpDm", "Dm Analysis"); RhoTuple * ntpDIndex = new RhoTuple("ntpDIndex", "D Index"); RhoTuple * ntpDpDm = new RhoTuple("ntpDpDm", "DpDmBar Analysis"); RhoTuple * ntpDpDmRAW = new RhoTuple("ntpDpDmRAW", "DpDmBar Analysis No Cuts"); if (!isDpmAnalysis) RhoTuple * nmc = new RhoTuple("nmc", "Mctruth info"); // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nevts == 0 || nevts > theAnalysis->GetEntries()) nevts = theAnalysis->GetEntries(); // *** RhoCandLists for the analysis RhoCandList kminus, piplus, piplus2, kplus, piminus, piminus2, dpluslist, dminuslist, all, mclist; // *** Mass selector double m0_dplus = TDatabasePDG::Instance()->GetParticle("D+")->Mass(); // Get nominal PDG mass of D+ std::cout << "mass d+ = " << m0_dplus << std::endl; RhoMassParticleSelector * dPlusMassSelector = new RhoMassParticleSelector("dp", m0_dplus, 0.3); // 0.3 is the full width of the window! RhoMassParticleSelector * dMinusMassSelector = new RhoMassParticleSelector("dm", m0_dplus, 0.3); // *** 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 beam double pbarmom = 6.5; 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); // TLorentzVector ini(0, 0, 6.231552, 7.240065); PndRhoTupleQA qa(theAnalysis, pbarmom); // *** // the event loop // *** while (theAnalysis->GetEvent() && i++ < nevts) { ntpDIndex->Column("evt", (Int_t) i); if ((i%100) == 0) cout << "evt " <FillList(mclist, "McTruth"); if (!isDpmAnalysis) qa.qaMcList("", mclist, nmc); if (!isDpmAnalysis) nmc->DumpData(); // *** Setup event shape object theAnalysis->FillList(all, "All", pidSelection); PndEventShape evsh(all, ini, 0.05, 0.1); // *** Select with no PID info theAnalysis->FillList(kplus, "KaonBestPlus", pidSelection); // Best = take only best probablity for Kaon theAnalysis->FillList(kminus, "KaonBestMinus", pidSelection); theAnalysis->FillList(piplus, "PionBestPlus", pidSelection); theAnalysis->FillList(piminus, "PionBestMinus", pidSelection); // ***** // *** combinatorics for D+ -> K- Pi+ Pi+ and the antiparticle's process // ***** // dpluslist.Cleanup(); dpluslist.Combine(kminus, piplus, piplus); dpluslist.Select(dPlusMassSelector); dpluslist.SetType(411); int nDplusMc = theAnalysis->McTruthMatch(dpluslist); // dminuslist.Cleanup(); dminuslist.Combine(kplus, piminus, piminus); dminuslist.Select(dMinusMassSelector); dminuslist.SetType(-411); int nDminusMc = theAnalysis->McTruthMatch(dminuslist); if (doPreAnalysis) { // Pre analysis for looking at the final state particles; only done sometimes, as this leads to large data files ntpPreParticles->Column("ev", (Float_t) i); ntpPreParticles->Column("nDplus", (Int_t) dpluslist.GetLength()); ntpPreParticles->Column("nDminus", (Int_t) dminuslist.GetLength()); ntpPreParticles->Column("nKplus", (Int_t) kplus.GetLength()); ntpPreParticles->Column("nKminus", (Int_t) kminus.GetLength()); ntpPreParticles->Column("nPiplus", (Int_t) piplus.GetLength()); ntpPreParticles->Column("nPiminus", (Int_t) piminus.GetLength()); ntpPreParticles->DumpData(); if (doPreDaughterAnalysis) { andi::qaFinalState(&qa, "Kplus", kplus, i, ntpPreKp); andi::qaFinalState(&qa, "Kminus", kminus, i, ntpPreKm); andi::qaFinalState(&qa, "Piplus", piplus, i, ntpPrePip); andi::qaFinalState(&qa, "Piminus", piminus, i, ntpPrePim); } } // *** // Following pre-analysis computes values to tag each composite D meson for its fit quality // Is it the best fit of all the reconstructed D mesons of the current event? Or is there a D meson candidate with a better fit? The following scheme is used (per event): // 1 - this is the best-fitting D meson candidate in the event // 2 - there's one better fitting D meson candidate; this is the 2nd best... // .. // nCand - this is the worst fitting candidate of all candidates // // negative number: same scheme as before, but the chi2 probability cut (> 0.01) is not passed // // Note: This is highly fragile; Andre adopted a more reliable scheme in which this macro-defined routine is extracted into a standalone class, which is compiled and with this much more secure // *** std::map dP_vtx_chi2ForIndex_good, dP_vtx_chi2ForIndex_bad; std::map dP_mass_chi2ForIndex_good, dP_mass_chi2ForIndex_bad; for (j = 0; j < dpluslist.GetLength(); ++j) { PndKinVtxFitter vertexFitter(dpluslist[j]); // vertexFitter.SetVerbose(); vertexFitter.Fit(); bool failVtxProb = TMath::IsNaN(vertexFitter.GetProb()); bool failVtxChi2 = TMath::IsNaN(vertexFitter.GetChi2()); if (!failVtxChi2 && !failVtxProb) { // std::cout << "Prob = " << vertexFitter.GetProb() << std::endl; // std::cout << ", Chi2 = " << vertexFitter.GetChi2() << std::endl; if (vertexFitter.GetProb() > 0.01) { dP_vtx_chi2ForIndex_good[vertexFitter.GetChi2()] = j; // std::cout << "i = " < 0.01 (" << vertexFitter.GetProb() << "). Chi2 = " << vertexFitter.GetChi2() << " saved " << j << std::endl; } else { // <= 0.01 dP_vtx_chi2ForIndex_bad[vertexFitter.GetChi2()] = j; // std::cout << "i = " <GetFit(); PndKinFitter massFitter(dplusfit); massFitter.AddMassConstraint(m0_dplus); massFitter.Fit(); bool failMassProb = TMath::IsNaN(massFitter.GetProb()); bool failMassChi2 = TMath::IsNaN(massFitter.GetChi2()); if (!failMassProb && !failMassChi2) { if (massFitter.GetProb() > 0.01) { dP_mass_chi2ForIndex_good[massFitter.GetChi2()] = j; } else { dP_mass_chi2ForIndex_bad[massFitter.GetChi2()] = j; } } } } std::map::iterator it_good, it_bad; std::map dP_vtx_indexOfBestFit; int runningI = 0; for (it_good = dP_vtx_chi2ForIndex_good.begin(); it_good != dP_vtx_chi2ForIndex_good.end(); it_good++, runningI++) { dP_vtx_indexOfBestFit[it_good->second] = runningI + 1; } runningI = 0; for (it_bad = dP_vtx_chi2ForIndex_bad.begin(); it_bad != dP_vtx_chi2ForIndex_bad.end(); it_bad++, runningI++) { // std::cout << "it->first = " << it_bad->first << ", it->second = " << it_bad->second << std::endl; dP_vtx_indexOfBestFit[it_bad->second] = - (runningI + 1); } std::map dP_mass_indexOfBestFit; runningI = 0; for (it_good = dP_mass_chi2ForIndex_good.begin(); it_good != dP_mass_chi2ForIndex_good.end(); it_good++, runningI++) { dP_mass_indexOfBestFit[it_good->second] = runningI + 1; } runningI = 0; for (it_bad = dP_mass_chi2ForIndex_bad.begin(); it_bad != dP_mass_chi2ForIndex_bad.end(); it_bad++, runningI++) { // std::cout << "it->first = " << it_bad->first << ", it->second = " << it_bad->second << std::endl; dP_mass_indexOfBestFit[it_bad->second] = - (runningI + 1); } // *** // *** do all kind of analysis and store in N-tuple // *** bool DPVtxOk = false, DPMassOk = false, DPBothOk = false; // bool values are used when assembling two D mesons of opposite charges into a D+D- system RhoCandidate * rcDPBestCandVtx, * rcDPBestCandMass; for (j = 0; j < dpluslist.GetLength(); ++j) { // some general info about event (actually written for each candidate!) ntpDp->Column("ev", (Float_t) i); ntpDp->Column("cand", (Float_t) j); ntpDp->Column("evI", (Int_t) i); ntpDp->Column("candI", (Int_t) j); ntpDp->Column("ncandDplus", (Float_t) dpluslist.GetLength()); ntpDp->Column("ncandDplusI", (Int_t) dpluslist.GetLength()); ntpDp->Column("nmctDplus", (Float_t) nDplusMc); ntpDp->Column("nmctDplusI", (Int_t) nDplusMc); // store info about initial 4-vector qa.qaP4("beam", ini, ntpDp); // store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA) qa.qaComp("Dplus", dpluslist[j], ntpDp); // store Reco information // if (doRecoAnalysis) qa.qaRecoFullTree("recoDplus", dpluslist[j], ntpDp); if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "Dplus", dpluslist[j], ntpDp); // store info about event shapes qa.qaEventShapeShort("es", &evsh, ntpDp); PndKinVtxFitter vertexFitter(dpluslist[j]); // kinfit.AddMassConstraint(m0_dplus); vertexFitter.Fit(); RhoCandidate * dplusfit = dpluslist[j]->GetFit(); andi::qaFitter(&qa, "fVtxDplus", dplusfit, ntpDp, &vertexFitter); ntpDp->Column("fVtxHowGood", (Int_t) dP_vtx_indexOfBestFit[j]); qa.qaMcDiff("fVtxDplus", dplusfit, ntpDp); if (dP_vtx_indexOfBestFit[j] == 1) { rcDPBestCandVtx = dplusfit; DPVtxOk = true; } if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "fVtxDplus", dplusfit, ntpDp); PndKinFitter massFitter(dplusfit); massFitter.AddMassConstraint(m0_dplus); massFitter.Fit(); RhoCandidate * dplusfit_mass = dplusfit->GetFit(); andi::qaFitter(&qa, "fMassDplus", dplusfit_mass, ntpDp, &massFitter); ntpDp->Column("fMassHowGood", (Int_t) dP_mass_indexOfBestFit[j]); if (dP_mass_indexOfBestFit[j] == 1) { DPMassOk = true; } if (dP_vtx_indexOfBestFit[j] == 1 && dP_mass_indexOfBestFit[j] == 1) { rcDPBestCandMass = dplusfit_mass; DPBothOk = true; } if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "fMassDplus", dplusfit_mass, ntpDp); qa.qaMcDiff("Dplus", dpluslist[j], ntpDp); // store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent) RhoCandidate *truth = dpluslist[j]->GetMcTruth(); TLorentzVector lv; if (0x0 != truth) lv = truth->P4(); qa.qaP4("trDplus", lv, ntpDp); if (truth != 0x) { // TVector3 tmp = truth->Daughter(0)->GetPosition(); qa.qaVtx("trDplus", truth, ntpDp); } // semi-inclusive stuff: RhoCandList isList, isNotList; andi::isPartOfComposite(* dpluslist[j], kminus, isList, isNotList); andi::isPartOfComposite(* dpluslist[j], piplus, isList, isNotList); andi::isPartOfComposite(* dpluslist[j], kplus, isList, isNotList); andi::isPartOfComposite(* dpluslist[j], piminus, isList, isNotList); std::vector distIs = andi::getMinDistZ(* dplusfit, isList); std::vector distIsNot = andi::getMinDistZ(* dplusfit, isNotList); ntpDp->Column("fVtxDplusMeanDistComposite", (Float_t) andi::getMean(distIs)); ntpDp->Column("fVtxDplusMeanAbsDistComposite", (Float_t) andi::getAbsMean(distIs)); ntpDp->Column("fVtxDplusLengthComposite", (Int_t) distIs.size()); ntpDp->Column("fVtxDplusMeanDistRest", (Float_t) andi::getMean(distIsNot)); ntpDp->Column("fVtxDplusMeanAbsDistRest", (Float_t) andi::getAbsMean(distIsNot)); ntpDp->Column("fVtxDplusLengthRest", (Float_t) distIsNot.size()); ntpDp->DumpData(); } // END loop over D+ ntpDIndex->Column("hasDplus", (Int_t) dpluslist.GetLength()); ntpDIndex->Column("fVtxDplus", (Int_t) DPVtxOk); ntpDIndex->Column("fMassDplus", (Int_t) DPMassOk); ntpDIndex->Column("fBothDplus", (Int_t) DPBothOk); // D- std::map dM_vtx_chi2ForIndex_good, dM_vtx_chi2ForIndex_bad; std::map dM_mass_chi2ForIndex_good, dM_mass_chi2ForIndex_bad; for (j = 0; j < dminuslist.GetLength(); ++j) { PndKinVtxFitter vertexFitter(dminuslist[j]); // vertexFitter.SetVerbose(); vertexFitter.Fit(); bool failVtxProb = TMath::IsNaN(vertexFitter.GetProb()); bool failVtxChi2 = TMath::IsNaN(vertexFitter.GetChi2()); if (!failVtxChi2 && !failVtxProb) { // std::cout << "Prob = " << vertexFitter.GetProb() << std::endl; // std::cout << ", Chi2 = " << vertexFitter.GetChi2() << std::endl; if (vertexFitter.GetProb() > 0.01) { dM_vtx_chi2ForIndex_good[vertexFitter.GetChi2()] = j; // std::cout << "i = " < 0.01 (" << vertexFitter.GetProb() << "). Chi2 = " << vertexFitter.GetChi2() << " saved " << j << std::endl; } else { // <= 0.01 dM_vtx_chi2ForIndex_bad[vertexFitter.GetChi2()] = j; // std::cout << "i = " <GetFit(); PndKinFitter massFitter(dminusfit); massFitter.AddMassConstraint(m0_dplus); massFitter.Fit(); bool failMassProb = TMath::IsNaN(massFitter.GetProb()); bool failMassChi2 = TMath::IsNaN(massFitter.GetChi2()); if (!failMassProb && !failMassChi2) { if (massFitter.GetProb() > 0.01) { dM_mass_chi2ForIndex_good[massFitter.GetChi2()] = j; } else { dM_mass_chi2ForIndex_bad[massFitter.GetChi2()] = j; } } } } // pre-loop over D- std::map::iterator it_good, it_bad; std::map dM_vtx_indexOfBestFit; int runningI = 0; for (it_good = dM_vtx_chi2ForIndex_good.begin(); it_good != dM_vtx_chi2ForIndex_good.end(); it_good++, runningI++) { dM_vtx_indexOfBestFit[it_good->second] = runningI + 1; } runningI = 0; for (it_bad = dM_vtx_chi2ForIndex_bad.begin(); it_bad != dM_vtx_chi2ForIndex_bad.end(); it_bad++, runningI++) { // std::cout << "it->first = " << it_bad->first << ", it->second = " << it_bad->second << std::endl; dM_vtx_indexOfBestFit[it_bad->second] = - (runningI + 1); } std::map dM_mass_indexOfBestFit; runningI = 0; for (it_good = dM_mass_chi2ForIndex_good.begin(); it_good != dM_mass_chi2ForIndex_good.end(); it_good++, runningI++) { dM_mass_indexOfBestFit[it_good->second] = runningI + 1; } runningI = 0; for (it_bad = dM_mass_chi2ForIndex_bad.begin(); it_bad != dM_mass_chi2ForIndex_bad.end(); it_bad++, runningI++) { // std::cout << "it->first = " << it_bad->first << ", it->second = " << it_bad->second << std::endl; dM_mass_indexOfBestFit[it_bad->second] = - (runningI + 1); } bool DMVtxOk = false, DMMassOk = false, DMBothOk = false; RhoCandidate * rcDMBestCandVtx, * rcDMBestCandMass; for (j = 0; j < dminuslist.GetLength(); ++j) { ntpDm->Column("ev", (Float_t) i); ntpDm->Column("cand", (Float_t) j); ntpDm->Column("evI", (Int_t) i); ntpDm->Column("candI", (Int_t) j); ntpDm->Column("ncandDminus", (Float_t) dminuslist.GetLength()); ntpDm->Column("ncandDminusI", (Int_t) dminuslist.GetLength()); ntpDm->Column("nmctDminus", (Float_t) nDminusMc); ntpDm->Column("nmctDminusI", (Int_t) nDminusMc); qa.qaP4("beam", ini, ntpDm); qa.qaComp("Dminus", dminuslist[j], ntpDm); if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "Dminus", dminuslist[j], ntpDm); qa.qaEventShapeShort("es", &evsh, ntpDm); // store Reco information // if (doRecoAnalysis) qa.qaRecoFullTree("recoDminus", dminuslist[j], ntpDm); PndKinVtxFitter vertexFitter(dminuslist[j]); // kinfit.AddMassConstraint(m0_dplus); vertexFitter.Fit(); // // std::cout << "Fit Chi2/NDF: " << kinfit.GetChi2() << "/" << kinfit.GetNdf() << ", Fit Prob: " << kinfit.GetProb() << std::endl; RhoCandidate * dminusfit = dminuslist[j]->GetFit(); andi::qaFitter(&qa, "fVtxDminus", dminusfit, ntpDm, &vertexFitter); ntpDm->Column("fVtxHowGood", (Int_t) dM_vtx_indexOfBestFit[j]); qa.qaMcDiff("fVtxDminus", dminusfit, ntpDm); if (dM_vtx_indexOfBestFit[j] == 1) { // std::cout << "event: " << i << ", j = " << j << ", dM_vtx_indexOfBestFit[j]: " << dM_vtx_indexOfBestFit[j] << std::endl; // std::cout << "dminusfit: " << dminusfit << std::endl; rcDMBestCandVtx = dminusfit; DMVtxOk = true; } if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "fVtxDminus", dminusfit, ntpDm); PndKinFitter massFitter(dminusfit); massFitter.AddMassConstraint(m0_dplus); massFitter.Fit(); RhoCandidate * dminusfit_mass = dminusfit->GetFit(); andi::qaFitter(&qa, "fMassDminus", dminusfit_mass, ntpDm, &massFitter); ntpDm->Column("fMassHowGood", (Int_t) dM_mass_indexOfBestFit[j]); if (dM_mass_indexOfBestFit[j] == 1) { DMMassOk = true; } if (dM_vtx_indexOfBestFit[j] == 1 && dM_mass_indexOfBestFit[j] == 1) { rcDMBestCandMass = dminusfit_mass; DMBothOk = true; } if (doPreDaughterAnalysis) andi::qaDaughters(&qa, "fMassDminus", dminusfit_mass, ntpDm); qa.qaMcDiff("Dminus", dminuslist[j], ntpDm); RhoCandidate *truth = dminuslist[j]->GetMcTruth(); TLorentzVector lv; if (truth) lv = truth->P4(); qa.qaP4("trDminus", lv, ntpDm); if (truth != 0x) qa.qaVtx("trDminus", truth, ntpDm); RhoCandList isList, isNotList; andi::isPartOfComposite(* dminuslist[j], kminus, isList, isNotList); andi::isPartOfComposite(* dminuslist[j], piplus, isList, isNotList); andi::isPartOfComposite(* dminuslist[j], kplus, isList, isNotList); andi::isPartOfComposite(* dminuslist[j], piminus, isList, isNotList); std::vector distIs = andi::getMinDistZ(* dminusfit, isList); std::vector distIsNot = andi::getMinDistZ(* dminusfit, isNotList); ntpDm->Column("fVtxDminusMeanDistComposite", (Float_t) andi::getMean(distIs)); ntpDm->Column("fVtxDminusMeanAbsDistComposite", (Float_t) andi::getAbsMean(distIs)); ntpDm->Column("fVtxDminusLengthComposite", (Int_t) distIs.size()); ntpDm->Column("fVtxDminusMeanDistRest", (Float_t) andi::getMean(distIsNot)); ntpDm->Column("fVtxDminusMeanAbsDistRest", (Float_t) andi::getAbsMean(distIsNot)); ntpDm->Column("fVtxDminusLengthRest", (Float_t) distIsNot.size()); ntpDm->DumpData(); } // END loop over D- list ntpDIndex->Column("hasDminus", (Int_t) dminuslist.GetLength()); ntpDIndex->Column("fVtxDminus", (Int_t) DMVtxOk); ntpDIndex->Column("fMassDminus", (Int_t) DMMassOk); ntpDIndex->Column("fBothDminus", (Int_t) DMBothOk); // combination: D+D- to pbarp RhoCandidate DpDmSystem, DpDmSystem2; RhoCandList DpDmSystemList; bool hasDD = false, hasGoodDD = false; if (dminuslist.GetLength() > 0 && dpluslist.GetLength() > 0) { if (DPBothOk && DMBothOk) { // DpDmSystem = rcDPBestCandVtx->Combine(rcDMBestCandVtx); DpDmSystem = rcDPBestCandVtx->Combine(rcDMBestCandVtx); DpDmSystem.SetType(88888); hasDD = true; ntpDpDm->Column("ev", (Float_t) i); qa.qaP4("beam", ini, ntpDpDm); qa.qaComp("DD", &DpDmSystem, ntpDpDm); Pnd4CFitter fourCFitter(&DpDmSystem, ini); fourCFitter.Fit(); RhoCandidate * DpDm4CFit = DpDmSystem.GetFit(); andi::qaFitter(&qa, "f4CDD", DpDm4CFit, ntpDpDm, &fourCFitter); if (fourCFitter.GetProb() > 0.01) hasGoodDD = true; for (int k = 0; k < DpDm4CFit->NDaughters(); k++) { RhoCandidate * truth = DpDm4CFit->Daughter(k)->GetMcTruth(); TLorentzVector lv; if (truth) lv = truth->P4(); TString prefix = "trf4CDD"; prefix += "d"; prefix += k; qa.qaP4(prefix, lv, ntpDpDm); if (truth != 0x0) qa.qaVtx(prefix, truth, ntpDpDm); qa.qaMcDiff(prefix, DpDm4CFit->Daughter(k), ntpDpDm); } ntpDpDm->DumpData(); } // The following doesnt go further after event 300 of my firstrun_2_ dataset... // DpDmSystemList.Cleanup(); // DpDmSystemList.Combine(dminuslist, dpluslist); // for (int k = 0; k < DpDmSystemList.GetLength(); k++) { // ntpDpDmRAW->Column("ev", (Float_t) i); // ntpDpDmRAW->Column("nCand", (Float_t) DpDmSystemList.GetLength()); // ntpDpDmRAW->Column("candI", (Float_t) k); // qa.qaComp("DD", DpDmSystemList[k], ntpDpDmRAW); // Pnd4CFitter fourCFitter(&DpDmSystemList[k], ini); // fourCFitter.Fit(); // RhoCandidate * DpDm4CFit = DpDmSystemList[k]->GetFit(); // andi::qaFitter(&qa, "f4CDD", DpDm4CFit, ntpDpDmRAW, &fourCFitter); // for (int k = 0; k < DpDm4CFit->NDaughters(); k++) { // RhoCandidate * truth = DpDm4CFit->Daughter(k)->GetMcTruth(); // TLorentzVector lv; // if (truth) lv = truth->P4(); // TString prefix = "trf4CDD"; // prefix += "d"; // prefix += k; // qa.qaP4(prefix, lv, ntpDpDmRAW); // } // ntpDpDmRAW->DumpData(); // } } ntpDIndex->Column("hasDD", (Int_t) hasDD); ntpDIndex->Column("f4CDD", (Int_t) hasGoodDD); ntpDIndex->DumpData(); } // loop over events // *** write out all the histos out->cd(); ntpDp->GetInternalTree()->Write(); ntpDm->GetInternalTree()->Write(); ntpDIndex->GetInternalTree()->Write(); if (!isDpmAnalysis) nmc->GetInternalTree()->Write(); ntpDpDm->GetInternalTree()->Write(); ntpDpDmRAW->GetInternalTree()->Write(); if (doPreAnalysis) ntpPreParticles->GetInternalTree()->Write(); if (doPreAnalysis && doPreDaughterAnalysis) { ntpPreKp->GetInternalTree()->Write(); ntpPreKm->GetInternalTree()->Write(); ntpPrePip->GetInternalTree()->Write(); ntpPrePim->GetInternalTree()->Write(); } out->Save(); }