diff --git a/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h b/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h index b051294c4f3..d4fddd28771 100644 --- a/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h +++ b/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h @@ -140,6 +140,39 @@ DECLARE_SOA_TABLE(K0sReducedCandidatesMcReco, "AOD", "K0SCANDMCRECO", lf_selection_k0s_reduced::Phi, lf_selection_k0s_reduced::InMassRegion); +namespace lf_selection_xi_reduced +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); + +DECLARE_SOA_COLUMN(M, m, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Y, y, float); +DECLARE_SOA_COLUMN(Phi, phi, float); + +DECLARE_SOA_DYNAMIC_COLUMN(InMassRegion, inMassRegion, + [](float m, float minM, float maxM) -> bool { + return (m >= minM && m <= maxM); + }); +} // namespace lf_selection_xi_reduced + +DECLARE_SOA_TABLE(XiReducedCandidatesData, "AOD", "XICANDDATA", + soa::Index<>, + lf_selection_xi_reduced::CollisionId, + lf_selection_xi_reduced::M, + lf_selection_xi_reduced::Pt, + lf_selection_xi_reduced::Y, + lf_selection_xi_reduced::Phi, + lf_selection_xi_reduced::InMassRegion); + +DECLARE_SOA_TABLE(XiReducedCandidatesMcReco, "AOD", "XICANDMCRECO", + soa::Index<>, + lf_selection_xi_reduced::CollisionId, + lf_selection_xi_reduced::M, + lf_selection_xi_reduced::Pt, + lf_selection_xi_reduced::Y, + lf_selection_xi_reduced::Phi, + lf_selection_xi_reduced::InMassRegion); + namespace lf_selection_pion_track { DECLARE_SOA_INDEX_COLUMN(Collision, collision); diff --git a/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx b/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx index 87baf4ea3c4..816c35150b5 100644 --- a/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx +++ b/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx @@ -13,6 +13,8 @@ /// \brief Table producer for Phi-Strangeness correlation analysis /// \author Stefano Cannito (stefano.cannito@cern.ch) +// o2-linter: disable=name/workflow-file (File contains multiple necessary structs) + #include "PWGLF/DataModel/LFPhiStrangeCorrelationTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/mcCentrality.h" @@ -152,21 +154,27 @@ struct PhiMesonCandProducer { template bool selectionTrackResonance(const T& track) { - if (trackConfigs.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + if (trackConfigs.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) { return false; - if (trackConfigs.cfgPVContributor && !track.isPVContributor()) + } + if (trackConfigs.cfgPVContributor && !track.isPVContributor()) { return false; + } - if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) { return false; + } - if (track.pt() < trackConfigs.cMinKaonPtcut) + if (track.pt() < trackConfigs.cMinKaonPtcut) { return false; + } - if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPhi->at(0) + (trackConfigs.cMaxDCArToPVPhi->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPhi->at(2)))) + if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPhi->at(0) + (trackConfigs.cMaxDCArToPVPhi->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPhi->at(2)))) { return false; - if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) + } + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) { return false; + } return true; } @@ -175,10 +183,12 @@ struct PhiMesonCandProducer { template bool selectionPIDKaonpTdependent(const T& track) { - if (track.pt() < trackConfigs.pTToUseTOF && std::abs(track.tpcNSigmaKa()) < trackConfigs.nSigmaCutTPCKa) + if (track.pt() < trackConfigs.pTToUseTOF && std::abs(track.tpcNSigmaKa()) < trackConfigs.nSigmaCutTPCKa) { return true; - if (track.pt() >= trackConfigs.pTToUseTOF && track.hasTOF() && (std::pow(track.tofNSigmaKa(), 2) + std::pow(track.tpcNSigmaKa(), 2)) < std::pow(trackConfigs.nSigmaCutCombinedKa, 2)) + } + if (track.pt() >= trackConfigs.pTToUseTOF && track.hasTOF() && (std::pow(track.tofNSigmaKa(), 2) + std::pow(track.tpcNSigmaKa(), 2)) < std::pow(trackConfigs.nSigmaCutCombinedKa, 2)) { return true; + } return false; } @@ -199,59 +209,82 @@ struct PhiMesonCandProducer { { // Compile-time selection of the track partition to use based on the type of the analysis (Data or MCReco) auto posThisColl = [&]() { - if constexpr (isMC) + if constexpr (isMC) { return posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - else + } else { return posTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + } }(); auto negThisColl = [&]() { - if constexpr (isMC) + if constexpr (isMC) { return negMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - else + } else { return negTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + } }(); // Loop over positive tracks for (const auto& track1 : posThisColl) { - if (!selectionTrackResonance(track1) || !selectionPIDKaonpTdependent(track1)) + if (!selectionTrackResonance(track1) || !selectionPIDKaonpTdependent(track1)) { continue; + } if constexpr (isMC) { - if (!track1.has_mcParticle()) + if (!track1.has_mcParticle()) { + continue; + } + if (!mcParticlesOpt) { continue; + } + const auto& mcParticles = mcParticlesOpt.value().get(); const auto track1McParticle = mcParticles.rawIteratorAt(track1.mcParticleId()); - if (track1McParticle.pdgCode() != PDG_t::kKPlus || !track1McParticle.isPhysicalPrimary()) + if (track1McParticle.pdgCode() != PDG_t::kKPlus || !track1McParticle.isPhysicalPrimary()) { continue; + } } // Loop over negative tracks for (const auto& track2 : negThisColl) { - if (!selectionTrackResonance(track2) || !selectionPIDKaonpTdependent(track2)) + if (!selectionTrackResonance(track2) || !selectionPIDKaonpTdependent(track2)) { continue; + } if constexpr (isMC) { - if (!track2.has_mcParticle()) + if (!track2.has_mcParticle()) { + continue; + } + if (!mcParticlesOpt) { continue; + } + const auto& mcParticles = mcParticlesOpt.value().get(); const auto track2McParticle = mcParticles.rawIteratorAt(track2.mcParticleId()); - if (track2McParticle.pdgCode() != PDG_t::kKMinus || !track2McParticle.isPhysicalPrimary()) + if (track2McParticle.pdgCode() != PDG_t::kKMinus || !track2McParticle.isPhysicalPrimary()) { continue; + } } // Kinematic reconstruction (common for data and MC) ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa); - if (recPhi.Pt() < phiConfigs.minPhiPt) + if (recPhi.Pt() < phiConfigs.minPhiPt) { continue; - if (recPhi.M() > phiConfigs.maxMPhi) + } + if (recPhi.M() > phiConfigs.maxMPhi) { continue; - if (std::abs(recPhi.Rapidity()) > phiConfigs.cfgYAcceptance) + } + if (std::abs(recPhi.Rapidity()) > phiConfigs.cfgYAcceptance) { continue; + } // PDG check and MC truth association for MCReco analysis if constexpr (isMC) { + if (!mcParticlesOpt) { + continue; + } + const auto& mcParticles = mcParticlesOpt.value().get(); const auto track1McParticle = mcParticles.rawIteratorAt(track1.mcParticleId()); @@ -263,9 +296,10 @@ struct PhiMesonCandProducer { auto genPhiMaybe = [&]() -> std::optional { for (const auto& mother1Index : track1mcPartMotherIndexes) { for (const auto& mother2Index : track2mcPartMotherIndexes) { - if (mother1Index != mother2Index) + if (mother1Index != mother2Index) { continue; - const auto motherMcParticle = mcParticles.rawIteratorAt(mother1Index); + } + auto motherMcParticle = mcParticles.rawIteratorAt(mother1Index); if (std::abs(motherMcParticle.pdgCode()) == o2::constants::physics::Pdg::kPhi) { return motherMcParticle; } @@ -274,9 +308,10 @@ struct PhiMesonCandProducer { return std::nullopt; }(); - if (!genPhiMaybe) + if (!genPhiMaybe) { continue; - const auto genPhi = *genPhiMaybe; + } + const auto& genPhi = *genPhiMaybe; phimesonCandidatesMcReco(collision.globalIndex(), recPhi.M(), genPhi.pt(), genPhi.y(), genPhi.phi()); } else { @@ -319,36 +354,45 @@ struct PhiMesonCandProducer { { if (trueGenPhi) { for (const auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.pdgCode()) != o2::constants::physics::Pdg::kPhi) + if (std::abs(mcParticle.pdgCode()) != o2::constants::physics::Pdg::kPhi) { continue; - if (mcParticle.pt() < phiConfigs.minPhiPt) + } + if (mcParticle.pt() < phiConfigs.minPhiPt) { continue; - if (std::abs(mcParticle.y()) > phiConfigs.cfgYAcceptance) + } + if (std::abs(mcParticle.y()) > phiConfigs.cfgYAcceptance) { continue; + } phimesonCandidatesMcGen(mcCollision.globalIndex(), 0, mcParticle.pt(), mcParticle.y(), mcParticle.phi()); } } else { for (const auto& mcParticle1 : mcParticles) { - if (!mcParticle1.isPhysicalPrimary() || std::abs(mcParticle1.eta()) > trackConfigs.etaMax) + if (!mcParticle1.isPhysicalPrimary() || std::abs(mcParticle1.eta()) > trackConfigs.etaMax) { continue; + } for (const auto& mcParticle2 : mcParticles) { - if (!mcParticle2.isPhysicalPrimary() || std::abs(mcParticle2.eta()) > trackConfigs.etaMax) + if (!mcParticle2.isPhysicalPrimary() || std::abs(mcParticle2.eta()) > trackConfigs.etaMax) { continue; + } - if (!(mcParticle1.pdgCode() == PDG_t::kKPlus && mcParticle2.pdgCode() == PDG_t::kKMinus) && - !(mcParticle1.pdgCode() == PDG_t::kKMinus && mcParticle2.pdgCode() == PDG_t::kKPlus)) + if ((mcParticle1.pdgCode() != PDG_t::kKPlus || mcParticle2.pdgCode() != PDG_t::kKMinus) && + (mcParticle1.pdgCode() != PDG_t::kKMinus || mcParticle2.pdgCode() != PDG_t::kKPlus)) { continue; + } ROOT::Math::PxPyPzMVector genKPair = recMother(mcParticle1, mcParticle2, massKa, massKa); - if (genKPair.Pt() < phiConfigs.minPhiPt) + if (genKPair.Pt() < phiConfigs.minPhiPt) { continue; - if (genKPair.M() > phiConfigs.maxMPhi) + } + if (genKPair.M() > phiConfigs.maxMPhi) { continue; - if (std::abs(genKPair.Rapidity()) > phiConfigs.cfgYAcceptance) + } + if (std::abs(genKPair.Rapidity()) > phiConfigs.cfgYAcceptance) { continue; + } phimesonCandidatesMcGen(mcCollision.globalIndex(), genKPair.M(), genKPair.Pt(), genKPair.Rapidity(), genKPair.Phi()); } @@ -382,7 +426,7 @@ struct K0sReducedCandProducer { // Configurables for V0 selection struct : ConfigurableGroup { Configurable v0SettingCosPA{"v0SettingCosPA", 0.98f, "V0 CosPA"}; - Configurable v0SettingRadius{"v0SettingRadius", 0.5f, "v0radius"}; + Configurable v0SettingRadius{"v0SettingRadius", 0.5f, "V0 decay radius"}; Configurable v0SettingDCAV0Dau{"v0SettingDCAV0Dau", 1.0f, "DCA V0 Daughters"}; Configurable v0SettingDCAPosToPV{"v0SettingDCAPosToPV", 0.1f, "DCA Pos To PV"}; Configurable v0SettingDCANegToPV{"v0SettingDCANegToPV", 0.1f, "DCA Neg To PV"}; @@ -402,15 +446,17 @@ struct K0sReducedCandProducer { Configurable> binspTK0S{"binspTK0S", {0.1, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S"}; // Constants - static constexpr double massK0S = o2::constants::physics::MassK0Short; - static constexpr double massLambda = o2::constants::physics::MassLambda0; + static constexpr double MassK0S = o2::constants::physics::MassK0Short; + static constexpr double MassLambda = o2::constants::physics::MassLambda0; // Filter on default selected collisions // Filter collisionFilter = aod::lf_selection_default_collision::defaultSel == true; Filter collisionFilter = aod::lf_selection_event::defaultSel == true; // Defining filters on V0s (cannot filter on dynamic columns) - Filter v0PreFilter = (nabs(aod::v0data::dcapostopv) > v0Configs.v0SettingDCAPosToPV && nabs(aod::v0data::dcanegtopv) > v0Configs.v0SettingDCANegToPV && aod::v0data::dcaV0daughters < v0Configs.v0SettingDCAV0Dau); + Filter v0PreFilter = (nabs(aod::v0data::dcapostopv) > v0Configs.v0SettingDCAPosToPV && + nabs(aod::v0data::dcanegtopv) > v0Configs.v0SettingDCANegToPV && + aod::v0data::dcaV0daughters < v0Configs.v0SettingDCAV0Dau); // Defining the type of the collisions for data and MC using SelCollisions = soa::Join; @@ -439,17 +485,22 @@ struct K0sReducedCandProducer { template bool selectionTrackStrangeness(const T& track) { - if (!track.hasTPC()) + if (!track.hasTPC()) { return false; - if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + } + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) { return false; - if (track.tpcNClsCrossedRows() < trackConfigs.minNCrossedRowsTPC) + } + if (track.tpcNClsCrossedRows() < trackConfigs.minNCrossedRowsTPC) { return false; - if (track.tpcChi2NCl() > trackConfigs.maxChi2TPC) + } + if (track.tpcChi2NCl() > trackConfigs.maxChi2TPC) { return false; + } - if (std::abs(track.eta()) > trackConfigs.etaMax) + if (std::abs(track.eta()) > trackConfigs.etaMax) { return false; + } return true; } @@ -462,34 +513,44 @@ struct K0sReducedCandProducer { const auto& posDaughterTrack = v0.template posTrack_as(); const auto& negDaughterTrack = v0.template negTrack_as(); - if (!selectionTrackStrangeness(posDaughterTrack) || !selectionTrackStrangeness(negDaughterTrack)) + if (!selectionTrackStrangeness(posDaughterTrack) || !selectionTrackStrangeness(negDaughterTrack)) { return false; + } if constexpr (!isMC) { - if (std::abs(posDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) + if (std::abs(posDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) { return false; - if (std::abs(negDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) + } + if (std::abs(negDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) { return false; + } } - if (v0.v0cosPA() < v0Configs.v0SettingCosPA) + if (v0.v0cosPA() < v0Configs.v0SettingCosPA) { return false; - if (v0.v0radius() < v0Configs.v0SettingRadius) + } + if (v0.v0radius() < v0Configs.v0SettingRadius) { return false; - if (v0.pt() < v0Configs.v0SettingMinPt) + } + if (v0.pt() < v0Configs.v0SettingMinPt) { return false; + } if (v0Configs.cfgFurtherV0Selection) { - if (v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massK0S > v0Configs.ctauK0s) + if (v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassK0S > v0Configs.ctauK0s) { return false; - if (v0.qtarm() < (v0Configs.paramArmenterosCut * std::abs(v0.alpha()))) + } + if (v0.qtarm() < (v0Configs.paramArmenterosCut * std::abs(v0.alpha()))) { return false; - if (std::abs(v0.mLambda() - massLambda) < v0Configs.v0rejK0s) + } + if (std::abs(v0.mLambda() - MassLambda) < v0Configs.v0rejK0s) { return false; + } } - if (std::abs(v0.yK0Short()) > v0Configs.cfgYAcceptance) + if (std::abs(v0.yK0Short()) > v0Configs.cfgYAcceptance) { return false; + } return true; } @@ -499,8 +560,9 @@ struct K0sReducedCandProducer { { for (const auto& v0 : V0s) { // Cut on V0 dynamic columns - if (!selectionV0(v0, collision)) + if (!selectionV0(v0, collision)) { continue; + } histos.fill(HIST("h3K0sCandidatesMass"), collision.centFT0M(), v0.pt(), v0.mK0Short()); @@ -514,14 +576,17 @@ struct K0sReducedCandProducer { void processMCReco(FilteredSimCollisions::iterator const& collision, FullMCV0s const& V0s, V0DauMCTracks const&, aod::McParticles const& mcParticles) { for (const auto& v0 : V0s) { - if (!selectionV0(v0, collision)) + if (!selectionV0(v0, collision)) { continue; - if (!v0.has_mcParticle()) + } + if (!v0.has_mcParticle()) { continue; + } const auto& v0McParticle = mcParticles.rawIteratorAt(v0.mcParticleId()); - if (std::abs(v0McParticle.pdgCode()) != PDG_t::kK0Short || !v0McParticle.isPhysicalPrimary()) + if (std::abs(v0McParticle.pdgCode()) != PDG_t::kK0Short || !v0McParticle.isPhysicalPrimary()) { continue; + } k0sReducedCandidatesMcReco(collision.globalIndex(), v0.mK0Short(), v0McParticle.pt(), v0McParticle.y(), v0McParticle.phi()); } @@ -530,6 +595,204 @@ struct K0sReducedCandProducer { PROCESS_SWITCH(K0sReducedCandProducer, processMCReco, "Process function to select reduced K0s candidates in MCReco w MC truth", false); }; +struct XiReducedCandProducer { + // Produce the table with the Xi candidates information + Produces xiReducedCandidatesData; + Produces xiReducedCandidatesMcReco; + + HistogramRegistry histos{"xiReducedCandidates", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // Configurable on multiplicity bins + Configurable> binsMult{"binsMult", {0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "Multiplicity bin limits"}; + + // Configurables for tracks selection + struct : ConfigurableGroup { + Configurable etaMax{"etaMax", 0.8f, "eta max"}; + Configurable nSigmaCutTPCSecPion{"nSigmaCutTPCSecPion", 4.0f, "Value of the TPC Nsigma cut for secondary Pions"}; + Configurable nSigmaCutTPCSecProton{"nSigmaCutTPCSecProton", 4.0f, "Value of the TPC Nsigma cut for secondary Protons"}; + + Configurable minTPCnClsFound{"minTPCnClsFound", 70, "min number of found TPC clusters"}; + Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 70, "min number of TPC crossed rows"}; + Configurable maxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; + } trackConfigs; + + // Configurables for Cascade selection + struct : ConfigurableGroup { + Configurable cascSettingV0CosPA{"cascSettingV0CosPA", 0.98f, "V0 CosPA"}; + Configurable cascSettingCosPA{"cascSettingCosPA", 0.98f, "Cascade CosPA"}; + Configurable cascSettingV0Radius{"cascSettingV0Radius", 0.5f, "V0 decay radius"}; + Configurable cascSettingRadius{"cascSettingRadius", 0.5f, "Cascade decay radius"}; + Configurable cascSettingDCAV0Dau{"cascSettingDCAV0Dau", 1.0f, "DCA V0 Daughters"}; + Configurable cascSettingDCACascDau{"cascSettingDCACascDau", 1.0f, "DCA Cascade Daughters"}; + Configurable cascSettingDCAPosToPV{"cascSettingDCAPosToPV", 0.1f, "DCA Pos To PV"}; + Configurable cascSettingDCANegToPV{"cascSettingDCANegToPV", 0.1f, "DCA Neg To PV"}; + Configurable cascSettingDCABachToPV{"cascSettingDCABachToPV", 0.1f, "DCA Bach To PV"}; + Configurable cascSettingMinPt{"cascSettingMinPt", 0.8f, "Cascade min pt"}; + + Configurable cfgYAcceptance{"cfgYAcceptance", 0.5f, "Rapidity acceptance"}; + } cascadeConfigs; + + // Configurables for Xi pT bins + Configurable> binspTXi{"binspTXi", {0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for Xi"}; + + // Constants + static constexpr double MassXi = o2::constants::physics::MassXiMinus; + + // Filter on default selected collisions + Filter collisionFilter = aod::lf_selection_event::defaultSel == true; + + // Defining filters on Cascades (cannot filter on dynamic columns) + Filter cascadePreFilter = (nabs(aod::cascdata::dcapostopv) > cascadeConfigs.cascSettingDCAPosToPV && + nabs(aod::cascdata::dcanegtopv) > cascadeConfigs.cascSettingDCANegToPV && + nabs(aod::cascdata::dcabachtopv) > cascadeConfigs.cascSettingDCABachToPV && + aod::cascdata::dcaV0daughters < cascadeConfigs.cascSettingDCAV0Dau && + aod::cascdata::dcacascdaughters < cascadeConfigs.cascSettingDCACascDau); + + // Defining the type of the collisions for data and MC + using SelCollisions = soa::Join; + using SimCollisions = soa::Join; + + using FilteredSelCollisions = soa::Filtered; + using FilteredSimCollisions = soa::Filtered; + + // Defining the type of the Cascades and corresponding daughter and bachelor tracks for data and MC + using FullCascades = soa::Filtered; + using FullMCCascades = soa::Filtered>; + + using DauTracks = soa::Join; + using DauMCTracks = soa::Join; + + void init(InitContext&) + { + AxisSpec binnedmultAxis{(std::vector)binsMult, "centFT0M"}; + AxisSpec binnedpTXiAxis{(std::vector)binspTXi, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec massXiAxis = {200, 1.2f, 1.4f, "#it{M}_{inv} [GeV/#it{c}^{2}]"}; + + histos.add("h3XiCandidatesMass", "#Xi candidate invariant mass", kTH3F, {binnedmultAxis, binnedpTXiAxis, massXiAxis}); + } + + // Single track selection for strangeness sector + template + bool selectionTrackStrangeness(const T& track) + { + if (!track.hasTPC()) { + return false; + } + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) { + return false; + } + if (track.tpcNClsCrossedRows() < trackConfigs.minNCrossedRowsTPC) { + return false; + } + if (track.tpcChi2NCl() > trackConfigs.maxChi2TPC) { + return false; + } + + if (std::abs(track.eta()) > trackConfigs.etaMax) { + return false; + } + return true; + } + + // Cascade selection + template + bool selectionCascade(const T1& cascade, const T2& collision) + { + using DauTrackType = std::conditional_t; + + const auto& posDaughterTrack = cascade.template posTrack_as(); + const auto& negDaughterTrack = cascade.template negTrack_as(); + const auto& bachDaughterTrack = cascade.template bachelor_as(); + + if (!selectionTrackStrangeness(posDaughterTrack) || !selectionTrackStrangeness(negDaughterTrack) || !selectionTrackStrangeness(bachDaughterTrack)) { + return false; + } + + if constexpr (!isMC) { + if (cascade.sign() < 0) { + if (std::abs(posDaughterTrack.tpcNSigmaPr()) > trackConfigs.nSigmaCutTPCSecProton) { + return false; + } + if (std::abs(negDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) { + return false; + } + } else if (cascade.sign() > 0) { + if (std::abs(posDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) { + return false; + } + if (std::abs(negDaughterTrack.tpcNSigmaPr()) > trackConfigs.nSigmaCutTPCSecProton) { + return false; + } + } + if (std::abs(bachDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) { + return false; + } + } + + const auto& pvx = collision.posX(); + const auto& pvy = collision.posY(); + const auto& pvz = collision.posZ(); + + if (cascade.v0cosPA(pvx, pvy, pvz) < cascadeConfigs.cascSettingV0CosPA) { + return false; + } + if (cascade.casccosPA(pvx, pvy, pvz) < cascadeConfigs.cascSettingCosPA) { + return false; + } + if (cascade.v0radius() < cascadeConfigs.cascSettingV0Radius) { + return false; + } + if (cascade.cascradius() < cascadeConfigs.cascSettingRadius) { + return false; + } + if (cascade.pt() < cascadeConfigs.cascSettingMinPt) { + return false; + } + + if (std::abs(cascade.yXi()) > cascadeConfigs.cfgYAcceptance) { + return false; + } + return true; + } + + void processData(FilteredSelCollisions::iterator const& collision, FullCascades const& cascades, DauTracks const&) + { + for (const auto& cascade : cascades) { + // Cut on cascade dynamic columns + if (!selectionCascade(cascade, collision)) { + continue; + } + + histos.fill(HIST("h3XiCandidatesMass"), collision.centFT0M(), cascade.pt(), cascade.mXi()); + + xiReducedCandidatesData(collision.globalIndex(), cascade.mXi(), cascade.pt(), cascade.yXi(), cascade.phi()); + } + } + + PROCESS_SWITCH(XiReducedCandProducer, processData, "Process function to select reduced Xi candidates in Data or in McReco (w/o McTruth) analysis", true); + + void processMCReco(FilteredSimCollisions::iterator const& collision, FullMCCascades const& cascades, DauMCTracks const&, aod::McParticles const& mcParticles) + { + for (const auto& cascade : cascades) { + if (!selectionCascade(cascade, collision)) { + continue; + } + if (!cascade.has_mcParticle()) { + continue; + } + + const auto& cascadeMcParticle = mcParticles.rawIteratorAt(cascade.mcParticleId()); + if (std::abs(cascadeMcParticle.pdgCode()) != PDG_t::kXiMinus || !cascadeMcParticle.isPhysicalPrimary()) { + continue; + } + + xiReducedCandidatesMcReco(collision.globalIndex(), cascade.mXi(), cascadeMcParticle.pt(), cascadeMcParticle.y(), cascadeMcParticle.phi()); + } + } + + PROCESS_SWITCH(XiReducedCandProducer, processMCReco, "Process function to select reduced Xi candidates in MCReco w MC truth", false); +}; + struct PionTrackProducer { // Produce the table with the pion tracks information Produces pionTracksData; @@ -574,7 +837,7 @@ struct PionTrackProducer { // Configurable on pion pT bins Configurable> binspTPi{"binspTPi", {0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0}, "pT bin limits for pions"}; - static constexpr double massPi = o2::constants::physics::MassPiPlus; + static constexpr double MassPi = o2::constants::physics::MassPiPlus; // Filter on default selected collisions // Filter collisionFilter = aod::lf_selection_default_collision::defaultSel == true; @@ -666,7 +929,7 @@ struct PionTrackProducer { } // Other hadron species rejection - for (auto pid : {o2::track::PID::Kaon, o2::track::PID::Proton}) { + for (const auto& pid : {o2::track::PID::Kaon, o2::track::PID::Proton}) { auto nSigmaTPC = aod::pidutils::tpcNSigma(pid, track); if (std::abs(nSigmaTPC) < trackConfigs.pidTPCMaxHadrons) { @@ -701,35 +964,44 @@ struct PionTrackProducer { template bool selectionPion(const T& track) { - if (!track.isGlobalTrackWoDCA()) + if (!track.isGlobalTrackWoDCA()) { return false; + } - if (track.itsNCls() < trackConfigs.minITSnCls) + if (track.itsNCls() < trackConfigs.minITSnCls) { return false; - if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + } + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) { return false; + } - if (track.pt() < trackConfigs.cMinPionPtcut) + if (track.pt() < trackConfigs.cMinPionPtcut) { return false; + } - if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPion->at(0) + (trackConfigs.cMaxDCArToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPion->at(2)))) + if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPion->at(0) + (trackConfigs.cMaxDCArToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPion->at(2)))) { return false; + } if (trackConfigs.cfgIsDCAzParameterized) { - if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVPion->at(0) + (trackConfigs.cMaxDCAzToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCAzToPVPion->at(2)))) + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVPion->at(0) + (trackConfigs.cMaxDCAzToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCAzToPVPion->at(2)))) { return false; + } } else { - if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) { return false; + } } - if (trackConfigs.cfgIsTOFChecked && track.pt() >= trackConfigs.tofPIDThreshold && !track.hasTOF()) + if (trackConfigs.cfgIsTOFChecked && track.pt() >= trackConfigs.tofPIDThreshold && !track.hasTOF()) { return false; + } /*if (analysisMode == 1 && !pidSelectionPion(track)) return false;*/ - if (!pidHypothesesRejection(track)) + if (!pidHypothesesRejection(track)) { return false; + } /* if (analysisMode == 1) { @@ -740,8 +1012,9 @@ struct PionTrackProducer { } */ - if (std::abs(track.rapidity(massPi)) > trackConfigs.cfgYAcceptance) + if (std::abs(track.rapidity(MassPi)) > trackConfigs.cfgYAcceptance) { return false; + } return true; } @@ -750,13 +1023,14 @@ struct PionTrackProducer { void processData(FilteredSelCollisions::iterator const& collision, FullTracks const& fullTracks) { for (const auto& track : fullTracks) { - if (!selectionPion(track)) + if (!selectionPion(track)) { continue; + } histos.fill(HIST("h3PionTPCnSigma"), collision.centFT0M(), track.pt(), track.tpcNSigmaPi()); histos.fill(HIST("h3PionTOFnSigma"), collision.centFT0M(), track.pt(), track.tofNSigmaPi()); - pionTracksData(collision.globalIndex(), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt(), track.rapidity(massPi), track.phi(), track.hasTOF()); + pionTracksData(collision.globalIndex(), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt(), track.rapidity(MassPi), track.phi(), track.hasTOF()); } } @@ -766,14 +1040,17 @@ struct PionTrackProducer { void processMCReco(FilteredSimCollisions::iterator const& collision, FullMCTracks const& fullTracks, aod::McParticles const& mcParticles) { for (const auto& track : fullTracks) { - if (!selectionPion(track)) + if (!selectionPion(track)) { continue; - if (!track.has_mcParticle()) + } + if (!track.has_mcParticle()) { continue; + } const auto trackMcParticle = mcParticles.rawIteratorAt(track.mcParticleId()); - if (std::abs(trackMcParticle.pdgCode()) != PDG_t::kPiPlus) + if (std::abs(trackMcParticle.pdgCode()) != PDG_t::kPiPlus) { continue; + } if (trackMcParticle.isPhysicalPrimary()) { histos.fill(HIST("h2RecMCDCAxyPrimPi"), track.pt(), track.dcaXY()); @@ -786,7 +1063,7 @@ struct PionTrackProducer { continue; } - pionTracksMcReco(collision.globalIndex(), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt(), track.rapidity(massPi), track.phi(), track.hasTOF()); + pionTracksMcReco(collision.globalIndex(), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt(), track.rapidity(MassPi), track.phi(), track.hasTOF()); } } @@ -836,7 +1113,7 @@ struct EventSelectionProducer { using FullMCTracks = soa::Join; // Necessary service to flag INEL>0 events in GenMC - Service pdgDB; + Service pdgDB{}; void init(InitContext&) { @@ -882,14 +1159,17 @@ struct EventSelectionProducer { if constexpr (!isMC) { // data event histos.fill(HIST("hEventSelectionData"), 0); // all collisions - if (!collision.sel8()) + if (!collision.sel8()) { return false; + } histos.fill(HIST("hEventSelectionData"), 1); // sel8 collisions - if (std::abs(collision.posZ()) >= cutZVertex) + if (std::abs(collision.posZ()) >= cutZVertex) { return false; + } histos.fill(HIST("hEventSelectionData"), 2); // vertex-Z selected - if (!collision.isInelGt0()) + if (!collision.isInelGt0()) { return false; + } histos.fill(HIST("hEventSelectionData"), 3); // INEL>0 collisions multPercentile = collision.centFT0M(); @@ -897,23 +1177,29 @@ struct EventSelectionProducer { static_assert(!std::is_same_v, "Need to set MC_T to MCCollisions for isMC = true"); histos.fill(HIST("hEventSelectionMC"), 0); // all collisions - if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) { return false; + } histos.fill(HIST("hEventSelectionMC"), 1); // kIsTriggerTVX collisions - if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { return false; + } histos.fill(HIST("hEventSelectionMC"), 2); // kNoTimeFrameBorder collisions - if (cfgiskNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + if (cfgiskNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { return false; + } histos.fill(HIST("hEventSelectionMC"), 3); // kNoITSROFrameBorder collisions (by default not requested by the selection) - if (std::abs(collision.posZ()) > cutZVertex) + if (std::abs(collision.posZ()) > cutZVertex) { return false; + } histos.fill(HIST("hEventSelectionMC"), 4); // vertex-Z selected - if (!collision.isInelGt0()) + if (!collision.isInelGt0()) { return false; + } histos.fill(HIST("hEventSelectionMC"), 5); // INEL>0 collisions - if (!collision.has_mcCollision()) + if (!collision.has_mcCollision()) { return false; + } histos.fill(HIST("hEventSelectionMC"), 6); // with at least a gen collision const auto& mcCollision = collision.template mcCollision_as(); @@ -932,14 +1218,16 @@ struct EventSelectionProducer { bool eventHasPhi(const T1& collision, const T2& phiCandidates) { if (withTrueGenPhi) { - if (phiCandidates.size() < 1) + if (phiCandidates.size() < 1) { return false; + } } else { uint16_t nPhi{0}; for (const auto& phiCand : phiCandidates) { - if (phiCand.inMassRegion(phiConfigs.minMPhiSignal, phiConfigs.maxMPhiSignal)) + if (phiCand.inMassRegion(phiConfigs.minMPhiSignal, phiConfigs.maxMPhiSignal)) { nPhi++; + } // histos.fill(HIST("hEta"), track1.eta()); // histos.fill(HIST("hNsigmaKaonTPC"), track1.tpcInnerParam(), track1.tpcNSigmaKa()); @@ -948,8 +1236,9 @@ struct EventSelectionProducer { // histos.fill(HIST("h2DauTracksPhiDCAz"), track1.pt(), track1.dcaZ()); } - if (nPhi == 0) + if (nPhi == 0) { return false; + } } float multPercentile{0.0f}; @@ -1035,6 +1324,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } diff --git a/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx index 14c2b3c3ad4..2cfef9a857b 100644 --- a/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx +++ b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx @@ -80,23 +80,45 @@ enum AssociatedParticleType { kAssocPartSize }; +enum PhiMassRegion { + kSignalRegion = 0, + kSidebandRegion, + kPhiMassRegions +}; + +template +constexpr int getPdgCode() +{ + if constexpr (PartType == kK0S) { + return PDG_t::kK0Short; + } else if constexpr (PartType == kXi) { + return PDG_t::kXiMinus; + } else if constexpr (PartType == kPion) { + return PDG_t::kPiPlus; + } else { + static_assert(PartType == kK0S || PartType == kXi || PartType == kPion, "Unsupported particle type in getPdgCode"); + return 0; + } +} + using EffMapPtr = std::variant, std::shared_ptr>; struct BoundEfficiencyMap { using CoordsTuple = std::tuple; - const EffMapPtr& effMap; + const EffMapPtr* effMap; CoordsTuple coords; - BoundEfficiencyMap(const EffMapPtr& effMap, float x, float y, float z) : effMap(effMap), coords(x, y, z) {} - BoundEfficiencyMap(const EffMapPtr& effMap, const CoordsTuple& coords) : effMap(effMap), coords(coords) {} + BoundEfficiencyMap(const EffMapPtr& effMap, float x, float y, float z) : effMap(&effMap), coords(x, y, z) {} + BoundEfficiencyMap(const EffMapPtr& effMap, const CoordsTuple& coords) : effMap(&effMap), coords(coords) {} - std::pair getBinEfficiencyAndError() const + [[nodiscard]] std::pair getBinEfficiencyAndError() const { return std::visit( [this](auto&& mapPtr) -> std::pair { - if (!mapPtr) + if (!mapPtr) { return {1.0f, 0.0f}; + } const auto& [x, y, z] = coords; @@ -112,15 +134,16 @@ struct BoundEfficiencyMap { return {mapPtr->GetBinContent(bin), mapPtr->GetBinError(bin)}; } }, - effMap); + *effMap); } - float interpolateEfficiency() const + [[nodiscard]] float interpolateEfficiency() const { return std::visit( [this](auto&& mapPtr) -> float { - if (!mapPtr) + if (!mapPtr) { return 1.0f; + } const auto& [x, y, z] = coords; @@ -132,11 +155,23 @@ struct BoundEfficiencyMap { return mapPtr->Interpolate(x, y, z); // Native 3D interpolation } }, - effMap); + *effMap); } }; -struct PhiStrangenessCorrelation { +template +struct AssocInput { + static constexpr AssociatedParticleType Type = PartType; + TCands const* candidates; +}; + +template +constexpr auto makeAssocInput(TCands const& cands) +{ + return AssocInput{&cands}; +} + +struct PhiStrangeCorrelation { HistogramRegistry histos{"phiStrangenessCorrelation", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // Configurable for selection type @@ -163,6 +198,12 @@ struct PhiStrangenessCorrelation { Configurable> rangeMK0sSignal{"rangeMK0sSignal", {0.47f, 0.53f}, "K0S mass range for signal extraction"}; } k0sConfigs; + // Configurables for Xi selection + struct : ConfigurableGroup { + Configurable selectXiInSigRegion{"selectXiInSigRegion", true, "Select Xi candidates in signal region"}; + Configurable> rangeMXiSignal{"rangeMXiSignal", {1.31f, 1.33f}, "Xi mass range for signal extraction"}; + } xiConfigs; + // Configurables for Pions selection struct : ConfigurableGroup { Configurable selectPionInSigRegion{"selectPionInSigRegion", true, "Select Pion candidates in signal region"}; @@ -179,15 +220,15 @@ struct PhiStrangenessCorrelation { Configurable> binspTK0S{"binspTK0S", {0.1, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S"}; Configurable> binspTK0SExt{"binspTK0SExt", {0.0, 0.1, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S extended for MC Gen"}; + // Configurable on Xi pT bins + Configurable> binspTXi{"binspTXi", {0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for Xi"}; + Configurable> binspTXiExt{"binspTXiExt", {0.0, 0.4, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for Xi extended for MC Gen"}; + // Configurable on pion pT bins Configurable> binspTPi{"binspTPi", {0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0}, "pT bin limits for Pions"}; Configurable> binspTPiExt{"binspTPiExt", {0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0}, "pT bin limits for Pions extended for MC Gen"}; - struct : ConfigurableGroup { - Configurable doK0SCorrelation{"doK0SCorrelation", true, "Enable Phi-K0S correlation"}; - Configurable doXiCorrelation{"doXiCorrelation", false, "Enable Phi-Xi correlation"}; - Configurable doPionCorrelation{"doPionCorrelation", true, "Enable Phi-Pion correlation"}; - } activeCorrelationConfigs; + Configurable> activeCorrelationTypes{"activeCorrelationTypes", {true, false, true}, "Enable correlation types: 0 - K0S, 1 - Xi, 2 - Pion"}; // Configurables for delta y selection struct : ConfigurableGroup { @@ -216,6 +257,7 @@ struct PhiStrangenessCorrelation { struct : ConfigurableGroup { Configurable minPhiPt{"minPhiPt", 0.4f, "Minimum pT for Phi candidates"}; Configurable v0SettingMinPt{"v0SettingMinPt", 0.1f, "V0 min pt"}; + Configurable cascadeSettingMinPt{"cascadeSettingMinPt", 0.8f, "Cascade min pt"}; Configurable cMinPionPtcut{"cMinPionPtcut", 0.2f, "Track minimum pt cut"}; Configurable bypassPtCut{"bypassPtCut", false, "Bypass the minimum pt cut at MCGen level"}; } minPtMcGenConfigs; @@ -240,6 +282,9 @@ struct PhiStrangenessCorrelation { Preslice k0sDataPerCollision = aod::v0::collisionId; Preslice k0sMcRecoPerCollision = aod::v0::collisionId; + Preslice xiDataPerCollision = aod::cascade::collisionId; + Preslice xiMcRecoPerCollision = aod::cascade::collisionId; + Preslice pionTrackDataPerCollision = aod::track::collisionId; Preslice pionTrackMcRecoPerCollision = aod::track::collisionId; @@ -247,12 +292,12 @@ struct PhiStrangenessCorrelation { } preslices; // Necessary service to retrieve efficiency maps from CCDB - Service ccdb; + Service ccdb{}; // std::shared_ptr effMapPhi{}; // std::array, kAssocPartSize> effMapsAssoc{}; - EffMapPtr effMapPhi{}; + EffMapPtr effMapPhi; std::array effMapsAssoc{}; // Binning policy and axes for mixed event @@ -262,18 +307,18 @@ struct PhiStrangenessCorrelation { using BinningTypeVertexCent = ColumnBinningPolicy; BinningTypeVertexCent binningOnVertexAndCent{{axisVertexMixing, axisCentralityMixing}, true}; - static constexpr std::array PhiMassRegionLabels{"Signal", "Sideband"}; + static constexpr std::array PhiMassRegionLabels{"Signal", "Sideband"}; static constexpr std::array AssocParticleLabels{"K0S", "Xi", "Pi"}; // Light structures to store only the necessary information for the correlation analysis at MCGen level struct MiniParticle { - float pt; - float y; - float phi; + float pt{}; + float y{}; + float phi{}; }; struct MiniEvent { - float multiplicity; + float multiplicity{}; std::vector phiParticles; std::vector k0sParticles; std::vector xiParticles; @@ -299,6 +344,10 @@ struct PhiStrangenessCorrelation { AxisSpec pTK0SAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"}; AxisSpec binnedpTK0SAxis{(std::vector)binspTK0S, "#it{p}_{T} (GeV/#it{c})"}; AxisSpec binnedpTK0SAxisExt{(std::vector)binspTK0SExt, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec massXiAxis = {200, 1.3f, 1.4f, "#it{M}_{inv} [GeV/#it{c}^{2}]"}; + AxisSpec pTXiAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec binnedpTXiAxis{(std::vector)binspTXi, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec binnedpTXiAxisExt{(std::vector)binspTXiExt, "#it{p}_{T} (GeV/#it{c})"}; AxisSpec nSigmaPiAxis = {100, -10.0f, 10.0f, "N#sigma #pi"}; AxisSpec pTPiAxis = {50, 0.0f, 5.0f, "#it{p}_{T} (GeV/#it{c})"}; AxisSpec binnedpTPiAxis{(std::vector)binspTPi, "#it{p}_{T} (GeV/#it{c})"}; @@ -307,18 +356,22 @@ struct PhiStrangenessCorrelation { histos.add("phi/h3PhiData", "Invariant mass of Phi in Data", kTH3F, {binnedmultAxis, binnedpTPhiAxis, massPhiAxis}); histos.add("phiK0S/h6PhiK0SData", "Invariant mass of Phi vs Invariant mass of K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, massPhiAxis, massK0SAxis}); + histos.add("phiXi/h6PhiXiData", "Invariant mass of Phi vs Invariant mass of Xi in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTXiAxis, deltayAxis, massPhiAxis, massXiAxis}); histos.add("phiPi/h6PhiPiTPCData", "Invariant mass of Phi vs nSigmaTPC Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); histos.add("phiPi/h6PhiPiTOFData", "Invariant mass of Phi vs nSigmaTOF Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); histos.add("phiK0S/h6PhiK0SDataME", "Invariant mass of Phi vs Invariant mass of K0Short in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, massPhiAxis, massK0SAxis}); + histos.add("phiXi/h6PhiXiDataME", "Invariant mass of Phi vs Invariant mass of Xi in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTXiAxis, deltayAxis, massPhiAxis, massXiAxis}); histos.add("phiPi/h6PhiPiTPCDataME", "Invariant mass of Phi vs nSigmaTPC Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); histos.add("phiPi/h6PhiPiTOFDataME", "Invariant mass of Phi vs nSigmaTOF Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, massPhiAxis, nSigmaPiAxis}); for (const auto& label : PhiMassRegionLabels) { histos.add(fmt::format("phiK0S/h5PhiK0SData{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis}); + histos.add(fmt::format("phiXi/h5PhiXiData{}", label).c_str(), "Deltay vs deltaphi for Phi and Xi in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTXiAxis, deltayAxis, deltaphiAxis}); histos.add(fmt::format("phiPi/h5PhiPiData{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis}); histos.add(fmt::format("phiK0S/h5PhiK0SDataME{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis}); + histos.add(fmt::format("phiXi/h5PhiXiDataME{}", label).c_str(), "Deltay vs deltaphi for Phi and Xi in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTXiAxis, deltayAxis, deltaphiAxis}); histos.add(fmt::format("phiPi/h5PhiPiDataME{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis}); } @@ -337,6 +390,10 @@ struct PhiStrangenessCorrelation { histos.add("k0s/h3K0SMCGen", "K0S in MC Gen", kTH3F, {binnedmultAxis, binnedpTK0SAxis, yAxis}); histos.add("k0s/h4K0SMCGenAssocReco", "K0S in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTK0SAxis, yAxis}); + histos.add("xi/h4XiMCReco", "Xi in MC Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTXiAxis, yAxis}); + histos.add("xi/h3XiMCGen", "Xi in MC Gen", kTH3F, {binnedmultAxis, binnedpTXiAxis, yAxis}); + histos.add("xi/h4XiMCGenAssocReco", "Xi in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTXiAxis, yAxis}); + histos.add("pi/h4PiMCReco", "Pion in MC Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPiAxis, yAxis}); histos.add("pi/h3PiMCGen", "Pion in MC Gen", kTH3F, {binnedmultAxis, binnedpTPiAxis, yAxis}); histos.add("pi/h4PiMCGenAssocReco", "Pion in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPiAxis, yAxis}); @@ -344,9 +401,11 @@ struct PhiStrangenessCorrelation { histos.add("phi/h3PhiMCClosureGen", "Phi in MC Gen for MC Closure Test", kTH3F, {binnedmultAxis, binnedpTPhiAxisExt, yAxis}); histos.add("phiK0S/h5PhiK0SClosureMCGen", "Deltay vs deltaphi for Phi and K0Short in MCGen", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTK0SAxisExt, deltayAxis, deltaphiAxis}); + histos.add("phiXi/h5PhiXiClosureMCGen", "Deltay vs deltaphi for Phi and Xi in MCGen", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTXiAxisExt, deltayAxis, deltaphiAxis}); histos.add("phiPi/h5PhiPiClosureMCGen", "Deltay vs deltaphi for Phi and Pion in MCGen", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTPiAxisExt, deltayAxis, deltaphiAxis}); histos.add("phiK0S/h5PhiK0SClosureMCGenME", "Deltay vs deltaphi for Phi and K0Short in MCGen ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTK0SAxisExt, deltayAxis, deltaphiAxis}); + histos.add("phiXi/h5PhiXiClosureMCGenME", "Deltay vs deltaphi for Phi and Xi in MCGen ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTXiAxisExt, deltayAxis, deltaphiAxis}); histos.add("phiPi/h5PhiPiClosureMCGenME", "Deltay vs deltaphi for Phi and Pion in MCGen ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxisExt, binnedpTPiAxisExt, deltayAxis, deltaphiAxis}); // Load efficiency maps from CCDB @@ -387,33 +446,25 @@ struct PhiStrangenessCorrelation { void loadEfficiencyMaps() { // Load the Trigger (Phi) map if requested by analysis method - if (efficiencyConfigs.applyPhiEfficiency) + if (efficiencyConfigs.applyPhiEfficiency) { fetchSingleEfficiencyMapFromCCDB(effMapPhi, "Phi"); - - // Map the user configurations for the associated particles - bool doAssocCorrelations[kAssocPartSize] = {activeCorrelationConfigs.doK0SCorrelation, activeCorrelationConfigs.doXiCorrelation, activeCorrelationConfigs.doPionCorrelation}; + } // Only load the associated maps that are explicitly enabled for (size_t i = 0; i < kAssocPartSize; ++i) { - if (doAssocCorrelations[i]) + if (activeCorrelationTypes->at(i)) { fetchSingleEfficiencyMapFromCCDB(effMapsAssoc[i], AssocParticleLabels[i]); + } } } - /*void loadEfficiencyMapFromCCDB(ParticleOfInterest poi) - { - effMaps[poi] = std::shared_ptr(ccdb->get(fmt::format("{}{}", ccdbEfficiencyPath.value, particleOfInterestLabels[poi]))); - if (!effMaps[poi]) - LOG(fatal) << "Could not load efficiency map for " << particleOfInterestLabels[poi] << "!"; - LOG(info) << "Efficiency map for " << particleOfInterestLabels[poi] << " loaded from CCDB"; - }*/ - // Compute weight based on efficiencies template std::pair computeWeightAndError(const BoundEffMaps&... boundEffMaps) { - if (!efficiencyConfigs.applyEfficiency) + if (!efficiencyConfigs.applyEfficiency) { return {1.0f, 0.0f}; + } float effTot = 1.0f; float relErrSqSum = 0.0f; @@ -436,8 +487,9 @@ struct PhiStrangenessCorrelation { (processMap(boundEffMaps), ...); // Fold expression to process all bound efficiency maps - if (effTot <= 0.0f) + if (effTot <= 0.0f) { return {1.0f, 0.0f}; + } float weight = 1.0f / effTot; float weightErr = weight * std::sqrt(relErrSqSum); // Propagate relative error to the weight @@ -445,6 +497,13 @@ struct PhiStrangenessCorrelation { return {weight, weightErr}; } + template + std::pair computeAssocWeight(float multiplicity, const auto& phiCand, const auto& assocCand) + { + return computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapsAssoc[PartType], multiplicity, assocCand.pt(), assocCand.y())); + } + float getDeltaPhi(float phiTrigger, float phiAssociated) { return RecoDecay::constrainAngle(phiTrigger - phiAssociated, -o2::constants::math::PIHalf); @@ -452,13 +511,44 @@ struct PhiStrangenessCorrelation { int getCentBin(float multiplicity) { - if (multiplicity < binsMult->front() || multiplicity >= binsMult->back()) + if (multiplicity < binsMult->front() || multiplicity >= binsMult->back()) { return -1; + } auto it = std::upper_bound(binsMult->begin(), binsMult->end(), multiplicity); return std::distance(binsMult->begin(), it) - 1; } + template + auto getPreslice() + { + if constexpr (PartType == kK0S) { + return preslices.k0sMcRecoPerCollision; + } else if constexpr (PartType == kXi) { + return preslices.xiMcRecoPerCollision; + } else if constexpr (PartType == kPion) { + return preslices.pionTrackMcRecoPerCollision; + } else { + static_assert(PartType == kK0S || PartType == kXi || PartType == kPion, "Unsupported particle type in getPreslice"); + return preslices.k0sMcRecoPerCollision; + } + } + + template + float getMinGenPt() const + { + if constexpr (PartType == kK0S) { + return minPtMcGenConfigs.v0SettingMinPt.value; + } else if constexpr (PartType == kXi) { + return minPtMcGenConfigs.cascadeSettingMinPt.value; + } else if constexpr (PartType == kPion) { + return minPtMcGenConfigs.cMinPionPtcut.value; + } else { + static_assert(PartType == kK0S || PartType == kXi || PartType == kPion, "Unsupported particle type in getMinGenPt"); + return 0.f; + } + } + template void customFillHist(auto histId, const std::pair& weightPair, Args... coords) { @@ -503,11 +593,11 @@ struct PhiStrangenessCorrelation { auto [w, wErr] = weightPair; // Create an array of floats for the coordinates to pass to GetBin and SetBinContent - double coordArray[] = {static_cast(coords)...}; + std::array coordArray{static_cast(coords)...}; // Find the bin number. // The 'true' flag is mandatory: it allocates the bin in memory if it doesn't exist yet. - auto bin = hist->GetBin(coordArray, true); + auto bin = hist->GetBin(coordArray.data(), true); // Retrieve the previous content and the squared error (variance) double prevContent = hist->GetBinContent(bin); @@ -525,16 +615,291 @@ struct PhiStrangenessCorrelation { hist->SetBinError2(bin, newErr2); // Note: SetBinError2 takes the squared error directly } - template - void processPhiK0SPionSE(TCollision const& collision, TPhiCands const& phiCandidates, TK0SCands const& k0sReduced, TPionCands const& pionTracks) + template + bool isAssocValid(TAssoc const& assoc) + { + if constexpr (PartType == kK0S) { + const bool applyMassCut = (analysisMode == kDeltaYvsDeltaPhi) && k0sConfigs.selectK0sInSigRegion; + const auto& [minMass, maxMass] = k0sConfigs.rangeMK0sSignal.value; + return (!efficiencyConfigs.applyEfficiency || assoc.pt() < binspTK0S->back()) && (!applyMassCut || assoc.inMassRegion(minMass, maxMass)); + } else if constexpr (PartType == kXi) { + const bool applyMassCut = (analysisMode == kDeltaYvsDeltaPhi) && xiConfigs.selectXiInSigRegion; + const auto& [minMass, maxMass] = xiConfigs.rangeMXiSignal.value; + return (!efficiencyConfigs.applyEfficiency || assoc.pt() < binspTXi->back()) && (!applyMassCut || assoc.inMassRegion(minMass, maxMass)); + } else if constexpr (PartType == kPion) { + const bool applyNSigmaCut = (analysisMode == kDeltaYvsDeltaPhi) && pionConfigs.selectPionInSigRegion; + return (!efficiencyConfigs.applyEfficiency || assoc.pt() < binspTPi->back()) && + (!applyNSigmaCut || assoc.inNSigmaRegion(pionConfigs.pidTPCMax, pionConfigs.pidTOFMax)); + } else { + static_assert(PartType == kK0S || PartType == kXi || PartType == kPion, "Unsupported particle type in isAssocValid"); + return false; + } + } + + template + bool isGenSpeciesValid(const auto& mcParticle, bool allowBypass) const + { + return mcParticle.isPhysicalPrimary() && (allowBypass || mcParticle.pt() >= getMinGenPt()); + } + + template + void fillPhiAssocCorrelation(TPhiCand const& phiCand, TAssoc const& assoc, float multiplicity, const std::pair& weight) + { + const std::array, kPhiMassRegions> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband}; + + if (analysisMode == kMassvsMass) { + auto fillMass = [&](auto histID, auto assocVal) { + customFillTHn(histID, weight, multiplicity, phiCand.pt(), assoc.pt(), phiCand.y() - assoc.y(), phiCand.m(), assocVal); + }; + + if constexpr (PartType == kK0S) { + if constexpr (IsME) { + fillMass(HIST("phiK0S/h6PhiK0SDataME"), assoc.m()); + } else { + fillMass(HIST("phiK0S/h6PhiK0SData"), assoc.m()); + } + } else if constexpr (PartType == kXi) { + if constexpr (IsME) { + fillMass(HIST("phiXi/h6PhiXiDataME"), assoc.m()); + } else { + fillMass(HIST("phiXi/h6PhiXiData"), assoc.m()); + } + } else if constexpr (PartType == kPion) { + if constexpr (IsME) { + fillMass(HIST("phiPi/h6PhiPiTPCDataME"), assoc.nSigmaTPC()); + fillMass(HIST("phiPi/h6PhiPiTOFDataME"), assoc.nSigmaTOF()); + } else { + fillMass(HIST("phiPi/h6PhiPiTPCData"), assoc.nSigmaTPC()); + fillMass(HIST("phiPi/h6PhiPiTOFData"), assoc.nSigmaTOF()); + } + } + } else if (analysisMode == kDeltaYvsDeltaPhi) { + static_for<0, kPhiMassRegions - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; + if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) { + return; + } + + auto fillDelta = [&](auto histID) { + customFillTHn(histID, weight, multiplicity, phiCand.pt(), assoc.pt(), phiCand.y() - assoc.y(), getDeltaPhi(phiCand.phi(), assoc.phi())); + }; + + if constexpr (PartType == kK0S) { + if constexpr (IsME) { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiK0S/h5PhiK0SDataMESignal")); + else + fillDelta(HIST("phiK0S/h5PhiK0SDataMESideband")); + } else { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiK0S/h5PhiK0SDataSignal")); + else + fillDelta(HIST("phiK0S/h5PhiK0SDataSideband")); + } + } else if constexpr (PartType == kXi) { + if constexpr (IsME) { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiXi/h5PhiXiDataMESignal")); + else + fillDelta(HIST("phiXi/h5PhiXiDataMESideband")); + } else { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiXi/h5PhiXiDataSignal")); + else + fillDelta(HIST("phiXi/h5PhiXiDataSideband")); + } + } else if constexpr (PartType == kPion) { + if constexpr (IsME) { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiPi/h5PhiPiDataMESignal")); + else + fillDelta(HIST("phiPi/h5PhiPiDataMESideband")); + } else { + if constexpr (Idx == kSignalRegion) + fillDelta(HIST("phiPi/h5PhiPiDataSignal")); + else + fillDelta(HIST("phiPi/h5PhiPiDataSideband")); + } + } + + /*if constexpr (PartType == kK0S) { + if constexpr (IsME) { + auto t = std::make_tuple(HIST("phiK0S/h5PhiK0SDataMESignal"), HIST("phiK0S/h5PhiK0SDataMESideband")); + fillDelta(std::get(t)); + } else { + auto t = std::make_tuple(HIST("phiK0S/h5PhiK0SDataSignal"), HIST("phiK0S/h5PhiK0SDataSideband")); + fillDelta(std::get(t)); + } + } else if constexpr (PartType == kXi) { + if constexpr (IsME) { + auto t = std::make_tuple(HIST("phiXi/h5PhiXiDataMESignal"), HIST("phiXi/h5PhiXiDataMESideband")); + fillDelta(std::get(t)); + } else { + auto t = std::make_tuple(HIST("phiXi/h5PhiXiDataSignal"), HIST("phiXi/h5PhiXiDataSideband")); + fillDelta(std::get(t)); + } + } else if constexpr (PartType == kPion) { + if constexpr (IsME) { + auto t = std::make_tuple(HIST("phiPi/h5PhiPiDataMESignal"), HIST("phiPi/h5PhiPiDataMESideband")); + fillDelta(std::get(t)); + } else { + auto t = std::make_tuple(HIST("phiPi/h5PhiPiDataSignal"), HIST("phiPi/h5PhiPiDataSideband")); + fillDelta(std::get(t)); + } + }*/ + }); + } + } + + template + void processOneAssocPair(TPhiCand const& phiCand, TAssoc const& assoc, float multiplicity) + { + if (!isAssocValid(assoc)) { + return; + } + + auto weight = computeAssocWeight(multiplicity, phiCand, assoc); + fillPhiAssocCorrelation(phiCand, assoc, multiplicity, weight); + } + + template + void processOneAssocSpeciesSE(TPhiCand const& phiCand, TAssocCands const& assocCands, float multiplicity) + { + if (!activeCorrelationTypes->at(PartType)) { + return; + } + + for (const auto& assoc : assocCands) { + processOneAssocPair(phiCand, assoc, multiplicity); + } + } + + template + void processPhiAssocSE(TCollision const& collision, TPhiCands const& phiCandidates, TAssocInputs const&... assocInputs) { float multiplicity = collision.centFT0M(); - const std::array, 2> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband}; + for (const auto& phiCand : phiCandidates) { + if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt() >= binspTPhi->back()) { + continue; + } + + auto weightPhi = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); + customFillHist(HIST("phi/h3PhiData"), weightPhi, multiplicity, phiCand.pt(), phiCand.m()); + + // Fold expression: elaborate the correlation for each associated particle type in the parameter pack + (processOneAssocSpeciesSE(phiCand, *assocInputs.candidates, multiplicity), ...); + } + } + + void processPhiAssocSEDataLike(SelCollisions::iterator const& collision, + aod::PhimesonCandidatesData const& phiCandidates, + aod::K0sReducedCandidatesData const& k0sReduced, + aod::XiReducedCandidatesData const& xiReduced, + aod::PionTracksData const& pionTracks) + { + processPhiAssocSE(collision, phiCandidates, + makeAssocInput(k0sReduced), + makeAssocInput(xiReduced), + makeAssocInput(pionTracks)); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiAssocSEDataLike, "Process function for Phi-Associated 2D Correlations in Data or MC w/o PDG SE", true); - bool doAssocCorrelations[kAssocPartSize] = {activeCorrelationConfigs.doK0SCorrelation, - activeCorrelationConfigs.doXiCorrelation, - activeCorrelationConfigs.doPionCorrelation}; + void processPhiAssocSEMCWithPDG(SimCollisions::iterator const& collision, + aod::PhimesonCandidatesMcReco const& phiCandidates, + aod::K0sReducedCandidatesMcReco const& k0sReduced, + aod::XiReducedCandidatesMcReco const& xiReduced, + aod::PionTracksMcReco const& pionTracks) + { + processPhiAssocSE(collision, phiCandidates, + makeAssocInput(k0sReduced), + makeAssocInput(xiReduced), + makeAssocInput(pionTracks)); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiAssocSEMCWithPDG, "Process function for Phi-Associated 2D Correlations in MC with PDG SE", true); + + template + void processPhiAssocME(TCollisions const& collisions, TPhiCands const& phiCandidates, TAssocCands const& assocCands) + { + auto tuplePhiAssoc = std::make_tuple(phiCandidates, assocCands); + Pair pairPhiAssoc{binningOnVertexAndCent, cfgNoMixedEvents, -1, collisions, tuplePhiAssoc, &cache}; + + for (const auto& [c1, phiCands, c2, assocRed] : pairPhiAssoc) { + float multiplicity = c1.centFT0M(); + + for (const auto& [phiCand, assoc] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(phiCands, assocRed))) { + if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt() >= binspTPhi->back()) { + continue; + } + + processOneAssocPair(phiCand, assoc, multiplicity); + } + } + } + + void processPhiK0SMEDataLike(SelCollisions const& collisions, + aod::PhimesonCandidatesData const& phiCandidates, + aod::K0sReducedCandidatesData const& k0sReduced) + { + processPhiAssocME(collisions, phiCandidates, k0sReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiK0SMEDataLike, "Process function for Phi-K0S 2D Correlations in Data or MC w/o PDG ME", true); + + void processPhiK0SMEMCWithPDG(SimCollisions const& collisions, + aod::PhimesonCandidatesMcReco const& phiCandidates, + aod::K0sReducedCandidatesMcReco const& k0sReduced) + { + processPhiAssocME(collisions, phiCandidates, k0sReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiK0SMEMCWithPDG, "Process function for Phi-K0S 2D Correlations in MC with PDG ME", true); + + void processPhiXiMEDataLike(SelCollisions const& collisions, + aod::PhimesonCandidatesData const& phiCandidates, + aod::XiReducedCandidatesData const& xiReduced) + { + processPhiAssocME(collisions, phiCandidates, xiReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiXiMEDataLike, "Process function for Phi-Xi 2D Correlations in Data or MC w/o PDG ME", true); + + void processPhiXiMEMCWithPDG(SimCollisions const& collisions, + aod::PhimesonCandidatesMcReco const& phiCandidates, + aod::XiReducedCandidatesMcReco const& xiReduced) + { + processPhiAssocME(collisions, phiCandidates, xiReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiXiMEMCWithPDG, "Process function for Phi-Xi 2D Correlations in MC with PDG ME", true); + + void processPhiPionMEDataLike(SelCollisions const& collisions, + aod::PhimesonCandidatesData const& phiCandidates, + aod::PionTracksData const& pionTracks) + { + processPhiAssocME(collisions, phiCandidates, pionTracks); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiPionMEDataLike, "Process function for Phi-Pion 2D Correlations in Data or MC w/o PDG ME", true); + + void processPhiPionMEMCWithPDG(SimCollisions const& collisions, + aod::PhimesonCandidatesMcReco const& phiCandidates, + aod::PionTracksMcReco const& pionTracks) + { + processPhiAssocME(collisions, phiCandidates, pionTracks); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiPionMEMCWithPDG, "Process function for Phi-Pion 2D Correlations in MC with PDG ME", true); + + /*template + void processPhiAssocSE(TCollision const& collision, TPhiCands const& phiCandidates, TK0SCands const& k0sReduced, TXiCands const& xiReduced, TPionCands const& pionTracks) + { + float multiplicity = collision.centFT0M(); + + const std::array, 2> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband}; const bool applyK0sMassCut = (analysisMode == kDeltaYvsDeltaPhi) && k0sConfigs.selectK0sInSigRegion; const auto& [minMassK0s, maxMassK0s] = k0sConfigs.rangeMK0sSignal.value; @@ -542,6 +907,12 @@ struct PhiStrangenessCorrelation { return (!efficiencyConfigs.applyEfficiency || k0s.pt() < binspTK0S->back()) && (!applyK0sMassCut || k0s.inMassRegion(minMassK0s, maxMassK0s)); }; + const bool applyXiMassCut = (analysisMode == kDeltaYvsDeltaPhi) && xiConfigs.selectXiInSigRegion; + const auto& [minMassXi, maxMassXi] = xiConfigs.rangeMXiSignal.value; + auto isXiValid = [&](const auto& xi) { + return (!efficiencyConfigs.applyEfficiency || xi.pt() < binspTXi->back()) && (!applyXiMassCut || xi.inMassRegion(minMassXi, maxMassXi)); + }; + const bool applyPionNSigmaCut = (analysisMode == kDeltaYvsDeltaPhi) && pionConfigs.selectPionInSigRegion; const float& pidTPCMax = pionConfigs.pidTPCMax; const float& pidTOFMax = pionConfigs.pidTOFMax; @@ -555,44 +926,55 @@ struct PhiStrangenessCorrelation { if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt() >= binspTPhi->back()) continue; - // float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); auto weightPhi = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); // histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi); customFillHist(HIST("phi/h3PhiData"), weightPhi, multiplicity, phiCand.pt(), phiCand.m()); - auto processCorrelations = [&](auto fillK0S, auto fillPion) { - if (doAssocCorrelations[kK0S]) { + auto processCorrelations = [&](auto fillK0S, auto fillXi, auto fillPion) { + if (activeCorrelationTypes->at(kK0S)) { // Loop over all reduced K0S candidates for (const auto& k0s : k0sReduced) { if (!isK0sValid(k0s)) continue; - // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); - auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); - fillK0S(k0s, weightPhiK0S); + + fillK0S(k0s, computeAssocWeight(multiplicity, phiCand, k0s)); } } - if (doAssocCorrelations[kPion]) { + if (activeCorrelationTypes->at(kXi)) { + // Loop over all reduced Xi candidates + for (const auto& xi : xiReduced) { + if (!isXiValid(xi)) + continue; + + // auto weightPhiXi = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kXi], multiplicity, xi.pt(), xi.y())); + + fillXi(xi, computeAssocWeight(multiplicity, phiCand, xi)); + } + } + + if (activeCorrelationTypes->at(kPion)) { // Loop over all primary pion candidates for (const auto& pionTrack : pionTracks) { if (!isPionValid(pionTrack)) continue; - // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y())); - auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y())); - fillPion(pionTrack, weightPhiPion); + + fillPion(pionTrack, computeAssocWeight(multiplicity, phiCand, pionTrack)); } } }; if (analysisMode == kMassvsMass) { auto k0sHistID = HIST("phiK0S/h6PhiK0SData"); + auto xiHistID = HIST("phiXi/h6PhiXiData"); auto piTPCHistID = HIST("phiPi/h6PhiPiTPCData"); auto piTOFHistID = HIST("phiPi/h6PhiPiTOFData"); @@ -602,6 +984,10 @@ struct PhiStrangenessCorrelation { [&](const auto& k0s, const std::pair& w) { customFillTHn(k0sHistID, w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m()); }, + //[&](const auto& xi, float w) { + [&](const auto& xi, const std::pair& w) { + customFillTHn(xiHistID, w, multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), phiCand.m(), xi.m()); + }, //[&](const auto& pion, float w) { // histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w); // histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w); @@ -611,9 +997,10 @@ struct PhiStrangenessCorrelation { }); } else if (analysisMode == kDeltaYvsDeltaPhi) { auto k0sHistID = std::make_tuple(HIST("phiK0S/h5PhiK0SDataSignal"), HIST("phiK0S/h5PhiK0SDataSideband")); + auto xiHistID = std::make_tuple(HIST("phiXi/h5PhiXiDataSignal"), HIST("phiXi/h5PhiXiDataSideband")); auto piHistID = std::make_tuple(HIST("phiPi/h5PhiPiDataSignal"), HIST("phiPi/h5PhiPiDataSideband")); - static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + static_for<0, kPhiMassRegions - 1>([&](auto i_idx) { constexpr unsigned int Idx = i_idx.value; const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; @@ -629,6 +1016,11 @@ struct PhiStrangenessCorrelation { [&](const auto& k0s, const std::pair& w) { customFillTHn(std::get(k0sHistID), w, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi())); }, + //[&](const auto& xi, float w) { + // histos.fill(std::get(xiHistID), multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), getDeltaPhi(phiCand.phi(), xi.phi()), w); + [&](const auto& xi, const std::pair& w) { + customFillTHn(std::get(xiHistID), w, multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), getDeltaPhi(phiCand.phi(), xi.phi())); + }, //[&](const auto& pion, float w) { // histos.fill(std::get(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w); [&](const auto& pion, const std::pair& w) { @@ -636,54 +1028,22 @@ struct PhiStrangenessCorrelation { }); }); } - - /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int Idx = i_idx.value; - - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; - if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) - return; - - // Loop over all reduced K0S candidates - for (const auto& k0s : k0sReduced) { - if (k0sConfigs.selectK0sInSigRegion) { - const auto& [minMassK0s, maxMassK0s] = k0sConfigs.rangeMK0sSignal.value; - if (!k0s.inMassRegion(minMassK0s, maxMassK0s)) - continue; - } - - float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMaps[K0S], multiplicity, k0s.pt(), k0s.y())); - - histos.fill(HIST("phiK0S/h5PhiK0SData") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); - } - - // Loop over all primary pion candidates - for (const auto& pionTrack : pionTracks) { - float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMaps[Pion], multiplicity, pionTrack.pt(), pionTrack.y())); - - histos.fill(HIST("phiPi/h5PhiPiData") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), pionTrack.pt(), phiCand.y() - pionTrack.y(), getDeltaPhi(phiCand.phi(), pionTrack.phi()), weightPhiPion); - } - });*/ } } - // PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionSE, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in SE", true); - - void processPhiK0SPionSEDataLike(SelCollisions::iterator const& collision, aod::PhimesonCandidatesData const& phiCandidates, aod::K0sReducedCandidatesData const& k0sReduced, aod::PionTracksData const& pionTracks) + void processPhiAssocSEDataLike(SelCollisions::iterator const& collision, aod::PhimesonCandidatesData const& phiCandidates, aod::K0sReducedCandidatesData const& k0sReduced, aod::XiReducedCandidatesData const& xiReduced, aod::PionTracksData const& pionTracks) { - processPhiK0SPionSE(collision, phiCandidates, k0sReduced, pionTracks); + processPhiAssocSE(collision, phiCandidates, k0sReduced, xiReduced, pionTracks); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionSEDataLike, "Process function for Phi-K0S and Phi-Pion 2D Correlations in Data or MC w/o PDG SE", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiAssocSEDataLike, "Process function for Phi-Associated 2D Correlations in Data or MC w/o PDG SE", true); - void processPhiK0SPionSEMCWithPDG(SimCollisions::iterator const& collision, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::K0sReducedCandidatesMcReco const& k0sReduced, aod::PionTracksMcReco const& pionTracks) + void processPhiAssocSEMCWithPDG(SimCollisions::iterator const& collision, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::K0sReducedCandidatesMcReco const& k0sReduced, aod::XiReducedCandidatesMcReco const& xiReduced, aod::PionTracksMcReco const& pionTracks) { - processPhiK0SPionSE(collision, phiCandidates, k0sReduced, pionTracks); + processPhiAssocSE(collision, phiCandidates, k0sReduced, xiReduced, pionTracks); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionSEMCWithPDG, "Process function for Phi-K0S and Phi-Pion 2D Correlations in MC with PDG SE", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiAssocSEMCWithPDG, "Process function for Phi-Associated 2D Correlations in MC with PDG SE", true); template void processPhiK0SME(TCollisions const& collisions, TPhiCands const& phiCandidates, TK0SCands const& k0sReduced) @@ -711,11 +1071,9 @@ struct PhiStrangenessCorrelation { continue; auto processCorrelations = [&](auto fillK0S) { - // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); - auto weightPhiK0S = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y())); - fillK0S(k0s, weightPhiK0S); + fillK0S(k0s, computeAssocWeight(multiplicity, phiCand, k0s)); }; if (analysisMode == kMassvsMass) { @@ -730,7 +1088,7 @@ struct PhiStrangenessCorrelation { } else if (analysisMode == kDeltaYvsDeltaPhi) { auto k0sHistID = std::make_tuple(HIST("phiK0S/h5PhiK0SDataMESignal"), HIST("phiK0S/h5PhiK0SDataMESideband")); - static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + static_for<0, kPhiMassRegions - 1>([&](auto i_idx) { constexpr unsigned int Idx = i_idx.value; const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; @@ -745,44 +1103,101 @@ struct PhiStrangenessCorrelation { }); }); } - - /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int Idx = i_idx.value; - - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; - if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) - return; - - if (k0sConfigs.selectK0sInSigRegion) { - const auto& [minMassK0s, maxMassK0s] = k0sConfigs.rangeMK0sSignal.value; - if (!k0s.inMassRegion(minMassK0s, maxMassK0s)) - return; - } - - float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMaps[K0S], multiplicity, k0s.pt(), k0s.y())); - - histos.fill(HIST("phiK0S/h5PhiK0SDataME") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), weightPhiK0S); - });*/ } } } - // PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SME, "Process function for Phi-K0S and Deltay and Deltaphi 2D Correlations in Data ME", false); + // PROCESS_SWITCH(PhiStrangeCorrelation, processPhiK0SME, "Process function for Phi-K0S and Deltay and Deltaphi 2D Correlations in Data ME", false); void processPhiK0SMEDataLike(SelCollisions const& collisions, aod::PhimesonCandidatesData const& phiCandidates, aod::K0sReducedCandidatesData const& k0sReduced) { processPhiK0SME(collisions, phiCandidates, k0sReduced); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SMEDataLike, "Process function for Phi-K0S 2D Correlations in Data or MC w/o PDG ME", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiK0SMEDataLike, "Process function for Phi-K0S 2D Correlations in Data or MC w/o PDG ME", true); void processPhiK0SMEMCWithPDG(SimCollisions const& collisions, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::K0sReducedCandidatesMcReco const& k0sReduced) { processPhiK0SME(collisions, phiCandidates, k0sReduced); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SMEMCWithPDG, "Process function for Phi-K0S 2D Correlations in MC with PDG ME", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiK0SMEMCWithPDG, "Process function for Phi-K0S 2D Correlations in MC with PDG ME", true); + + template + void processPhiXiME(TCollisions const& collisions, TPhiCands const& phiCandidates, TXiCands const& xiReduced) + { + const std::array, 2> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband}; + + const bool applyXiMassCut = (analysisMode == kDeltaYvsDeltaPhi) && xiConfigs.selectXiInSigRegion; + const auto& [minMassXi, maxMassXi] = xiConfigs.rangeMXiSignal.value; + + auto isXiValid = [&](const auto& xi) { + return (!efficiencyConfigs.applyEfficiency || xi.pt() < binspTXi->back()) && (!applyXiMassCut || xi.inMassRegion(minMassXi, maxMassXi)); + }; + + auto tuplePhiXi = std::make_tuple(phiCandidates, xiReduced); + Pair pairPhiXi{binningOnVertexAndCent, cfgNoMixedEvents, -1, collisions, tuplePhiXi, &cache}; + + for (const auto& [c1, phiCands, c2, xiRed] : pairPhiXi) { + + float multiplicity = c1.centFT0M(); + + for (const auto& [phiCand, xi] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(phiCands, xiRed))) { + if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt() >= binspTPhi->back()) + continue; + if (!isXiValid(xi)) + continue; + + auto processCorrelations = [&](auto fillXi) { + // auto weightPhiXi = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // BoundEfficiencyMap(effMapsAssoc[kXi], multiplicity, xi.pt(), xi.y())); + fillXi(xi, computeAssocWeight(multiplicity, phiCand, xi)); + }; + + if (analysisMode == kMassvsMass) { + auto xiHistID = HIST("phiXi/h6PhiXiDataME"); + + processCorrelations( + //[&](const auto& xi, float w) { + // histos.fill(xiHistID, multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), phiCand.m(), xi.m(), w); + [&](const auto& xi, const std::pair& w) { + customFillTHn(xiHistID, w, multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), phiCand.m(), xi.m()); + }); + } else if (analysisMode == kDeltaYvsDeltaPhi) { + auto xiHistID = std::make_tuple(HIST("phiXi/h5PhiXiDataMESignal"), HIST("phiXi/h5PhiXiDataMESideband")); + + static_for<0, kPhiMassRegions - 1>([&](auto i_idx) { + constexpr unsigned int Idx = i_idx.value; + + const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; + if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) + return; + + processCorrelations( + //[&](const auto& xi, float w) { + // histos.fill(std::get(xiHistID), multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), getDeltaPhi(phiCand.phi(), xi.phi()), w); + [&](const auto& xi, const std::pair& w) { + customFillTHn(std::get(xiHistID), w, multiplicity, phiCand.pt(), xi.pt(), phiCand.y() - xi.y(), getDeltaPhi(phiCand.phi(), xi.phi())); + }); + }); + } + } + } + } + + void processPhiXiMEDataLike(SelCollisions const& collisions, aod::PhimesonCandidatesData const& phiCandidates, aod::XiReducedCandidatesData const& xiReduced) + { + processPhiXiME(collisions, phiCandidates, xiReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiXiMEDataLike, "Process function for Phi-Xi 2D Correlations in Data or MC w/o PDG ME", true); + + void processPhiXiMEMCWithPDG(SimCollisions const& collisions, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::XiReducedCandidatesMcReco const& xiReduced) + { + processPhiXiME(collisions, phiCandidates, xiReduced); + } + + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiXiMEMCWithPDG, "Process function for Phi-Xi 2D Correlations in MC with PDG ME", true); template void processPhiPionME(TCollisions const& collisions, TPhiCands const& phiCandidates, TPionCands const& pionTracks) @@ -812,11 +1227,9 @@ struct PhiStrangenessCorrelation { continue; auto processCorrelations = [&](auto fillPion) { - // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + // auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y())); - auto weightPhiPion = computeWeightAndError(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y())); - fillPion(piTrack, weightPhiPion); + fillPion(piTrack, computeAssocWeight(multiplicity, phiCand, piTrack)); }; if (analysisMode == kMassvsMass) { @@ -834,7 +1247,7 @@ struct PhiStrangenessCorrelation { } else if (analysisMode == kDeltaYvsDeltaPhi) { auto piHistID = std::make_tuple(HIST("phiPi/h5PhiPiDataMESignal"), HIST("phiPi/h5PhiPiDataMESideband")); - static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { + static_for<0, kPhiMassRegions - 1>([&](auto i_idx) { constexpr unsigned int Idx = i_idx.value; const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; @@ -849,47 +1262,45 @@ struct PhiStrangenessCorrelation { }); }); } - - /*static_for<0, PhiMassRegionLabels.size() - 1>([&](auto i_idx) { - constexpr unsigned int Idx = i_idx.value; - - const auto& [minMassPhi, maxMassPhi] = phiMassRegions[Idx]; - if (!phiCand.inMassRegion(minMassPhi, maxMassPhi)) - return; - - float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), - BoundEfficiencyMap(effMaps[Pion], multiplicity, piTrack.pt(), piTrack.y())); - - histos.fill(HIST("phiPi/h5PhiPiDataME") + HIST(PhiMassRegionLabels[Idx]), multiplicity, phiCand.pt(), piTrack.pt(), phiCand.y() - piTrack.y(), getDeltaPhi(phiCand.phi(), piTrack.phi()), weightPhiPion); - });*/ } } } - // PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiPionME, "Process function for Phi-Pion Deltay and Deltaphi 2D Correlations in Data ME", false); + // PROCESS_SWITCH(PhiStrangeCorrelation, processPhiPionME, "Process function for Phi-Pion Deltay and Deltaphi 2D Correlations in Data ME", false); void processPhiPionMEDataLike(SelCollisions const& collisions, aod::PhimesonCandidatesData const& phiCandidates, aod::PionTracksData const& pionTracks) { processPhiPionME(collisions, phiCandidates, pionTracks); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiPionMEDataLike, "Process function for Phi-Pion 2D Correlations in Data or MC w/o PDG ME", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiPionMEDataLike, "Process function for Phi-Pion 2D Correlations in Data or MC w/o PDG ME", true); void processPhiPionMEMCWithPDG(SimCollisions const& collisions, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::PionTracksMcReco const& pionTracks) { processPhiPionME(collisions, phiCandidates, pionTracks); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiPionMEMCWithPDG, "Process function for Phi-Pion 2D Correlations in MC with PDG ME", true); + PROCESS_SWITCH(PhiStrangeCorrelation, processPhiPionMEMCWithPDG, "Process function for Phi-Pion 2D Correlations in MC with PDG ME", true);*/ - void processParticleEfficiency(MCCollisions const& mcCollisions, SimCollisions const& collisions, aod::PhimesonCandidatesMcReco const& phiCandidates, aod::K0sReducedCandidatesMcReco const& k0sReduced, aod::PionTracksMcReco const& pionTracks, aod::McParticles const& mcParticles) + void processParticleEfficiency(MCCollisions const& mcCollisions, + SimCollisions const& collisions, + aod::PhimesonCandidatesMcReco const& phiCandidates, + aod::K0sReducedCandidatesMcReco const& k0sReduced, + aod::XiReducedCandidatesMcReco const& xiReduced, + aod::PionTracksMcReco const& pionTracks, aod::McParticles const& mcParticles) { - const bool applyK0sMassCut = (analysisMode == kDeltaYvsDeltaPhi) && k0sConfigs.selectK0sInSigRegion; + /*const bool applyK0sMassCut = (analysisMode == kDeltaYvsDeltaPhi) && k0sConfigs.selectK0sInSigRegion; const auto& [minMassK0s, maxMassK0s] = k0sConfigs.rangeMK0sSignal.value; auto isK0sValid = [&](const auto& k0s) { return !applyK0sMassCut || k0s.inMassRegion(minMassK0s, maxMassK0s); }; + const bool applyXiMassCut = (analysisMode == kDeltaYvsDeltaPhi) && xiConfigs.selectXiInSigRegion; + const auto& [minMassXi, maxMassXi] = xiConfigs.rangeMXiSignal.value; + auto isXiValid = [&](const auto& xi) { + return !applyXiMassCut || xi.inMassRegion(minMassXi, maxMassXi); + }; + const bool applyPionNSigmaCut = (analysisMode == kDeltaYvsDeltaPhi) && pionConfigs.selectPionInSigRegion; const float& pidTPCMax = pionConfigs.pidTPCMax; const float& pidTOFMax = pionConfigs.pidTOFMax; @@ -897,14 +1308,15 @@ struct PhiStrangenessCorrelation { auto isPionValid = [&](const auto& pion) { return !applyPionNSigmaCut || pion.inNSigmaRegion(pidTPCMax, pidTOFMax); - }; + };*/ std::unordered_map> collsGrouped; collsGrouped.reserve(mcCollisions.size()); for (const auto& collision : collisions) { - if (!collision.has_mcCollision()) + if (!collision.has_mcCollision()) { continue; + } const auto& mcCollision = collision.mcCollision_as(); collsGrouped[mcCollision.globalIndex()].push_back(collision.globalIndex()); } @@ -936,22 +1348,50 @@ struct PhiStrangenessCorrelation { } } - const auto k0sThisColl = k0sReduced.sliceBy(preslices.k0sMcRecoPerCollision, collision.globalIndex()); + auto fillRecoAssocSpecies = [&](TAssocInputs const& assocCands, auto histoKey) { + if (!activeCorrelationTypes->at(TAssocInputs::Type)) { + return; + } + + const auto assocThisColl = assocCands.candidates->sliceBy(getPreslice(), collision.globalIndex()); + + for (const auto& assoc : assocThisColl) { + if (!isAssocValid(assoc)) { + continue; + } + + histos.fill(histoKey, collision.posZ(), mcCollision.centFT0M(), assoc.pt(), assoc.y()); + } + }; + + fillRecoAssocSpecies(makeAssocInput(k0sReduced), HIST("k0s/h4K0SMCReco")); + fillRecoAssocSpecies(makeAssocInput(xiReduced), HIST("xi/h4XiMCReco")); + fillRecoAssocSpecies(makeAssocInput(pionTracks), HIST("pi/h4PiMCReco")); + + /*const auto k0sThisColl = k0sReduced.sliceBy(preslices.k0sMcRecoPerCollision, collision.globalIndex()); + const auto xiThisColl = xiReduced.sliceBy(preslices.xiMcRecoPerCollision, collision.globalIndex()); const auto pionTracksThisColl = pionTracks.sliceBy(preslices.pionTrackMcRecoPerCollision, collision.globalIndex()); for (const auto& k0s : k0sThisColl) { - if (!isK0sValid(k0s)) + if (!isAssocValid(k0s)) continue; histos.fill(HIST("k0s/h4K0SMCReco"), collision.posZ(), mcCollision.centFT0M(), k0s.pt(), k0s.y()); } + for (const auto& xi : xiThisColl) { + if (!isAssocValid(xi)) + continue; + + histos.fill(HIST("xi/h4XiMCReco"), collision.posZ(), mcCollision.centFT0M(), xi.pt(), xi.y()); + } + for (const auto& pionTrack : pionTracksThisColl) { - if (!isPionValid(pionTrack)) + if (!isAssocValid(pionTrack)) continue; histos.fill(HIST("pi/h4PiMCReco"), collision.posZ(), mcCollision.centFT0M(), pionTrack.pt(), pionTrack.y()); - } + }*/ numberAssocColls++; } @@ -962,7 +1402,6 @@ struct PhiStrangenessCorrelation { const bool hasAssoc = (numberAssocColls > 0); const float zVtxRef = hasAssoc ? zVtxs[0] : 0.0f; - //////TOBECHANGED////// if (hasAssoc) { if (zVtxs.size() > 1) { for (size_t i = 1; i < zVtxs.size(); ++i) { @@ -973,36 +1412,59 @@ struct PhiStrangenessCorrelation { histos.fill(HIST("event/hGenMCAssocRecoMultiplicityPercent"), mcCollision.centFT0M()); histos.fill(HIST("event/h2GenMCAssocRecoVertexZvsMult"), zVtxRef, mcCollision.centFT0M()); } - /////////////////////// const auto mcParticlesThisMcColl = mcParticles.sliceBy(preslices.mcPartPerMcCollision, mcCollision.globalIndex()); - auto inYAcceptance = [&](const auto& mcParticle) { - return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance; - }; + for (const auto& mcParticle : mcParticlesThisMcColl) { + auto inYAcceptance = [&]() { + return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance; + }; - auto fillGenHistos = [&](auto h3Key, auto h4Key, const auto& mcParticle) { - histos.fill(h3Key, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); - if (hasAssoc) - histos.fill(h4Key, zVtxRef, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); - }; + auto fillGenHistos = [&](auto h3Key, auto h4Key) { + histos.fill(h3Key, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); + if (hasAssoc) { + histos.fill(h4Key, zVtxRef, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); + } + }; + + auto fillGenAssocSpecies = [&](auto h3Key, auto h4Key) { + if (!activeCorrelationTypes->at(PartType) || !isGenSpeciesValid(mcParticle, false)) { + return; + } - for (const auto& mcParticle : mcParticlesThisMcColl /*| std::views::filter(inYAcceptance)*/) { - if (!inYAcceptance(mcParticle)) + fillGenHistos(h3Key, h4Key); + }; + + if (!inYAcceptance()) { continue; + } switch (std::abs(mcParticle.pdgCode())) { case o2::constants::physics::Pdg::kPhi: - if (eventSelectionType == 0 && mcParticle.pt() >= minPtMcGenConfigs.minPhiPt) - fillGenHistos(HIST("phi/h3PhiMCGen"), HIST("phi/h4PhiMCGenAssocReco"), mcParticle); + if (eventSelectionType == 0 && mcParticle.pt() >= minPtMcGenConfigs.minPhiPt) { + fillGenHistos(HIST("phi/h3PhiMCGen"), HIST("phi/h4PhiMCGenAssocReco")); + } break; - case PDG_t::kK0Short: + /*case PDG_t::kK0Short: if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.v0SettingMinPt) - fillGenHistos(HIST("k0s/h3K0SMCGen"), HIST("k0s/h4K0SMCGenAssocReco"), mcParticle); + fillGenHistos(HIST("k0s/h3K0SMCGen"), HIST("k0s/h4K0SMCGenAssocReco")); + break; + case PDG_t::kXiMinus: + if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.cascadeSettingMinPt) + fillGenHistos(HIST("xi/h3XiMCGen"), HIST("xi/h4XiMCGenAssocReco")); break; case PDG_t::kPiPlus: if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.cMinPionPtcut) - fillGenHistos(HIST("pi/h3PiMCGen"), HIST("pi/h4PiMCGenAssocReco"), mcParticle); + fillGenHistos(HIST("pi/h3PiMCGen"), HIST("pi/h4PiMCGenAssocReco")); + break;*/ + case getPdgCode(): + fillGenAssocSpecies.template operator()(HIST("k0s/h3K0SMCGen"), HIST("k0s/h4K0SMCGenAssocReco")); + break; + case getPdgCode(): + fillGenAssocSpecies.template operator()(HIST("xi/h3XiMCGen"), HIST("xi/h4XiMCGenAssocReco")); + break; + case getPdgCode(): + fillGenAssocSpecies.template operator()(HIST("pi/h3PiMCGen"), HIST("pi/h4PiMCGenAssocReco")); break; default: break; @@ -1011,7 +1473,7 @@ struct PhiStrangenessCorrelation { } } - PROCESS_SWITCH(PhiStrangenessCorrelation, processParticleEfficiency, "Process function for Efficiency Computation for Particles of Interest", false); + PROCESS_SWITCH(PhiStrangeCorrelation, processParticleEfficiency, "Process function for Efficiency Computation for Particles of Interest", false); /*void processMCGenClosureSE(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) { @@ -1066,7 +1528,7 @@ struct PhiStrangenessCorrelation { } } - PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosureSE, "Process function for MC Gen Closure Test in SE", false); + PROCESS_SWITCH(PhiStrangeCorrelation, processMCGenClosureSE, "Process function for MC Gen Closure Test in SE", false); void processMCGenClosureME(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) { @@ -1131,7 +1593,7 @@ struct PhiStrangenessCorrelation { eventBuffer[multBin].pop_back(); } - PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosureME, "Process function for MC Gen Closure Test in ME", false);*/ + PROCESS_SWITCH(PhiStrangeCorrelation, processMCGenClosureME, "Process function for MC Gen Closure Test in ME", false);*/ void processMCGenClosure(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) { @@ -1142,54 +1604,80 @@ struct PhiStrangenessCorrelation { std::vector xiParticles; std::vector pionParticles; - auto inYAcceptance = [&](const auto& mcParticle) { - return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance; - }; - // Preliminary loop to fill vectors of particles of interest for the current event, applying pt and y cuts for (const auto& mcParticle : mcParticles) { - if (!inYAcceptance(mcParticle)) + auto inYAcceptance = [&]() { + return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance; + }; + + auto fillPartCollection = [&](auto& collection) { + if (!activeCorrelationTypes->at(PartType) || !isGenSpeciesValid(mcParticle, minPtMcGenConfigs.bypassPtCut.value)) { + return; + } + + collection.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi()); + }; + + if (!inYAcceptance()) { continue; + } switch (std::abs(mcParticle.pdgCode())) { case o2::constants::physics::Pdg::kPhi: - if (eventSelectionType == 0 && (minPtMcGenConfigs.bypassPtCut || mcParticle.pt() >= minPtMcGenConfigs.minPhiPt)) + if (eventSelectionType == 0 && (minPtMcGenConfigs.bypassPtCut || mcParticle.pt() >= minPtMcGenConfigs.minPhiPt)) { phiParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi()); + } break; - case PDG_t::kK0Short: + /*case PDG_t::kK0Short: if (mcParticle.isPhysicalPrimary() && (minPtMcGenConfigs.bypassPtCut || mcParticle.pt() >= minPtMcGenConfigs.v0SettingMinPt)) k0sParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi()); break; + case PDG_t::kXiMinus: + if (mcParticle.isPhysicalPrimary() && (minPtMcGenConfigs.bypassPtCut || mcParticle.pt() >= minPtMcGenConfigs.cascadeSettingMinPt)) + xiParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi()); + break; case PDG_t::kPiPlus: if (mcParticle.isPhysicalPrimary() && (minPtMcGenConfigs.bypassPtCut || mcParticle.pt() >= minPtMcGenConfigs.cMinPionPtcut)) pionParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi()); + break;*/ + case getPdgCode(): + fillPartCollection.template operator()(k0sParticles); + break; + case getPdgCode(): + fillPartCollection.template operator()(xiParticles); + break; + case getPdgCode(): + fillPartCollection.template operator()(pionParticles); break; default: break; } } - if (phiParticles.empty() && k0sParticles.empty() && pionParticles.empty()) + const bool skipK0s = !activeCorrelationTypes->at(kK0S) || k0sParticles.empty(); + const bool skipXi = !activeCorrelationTypes->at(kXi) || xiParticles.empty(); + const bool skipPion = !activeCorrelationTypes->at(kPion) || pionParticles.empty(); + + if (phiParticles.empty() && skipK0s && skipXi && skipPion) { return; + } int multBin = getCentBin(multiplicity); - if (multBin < 0) + if (multBin < 0) { return; - - bool doAssocCorrelations[kAssocPartSize] = {activeCorrelationConfigs.doK0SCorrelation, - activeCorrelationConfigs.doXiCorrelation, - activeCorrelationConfigs.doPionCorrelation}; + } // Same Event Correlations - std::vector* currentAssocParticles[] = {&k0sParticles, &xiParticles, &pionParticles}; + std::array*, kAssocPartSize> currentAssocParticles{&k0sParticles, &xiParticles, &pionParticles}; for (const auto& phiParticle : phiParticles) { histos.fill(HIST("phi/h3PhiMCClosureGen"), multiplicity, phiParticle.pt, phiParticle.y); static_for<0, AssocParticleLabels.size() - 1>([&](auto i_idx) { constexpr unsigned int Idx = i_idx.value; - if (!doAssocCorrelations[Idx]) + if (!activeCorrelationTypes->at(Idx)) { return; + } for (const auto& assocParticle : *(currentAssocParticles[Idx])) { histos.fill(HIST("phi") + HIST(AssocParticleLabels[Idx]) + HIST("/h5Phi") + HIST(AssocParticleLabels[Idx]) + HIST("ClosureMCGen"), @@ -1202,14 +1690,15 @@ struct PhiStrangenessCorrelation { // Mixed Event Correlations for (const auto& pastEvent : eventBuffer[multBin]) { - const std::vector* pastAssocParticles[] = {&pastEvent.k0sParticles, &pastEvent.xiParticles, &pastEvent.pionParticles}; + std::array*, kAssocPartSize> pastAssocParticles{&pastEvent.k0sParticles, &pastEvent.xiParticles, &pastEvent.pionParticles}; // Loop over past events in the same multiplicity bin and fill histograms with all combinations of current phi particles and past associated particles for (const auto& phiParticle : phiParticles) { static_for<0, AssocParticleLabels.size() - 1>([&](auto i_idx) { constexpr unsigned int Idx = i_idx.value; - if (!doAssocCorrelations[Idx]) + if (!activeCorrelationTypes->at(Idx)) { return; + } for (const auto& assocParticle : *(pastAssocParticles[Idx])) { histos.fill(HIST("phi") + HIST(AssocParticleLabels[Idx]) + HIST("/h5Phi") + HIST(AssocParticleLabels[Idx]) + HIST("ClosureMCGenME"), @@ -1229,14 +1718,15 @@ struct PhiStrangenessCorrelation { currentEvent.pionParticles = std::move(pionParticles); eventBuffer[multBin].push_front(std::move(currentEvent)); - if (eventBuffer[multBin].size() > static_cast(cfgNoMixedEvents.value)) + if (eventBuffer[multBin].size() > static_cast(cfgNoMixedEvents.value)) { eventBuffer[multBin].pop_back(); + } } - PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosure, "Process function for MC Gen Closure Test in SE and ME", false); + PROCESS_SWITCH(PhiStrangeCorrelation, processMCGenClosure, "Process function for MC Gen Closure Test in SE and ME", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }