From 4f61a014dc05cef40b4e52afcf2da09d170e92f5 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Wed, 12 Aug 2026 10:10:48 +0200 Subject: [PATCH 1/4] remove unused configurables --- PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx index 654d320f3e8..5a564605877 100644 --- a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx @@ -206,10 +206,8 @@ struct FlowGenericFramework { Configurable> resonanceSwitches{"resonanceSwitches", {LongArrayInt.front().data(), 8, 3, {"UseParticle", "UseCosPA", "NMassBins", "UseDCAxDaughters", "UseProperLifetime", "UseV0Radius", "UseArmPodCut", "UseCompetingMassRejection"}, {"K0", "Lambda", "Phi"}}, "Labeled array (int) for various cuts on resonances"}; struct : ConfigurableGroup { - O2_DEFINE_CONFIGURABLE(cfgUseLsPhi, bool, true, "Use LikeSign for Phi v2") O2_DEFINE_CONFIGURABLE(cfgUseOnlyTPC, bool, true, "Use only TPC PID for daughter selection") O2_DEFINE_CONFIGURABLE(cfgDaughterPIDRejection, bool, true, "Reject daughters if not consistent with expected K0/Lambda decay products") - O2_DEFINE_CONFIGURABLE(cfgFakeKaonCut, float, 0.1f, "Maximum difference in measured momentum and TPC inner ring momentum of particle") O2_DEFINE_CONFIGURABLE(cfgUseAsymmetricPID, bool, false, "Use asymmetric PID cuts") O2_DEFINE_CONFIGURABLE(cfgTPCNsigmaCut, float, 3.0f, "TPC N-sigma cut for pions, kaons, protons") O2_DEFINE_CONFIGURABLE(cfgUseStrictPID, bool, true, "Use strict PID cuts for TPC") From 41a08d24ac3835721b87adf333acbadac065940a Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Wed, 12 Aug 2026 13:22:02 +0200 Subject: [PATCH 2/4] add MC validation of estimator with truth nMPI --- PWGCF/DataModel/CorrelationsDerived.h | 10 + PWGCF/TableProducer/filterCorrelations.cxx | 16 +- .../Tasks/twoParticleCorrelationsMpi.cxx | 198 ++++++++++++++++-- 3 files changed, 200 insertions(+), 24 deletions(-) diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index 09f41d5ebfc..ec0f2aa11c2 100644 --- a/PWGCF/DataModel/CorrelationsDerived.h +++ b/PWGCF/DataModel/CorrelationsDerived.h @@ -28,6 +28,16 @@ DECLARE_SOA_TABLE(CFMcCollisions, "AOD", "CFMCCOLLISION", //! Reduced MC collisi mccollision::PosZ, cfmccollision::Multiplicity); using CFMcCollision = CFMcCollisions::iterator; +namespace cfmccollisionextra +{ +DECLARE_SOA_COLUMN(NMPI, nMPI, int); //! Number of multi-parton interactions from HepMC +} // namespace cfmccollisionextra +DECLARE_SOA_TABLE(CFMcCollisionExtras, "AOD", "CFMCCOLLEXTRA", //! Row-aligned extension of CFMcCollisions + cfmccollisionextra::NMPI); +using CFMcCollisionExtra = CFMcCollisionExtras::iterator; +using CFMcCollisionsWithExtra = soa::Join; +using CFMcCollisionWithExtra = CFMcCollisionsWithExtra::iterator; + namespace cfmcparticle { DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index bdc404ea7f0..3063ff72605 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -125,6 +125,7 @@ struct FilterCF { Produces outputTrackLabels; Produces outputMcCollisions; + Produces outputMcCollisionExtras; Produces outputMcParticles; Produces outputCollRefs; @@ -385,8 +386,8 @@ struct FilterCF { /// event selections /// \param tracks The collection of tracks, filtered by selection criteria /// \param bcs The collection of bunch crossings with timestamps - template - void processMCT(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + template + void processMCT(MCs const& mcCollisions, aod::McParticles const& allParticles, C1 const& allCollisions, T1 const& tracks, aod::BCsWithTimestamps const&) @@ -459,6 +460,7 @@ struct FilterCF { } } outputMcCollisions(mcCollision.posZ(), multiplicity); + outputMcCollisionExtras(mcCollision.nMPI()); } // PASS 2 on collisions: store collisions and tracks @@ -517,7 +519,8 @@ struct FilterCF { // NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions Preslice perMcCollision = aod::mcparticle::mcCollisionId; Preslice perCollision = aod::track::collisionId; - void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + using McCollisionsWithHepMC = soa::Join; + void processMC(McCollisionsWithHepMC const& mcCollisions, aod::McParticles const& allParticles, soa::Join const& allCollisions, soa::Filtered> const& tracks, aod::BCsWithTimestamps const& bcs) @@ -527,7 +530,7 @@ struct FilterCF { PROCESS_SWITCH(FilterCF, processMC, "Process MC", false); // NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions - void processMCPid(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + void processMCPid(McCollisionsWithHepMC const& mcCollisions, aod::McParticles const& allParticles, soa::Join const& allCollisions, soa::Filtered> const& tracks, aod::BCsWithTimestamps const& bcs) @@ -536,7 +539,7 @@ struct FilterCF { } PROCESS_SWITCH(FilterCF, processMCPid, "Process MC with PID", false); - void processMCMults(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + void processMCMults(McCollisionsWithHepMC const& mcCollisions, aod::McParticles const& allParticles, soa::Join const& allCollisions, soa::Filtered> const& tracks, aod::BCsWithTimestamps const& bcs) @@ -546,7 +549,7 @@ struct FilterCF { PROCESS_SWITCH(FilterCF, processMCMults, "Process MC with multiplicity sets", false); - void processMCGen(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& particles) + void processMCGen(McCollisionsWithHepMC::iterator const& mcCollision, aod::McParticles const& particles) { float multiplicity = 0.0f; for (auto& particle : particles) { @@ -562,6 +565,7 @@ struct FilterCF { sign, particle.pdgCode(), particle.flags()); } outputMcCollisions(mcCollision.posZ(), multiplicity); + outputMcCollisionExtras(mcCollision.nMPI()); } PROCESS_SWITCH(FilterCF, processMCGen, "Process MCGen", false); }; diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index abc1a16adc3..a4071c3d2ac 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -38,6 +38,7 @@ #include #include #include +#include #include #include @@ -144,6 +145,7 @@ struct TwoParticleCorrelationsMpi { // MC filters Filter cfMCCollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex; Filter cfMCParticleFilter = (nabs(aod::cfmcparticle::eta) < cfgCutEta) && (aod::cfmcparticle::pt > cfgCutPt); // && (aod::cfmcparticle::sign != 0); //check the sign manually, some specials may be neutral + Filter mcParticleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPt); // Output definitions OutputObj same{"sameEvent"}; @@ -206,15 +208,21 @@ struct TwoParticleCorrelationsMpi { PairCuts mPairCuts; Service ccdb{}; + Service pdg{}; using AodCollisions = soa::Filtered>; using AodTracks = soa::Filtered>; + using FilteredMcParticles = soa::Filtered; + using McCollisionsWithHepMC = soa::Join; using DerivedCollisions = soa::Filtered; using DerivedTracks = soa::Filtered; void init(o2::framework::InitContext&) { + if (doprocessMCSameDerived && (doprocessSameDerived || doprocessSameDerivedMultSet)) { + LOGF(fatal, "processMCSameDerived is mutually exclusive with the reconstructed derived same-event processes because it also fills those outputs"); + } registry.add("yields", "multiplicity/centrality vs pT vs eta", {HistType::kTH3F, {{100, 0, 100, "/multiplicity/centrality"}, {40, 0, 20, "p_{T}"}, {100, -2, 2, "#eta"}}}); registry.add("etaphi", "multiplicity/centrality vs eta vs phi", {HistType::kTH3F, {{100, 0, 100, "multiplicity/centrality"}, {100, -2, 2, "#eta"}, {200, 0, o2::constants::math::TwoPI, "#varphi"}}}); if (doprocessSameDerivedMultSet) { @@ -260,6 +268,33 @@ struct TwoParticleCorrelationsMpi { estimatorStatus->GetXaxis()->SetBinLabel(2, "no template-covered triggers"); estimatorStatus->GetXaxis()->SetBinLabel(3, "no template-covered pairs"); estimatorStatus->GetXaxis()->SetBinLabel(4, "valid estimate"); + registry.add("mcValidation/estimatedSeedsVsTrueNMPI", "template estimator response;N_{MPI}^{true};N_{seed}^{estimated}", {HistType::kTH2F, {{101, -0.5, 100.5}, {202, -0.5, 100.5}}}); + registry.add("mcValidation/profileEstimatedSeedsVsTrueNMPI", "mean template estimate;N_{MPI}^{true};#LT N_{seed}^{estimated} #GT", {HistType::kTProfile, {{101, -0.5, 100.5}}}); + registry.add("mcValidation/profileBiasVsTrueNMPI", "mean estimator bias;N_{MPI}^{true};#LT N_{seed}^{estimated} - N_{MPI}^{true} #GT", {HistType::kTProfile, {{101, -0.5, 100.5}}}); + registry.add("mcValidation/relativeResidualVsTrueNMPI", "relative estimator residual;N_{MPI}^{true};(N_{seed}^{estimated} - N_{MPI}^{true}) / N_{MPI}^{true}", {HistType::kTH2F, {{101, -0.5, 100.5}, {240, -3., 3.}}}); + registry.add("mcValidation/trueNMPIVsMultiplicity", "true MPI count versus reconstructed multiplicity;multiplicity;N_{MPI}^{true}", {HistType::kTH2F, {{100, 0., 100.}, {101, -0.5, 100.5}}}); + registry.add("mcValidation/templateCoverageVsTrueNMPI", "template pair coverage versus true MPI count;N_{MPI}^{true};matched / candidate pairs", {HistType::kTH2F, {{101, -0.5, 100.5}, {102, -0.01, 1.01}}}); + registry.add("mcValidation/status", "MC template-estimator validation status;status;events", {HistType::kTH1F, {{6, -0.5, 5.5}}}); + auto* mcValidationStatus = registry.get(HIST("mcValidation/status")).get(); + mcValidationStatus->GetXaxis()->SetBinLabel(1, "no selected reconstructed collision"); + mcValidationStatus->GetXaxis()->SetBinLabel(2, "invalid N MPI"); + mcValidationStatus->GetXaxis()->SetBinLabel(3, "outside template multiplicity"); + mcValidationStatus->GetXaxis()->SetBinLabel(4, "no template-covered triggers"); + mcValidationStatus->GetXaxis()->SetBinLabel(5, "no template-covered pairs"); + mcValidationStatus->GetXaxis()->SetBinLabel(6, "valid response"); + registry.add("mcValidation/generated/estimatedSeedsVsTrueNMPI", "generated-level template estimator response;N_{MPI}^{true};N_{seed,gen}^{estimated}", {HistType::kTH2F, {{101, -0.5, 100.5}, {202, -0.5, 100.5}}}); + registry.add("mcValidation/generated/profileEstimatedSeedsVsTrueNMPI", "mean generated-level template estimate;N_{MPI}^{true};#LT N_{seed,gen}^{estimated} #GT", {HistType::kTProfile, {{101, -0.5, 100.5}}}); + registry.add("mcValidation/generated/profileBiasVsTrueNMPI", "mean generated-level estimator bias;N_{MPI}^{true};#LT N_{seed,gen}^{estimated} - N_{MPI}^{true} #GT", {HistType::kTProfile, {{101, -0.5, 100.5}}}); + registry.add("mcValidation/generated/relativeResidualVsTrueNMPI", "generated-level relative estimator residual;N_{MPI}^{true};(N_{seed,gen}^{estimated} - N_{MPI}^{true}) / N_{MPI}^{true}", {HistType::kTH2F, {{101, -0.5, 100.5}, {240, -3., 3.}}}); + registry.add("mcValidation/generated/trueNMPIVsMultiplicity", "true MPI count versus generated charged multiplicity;N_{ch}^{gen};N_{MPI}^{true}", {HistType::kTH2F, {{101, -0.5, 100.5}, {101, -0.5, 100.5}}}); + registry.add("mcValidation/generated/templateCoverageVsTrueNMPI", "generated-level template pair coverage versus true MPI count;N_{MPI}^{true};matched / candidate pairs", {HistType::kTH2F, {{101, -0.5, 100.5}, {102, -0.01, 1.01}}}); + registry.add("mcValidation/generated/status", "generated-level MC template-estimator validation status;status;events", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + auto* generatedValidationStatus = registry.get(HIST("mcValidation/generated/status")).get(); + generatedValidationStatus->GetXaxis()->SetBinLabel(1, "invalid N MPI"); + generatedValidationStatus->GetXaxis()->SetBinLabel(2, "outside template multiplicity"); + generatedValidationStatus->GetXaxis()->SetBinLabel(3, "no template-covered triggers"); + generatedValidationStatus->GetXaxis()->SetBinLabel(4, "no template-covered pairs"); + generatedValidationStatus->GetXaxis()->SetBinLabel(5, "valid response"); registry.add("yvspt", "y vs pT", {HistType::kTH2F, {{100, -1, 1, "y"}, {100, 0, 20, "p_{T}"}}}); // y vs pT for all tracks (control histogram) const int maxMixBin = AxisSpec(axisMultiplicity).getNbins() * AxisSpec(axisVertex).getNbins(); @@ -472,6 +507,20 @@ struct TwoParticleCorrelationsMpi { template using HasPartDaugh1Id = decltype(std::declval().cfParticleDaugh1Id()); + template + int getParticleSign(const TParticle& particle) + { + if constexpr (std::experimental::is_detected::value) { + return particle.sign(); + } else if constexpr (std::experimental::is_detected::value) { + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle) { + return (pdgParticle->Charge() > 0.0) - (pdgParticle->Charge() < 0.0); + } + } + return 0; + } + template bool passOutlier(CollType const& collision) { @@ -718,6 +767,72 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("profileEventNuncSeeds"), multiplicity, estimate.nuncSeeds()); } + void fillMCValidation(double multiplicity, const EventSeedEstimate& estimate, int trueNMPI) + { + if (trueNMPI < 0) { + registry.fill(HIST("mcValidation/status"), 1.0); + return; + } + registry.fill(HIST("mcValidation/trueNMPIVsMultiplicity"), multiplicity, trueNMPI); + if (!hasMultiplicityTemplate(multiplicity)) { + registry.fill(HIST("mcValidation/status"), 2.0); + return; + } + if (!estimate.isValid()) { + registry.fill(HIST("mcValidation/status"), 3.0); + return; + } + const double coverage = estimate.nCandidatePairs > 0 ? static_cast(estimate.nPairs) / estimate.nCandidatePairs : 0.0; + registry.fill(HIST("mcValidation/templateCoverageVsTrueNMPI"), trueNMPI, coverage); + if (estimate.nPairs == 0) { + registry.fill(HIST("mcValidation/status"), 4.0); + return; + } + + const double estimatedSeeds = estimate.nuncSeeds(); + const double bias = estimatedSeeds - trueNMPI; + registry.fill(HIST("mcValidation/status"), 5.0); + registry.fill(HIST("mcValidation/estimatedSeedsVsTrueNMPI"), trueNMPI, estimatedSeeds); + registry.fill(HIST("mcValidation/profileEstimatedSeedsVsTrueNMPI"), trueNMPI, estimatedSeeds); + registry.fill(HIST("mcValidation/profileBiasVsTrueNMPI"), trueNMPI, bias); + if (trueNMPI > 0) { + registry.fill(HIST("mcValidation/relativeResidualVsTrueNMPI"), trueNMPI, bias / trueNMPI); + } + } + + void fillGeneratedMCValidation(double multiplicity, const EventSeedEstimate& estimate, int trueNMPI) + { + if (trueNMPI < 0) { + registry.fill(HIST("mcValidation/generated/status"), 0.0); + return; + } + registry.fill(HIST("mcValidation/generated/trueNMPIVsMultiplicity"), multiplicity, trueNMPI); + if (!hasMultiplicityTemplate(multiplicity)) { + registry.fill(HIST("mcValidation/generated/status"), 1.0); + return; + } + if (!estimate.isValid()) { + registry.fill(HIST("mcValidation/generated/status"), 2.0); + return; + } + const double coverage = estimate.nCandidatePairs > 0 ? static_cast(estimate.nPairs) / estimate.nCandidatePairs : 0.0; + registry.fill(HIST("mcValidation/generated/templateCoverageVsTrueNMPI"), trueNMPI, coverage); + if (estimate.nPairs == 0) { + registry.fill(HIST("mcValidation/generated/status"), 3.0); + return; + } + + const double estimatedSeeds = estimate.nuncSeeds(); + const double bias = estimatedSeeds - trueNMPI; + registry.fill(HIST("mcValidation/generated/status"), 4.0); + registry.fill(HIST("mcValidation/generated/estimatedSeedsVsTrueNMPI"), trueNMPI, estimatedSeeds); + registry.fill(HIST("mcValidation/generated/profileEstimatedSeedsVsTrueNMPI"), trueNMPI, estimatedSeeds); + registry.fill(HIST("mcValidation/generated/profileBiasVsTrueNMPI"), trueNMPI, bias); + if (trueNMPI > 0) { + registry.fill(HIST("mcValidation/generated/relativeResidualVsTrueNMPI"), trueNMPI, bias / trueNMPI); + } + } + template void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr) { @@ -754,11 +869,12 @@ struct TwoParticleCorrelationsMpi { continue; } } else { // otherwise check the sign against the configuration + const int sign = getParticleSign(track1); if (cfgTriggerCharge != 0) { - if (cfgTriggerCharge * track1.sign() < 0) { + if (cfgTriggerCharge * sign < 0) { continue; } - } else if (track1.sign() == 0) { + } else if (sign == 0) { continue; // reject neutral MC particles } } @@ -865,18 +981,20 @@ struct TwoParticleCorrelationsMpi { continue; } - if constexpr (std::experimental::is_detected::value) { + if constexpr (std::experimental::is_detected::value || std::experimental::is_detected::value) { + const int associatedSign = getParticleSign(track2); if (cfgAssociatedCharge != 0) { - if (cfgAssociatedCharge * track2.sign() < 0) { + if (cfgAssociatedCharge * associatedSign < 0) { continue; } - } else if (track2.sign() == 0) { // mc particles come in neutrals, need to check explicitly + } else if (associatedSign == 0) { // mc particles come in neutrals, need to check explicitly continue; } } - if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { - if (cfgPairCharge != 0 && cfgPairCharge * track1.sign() * track2.sign() < 0) { + if constexpr ((std::experimental::is_detected::value || std::experimental::is_detected::value) && + (std::experimental::is_detected::value || std::experimental::is_detected::value)) { + if (cfgPairCharge != 0 && cfgPairCharge * getParticleSign(track1) * getParticleSign(track2) < 0) { continue; } } @@ -975,8 +1093,8 @@ struct TwoParticleCorrelationsMpi { return multiplicity; } - // Version with explicit nested loop - void processSameAOD(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracks const& tracks) + template + void processSameAODT(TCollision const& collision, TTracks const& tracks, const int* trueNMPI = nullptr) { // NOTE legacy function for O2 integration tests. Full version needs derived data @@ -985,7 +1103,7 @@ struct TwoParticleCorrelationsMpi { } // TODO will go to CCDBConfigurable - auto bc = collision.bc_as(); + auto bc = collision.template bc_as(); loadEfficiency(bc.timestamp()); loadCcdbYieldTemplates(bc.timestamp()); @@ -999,11 +1117,46 @@ struct TwoParticleCorrelationsMpi { EventSeedEstimate seedEstimate; fillCorrelations(same, tracks, tracks, multiplicity, collision.posZ(), getMagneticField(bc.timestamp()), 1.0f, &seedEstimate); fillEventSeedEstimatorQA(multiplicity, seedEstimate); + if (trueNMPI) { + fillMCValidation(multiplicity, seedEstimate, *trueNMPI); + } + } + + // Version with explicit nested loop + void processSameAOD(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracks const& tracks) + { + processSameAODT(collision, tracks); } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameAOD, "Process same event on AOD", true); + void processSameGenMC(McCollisionsWithHepMC::iterator const& mcCollision, FilteredMcParticles const& mcParticles, aod::BCsWithTimestamps const&) + { + if (std::abs(mcCollision.posZ()) >= cfgCutVertex) { + return; + } + const auto bc = mcCollision.bc_as(); + loadCcdbYieldTemplates(bc.timestamp()); + + int generatedMultiplicity = 0; + for (const auto& particle : mcParticles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle && pdgParticle->Charge() != 0.0) { + ++generatedMultiplicity; + } + } + + fillContainerEvent(same, generatedMultiplicity, CorrelationContainer::kCFStepAll); + EventSeedEstimate seedEstimate; + fillCorrelations(same, mcParticles, mcParticles, generatedMultiplicity, mcCollision.posZ(), 0, 1.0f, &seedEstimate); + fillGeneratedMCValidation(generatedMultiplicity, seedEstimate, mcCollision.nMPI()); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameGenMC, "Process generated MC events and validate the template estimator against HepMC N MPI", false); + template - void processSameDerivedT(CollType const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2) + void processSameDerivedT(CollType const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2, const int* trueNMPI = nullptr) { using BinningTypeDerived = ColumnBinningPolicy; BinningTypeDerived configurableBinningDerived{{axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1. @@ -1042,6 +1195,9 @@ struct TwoParticleCorrelationsMpi { fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, fillReco ? nullptr : &seedEstimate); } fillEventSeedEstimatorQA(multiplicity, seedEstimate); + if (trueNMPI) { + fillMCValidation(multiplicity, seedEstimate, *trueNMPI); + } } void processSameDerived(DerivedCollisions::iterator const& collision, soa::Filtered const& tracks) @@ -1253,8 +1409,8 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMCEfficiency, "MC: Extract efficiencies", false); - template - void processMCSameDerivedT(soa::Filtered::iterator const& mcCollision, Particles1 const& mcParticles1, Particles2 const& mcParticles2, soa::SmallGroups const& collisions) + template + void processMCSameDerivedT(McCollision const& mcCollision, Particles1 const& mcParticles1, Particles2 const& mcParticles2, soa::SmallGroups const& collisions) { if (cfgVerbosity > 0) { LOGF(info, "processMCSameDerivedT. MC collision: %d, particles1: %d, particles2: %d, collisions: %d", mcCollision.globalIndex(), mcParticles1.size(), mcParticles2.size(), collisions.size()); @@ -1270,7 +1426,7 @@ struct TwoParticleCorrelationsMpi { } } - if (!(doprocessSameDerived || doprocessSameDerivedMultSet)) { + if (!(doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedMultSet)) { if constexpr (std::experimental::is_detected::value) { fillQA(mcCollision, multiplicity, mcCollision.posZ(), mcParticles1, mcParticles2); } else { @@ -1294,16 +1450,22 @@ struct TwoParticleCorrelationsMpi { fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepTracked); fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); - // NOTE kCFStepReconstructed and kCFStepCorrected are filled in processSameDerived - // This also means that if a MC collision had several reconstructed vertices (collisions), all of them are filled + // kCFStepReconstructed and kCFStepCorrected are filled below for every + // reconstructed collision associated with this MC collision. } // NOTE SmallGroups includes soa::Filtered always - void processMCSameDerived(soa::Filtered::iterator const& mcCollision, soa::Filtered const& mcParticles, soa::SmallGroups const& collisions) // TODO. For mixed no need to check the daughters since the events are different + Preslice derivedTracksPerCollision = aod::cftrack::cfCollisionId; + void processMCSameDerived(soa::Filtered::iterator const& mcCollision, soa::Filtered const& mcParticles, soa::SmallGroups const& collisions, soa::Filtered const& tracks) // TODO. For mixed no need to check the daughters since the events are different { processMCSameDerivedT(mcCollision, mcParticles, mcParticles, collisions); + const int trueNMPI = mcCollision.nMPI(); + for (const auto& collision : collisions) { + auto collisionTracks = tracks.sliceBy(derivedTracksPerCollision, collision.globalIndex()); + processSameDerivedT(collision, collisionTracks, collisionTracks, &trueNMPI); + } } - PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMCSameDerived, "Process MC same event on derived data", false); + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMCSameDerived, "Process generated and reconstructed MC same events on derived data and validate against HepMC N MPI", false); PresliceUnsorted collisionPerMCCollision = aod::cfcollision::cfMcCollisionId; template From d8ce1c52638d0b4d9c2b0c26d01ff331756aa4c0 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Wed, 12 Aug 2026 20:31:05 +0200 Subject: [PATCH 3/4] linter for twoParticleCorrelationsMpi.cxx --- .../Tasks/twoParticleCorrelationsMpi.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index a4071c3d2ac..cb744db3a3b 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -602,7 +602,7 @@ struct TwoParticleCorrelationsMpi { yieldTemplates.clear(); yieldTemplates.reserve(tree->GetEntries()); - for (Long64_t iEntry = 0; iEntry < tree->GetEntries(); ++iEntry) { + for (int64_t iEntry = 0; iEntry < tree->GetEntries(); ++iEntry) { tree->GetEntry(iEntry); if (fitStatus != 0) { LOGF(warning, "Skipping failed yield template (%d, %d, %d), fit status %d", value.trigBin, value.assocBin, value.multBin, fitStatus); @@ -614,11 +614,11 @@ struct TwoParticleCorrelationsMpi { continue; } const auto intervalMatchesAxis = [](const AxisSpec& axis, double low, double high) { - constexpr double tolerance = 1e-6; + constexpr double Tolerance = 1e-6; const auto& edges = axis.binEdges; return std::any_of(edges.begin(), edges.end() - 1, [&](const auto& edge) { const auto index = static_cast(&edge - edges.data()); - return std::abs(edge - low) < tolerance && std::abs(edges[index + 1] - high) < tolerance; + return std::abs(edge - low) < Tolerance && std::abs(edges[index + 1] - high) < Tolerance; }); }; if (!intervalMatchesAxis(AxisSpec(axisMultiplicity), value.nchLow, value.nchHigh) || From 6298cf050fc470dd7b8b0171ddc961e7a25d6a59 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Thu, 13 Aug 2026 00:14:41 +0200 Subject: [PATCH 4/4] code check for twoParticleCorrelationsMpi.cxx --- .../Tasks/twoParticleCorrelationsMpi.cxx | 36 ++++++++++--------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index cb744db3a3b..04f6d6622b9 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -48,6 +48,7 @@ #include #include #include +#include #include #include @@ -58,6 +59,7 @@ #include #include #include +#include #include #include #include @@ -186,10 +188,10 @@ struct TwoParticleCorrelationsMpi { double awayPairs = 0.0; double baselinePairs = 0.0; - bool isValid() const { return nTriggers > 0; } - double nearYield() const { return isValid() ? nearPairs / nTriggers : 0.0; } - double awayYield() const { return isValid() ? awayPairs / nTriggers : 0.0; } - double nuncSeeds() const + [[nodiscard]] bool isValid() const { return nTriggers > 0; } + [[nodiscard]] double nearYield() const { return isValid() ? nearPairs / nTriggers : 0.0; } + [[nodiscard]] double awayYield() const { return isValid() ? awayPairs / nTriggers : 0.0; } + [[nodiscard]] double nuncSeeds() const { const double denominator = 1.0 + nearYield() + awayYield(); return isValid() && denominator > 0.0 ? nTriggers / denominator : -1.0; @@ -576,18 +578,18 @@ struct TwoParticleCorrelationsMpi { LOGF(fatal, "Missing ensembleYieldTemplates in %s", source.c_str()); return; } - if (!schemaVersion || TString(schemaVersion->GetTitle()) != "1") { + if (schemaVersion == nullptr || TString(schemaVersion->GetTitle()) != "1") { LOGF(fatal, "Unsupported or missing ensemble-yield template schema version in %s", source.c_str()); return; } - if (!correlationStep || TString(correlationStep->GetTitle()) != "kCFStepReconstructed") { + if (correlationStep == nullptr || TString(correlationStep->GetTitle()) != "kCFStepReconstructed") { LOGF(fatal, "Ensemble-yield templates must be derived at kCFStepReconstructed"); return; } YieldTemplate value; int fitStatus = -1; - double parameters[10] = {}; + std::array parameters{}; tree->SetBranchAddress("trigBin", &value.trigBin); tree->SetBranchAddress("assocBin", &value.assocBin); tree->SetBranchAddress("multBin", &value.multBin); @@ -597,7 +599,7 @@ struct TwoParticleCorrelationsMpi { tree->SetBranchAddress("trigPtHigh", &value.trigPtHigh); tree->SetBranchAddress("assocPtLow", &value.assocPtLow); tree->SetBranchAddress("assocPtHigh", &value.assocPtHigh); - tree->SetBranchAddress("parameters", parameters); + tree->SetBranchAddress("parameters", parameters.data()); tree->SetBranchAddress("fitStatus", &fitStatus); yieldTemplates.clear(); @@ -608,7 +610,7 @@ struct TwoParticleCorrelationsMpi { LOGF(warning, "Skipping failed yield template (%d, %d, %d), fit status %d", value.trigBin, value.assocBin, value.multBin, fitStatus); continue; } - std::copy_n(parameters, value.parameters.size(), value.parameters.begin()); + value.parameters = parameters; if (value.parameters[2] <= 0.0 || value.parameters[5] <= 0.0 || value.parameters[8] <= 0.0) { LOGF(warning, "Skipping yield template (%d, %d, %d) with non-positive Gaussian width", value.trigBin, value.assocBin, value.multBin); continue; @@ -616,10 +618,12 @@ struct TwoParticleCorrelationsMpi { const auto intervalMatchesAxis = [](const AxisSpec& axis, double low, double high) { constexpr double Tolerance = 1e-6; const auto& edges = axis.binEdges; - return std::any_of(edges.begin(), edges.end() - 1, [&](const auto& edge) { - const auto index = static_cast(&edge - edges.data()); - return std::abs(edge - low) < Tolerance && std::abs(edges[index + 1] - high) < Tolerance; - }); + for (std::size_t index = 0; index + 1 < edges.size(); ++index) { + if (std::abs(edges[index] - low) < Tolerance && std::abs(edges[index + 1] - high) < Tolerance) { + return true; + } + } + return false; }; if (!intervalMatchesAxis(AxisSpec(axisMultiplicity), value.nchLow, value.nchHigh) || !intervalMatchesAxis(AxisSpec(axisPtTrigger), value.trigPtLow, value.trigPtHigh) || @@ -651,7 +655,7 @@ struct TwoParticleCorrelationsMpi { return; } - TList* calibration = dynamic_cast(input->Get("ccdb_object")); + auto* calibration = dynamic_cast(input->Get("ccdb_object")); auto findObject = [&](const char* name) -> TObject* { if (auto* object = input->Get(name)) { return object; @@ -724,8 +728,8 @@ struct TwoParticleCorrelationsMpi { void addPairProbabilities(EventSeedEstimate& estimate, const YieldTemplate& yieldTemplate, double deltaPhi) const { const auto& parameters = yieldTemplate.parameters; - const double near = std::max(0.0, evaluateGaussian(deltaPhi, ¶meters[0]) + evaluateGaussian(deltaPhi, ¶meters[3])); - const double away = std::max(0.0, evaluateGaussian(deltaPhi, ¶meters[6])); + const double near = std::max(0.0, evaluateGaussian(deltaPhi, parameters.data()) + evaluateGaussian(deltaPhi, parameters.data() + 3)); + const double away = std::max(0.0, evaluateGaussian(deltaPhi, parameters.data() + 6)); const double baseline = std::max(0.0, parameters[9]); const double total = baseline + near + away; if (total <= 0.0) {