diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index 3c5b8cf86fd..3b8e68cdbd8 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -57,6 +57,12 @@ DECLARE_SOA_COLUMN(PtD, ptD, float); //! Transverse mom DECLARE_SOA_COLUMN(PtHadron, ptHadron, float); //! Transverse momentum of Hadron DECLARE_SOA_COLUMN(MD, mD, float); //! Invariant mass of D0 DECLARE_SOA_COLUMN(MDbar, mDbar, float); //! Invariant mass of D0bar +DECLARE_SOA_COLUMN(MlScoreBkgD0, mlScoreBkgD0, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0, mlScoreNonPromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0, mlScorePromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScoreBkgD0bar, mlScoreBkgD0bar, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0bar, mlScoreNonPromptD0bar, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0bar, mlScorePromptD0bar, float); //! ML prompt score for D0 selection DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); //! Tag for D0,D0bar DECLARE_SOA_COLUMN(PoolBin, poolBin, int); //! Pool Bin for the MixedEvent DECLARE_SOA_COLUMN(IsAutoCorrelated, isAutoCorrelated, bool); //! Correlation Status @@ -81,7 +87,7 @@ enum ParticleTypeMcRec { }; } // namespace hf_correlation_d0_hadron -DECLARE_SOA_TABLE(DHadronPair, "AOD", "DHADRONPAIR", //! D0-Hadrons pairs Informations +DECLARE_SOA_TABLE(D0HadronPair, "AOD", "D0HPAIR", //! D0-Hadrons pairs Informations aod::hf_correlation_d0_hadron::DeltaPhi, aod::hf_correlation_d0_hadron::DeltaEta, aod::hf_correlation_d0_hadron::PtD, @@ -89,11 +95,28 @@ DECLARE_SOA_TABLE(DHadronPair, "AOD", "DHADRONPAIR", //! D0-Hadrons pairs Inform aod::hf_correlation_d0_hadron::PoolBin, aod::hf_correlation_d0_hadron::IsAutoCorrelated); -DECLARE_SOA_TABLE(DHadronRecoInfo, "AOD", "DHADRONRECOINFO", //! D0-Hadrons pairs Reconstructed Informations +DECLARE_SOA_TABLE(D0HadronRecoInfo, "AOD", "D0HRECOINFO", //! D0-Hadrons pairs Reconstructed Informations aod::hf_correlation_d0_hadron::MD, aod::hf_correlation_d0_hadron::MDbar, aod::hf_correlation_d0_hadron::SignalStatus); +DECLARE_SOA_TABLE(D0HadronMlInfo, "AOD", "D0HMLINFO", //! D0-Hadrons pairs Machine Learning Information + aod::hf_correlation_d0_hadron::MlScoreBkgD0, + aod::hf_correlation_d0_hadron::MlScoreNonPromptD0, + aod::hf_correlation_d0_hadron::MlScorePromptD0, + aod::hf_correlation_d0_hadron::MlScoreBkgD0bar, + aod::hf_correlation_d0_hadron::MlScoreNonPromptD0bar, + aod::hf_correlation_d0_hadron::MlScorePromptD0bar); + +DECLARE_SOA_TABLE(D0CandRecoInfo, "AOD", "D0CANDRECOINFO", //! Ds candidates Reconstructed Information + aod::hf_correlation_d0_hadron::MD, + aod::hf_correlation_d0_hadron::MDbar, + aod::hf_correlation_d0_hadron::PtD, + aod::hf_correlation_d0_hadron::MlScoreBkgD0, + aod::hf_correlation_d0_hadron::MlScorePromptD0, + aod::hf_correlation_d0_hadron::MlScoreBkgD0bar, + aod::hf_correlation_d0_hadron::MlScorePromptD0bar); + // Note: definition of columns and tables for Lc-Hadron correlation pairs namespace hf_correlation_lc_hadron { diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 9ec7e5fc7c2..3500b09dc92 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -72,6 +72,7 @@ using BinningType = ColumnBinningPolicy>; using SelectedTracks = soa::Filtered>; using SelectedCandidatesData = soa::Filtered>; +using SelectedCandidatesDataMl = soa::Filtered>; using SelectedTracksMcRec = soa::Filtered>; using SelectedCandidatesMcRec = soa::Filtered>; using SelectedCollisionsMcGen = soa::Filtered>; @@ -79,8 +80,6 @@ using SelectedParticlesMcGen = soa::Join d0Sel; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; @@ -89,8 +88,9 @@ struct HfCorrelatorD0HadronsSelection { Configurable ptCandMin{"ptCandMin", -1., "min. cand. pT"}; HfHelper hfHelper; - + SliceCache cache; Preslice perCol = aod::hf_cand::collisionId; + Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; @@ -161,10 +161,10 @@ struct HfCorrelatorD0HadronsSelection { /// D0-Hadron correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) struct HfCorrelatorD0Hadrons { - SliceCache cache; - - Produces entryD0HadronPair; - Produces entryD0HadronRecoInfo; + Produces entryD0HadronPair; + Produces entryD0HadronRecoInfo; + Produces entryD0HadronMlInfo; + Produces entryD0CandRecoInfo; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; @@ -175,21 +175,16 @@ struct HfCorrelatorD0Hadrons { Configurable ptCandMin{"ptCandMin", 1., "min. cand. pT"}; Configurable ptTrackMin{"ptTrackMin", 0.3, "min. track pT"}; Configurable ptTrackMax{"ptTrackMax", 99., "max. track pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; + Configurable> bins{"bins", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D0 meson"}; - Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying D-meson efficiency weights"}; Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; Configurable ptSoftPionMax{"ptSoftPionMax", 3.f * 800.f * std::pow(10.f, -6.f), "max. pT cut for soft pion identification"}; + Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; Configurable correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"}; Configurable storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"}; Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; - ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; - ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; - ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks - - HfHelper hfHelper; - BinningType corrBinning{{zPoolBins, multPoolBins}, true}; int leadingIndex = 0; double massD0{0.}; @@ -197,10 +192,8 @@ struct HfCorrelatorD0Hadrons { double massK{0.}; double softPiMass = 0.14543; // pion mass + Q-value of the D*->D0pi decay - Preslice perCol = aod::hf_cand::collisionId; - - Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + HfHelper hfHelper; + SliceCache cache; Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); @@ -208,6 +201,14 @@ struct HfCorrelatorD0Hadrons { Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); + Preslice perCol = aod::hf_cand::collisionId; + + ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; + ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; + ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks + + BinningType corrBinning{{zPoolBins, multPoolBins}, true}; + HistogramRegistry registry{ "registry", // NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes @@ -273,12 +274,8 @@ struct HfCorrelatorD0Hadrons { /// D0-h correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) void processData(SelectedCollisions::iterator const& collision, SelectedTracks const& tracks, - soa::Join const&) + SelectedCandidatesDataMl const& candidates) { - // protection against empty tables to be sliced - if (selectedD0Candidates.size() == 0) { - return; - } // find leading particle if (correlateD0WithLeadingParticle) { leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); @@ -304,46 +301,54 @@ struct HfCorrelatorD0Hadrons { return; } registry.fill(HIST("hMultiplicity"), nTracks); + std::vector outputMlD0 = {-1., -1., -1.}; + std::vector outputMlD0bar = {-1., -1., -1.}; - auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - - for (const auto& candidate1 : selectedD0CandidatesGrouped) { - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { + for (const auto& candidate : candidates) { + if (std::abs(hfHelper.yD0(candidate)) >= yCandMax || candidate.pt() <= ptCandMin || candidate.pt() >= ptTrackMax) { continue; } - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { continue; } // ========================== Define parameters for soft pion removal ================================ - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + auto ePiK = RecoDecay::e(candidate.pVectorProng0(), massPi) + RecoDecay::e(candidate.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate.pVectorProng0(), massK) + RecoDecay::e(candidate.pVectorProng1(), massPi); // ========================== trigger efficiency ================================ double efficiencyWeight = 1.; if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate.pt())); } // ========================== Fill mass histo ================================ - if (candidate1.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); + if (candidate.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate), efficiencyWeight); + registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate), efficiencyWeight); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0[iclass] = candidate.mlProbD0()[classMl->at(iclass)]; + } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate), efficiencyWeight); + registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate), efficiencyWeight); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0bar[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)]; + } } + entryD0CandRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), outputMlD0[0], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[2]); + // ========================== Fill general histos ================================ - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + registry.fill(HIST("hPtCand"), candidate.pt()); + registry.fill(HIST("hPtProng0"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate.ptProng1()); + registry.fill(HIST("hEta"), candidate.eta()); + registry.fill(HIST("hPhi"), candidate.phi()); + registry.fill(HIST("hY"), hfHelper.yD0(candidate)); + registry.fill(HIST("hSelectionStatus"), candidate.isSelD0bar() + (candidate.isSelD0() * 2)); registry.fill(HIST("hD0Bin"), poolBin); // ============ D-h correlation dedicated section ================================== @@ -353,7 +358,7 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hTrackCounter"), 1); // fill total no. of tracks // Remove D0 daughters by checking track indices bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex())) { if (!storeAutoCorrelationFlag) { continue; } @@ -365,20 +370,20 @@ struct HfCorrelatorD0Hadrons { // ========== soft pion removal =================================================== double invMassDstar1 = 0., invMassDstar2 = 0.; bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto pSum2 = RecoDecay::p2(candidate.pVector(), track.pVector()); auto ePion = track.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; continue; } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; continue; } @@ -386,10 +391,10 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hTrackCounter"), 3); // fill no. of tracks after soft pion removal int signalStatus = 0; - if ((candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if ((candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; } - if ((candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if ((candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; } @@ -399,13 +404,14 @@ struct HfCorrelatorD0Hadrons { } registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle } - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), + entryD0HadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), track.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); + entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]); } // end inner loop (tracks) @@ -417,12 +423,8 @@ struct HfCorrelatorD0Hadrons { void processMcRec(SelectedCollisions::iterator const& collision, SelectedTracksMcRec const& tracks, - soa::Join const&) + SelectedCandidatesMcRec const& candidates) { - // protection against empty tables to be sliced - if (selectedD0candidatesMc.size() == 0) { - return; - } // find leading particle if (correlateD0WithLeadingParticle) { leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); @@ -446,63 +448,62 @@ struct HfCorrelatorD0Hadrons { } registry.fill(HIST("hMultiplicity"), nTracks); - auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); // MC reco level bool flagD0 = false; bool flagD0bar = false; - for (const auto& candidate1 : selectedD0CandidatesGroupedMc) { - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + for (const auto& candidate : candidates) { + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { continue; } - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { + if (std::abs(hfHelper.yD0(candidate)) >= yCandMax || candidate.pt() <= ptCandMin || candidate.pt() >= ptTrackMax) { continue; } double efficiencyWeight = 1.; if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate.pt())); } - if (std::abs(candidate1.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + if (std::abs(candidate.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // fill per-candidate distributions from D0/D0bar true candidates - registry.fill(HIST("hPtCandRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0Rec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1Rec"), candidate1.ptProng1()); - registry.fill(HIST("hEtaRec"), candidate1.eta()); - registry.fill(HIST("hPhiRec"), candidate1.phi()); - registry.fill(HIST("hYRec"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatusRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + registry.fill(HIST("hPtCandRec"), candidate.pt()); + registry.fill(HIST("hPtProng0Rec"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1Rec"), candidate.ptProng1()); + registry.fill(HIST("hEtaRec"), candidate.eta()); + registry.fill(HIST("hPhiRec"), candidate.phi()); + registry.fill(HIST("hYRec"), hfHelper.yD0(candidate)); + registry.fill(HIST("hSelectionStatusRec"), candidate.isSelD0bar() + (candidate.isSelD0() * 2)); } // fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 - if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 - registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { - registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + if (candidate.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 + registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); + } else if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { + registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); } else { - registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar - if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar - registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { - registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + if (candidate.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar + registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); + } else if (candidate.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); } else { - registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); } } // ===================== Define parameters for soft pion removal ======================== - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + auto ePiK = RecoDecay::e(candidate.pVectorProng0(), massPi) + RecoDecay::e(candidate.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate.pVectorProng0(), massK) + RecoDecay::e(candidate.pVectorProng1(), massPi); // ============== D-h correlation dedicated section ==================================== - flagD0 = candidate1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0 = candidate.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate matched to D0 (particle) + flagD0bar = candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate, selected as D0 (particle), is matched to D0bar (antiparticle) // ========== track loop starts here ======================== @@ -511,7 +512,7 @@ struct HfCorrelatorD0Hadrons { // Removing D0 daughters by checking track indices bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex())) { if (!storeAutoCorrelationFlag) { continue; } @@ -522,20 +523,20 @@ struct HfCorrelatorD0Hadrons { // ===== soft pion removal =================================================== double invMassDstar1 = 0, invMassDstar2 = 0; bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto pSum2 = RecoDecay::p2(candidate.pVector(), track.pVector()); auto ePion = track.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; continue; } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; continue; } @@ -551,35 +552,35 @@ struct HfCorrelatorD0Hadrons { } int signalStatus = 0; - if (flagD0 && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (flagD0 && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); } // signal case D0 - if (flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (flagD0bar && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); } // reflection case D0 - if (!flagD0 && !flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (!flagD0 && !flagD0bar && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); } // background case D0 - if (flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (flagD0bar && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); } // signal case D0bar - if (flagD0 && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (flagD0 && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); } // reflection case D0bar - if (!flagD0 && !flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (!flagD0 && !flagD0bar && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); } // background case D0bar - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), + entryD0HadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), track.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); } // end inner loop (Tracks) - } // end of outer loop (D0) + } // end of outer loop (D0) registry.fill(HIST("hZvtx"), collision.posZ()); registry.fill(HIST("hMultV0M"), collision.multFT0M()); } @@ -684,7 +685,7 @@ struct HfCorrelatorD0Hadrons { // ====================== Implement Event mixing on Data =================================== void processDataMixedEvent(SelectedCollisions const& collisions, - SelectedCandidatesData const& candidates, + SelectedCandidatesDataMl const& candidates, SelectedTracks const& tracks) { for (const auto& collision : collisions) { @@ -693,7 +694,7 @@ struct HfCorrelatorD0Hadrons { } auto tracksTuple = std::make_tuple(candidates, tracks); - Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; for (const auto& [c1, tracks1, c2, tracks2] : pairData) { // LOGF(info, "Mixed event collisions: Index = (%d, %d), tracks Size: (%d, %d), Z Vertex: (%f, %f), Pool Bin: (%d, %d)", c1.globalIndex(), c2.globalIndex(), tracks1.size(), tracks2.size(), c1.posZ(), c2.posZ(), corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())),corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()))); // For debug @@ -713,16 +714,24 @@ struct HfCorrelatorD0Hadrons { auto ePion = t2.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); + std::vector outputMlD0 = {-1., -1., -1.}; + std::vector outputMlD0bar = {-1., -1., -1.}; if (t1.isSelD0() >= selectionFlagD0) { if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; } + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0[iclass] = t1.mlProbD0()[classMl->at(iclass)]; + } } if (t1.isSelD0bar() >= selectionFlagD0bar) { if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; } + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0bar[iclass] = t1.mlProbD0bar()[classMl->at(iclass)]; + } } int signalStatus = 0; @@ -743,6 +752,7 @@ struct HfCorrelatorD0Hadrons { bool correlationStatus = false; entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); + entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]); } } } @@ -789,8 +799,8 @@ struct HfCorrelatorD0Hadrons { } } - flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate matched to D0 (particle) + flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate, selected as D0 (particle), is matched to D0bar (antiparticle) int signalStatus = 0; if (flagD0 && (t1.isSelD0() >= selectionFlagD0)) { diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index 416460e4689..2654984c45a 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -36,7 +36,8 @@ using namespace o2::analysis::hf_correlations; namespace o2::aod { -using DHadronPairFull = soa::Join; +using D0HadronPairFull = soa::Join; +using D0HadronPairFullMl = soa::Join; } // namespace o2::aod // string definitions, used for histogram axis labels @@ -58,6 +59,8 @@ AxisSpec axisPtHadron = {11, 0., 11., ""}; AxisSpec axisPoolBin = {9, 0., 9., ""}; AxisSpec axisCorrelationState = {2, 0., 2., ""}; ConfigurableAxis axisMass{"axisMass", {250, 1.65f, 2.15f}, ""}; +ConfigurableAxis binsBdtScore{"binsBdtScore", {100, 0., 1.}, "Bdt output scores"}; +AxisSpec axisBdtScore = {binsBdtScore, "Bdt score"}; // definition of vectors for standard ptbin and invariant mass configurables const int nPtBinsCorrelations = 12; @@ -83,9 +86,14 @@ const double ptHadronMax = 10.0; struct HfTaskCorrelationD0Hadrons { // pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_to_pi_k (i.e. the mass pT bins), but can be redefined via json files - Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{vecPtBinsCorrelations}, "pT bin limits for correlation plots"}; + Configurable> binsCorrelations{"binsCorrelations", std::vector{vecPtBinsCorrelations}, "pT bin limits for correlation plots"}; // pT bins for effiencies: same as above - Configurable> binsEfficiency{"ptBinsForEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; + Configurable> binsEfficiency{"binsEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; + Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; + Configurable> mlOutputPromptD0{"mlOutputPromptD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0{"mlOutputBkgD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; + Configurable> mlOutputPromptD0bar{"mlOutputPromptD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0bar{"mlOutputBkgD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; // signal and sideband region edges, to be defined via json file (initialised to empty) Configurable> signalRegionLeft{"signalRegionLeft", std::vector{vecsignalRegionLeft}, "Inner values of signal region vs pT"}; Configurable> signalRegionRight{"signalRegionRight", std::vector{vecsignalRegionRight}, "Outer values of signal region vs pT"}; @@ -96,7 +104,7 @@ struct HfTaskCorrelationD0Hadrons { Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D meson specie under study"}; Configurable isTowardTransverseAway{"isTowardTransverseAway", false, "Divide into three regions: toward, transverse, and away"}; Configurable leadingParticlePtMin{"leadingParticlePtMin", 0., "Min for leading particle pt"}; - Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying efficiency weights"}; + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying efficiency weights"}; HistogramRegistry registry{ "registry", @@ -191,6 +199,10 @@ struct HfTaskCorrelationD0Hadrons { { int nBinsPtAxis = binsCorrelations->size() - 1; const double* valuesPtAxis = binsCorrelations->data(); + registry.add("hBdtScorePromptD0", "D0 BDT prompt score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScoreBkgD0", "D0 BDT bkg score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScorePromptD0bar", "D0bar BDT prompt score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScoreBkgD0bar", "D0bar BDT bkg score", {HistType::kTH1F, {axisBdtScore}}); registry.get(HIST("hCorrel2DVsPtSignalRegion"))->GetAxis(2)->Set(nBinsPtAxis, valuesPtAxis); registry.get(HIST("hCorrel2DVsPtSidebands"))->GetAxis(2)->Set(nBinsPtAxis, valuesPtAxis); registry.get(HIST("hCorrel2DVsPtSignalRegion"))->Sumw2(); @@ -233,8 +245,26 @@ struct HfTaskCorrelationD0Hadrons { /// D-h correlation pair filling task, from pair tables - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) /// Works on both USL and LS analyses pair tables - void processData(aod::DHadronPairFull const& pairEntries) + void processData(aod::D0HadronPairFullMl const& pairEntries, + aod::D0CandRecoInfo const& candidates) { + for (const auto& candidate : candidates) { + float ptD = candidate.ptD(); + float bdtScorePromptD0 = candidate.mlScorePromptD0(); + float bdtScoreBkgD0 = candidate.mlScoreBkgD0(); + float bdtScorePromptD0bar = candidate.mlScorePromptD0bar(); + float bdtScoreBkgD0bar = candidate.mlScoreBkgD0bar(); + int effBinD = o2::analysis::findBin(binsEfficiency, ptD); + if (bdtScorePromptD0 < mlOutputPromptD0->at(effBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(effBinD) || + bdtScorePromptD0bar < mlOutputPromptD0bar->at(effBinD) || bdtScoreBkgD0bar > mlOutputBkgD0bar->at(effBinD)) { + continue; + } + registry.fill(HIST("hBdtScorePromptD0"), bdtScorePromptD0); + registry.fill(HIST("hBdtScoreBkgD0"), bdtScoreBkgD0); + registry.fill(HIST("hBdtScorePromptD0bar"), bdtScorePromptD0bar); + registry.fill(HIST("hBdtScoreBkgD0bar"), bdtScoreBkgD0bar); + } + for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities double deltaPhi = pairEntry.deltaPhi(); @@ -248,6 +278,10 @@ struct HfTaskCorrelationD0Hadrons { int ptBinD = o2::analysis::findBin(binsCorrelations, ptD); int poolBin = pairEntry.poolBin(); bool isAutoCorrelated = pairEntry.isAutoCorrelated(); + float bdtScorePromptD0 = pairEntry.mlScorePromptD0(); + float bdtScoreBkgD0 = pairEntry.mlScoreBkgD0(); + float bdtScorePromptD0bar = pairEntry.mlScorePromptD0bar(); + float bdtScoreBkgD0bar = pairEntry.mlScoreBkgD0bar(); // reject entries outside pT ranges of interest if (ptBinD < 0 || effBinD < 0) { continue; @@ -255,7 +289,10 @@ struct HfTaskCorrelationD0Hadrons { if (ptHadron > ptHadronMax) { ptHadron = ptHadronMax + 0.5; } - + if (bdtScorePromptD0 < mlOutputPromptD0->at(ptBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(ptBinD) || + bdtScorePromptD0bar < mlOutputPromptD0bar->at(ptBinD) || bdtScoreBkgD0bar > mlOutputBkgD0bar->at(ptBinD)) { + continue; + } double efficiencyWeight = 1.; if (applyEfficiency) { efficiencyWeight = 1. / (efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, ptD))); // ***** track efficiency to be implemented ***** @@ -362,7 +399,7 @@ struct HfTaskCorrelationD0Hadrons { } PROCESS_SWITCH(HfTaskCorrelationD0Hadrons, processData, "Process data", false); - void processMcRec(aod::DHadronPairFull const& pairEntries) + void processMcRec(aod::D0HadronPairFull const& pairEntries) { for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities @@ -580,7 +617,7 @@ struct HfTaskCorrelationD0Hadrons { PROCESS_SWITCH(HfTaskCorrelationD0Hadrons, processMcRec, "Process MC Reco mode", true); /// D-Hadron correlation pair filling task, from pair tables - for MC gen-level analysis (no filter/selection, only true signal) - void processMcGen(aod::DHadronPairFull const& pairEntries) + void processMcGen(aod::D0HadronPairFull const& pairEntries) { for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities