From 75219df023913805d18beda122cf78552824e0af Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 14 Sep 2026 20:20:20 +0200 Subject: [PATCH 1/5] Added GFW cumulants with MPI proxy classification --- .../Tasks/CMakeLists.txt | 2 +- .../Tasks/twoParticleCorrelationsMpi.cxx | 411 ++++++++++++++++++ 2 files changed, 412 insertions(+), 1 deletion(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt index ea65ecf1162..284098d88a7 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt +++ b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt @@ -120,7 +120,7 @@ o2physics_add_dpl_workflow(casc-di-hadron-corr o2physics_add_dpl_workflow(two-particle-correlations-mpi SOURCES twoParticleCorrelationsMpi.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore O2Physics::GFWCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(charge-balance-function diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index c17eb1de1b7..8fa85114fdd 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -16,6 +16,11 @@ #include "PWGCF/Core/CorrelationContainer.h" #include "PWGCF/Core/PairCuts.h" #include "PWGCF/DataModel/CorrelationsDerived.h" +#include "PWGCF/GenericFramework/Core/FlowContainer.h" +#include "PWGCF/GenericFramework/Core/FlowPtContainer.h" +#include "PWGCF/GenericFramework/Core/GFW.h" +#include "PWGCF/GenericFramework/Core/GFWConfig.h" +#include "PWGCF/GenericFramework/Core/GFWWeights.h" #include "Common/CCDB/TriggerAliases.h" #include "Common/Core/RecoDecay.h" @@ -41,6 +46,7 @@ #include #include +#include #include #include #include @@ -49,7 +55,9 @@ #include #include #include +#include #include +#include #include #include @@ -65,6 +73,7 @@ #include #include #include +#include #include #include #include @@ -76,8 +85,15 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace constants::math; +using namespace o2::analysis::genericframework; static constexpr std::array, 1> CfgPairCutDefaults = {{{-1, -1, -1, -1, -1}}}; +static constexpr std::array, 4> LongArrayDouble = {{{{-0.8, -0.5}}, {{0.5, 0.8}}, {{-2, -2}}, {{-2, -2}}}}; +static constexpr int DefaultAcceptancePhiBins = 72; +static constexpr int DefaultAcceptanceEtaBins = 32; +static constexpr int InvalidRunNumber = -1; +static constexpr int FlowNBootstrap = 10; +static constexpr int MaxPtCorrelationOrder = 4; struct TwoParticleCorrelationsMpi { SliceCache cache; @@ -108,6 +124,11 @@ struct TwoParticleCorrelationsMpi { Configurable cfgEfficiencyTrigger{"cfgEfficiencyTrigger", "", "CCDB path to efficiency object for trigger particles"}; Configurable cfgEfficiencyAssociated{"cfgEfficiencyAssociated", "", "CCDB path to efficiency object for associated particles"}; + Configurable cfgAcceptance{"cfgAcceptance", "", "CCDB path to the GFW acceptance object"}; + Configurable cfgFillAcceptanceWeights{"cfgFillAcceptanceWeights", false, "Fill acceptance maps instead of physics outputs"}; + Configurable cfgAcceptanceRunByRun{"cfgAcceptanceRunByRun", false, "Fill or load run-by-run acceptance weights"}; + Configurable cfgAcceptancePhiBins{"cfgAcceptancePhiBins", int{DefaultAcceptancePhiBins}, "Number of phi bins in acceptance maps"}; + Configurable cfgAcceptanceEtaBins{"cfgAcceptanceEtaBins", int{DefaultAcceptanceEtaBins}, "Number of eta bins in acceptance maps"}; Configurable cfgNuncSeedsTemplateFile{"cfgNuncSeedsTemplateFile", "", "Local ROOT file containing ensembleYieldTemplates"}; Configurable cfgNuncSeedsTemplate{"cfgNuncSeedsTemplate", "", "CCDB path to the ensemble-yield template ccdb_object"}; Configurable cfgMinPairAcceptance{"cfgMinPairAcceptance", 0.05f, "Minimum pair acceptance used by the event seed estimator"}; @@ -145,6 +166,15 @@ struct TwoParticleCorrelationsMpi { ConfigurableAxis axisMultCorrV0{"axisMultCorrV0", {1000, 0, 100000}, "multiplicity correlation axis for V0 multiplicities"}; ConfigurableAxis axisMultCorrMult{"axisMultCorrMult", {1000, 0, 1000}, "multiplicity correlation axis for track multiplicities"}; + Configurable cfgRegions{"cfgRegions", {{"refN", "refP", "refFull"}, {-0.8, 0.4, -0.8}, {-0.4, 0.8, 0.8}, {0, 0, 0}, {1, 1, 1}}, "Configurations for GFW regions"}; + Configurable cfgCorrConfig{"cfgCorrConfig", {{"refP {2} refN {-2}", "refP {3} refN {-3}", "refP {4} refN {-4}", "refFull {2 -2}", "refFull {2 2 -2 -2}"}, {"ChGap22", "ChGap32", "ChGap42", "ChFull22", "ChFull24"}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}}, "Configurations for each correlation to calculate"}; + Configurable> cfgPtPtGaps{"cfgPtPtGaps", {LongArrayDouble.front().data(), 4, 2, {"subevent 1", "subevent 2", "subevent 3", "subevent 4"}, {"etamin", "etamax"}}, "{etamin,etamax} for all ptpt-subevents"}; + struct : ConfigurableGroup { + Configurable cfgUsePtCorrWeights{"cfgUsePtCorrWeights", true, "Use multiplicity-based event weighting for pt-pt correlations"}; + Configurable cfgUseMultiplicityFlowWeights{"cfgUseMultiplicityFlowWeights", true, "Use multiplicity-based event weighting for azimuthal correlations"}; + Configurable cfgUseMultiplicityFractionWeights{"cfgUseMultiplicityFractionWeights", false, "Use multiplicity-based event weighting for the spectral fraction"}; + } cfgEventWeight; + // This filter is applied to AOD and derived data (column names are identical) Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; // This filter is only applied to AOD @@ -173,7 +203,10 @@ struct TwoParticleCorrelationsMpi { bool mPairCuts = false; THn* mEfficiencyTrigger = nullptr; THn* mEfficiencyAssociated = nullptr; + GFWWeights* mAcceptance = nullptr; bool efficiencyLoaded = false; + bool acceptanceLoaded = false; + int acceptanceRunNumber = InvalidRunNumber; } cfg; struct YieldTemplate { @@ -276,6 +309,18 @@ struct TwoParticleCorrelationsMpi { HistogramRegistry registry{"registry"}; PairCuts mPairCuts; + // Generic Framework outputs + OutputObj fFC{"FlowContainer"}; + OutputObj fFCpt{"FlowPtContainer"}; + OutputObj fFCGen{"FlowContainer_gen"}; + OutputObj fFCptGen{"FlowPtContainer_gen"}; + std::unique_ptr fGFW{std::make_unique()}; + std::unique_ptr fPtAxis{}; + std::unique_ptr fRndm{std::make_unique()}; + std::vector mCorrConfigs; + std::shared_ptr mAcceptanceWeights; + std::map> mAcceptanceWeightsByRun; + Service ccdb{}; Service pdg{}; @@ -310,6 +355,15 @@ struct TwoParticleCorrelationsMpi { if (doprocessSameGenMC && doprocessMCSameDerived) { LOGF(fatal, "processSameGenMC and processMCSameDerived are mutually exclusive because both fill the generated same-event outputs"); } + if (cfgFillAcceptanceWeights && !cfgAcceptance.value.empty()) { + LOGF(fatal, "cfgFillAcceptanceWeights and cfgAcceptance are mutually exclusive: produce and apply acceptance weights in separate jobs"); + } + if (cfgFillAcceptanceWeights && !(doprocessSameDerived || doprocessSameDerivedMultSet || doprocessMCSameDerived)) { + LOGF(fatal, "cfgFillAcceptanceWeights requires a reconstructed derived same-event process"); + } + if (cfgAcceptancePhiBins < 1 || cfgAcceptanceEtaBins < 1) { + LOGF(fatal, "Acceptance-map bin counts must be positive"); + } if (!cfgNuncSeedsTemplateFile.value.empty() && !cfgNuncSeedsTemplate.value.empty()) { LOGF(fatal, "Configure only one template source: cfgNuncSeedsTemplateFile or cfgNuncSeedsTemplate"); } @@ -468,6 +522,11 @@ struct TwoParticleCorrelationsMpi { registry.add("eventcount_mixed", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}}); registry.add("trackcount_same", "bin", {HistType::kTH2F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {10, -0.5, 9.5}}}); registry.add("trackcount_mixed", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {10, -0.5, 9.5}, {10, -0.5, 9.5}}}); + if (cfgFillAcceptanceWeights && !cfgAcceptanceRunByRun) { + AxisSpec phiAxis{cfgAcceptancePhiBins, 0., o2::constants::math::TwoPI, "#varphi"}; + AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgCutEta.value, cfgCutEta.value, "#eta"}; + mAcceptanceWeights = registry.add("phi_eta_vtxz_ref", "Reference-track acceptance;#varphi;#eta;z_{vtx} (cm)", {HistType::kTH3D, {phiAxis, etaAxis, axisVertex}}); + } mPairCuts.SetHistogramRegistry(®istry); @@ -546,6 +605,276 @@ struct TwoParticleCorrelationsMpi { ccdb->setCreatedNotAfter(now); // TODO must become global parameter from the train creation time loadLocalYieldTemplates(); + + // Generic Framework init + AxisSpec ptAxisSpec = axisPtTrigger; + const int nPtBins = static_cast(ptAxisSpec.binEdges.size()) - 1; + fPtAxis = std::make_unique(nPtBins, ptAxisSpec.binEdges.data()); + if (cfgRegions->GetSize() < 0) { + LOGF(error, "Configuration contains vectors of different size - check the GFWRegions configurable"); + } + for (auto i(0); i < cfgRegions->GetSize(); ++i) { + fGFW->AddRegion(cfgRegions->GetNames()[i], cfgRegions->GetEtaMin()[i], cfgRegions->GetEtaMax()[i], (cfgRegions->GetpTDifs()[i] != 0) ? nPtBins + 1 : 1, cfgRegions->GetBitmasks()[i]); + } + for (auto i = 0; i < cfgCorrConfig->GetSize(); ++i) { + mCorrConfigs.push_back(fGFW->GetCorrelatorConfig(cfgCorrConfig->GetCorrs()[i], cfgCorrConfig->GetHeads()[i], cfgCorrConfig->GetpTDifs()[i] != 0)); + } + if (mCorrConfigs.empty()) { + LOGF(error, "Configuration contains vectors of different size - check the GFWCorrConfig configurable"); + } + fGFW->CreateRegions(); + int nPtPtSubevents = 0; + for (uint32_t i = 0; i < cfgPtPtGaps->rows(); ++i) { + if (cfgPtPtGaps->get(i, uint32_t{0}) < -1. || cfgPtPtGaps->get(i, uint32_t{1}) < -1.) { + continue; + } + ++nPtPtSubevents; + } + + if (!cfgFillAcceptanceWeights && (doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedMultSet)) { + auto recoProfiles = std::make_unique(); + addConfigObjectsToObjArray(recoProfiles.get(), mCorrConfigs, eventSeedPercentileAxisEnabled); + fFC.setObject(new FlowContainer("FlowContainer")); + fFC->SetXAxis(fPtAxis.get()); + fFC->Initialize(recoProfiles.get(), axisMultiplicity, FlowNBootstrap); + + fFCpt.setObject(new FlowPtContainer("FlowPtContainer")); + fFCpt->setEventWeight(cfgEventWeight.cfgUsePtCorrWeights ? eventweight::TupleWeight : eventweight::UnityWeight); + fFCpt->setUseCentralMoments(true); + fFCpt->setUseGapMethod(true); + fFCpt->initialise(axisMultiplicity, MaxPtCorrelationOrder, cfgCorrConfig.value, FlowNBootstrap); + fFCpt->initialiseSubevent(axisMultiplicity, MaxPtCorrelationOrder, nPtPtSubevents, FlowNBootstrap); + } + if (!cfgFillAcceptanceWeights && (doprocessMCSameDerived || doprocessSameGenMC)) { + auto generatedProfiles = std::make_unique(); + addConfigObjectsToObjArray(generatedProfiles.get(), mCorrConfigs, false); + fFCGen.setObject(new FlowContainer("FlowContainer_gen")); + fFCGen->SetXAxis(fPtAxis.get()); + fFCGen->Initialize(generatedProfiles.get(), axisMultiplicity, FlowNBootstrap); + + fFCptGen.setObject(new FlowPtContainer("FlowPtContainer_gen")); + fFCptGen->setEventWeight(cfgEventWeight.cfgUsePtCorrWeights ? eventweight::TupleWeight : eventweight::UnityWeight); + fFCptGen->setUseCentralMoments(true); + fFCptGen->setUseGapMethod(true); + fFCptGen->initialise(axisMultiplicity, MaxPtCorrelationOrder, cfgCorrConfig.value, FlowNBootstrap); + fFCptGen->initialiseSubevent(axisMultiplicity, MaxPtCorrelationOrder, nPtPtSubevents, FlowNBootstrap); + } + } + + enum DataType { + Reco, + Gen + }; + enum GFWRegionMask { + ReferenceMask = 1, + PoiMask = 2, + OverlapMask = 4 + }; + enum SeedClass { + InvalidSeedClass = -1, + LowSeedClass, + MiddleSeedClass, + HighSeedClass, + SeedClassCount + }; + + const char* getSeedClassSuffix(const SeedClass seedClass) const + { + switch (seedClass) { + case LowSeedClass: + return "__seedLow"; + case MiddleSeedClass: + return "__seedMiddle"; + case HighSeedClass: + return "__seedHigh"; + default: + return ""; + } + } + + std::string getSeedClassProfileName(const std::string& baseName, const SeedClass seedClass) const + { + return baseName + getSeedClassSuffix(seedClass); + } + + void addFlowProfileDefinition(TObjArray* profiles, const std::string& name, const std::string& title, const bool addSeedClasses) + { + profiles->Add(new TNamed(name.c_str(), title.c_str())); + if (!addSeedClasses) { + return; + } + for (int seedClassIndex = LowSeedClass; seedClassIndex < SeedClassCount; ++seedClassIndex) { + const auto seedClass = static_cast(seedClassIndex); + const auto classifiedName = getSeedClassProfileName(name, seedClass); + const auto classifiedTitle = getSeedClassProfileName(title, seedClass); + profiles->Add(new TNamed(classifiedName.c_str(), classifiedTitle.c_str())); + } + } + + void addConfigObjectsToObjArray(TObjArray* profiles, const std::vector& configs, const bool addSeedClasses) + { + profiles->SetOwner(kTRUE); + for (const auto& config : configs) { + if (config.pTDif) { + const std::string suffix = "_ptDiff"; + for (auto i = 0; i < fPtAxis->GetNbins(); ++i) { + const std::string index = Form("_pt_%i", i + 1); + const std::string name = config.Head + index; + const std::string title = config.Head + suffix; + addFlowProfileDefinition(profiles, name, title, addSeedClasses); + } + } else { + addFlowProfileDefinition(profiles, config.Head, config.Head, addSeedClasses); + } + } + } + + template + void fillOutputContainers(const float& centMult, const double& randomNumber, const SeedClass seedClass = InvalidSeedClass) + { + auto& flowContainer = (dt == Gen) ? fFCGen : fFC; + auto& flowPtContainer = (dt == Gen) ? fFCptGen : fFCpt; + const auto fillFlowProfile = [&](const std::string& profileName, const double value, const double weight) { + flowContainer->FillProfile(profileName.c_str(), centMult, value, weight, randomNumber); + if constexpr (dt == Reco) { + if (seedClass != InvalidSeedClass) { + const auto classifiedName = getSeedClassProfileName(profileName, seedClass); + flowContainer->FillProfile(classifiedName.c_str(), centMult, value, weight, randomNumber); + } + } + }; + + flowPtContainer->calculateCorrelations(); + flowPtContainer->calculateSubeventCorrelations(); + flowPtContainer->fillPtProfiles(centMult, randomNumber); + flowPtContainer->fillSubeventPtProfiles(centMult, randomNumber); + flowPtContainer->fillCMProfiles(centMult, randomNumber); + flowPtContainer->fillCMSubeventProfiles(centMult, randomNumber); + + for (std::size_t configIndex = 0; configIndex < mCorrConfigs.size(); ++configIndex) { + if (!mCorrConfigs.at(configIndex).pTDif) { + const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, kTRUE).real(); + if (dnx == 0) { + continue; + } + const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, kFALSE).real() / dnx; + if (std::abs(val) < 1) { + fillFlowProfile(mCorrConfigs.at(configIndex).Head, val, cfgEventWeight.cfgUseMultiplicityFlowWeights ? dnx : 1.0); + flowPtContainer->fillVnPtProfiles(centMult, val, dnx, randomNumber, cfgCorrConfig->GetpTCorrMasks()[configIndex]); + } + continue; + } + for (int i = 1; i <= fPtAxis->GetNbins(); i++) { + const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, kTRUE).real(); + if (dnx == 0) { + continue; + } + const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, kFALSE).real() / dnx; + if (std::abs(val) < 1) { + const std::string profileName = Form("%s_pt_%i", mCorrConfigs.at(configIndex).Head.c_str(), i); + fillFlowProfile(profileName, val, cfgEventWeight.cfgUseMultiplicityFlowWeights ? dnx : 1.0); + } + } + } + } + + template + inline void fillPtSums(const TTrack& track, const float& centMult, const double& posZ) + { + auto& flowPtContainer = (dt == Gen) ? fFCptGen : fFCpt; + double weff = 1.; + if constexpr (dt == Reco) { + if (cfg.mEfficiencyTrigger) { + weff = getEfficiencyCorrection(cfg.mEfficiencyTrigger, track.eta(), track.pt(), centMult, posZ); + } + } + if (weff < 0) { + return; + } + + // Fill the nominal sums + if (std::abs(track.eta()) < cfgCutEta) { + flowPtContainer->fill(weff, track.pt()); + } + + // Fill the subevent sums + std::size_t index = 0; + for (uint32_t i = 0; i < cfgPtPtGaps->rows(); ++i) { + const auto etamin = cfgPtPtGaps->get(i, uint32_t{0}); + const auto etamax = cfgPtPtGaps->get(i, uint32_t{1}); + if (etamin < -1. || etamax < -1.) { + continue; + } + if (etamin < track.eta() && track.eta() < etamax) { + flowPtContainer->fillSub(weff, track.pt(), index); + } + ++index; + } + } + + template + double getAcceptance(const TTrack& track, const double& posZ) + { + if (!cfg.mAcceptance) { + return 1.; + } + return cfg.mAcceptance->getNUA(track.phi(), track.eta(), posZ); + } + + std::shared_ptr getAcceptanceWeightsForRun(const int runNumber) + { + const auto existing = mAcceptanceWeightsByRun.find(runNumber); + if (existing != mAcceptanceWeightsByRun.end()) { + return existing->second; + } + + AxisSpec phiAxis{cfgAcceptancePhiBins, 0., o2::constants::math::TwoPI, "#varphi"}; + AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgCutEta.value, cfgCutEta.value, "#eta"}; + const std::string histogramName = Form("%d/phi_eta_vtxz_ref", runNumber); + const auto histogram = registry.add(histogramName.c_str(), "Reference-track acceptance;#varphi;#eta;z_{vtx} (cm)", {HistType::kTH3D, {phiAxis, etaAxis, axisVertex}}); + mAcceptanceWeightsByRun.emplace(runNumber, histogram); + return histogram; + } + + template + void fillAcceptanceWeights(const TTrack& track, const double& posZ, const int runNumber) + { + if (track.pt() <= fPtAxis->GetXmin() || track.pt() >= fPtAxis->GetXmax()) { + return; + } + const auto histogram = cfgAcceptanceRunByRun ? getAcceptanceWeightsForRun(runNumber) : mAcceptanceWeights; + histogram->Fill(track.phi(), track.eta(), posZ); + } + + template + inline void fillGFW(const TTrack& track, const float& centMult, const double& posZ) + { + const bool withinPtRef = track.pt() > fPtAxis->GetXmin() && track.pt() < fPtAxis->GetXmax(); + const bool withinPtPOI = withinPtRef; + if (!withinPtPOI && !withinPtRef) { + return; + } + double wacc = 1.; + double weff = 1.; + if constexpr (dt == Reco) { + wacc = getAcceptance(track, posZ); + if (cfg.mEfficiencyAssociated) { + weff = getEfficiencyCorrection(cfg.mEfficiencyAssociated, track.eta(), track.pt(), centMult, posZ); + } + } + if (weff < 0 || !std::isfinite(wacc) || wacc <= 0.) { + return; + } + const double weight = weff * wacc; + if (withinPtRef) { + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, ReferenceMask); + } + if (withinPtPOI) { + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, PoiMask); + } + if (withinPtRef && withinPtPOI) { + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, OverlapMask); + } } int getMagneticField(uint64_t timestamp) @@ -1008,6 +1337,20 @@ struct TwoParticleCorrelationsMpi { return 1.0; } + SeedClass getEventSeedClass(const double multiplicity, const double eventSeed) const + { + if (!eventSeedPercentileAxisEnabled) { + return InvalidSeedClass; + } + const double userAxisValue = getEventSeedUserAxisValue(multiplicity, eventSeed); + if (!std::isfinite(userAxisValue) || + userAxisValue < static_cast(LowSeedClass) || + userAxisValue >= static_cast(SeedClassCount)) { + return InvalidSeedClass; + } + return static_cast(userAxisValue); + } + double getPairAcceptance(double multiplicity, double deltaPhi, double deltaEta, double posZ) const { const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); @@ -1596,6 +1939,33 @@ struct TwoParticleCorrelationsMpi { cfg.efficiencyLoaded = true; } + void loadAcceptance(const uint64_t timestamp, const int runNumber) + { + if (cfgFillAcceptanceWeights || cfgAcceptance.value.empty()) { + cfg.mAcceptance = nullptr; + cfg.acceptanceLoaded = true; + return; + } + if (cfg.acceptanceLoaded && (!cfgAcceptanceRunByRun || cfg.acceptanceRunNumber == runNumber)) { + return; + } + + std::string path = cfgAcceptance.value; + if (cfgAcceptanceRunByRun) { + if (path.back() != '/') { + path += '/'; + } + path += "RunByRun/"; + } + cfg.mAcceptance = ccdb->getForTimeStamp(path, timestamp); + if (!cfg.mAcceptance || !cfg.mAcceptance->isDataFilled()) { + LOGF(fatal, "Could not load a populated GFW acceptance object from %s for run %d at timestamp %llu", path.c_str(), runNumber, timestamp); + } + cfg.acceptanceLoaded = true; + cfg.acceptanceRunNumber = runNumber; + LOGF(info, "Loaded GFW acceptance object from %s for run %d", path.c_str(), runNumber); + } + double getEfficiencyCorrection(THn* eff, float eta, float pt, float multiplicity, float posZ) { std::array effVars{}; @@ -1658,6 +2028,9 @@ struct TwoParticleCorrelationsMpi { soa::Filtered const& mcParticles, soa::SmallGroups const& collisions) { + if (cfgFillAcceptanceWeights) { + return; + } if (!cfgNuncSeedsTemplate.value.empty()) { for (const auto& collision : collisions) { loadCcdbYieldTemplates(collision.timestamp()); @@ -1697,7 +2070,14 @@ struct TwoParticleCorrelationsMpi { if (cfgVerbosity > 0) { LOGF(info, "processSameDerivedT: Tracks for collision: %d/%d | Vertex: %.1f | Multiplicity/Centrality: %.1f", tracks1.size(), tracks2.size(), collision.posZ(), multiplicity); } + if (cfgFillAcceptanceWeights) { + for (const auto& track : tracks2) { + fillAcceptanceWeights(track, collision.posZ(), collision.runNumber()); + } + return; + } loadEfficiency(collision.timestamp()); + loadAcceptance(collision.timestamp(), collision.runNumber()); loadCcdbYieldTemplates(collision.timestamp()); int field = 0; @@ -1705,6 +2085,13 @@ struct TwoParticleCorrelationsMpi { field = getMagneticField(collision.timestamp()); } + fGFW->Clear(); + fFCpt->clearVector(); + for (const auto& track : tracks2) { + fillPtSums(track, multiplicity, collision.posZ()); + fillGFW(track, multiplicity, collision.posZ()); + } + int bin = configurableBinningDerived.getBin(std::tuple(collision.posZ(), multiplicity)); registry.fill(HIST("eventcount_same"), bin); registry.fill(HIST("trackcount_same"), bin, tracks1.size()); @@ -1740,6 +2127,8 @@ struct TwoParticleCorrelationsMpi { fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepCorrected); fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, nullptr, nullptr, nullptr, getEventSeedUserAxisValue(multiplicity, seedEstimate.nuncSeeds())); } + const auto seedClass = getEventSeedClass(multiplicity, seedEstimate.nuncSeeds()); + fillOutputContainers(multiplicity, fRndm->Rndm(), seedClass); fillEventSeedEstimatorQA(multiplicity, seedEstimate); if (trueNMPI) { fillMCValidation(multiplicity, seedEstimate, *trueNMPI); @@ -1756,6 +2145,7 @@ struct TwoParticleCorrelationsMpi { fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, fillReco ? nullptr : &seedEstimate); } finalizeEventSeedEstimate(seedEstimate); + fillOutputContainers(multiplicity, fRndm->Rndm()); fillEventSeedEstimatorQA(multiplicity, seedEstimate); if (trueNMPI) { fillMCValidation(multiplicity, seedEstimate, *trueNMPI); @@ -1945,6 +2335,9 @@ struct TwoParticleCorrelationsMpi { template void processMCSameDerivedT(McCollision const& mcCollision, Particles1 const& mcParticles1, Particles2 const& mcParticles2, soa::SmallGroups const& collisions) { + if (cfgFillAcceptanceWeights) { + return; + } if (cfgVerbosity > 0) { LOGF(info, "processMCSameDerivedT. MC collision: %d, particles1: %d, particles2: %d, collisions: %d", mcCollision.globalIndex(), mcParticles1.size(), mcParticles2.size(), collisions.size()); } @@ -1967,6 +2360,24 @@ struct TwoParticleCorrelationsMpi { } } + fGFW->Clear(); + fFCptGen->clearVector(); + for (const auto& mcParticle : mcParticles2) { + if (!mcParticle.isPhysicalPrimary()) { + continue; + } + const int sign = getParticleSign(mcParticle); + if (sign == 0 || (cfgAssociatedCharge != 0 && cfgAssociatedCharge * sign < 0)) { + continue; + } + if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) { + continue; + } + fillPtSums(mcParticle, multiplicity, mcCollision.posZ()); + fillGFW(mcParticle, multiplicity, mcCollision.posZ()); + } + fillOutputContainers(multiplicity, fRndm->Rndm()); + fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepAll); fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); From 734718fa1e573d4b352ee626ef043e5f2dd31b06 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 14 Sep 2026 20:47:34 +0200 Subject: [PATCH 2/5] megalinter --- .../Tasks/twoParticleCorrelationsMpi.cxx | 26 ++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 8fa85114fdd..14ca38f33f6 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -102,6 +102,8 @@ struct TwoParticleCorrelationsMpi { Configurable cfgCutVertex{"cfgCutVertex", 7.0f, "Accepted z-vertex range"}; Configurable cfgCutPt{"cfgCutPt", 0.5f, "Minimal pT for tracks"}; Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable> cfgFlowPtRef{"cfgFlowPtRef", {0.5f, 3.0f}, "Minimum and maximum pT for reference flow particles"}; + Configurable> cfgFlowPtPOI{"cfgFlowPtPOI", {0.5f, 10.0f}, "Minimum and maximum pT for particles of interest in flow correlations"}; Configurable cfgPtOrder{"cfgPtOrder", 1, "Only consider pairs for which pT,1 < pT,2 (0 = OFF, 1 = ON)"}; Configurable cfgTriggerCharge{"cfgTriggerCharge", 0, "Select on charge of trigger particle: 0 = all; 1 = positive; -1 = negative"}; @@ -608,6 +610,17 @@ struct TwoParticleCorrelationsMpi { // Generic Framework init AxisSpec ptAxisSpec = axisPtTrigger; + const auto validateFlowPtRange = [&ptAxisSpec](const std::pair& range, const char* name) { + if (!std::isfinite(range.first) || !std::isfinite(range.second) || range.first >= range.second) { + LOGF(fatal, "%s must define a finite, increasing pT range; received {%g, %g}", name, range.first, range.second); + } + if (range.first < ptAxisSpec.binEdges.front() || range.second > ptAxisSpec.binEdges.back()) { + LOGF(fatal, "%s={%g, %g} must be contained in axisPtTrigger={%g, %g}", + name, range.first, range.second, ptAxisSpec.binEdges.front(), ptAxisSpec.binEdges.back()); + } + }; + validateFlowPtRange(cfgFlowPtRef.value, "cfgFlowPtRef"); + validateFlowPtRange(cfgFlowPtPOI.value, "cfgFlowPtPOI"); const int nPtBins = static_cast(ptAxisSpec.binEdges.size()) - 1; fPtAxis = std::make_unique(nPtBins, ptAxisSpec.binEdges.data()); if (cfgRegions->GetSize() < 0) { @@ -849,9 +862,9 @@ struct TwoParticleCorrelationsMpi { template inline void fillGFW(const TTrack& track, const float& centMult, const double& posZ) { - const bool withinPtRef = track.pt() > fPtAxis->GetXmin() && track.pt() < fPtAxis->GetXmax(); - const bool withinPtPOI = withinPtRef; - if (!withinPtPOI && !withinPtRef) { + const bool withinPtRef = track.pt() > cfgFlowPtRef->first && track.pt() < cfgFlowPtRef->second; + const bool withinPtPOI = track.pt() > cfgFlowPtPOI->first && track.pt() < cfgFlowPtPOI->second; + if (!withinPtRef && !withinPtPOI) { return; } double wacc = 1.; @@ -866,14 +879,15 @@ struct TwoParticleCorrelationsMpi { return; } const double weight = weff * wacc; + const int ptBin = fPtAxis->FindBin(track.pt()) - 1; if (withinPtRef) { - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, ReferenceMask); + fGFW->Fill(track.eta(), ptBin, track.phi(), weight, ReferenceMask); } if (withinPtPOI) { - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, PoiMask); + fGFW->Fill(track.eta(), ptBin, track.phi(), weight, PoiMask); } if (withinPtRef && withinPtPOI) { - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), weight, OverlapMask); + fGFW->Fill(track.eta(), ptBin, track.phi(), weight, OverlapMask); } } From 97eb962701867bdb6d4324db86824c850c9a6890 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Tue, 15 Sep 2026 08:35:44 +0200 Subject: [PATCH 3/5] Reduce task members --- .../Tasks/twoParticleCorrelationsMpi.cxx | 113 +++++++++--------- 1 file changed, 56 insertions(+), 57 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 14ca38f33f6..ccecdc2eefc 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -99,24 +99,23 @@ struct TwoParticleCorrelationsMpi { SliceCache cache; // Configuration - Configurable cfgCutVertex{"cfgCutVertex", 7.0f, "Accepted z-vertex range"}; - Configurable cfgCutPt{"cfgCutPt", 0.5f, "Minimal pT for tracks"}; - Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; - Configurable> cfgFlowPtRef{"cfgFlowPtRef", {0.5f, 3.0f}, "Minimum and maximum pT for reference flow particles"}; - Configurable> cfgFlowPtPOI{"cfgFlowPtPOI", {0.5f, 10.0f}, "Minimum and maximum pT for particles of interest in flow correlations"}; - - Configurable cfgPtOrder{"cfgPtOrder", 1, "Only consider pairs for which pT,1 < pT,2 (0 = OFF, 1 = ON)"}; - Configurable cfgTriggerCharge{"cfgTriggerCharge", 0, "Select on charge of trigger particle: 0 = all; 1 = positive; -1 = negative"}; - Configurable cfgAssociatedCharge{"cfgAssociatedCharge", 0, "Select on charge of associated particle: 0 = all charged; 1 = positive; -1 = negative"}; - Configurable cfgPairCharge{"cfgPairCharge", 0, "Select on charge of particle pair: 0 = all; 1 = like sign; -1 = unlike sign"}; - Configurable cfgCorrelationMethod{"cfgCorrelationMethod", 0, "Correlation method, 0 = all, 1 = dd, 2 = ddbar"}; - - Configurable cfgTwoTrackCut{"cfgTwoTrackCut", -1, "Two track cut: -1 = off; >0 otherwise distance value (suggested: 0.02)"}; - Configurable cfgTwoTrackCutMinRadius{"cfgTwoTrackCutMinRadius", 0.8f, "Two track cut: radius in m from which two track cuts are applied"}; - ; - Configurable cfgLocalEfficiency{"cfgLocalEfficiency", 0, "0 = OFF and 1 = ON for local efficiency"}; - Configurable cfgDropStepRECO{"cfgDropStepRECO", false, "choice to drop step RECO if efficiency correction is used"}; - Configurable cfgCentBinsForMC{"cfgCentBinsForMC", 0, "0 = generated multiplicity; 1 = reconstructed multiplicity and all associated collisions"}; + struct : ConfigurableGroup { + Configurable cfgCutVertex{"cfgCutVertex", 7.0f, "Accepted z-vertex range"}; + Configurable cfgCutPt{"cfgCutPt", 0.5f, "Minimal pT for tracks"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable> cfgFlowPtRef{"cfgFlowPtRef", {0.5f, 3.0f}, "Minimum and maximum pT for reference flow particles"}; + Configurable> cfgFlowPtPOI{"cfgFlowPtPOI", {0.5f, 10.0f}, "Minimum and maximum pT for particles of interest in flow correlations"}; + Configurable cfgPtOrder{"cfgPtOrder", 1, "Only consider pairs for which pT,1 < pT,2 (0 = OFF, 1 = ON)"}; + Configurable cfgTriggerCharge{"cfgTriggerCharge", 0, "Select on charge of trigger particle: 0 = all; 1 = positive; -1 = negative"}; + Configurable cfgAssociatedCharge{"cfgAssociatedCharge", 0, "Select on charge of associated particle: 0 = all charged; 1 = positive; -1 = negative"}; + Configurable cfgPairCharge{"cfgPairCharge", 0, "Select on charge of particle pair: 0 = all; 1 = like sign; -1 = unlike sign"}; + Configurable cfgCorrelationMethod{"cfgCorrelationMethod", 0, "Correlation method, 0 = all, 1 = dd, 2 = ddbar"}; + Configurable cfgTwoTrackCut{"cfgTwoTrackCut", -1, "Two track cut: -1 = off; >0 otherwise distance value (suggested: 0.02)"}; + Configurable cfgTwoTrackCutMinRadius{"cfgTwoTrackCutMinRadius", 0.8f, "Two track cut: radius in m from which two track cuts are applied"}; + Configurable cfgLocalEfficiency{"cfgLocalEfficiency", 0, "0 = OFF and 1 = ON for local efficiency"}; + Configurable cfgDropStepRECO{"cfgDropStepRECO", false, "choice to drop step RECO if efficiency correction is used"}; + Configurable cfgCentBinsForMC{"cfgCentBinsForMC", 0, "0 = generated multiplicity; 1 = reconstructed multiplicity and all associated collisions"}; + } cfgGeneral; Configurable cfgTrackBitMask{"cfgTrackBitMask", 0, "BitMask for track selection systematics; refer to the enum TrackSelectionCuts in filtering task"}; Configurable cfgMultCorrelationsMask{"cfgMultCorrelationsMask", 0, "Selection bitmask for the multiplicity correlations. This should match the filter selection cfgEstimatorBitMask."}; Configurable cfgMultCutFormula{"cfgMultCutFormula", "", "Multiplicity correlations cut formula. A result greater than zero results in accepted event. Parameters: [cFT0C] FT0C centrality, [mFV0A] V0A multiplicity, [mGlob] global track multiplicity, [mPV] PV track multiplicity, [cFT0M] FT0M centrality"}; @@ -178,18 +177,18 @@ struct TwoParticleCorrelationsMpi { } cfgEventWeight; // This filter is applied to AOD and derived data (column names are identical) - Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgGeneral.cfgCutVertex; // This filter is only applied to AOD Filter collisionVertexTypeFilter = (aod::collision::flags & static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks)) == static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks); // Track filters - Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)); - Filter cfTrackFilter = (nabs(aod::cftrack::eta) < cfgCutEta) && (aod::cftrack::pt > cfgCutPt) && ncheckbit(aod::track::trackType, as(cfgTrackBitMask)); + Filter trackFilter = (nabs(aod::track::eta) < cfgGeneral.cfgCutEta) && (aod::track::pt > cfgGeneral.cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)); + Filter cfTrackFilter = (nabs(aod::cftrack::eta) < cfgGeneral.cfgCutEta) && (aod::cftrack::pt > cfgGeneral.cfgCutPt) && ncheckbit(aod::track::trackType, as(cfgTrackBitMask)); // 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); + Filter cfMCCollisionFilter = nabs(aod::mccollision::posZ) < cfgGeneral.cfgCutVertex; + Filter cfMCParticleFilter = (nabs(aod::cfmcparticle::eta) < cfgGeneral.cfgCutEta) && (aod::cfmcparticle::pt > cfgGeneral.cfgCutPt); // && (aod::cfmcparticle::sign != 0); //check the sign manually, some specials may be neutral + Filter mcParticleFilter = (nabs(aod::mcparticle::eta) < cfgGeneral.cfgCutEta) && (aod::mcparticle::pt > cfgGeneral.cfgCutPt); // Output definitions OutputObj same{"sameEvent"}; @@ -340,8 +339,8 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis < NoUserAxis || cfgUserAxis > EventSeedAxis) { LOGF(fatal, "Unsupported cfgUserAxis=%d; use 0 (off), 1 (invariant mass), or 2 (event seed)", cfgUserAxis.value); } - if (cfgCentBinsForMC < 0 || cfgCentBinsForMC > 1) { - LOGF(fatal, "Unsupported cfgCentBinsForMC=%d; use 0 (generated multiplicity), 1 (reconstructed multiplicity and all associated collisions)", cfgCentBinsForMC.value); + if (cfgGeneral.cfgCentBinsForMC < 0 || cfgGeneral.cfgCentBinsForMC > 1) { + LOGF(fatal, "Unsupported cfgCentBinsForMC=%d; use 0 (generated multiplicity), 1 (reconstructed multiplicity and all associated collisions)", cfgGeneral.cfgCentBinsForMC.value); } const int enabledDerivedSameProcesses = static_cast(doprocessSameDerived) + static_cast(doprocessSameDerivedCorrected) + static_cast(doprocessSameDerivedMultSet) + static_cast(doprocessSameDerivedMultSetCorrected); if (enabledDerivedSameProcesses > 1) { @@ -526,7 +525,7 @@ struct TwoParticleCorrelationsMpi { registry.add("trackcount_mixed", "bin", {HistType::kTH3F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}, {10, -0.5, 9.5}, {10, -0.5, 9.5}}}); if (cfgFillAcceptanceWeights && !cfgAcceptanceRunByRun) { AxisSpec phiAxis{cfgAcceptancePhiBins, 0., o2::constants::math::TwoPI, "#varphi"}; - AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgCutEta.value, cfgCutEta.value, "#eta"}; + AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgGeneral.cfgCutEta.value, cfgGeneral.cfgCutEta.value, "#eta"}; mAcceptanceWeights = registry.add("phi_eta_vtxz_ref", "Reference-track acceptance;#varphi;#eta;z_{vtx} (cm)", {HistType::kTH3D, {phiAxis, etaAxis, axisVertex}}); } @@ -541,8 +540,8 @@ struct TwoParticleCorrelationsMpi { cfg.mPairCuts = true; } - if (cfgTwoTrackCut > 0) { - mPairCuts.SetTwoTrackCuts(cfgTwoTrackCut, cfgTwoTrackCutMinRadius); + if (cfgGeneral.cfgTwoTrackCut > 0) { + mPairCuts.SetTwoTrackCuts(cfgGeneral.cfgTwoTrackCut, cfgGeneral.cfgTwoTrackCutMinRadius); } // --- OBJECT INIT --- @@ -592,8 +591,8 @@ struct TwoParticleCorrelationsMpi { same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userMixingAxis)); - same->setTrackEtaCut(cfgCutEta); - mixed->setTrackEtaCut(cfgCutEta); + same->setTrackEtaCut(cfgGeneral.cfgCutEta); + mixed->setTrackEtaCut(cfgGeneral.cfgCutEta); if (!cfgEfficiencyAssociated.value.empty()) { efficiencyAssociatedCache.reserve(512); @@ -619,8 +618,8 @@ struct TwoParticleCorrelationsMpi { name, range.first, range.second, ptAxisSpec.binEdges.front(), ptAxisSpec.binEdges.back()); } }; - validateFlowPtRange(cfgFlowPtRef.value, "cfgFlowPtRef"); - validateFlowPtRange(cfgFlowPtPOI.value, "cfgFlowPtPOI"); + validateFlowPtRange(cfgGeneral.cfgFlowPtRef.value, "cfgFlowPtRef"); + validateFlowPtRange(cfgGeneral.cfgFlowPtPOI.value, "cfgFlowPtPOI"); const int nPtBins = static_cast(ptAxisSpec.binEdges.size()) - 1; fPtAxis = std::make_unique(nPtBins, ptAxisSpec.binEdges.data()); if (cfgRegions->GetSize() < 0) { @@ -806,7 +805,7 @@ struct TwoParticleCorrelationsMpi { } // Fill the nominal sums - if (std::abs(track.eta()) < cfgCutEta) { + if (std::abs(track.eta()) < cfgGeneral.cfgCutEta) { flowPtContainer->fill(weff, track.pt()); } @@ -842,7 +841,7 @@ struct TwoParticleCorrelationsMpi { } AxisSpec phiAxis{cfgAcceptancePhiBins, 0., o2::constants::math::TwoPI, "#varphi"}; - AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgCutEta.value, cfgCutEta.value, "#eta"}; + AxisSpec etaAxis{cfgAcceptanceEtaBins, -cfgGeneral.cfgCutEta.value, cfgGeneral.cfgCutEta.value, "#eta"}; const std::string histogramName = Form("%d/phi_eta_vtxz_ref", runNumber); const auto histogram = registry.add(histogramName.c_str(), "Reference-track acceptance;#varphi;#eta;z_{vtx} (cm)", {HistType::kTH3D, {phiAxis, etaAxis, axisVertex}}); mAcceptanceWeightsByRun.emplace(runNumber, histogram); @@ -862,8 +861,8 @@ struct TwoParticleCorrelationsMpi { template inline void fillGFW(const TTrack& track, const float& centMult, const double& posZ) { - const bool withinPtRef = track.pt() > cfgFlowPtRef->first && track.pt() < cfgFlowPtRef->second; - const bool withinPtPOI = track.pt() > cfgFlowPtPOI->first && track.pt() < cfgFlowPtPOI->second; + const bool withinPtRef = track.pt() > cfgGeneral.cfgFlowPtRef->first && track.pt() < cfgGeneral.cfgFlowPtRef->second; + const bool withinPtPOI = track.pt() > cfgGeneral.cfgFlowPtPOI->first && track.pt() < cfgGeneral.cfgFlowPtPOI->second; if (!withinPtRef && !withinPtPOI) { return; } @@ -1727,8 +1726,8 @@ struct TwoParticleCorrelationsMpi { } } else { // otherwise check the sign against the configuration const int sign = getParticleSign(track1); - if (cfgTriggerCharge != 0) { - if (cfgTriggerCharge * sign < 0) { + if (cfgGeneral.cfgTriggerCharge != 0) { + if (cfgGeneral.cfgTriggerCharge * sign < 0) { continue; } } else if (sign == 0) { @@ -1737,7 +1736,7 @@ struct TwoParticleCorrelationsMpi { } } else if constexpr (std::experimental::is_detected::value) { // Check reco objects that have the sign attribute. There are no neutrals to deal with. - if (cfgTriggerCharge != 0 && cfgTriggerCharge * track1.sign() < 0) { + if (cfgGeneral.cfgTriggerCharge != 0 && cfgGeneral.cfgTriggerCharge * track1.sign() < 0) { continue; } } @@ -1837,22 +1836,22 @@ struct TwoParticleCorrelationsMpi { } if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { - if (cfgCorrelationMethod == CorrelationMethod::Dd && track1.decay() != track2.decay()) { + if (cfgGeneral.cfgCorrelationMethod == CorrelationMethod::Dd && track1.decay() != track2.decay()) { continue; } - if (cfgCorrelationMethod == CorrelationMethod::Ddbar && track1.decay() == track2.decay()) { + if (cfgGeneral.cfgCorrelationMethod == CorrelationMethod::Ddbar && track1.decay() == track2.decay()) { continue; } } - if (cfgPtOrder != 0 && track2.pt() >= track1.pt()) { + if (cfgGeneral.cfgPtOrder != 0 && track2.pt() >= track1.pt()) { continue; } if constexpr (std::experimental::is_detected::value || std::experimental::is_detected::value) { const int associatedSign = getParticleSign(track2); - if (cfgAssociatedCharge != 0) { - if (cfgAssociatedCharge * associatedSign < 0) { + if (cfgGeneral.cfgAssociatedCharge != 0) { + if (cfgGeneral.cfgAssociatedCharge * associatedSign < 0) { continue; } } else if (associatedSign == 0) { // mc particles come in neutrals, need to check explicitly @@ -1862,7 +1861,7 @@ struct TwoParticleCorrelationsMpi { 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) { + if (cfgGeneral.cfgPairCharge != 0 && cfgGeneral.cfgPairCharge * getParticleSign(track1) * getParticleSign(track2) < 0) { continue; } } @@ -1873,7 +1872,7 @@ struct TwoParticleCorrelationsMpi { if (cfg.mPairCuts && mPairCuts.conversionCuts(track1, track2)) { continue; } - if (cfgTwoTrackCut > 0 && mPairCuts.twoTrackCut(track1, track2, magField)) { + if (cfgGeneral.cfgTwoTrackCut > 0 && mPairCuts.twoTrackCut(track1, track2, magField)) { continue; } } @@ -1927,7 +1926,7 @@ struct TwoParticleCorrelationsMpi { return; } if (!cfgEfficiencyTrigger.value.empty()) { - if (cfgLocalEfficiency > 0) { + if (cfgGeneral.cfgLocalEfficiency > 0) { TFile* fEfficiencyTrigger = TFile::Open(cfgEfficiencyTrigger.value.c_str(), "READ"); cfg.mEfficiencyTrigger = dynamic_cast(fEfficiencyTrigger->Get("ccdb_object")); } else { @@ -1939,7 +1938,7 @@ struct TwoParticleCorrelationsMpi { LOGF(info, "Loaded efficiency histogram for trigger particles from %s", cfgEfficiencyTrigger.value.c_str()); } if (!cfgEfficiencyAssociated.value.empty()) { - if (cfgLocalEfficiency > 0) { + if (cfgGeneral.cfgLocalEfficiency > 0) { TFile* fEfficiencyAssociated = TFile::Open(cfgEfficiencyAssociated.value.c_str(), "READ"); cfg.mEfficiencyAssociated = dynamic_cast(fEfficiencyAssociated->Get("ccdb_object")); } else { @@ -2095,7 +2094,7 @@ struct TwoParticleCorrelationsMpi { loadCcdbYieldTemplates(collision.timestamp()); int field = 0; - if (cfgTwoTrackCut > 0) { + if (cfgGeneral.cfgTwoTrackCut > 0) { field = getMagneticField(collision.timestamp()); } @@ -2116,7 +2115,7 @@ struct TwoParticleCorrelationsMpi { } const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); - const bool fillReco = !(cfgDropStepRECO && hasEfficiency); + const bool fillReco = !(cfgGeneral.cfgDropStepRECO && hasEfficiency); EventSeedEstimate seedEstimate; if (cfgUserAxis == EventSeedAxis) { @@ -2275,7 +2274,7 @@ struct TwoParticleCorrelationsMpi { int bin = configurableBinningDerived.getBin(std::tuple(collision1.posZ(), multiplicity)); float eventWeight = 1.0f / it.currentWindowNeighbours(); int field = 0; - if (cfgTwoTrackCut > 0) { + if (cfgGeneral.cfgTwoTrackCut > 0) { field = getMagneticField(collision1.timestamp()); } @@ -2284,12 +2283,12 @@ struct TwoParticleCorrelationsMpi { } bool hasEfficiencyMixed = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); - bool fillRecoMixed = !(cfgDropStepRECO && hasEfficiencyMixed); + bool fillRecoMixed = !(cfgGeneral.cfgDropStepRECO && hasEfficiencyMixed); if (it.isNewWindow()) { loadEfficiency(collision1.timestamp()); hasEfficiencyMixed = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); - fillRecoMixed = !(cfgDropStepRECO && hasEfficiencyMixed); + fillRecoMixed = !(cfgGeneral.cfgDropStepRECO && hasEfficiencyMixed); if (cfgUserAxis == EventSeedAxis) { loadCcdbYieldTemplates(collision1.timestamp()); @@ -2357,7 +2356,7 @@ struct TwoParticleCorrelationsMpi { } auto multiplicity = mcCollision.multiplicity(); - if (cfgCentBinsForMC > 0) { + if (cfgGeneral.cfgCentBinsForMC > 0) { if (collisions.size() == 0) { return; } @@ -2381,7 +2380,7 @@ struct TwoParticleCorrelationsMpi { continue; } const int sign = getParticleSign(mcParticle); - if (sign == 0 || (cfgAssociatedCharge != 0 && cfgAssociatedCharge * sign < 0)) { + if (sign == 0 || (cfgGeneral.cfgAssociatedCharge != 0 && cfgGeneral.cfgAssociatedCharge * sign < 0)) { continue; } if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) { @@ -2432,7 +2431,7 @@ struct TwoParticleCorrelationsMpi { template void processMCMixedDerivedT(soa::Filtered const& mcCollisions, soa::Filtered const& collisions, ParticleTypes&&... particles) { - bool useMCMultiplicity = (cfgCentBinsForMC == 0); + bool useMCMultiplicity = (cfgGeneral.cfgCentBinsForMC == 0); auto getMultiplicity = [&collisions, &useMCMultiplicity, this](auto& col) { if (useMCMultiplicity) { From efe475a9105ff51ac7fda7c222b72948214fc819 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Tue, 15 Sep 2026 09:31:18 +0200 Subject: [PATCH 4/5] codecheck --- .../Tasks/twoParticleCorrelationsMpi.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index ccecdc2eefc..2a27e6df227 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -316,7 +316,7 @@ struct TwoParticleCorrelationsMpi { OutputObj fFCGen{"FlowContainer_gen"}; OutputObj fFCptGen{"FlowPtContainer_gen"}; std::unique_ptr fGFW{std::make_unique()}; - std::unique_ptr fPtAxis{}; + std::unique_ptr fPtAxis; std::unique_ptr fRndm{std::make_unique()}; std::vector mCorrConfigs; std::shared_ptr mAcceptanceWeights; @@ -725,7 +725,7 @@ struct TwoParticleCorrelationsMpi { void addConfigObjectsToObjArray(TObjArray* profiles, const std::vector& configs, const bool addSeedClasses) { - profiles->SetOwner(kTRUE); + profiles->SetOwner(true); for (const auto& config : configs) { if (config.pTDif) { const std::string suffix = "_ptDiff"; @@ -765,11 +765,11 @@ struct TwoParticleCorrelationsMpi { for (std::size_t configIndex = 0; configIndex < mCorrConfigs.size(); ++configIndex) { if (!mCorrConfigs.at(configIndex).pTDif) { - const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, kTRUE).real(); + const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, true).real(); if (dnx == 0) { continue; } - const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, kFALSE).real() / dnx; + const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), 0, false).real() / dnx; if (std::abs(val) < 1) { fillFlowProfile(mCorrConfigs.at(configIndex).Head, val, cfgEventWeight.cfgUseMultiplicityFlowWeights ? dnx : 1.0); flowPtContainer->fillVnPtProfiles(centMult, val, dnx, randomNumber, cfgCorrConfig->GetpTCorrMasks()[configIndex]); @@ -777,11 +777,11 @@ struct TwoParticleCorrelationsMpi { continue; } for (int i = 1; i <= fPtAxis->GetNbins(); i++) { - const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, kTRUE).real(); + const auto dnx = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, true).real(); if (dnx == 0) { continue; } - const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, kFALSE).real() / dnx; + const auto val = fGFW->Calculate(mCorrConfigs.at(configIndex), i - 1, false).real() / dnx; if (std::abs(val) < 1) { const std::string profileName = Form("%s_pt_%i", mCorrConfigs.at(configIndex).Head.c_str(), i); fillFlowProfile(profileName, val, cfgEventWeight.cfgUseMultiplicityFlowWeights ? dnx : 1.0); From a8082c71a237a66c717d41127e2591ae03536daf Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Tue, 15 Sep 2026 12:07:58 +0200 Subject: [PATCH 5/5] Added nMPI quantile axis for MCGen --- .../Tasks/twoParticleCorrelationsMpi.cxx | 187 +++++++++++++----- 1 file changed, 143 insertions(+), 44 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 2a27e6df227..39988d92474 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -138,8 +138,8 @@ struct TwoParticleCorrelationsMpi { Configurable cfgEventSeedPriorExposure{"cfgEventSeedPriorExposure", 1.f, "Gamma-prior strength in equivalent event exposures for MAP-EM"}; Configurable cfgEventSeedEMMaxIterations{"cfgEventSeedEMMaxIterations", 10, "Maximum number of weighted MAP-EM iterations per event"}; Configurable cfgEventSeedEMTolerance{"cfgEventSeedEMTolerance", 1.e-5f, "Relative component-yield convergence tolerance for MAP-EM"}; - Configurable> cfgEventSeedPercentileLower{"cfgEventSeedPercentileLower", {}, "Lower event-seed percentile boundary in each axisMultiplicity bin; empty keeps the raw Seed user axis"}; - Configurable> cfgEventSeedPercentileUpper{"cfgEventSeedPercentileUpper", {}, "Upper event-seed percentile boundary in each axisMultiplicity bin; empty keeps the raw Seed user axis"}; + Configurable> cfgEventClassifierPercentileLower{"cfgEventClassifierPercentileLower", {}, "Lower event-classifier percentile boundary in each axisMultiplicity bin; empty keeps the raw classifier user axis"}; + Configurable> cfgEventClassifierPercentileUpper{"cfgEventClassifierPercentileUpper", {}, "Upper event-classifier percentile boundary in each axisMultiplicity bin; empty keeps the raw classifier user axis"}; Configurable cfgNumMixedEvents{"cfgNumMixedEvents", 5, "Number of mixed events per event"}; @@ -147,7 +147,7 @@ struct TwoParticleCorrelationsMpi { Configurable cfgDecayParticleMask{"cfgDecayParticleMask", 0, "Selection bitmask for the decay particles: 0 = no selection"}; Configurable cfgV0RapidityMax{"cfgV0RapidityMax", 0.8, "Maximum rapidity for the decay particles (0 = no selection)"}; - Configurable cfgUserAxis{"cfgUserAxis", 0, "Additional user axis: 0 = OFF, 1 = invariant mass, 2 = event seed"}; + Configurable cfgUserAxis{"cfgUserAxis", 0, "Additional user axis: 0 = OFF, 1 = invariant mass, 2 = event seed, 3 = nMPI (MCGen only)"}; Configurable> cfgMcTriggerPDGs{"cfgMcTriggerPDGs", {}, "MC PDG codes to use exclusively as trigger particles and exclude from associated particles. Empty = no selection."}; ConfigurableAxis axisVertex{"axisVertex", {7, -7, 7}, "vertex axis for histograms"}; @@ -295,7 +295,7 @@ struct TwoParticleCorrelationsMpi { int pairAcceptanceSchemaVersion = 0; const TList* loadedCcdbYieldTemplateObject = nullptr; bool eventSeedEstimatorEnabled = false; - bool eventSeedPercentileAxisEnabled = false; + bool eventClassifierPercentileAxisEnabled = false; enum CorrelationMethod { All = 0, @@ -305,7 +305,8 @@ struct TwoParticleCorrelationsMpi { enum UserAxisMode { NoUserAxis = 0, InvariantMassAxis, - EventSeedAxis + EventSeedAxis, + NMPIAxis }; HistogramRegistry registry{"registry"}; PairCuts mPairCuts; @@ -336,8 +337,8 @@ struct TwoParticleCorrelationsMpi { void init(o2::framework::InitContext&) { - if (cfgUserAxis < NoUserAxis || cfgUserAxis > EventSeedAxis) { - LOGF(fatal, "Unsupported cfgUserAxis=%d; use 0 (off), 1 (invariant mass), or 2 (event seed)", cfgUserAxis.value); + if (cfgUserAxis < NoUserAxis || cfgUserAxis > NMPIAxis) { + LOGF(fatal, "Unsupported cfgUserAxis=%d; use 0 (off), 1 (invariant mass), 2 (event seed), or 3 (nMPI)", cfgUserAxis.value); } if (cfgGeneral.cfgCentBinsForMC < 0 || cfgGeneral.cfgCentBinsForMC > 1) { LOGF(fatal, "Unsupported cfgCentBinsForMC=%d; use 0 (generated multiplicity), 1 (reconstructed multiplicity and all associated collisions)", cfgGeneral.cfgCentBinsForMC.value); @@ -356,6 +357,13 @@ struct TwoParticleCorrelationsMpi { if (doprocessSameGenMC && doprocessMCSameDerived) { LOGF(fatal, "processSameGenMC and processMCSameDerived are mutually exclusive because both fill the generated same-event outputs"); } + if (cfgUserAxis == NMPIAxis && + (!(doprocessSameGenMC || doprocessMCSameDerived) || + doprocessSameAOD || enabledDerivedSameProcesses > 0 || + doprocessMixedAOD || enabledDerivedMixedProcesses > 0 || doprocessMCMixedDerived || + cfgGeneral.cfgCentBinsForMC != 0 || cfgFillAcceptanceWeights)) { + LOGF(fatal, "cfgUserAxis=3 supports generated MC same events only: enable processSameGenMC or processMCSameDerived and disable all other same/mixed processes"); + } if (cfgFillAcceptanceWeights && !cfgAcceptance.value.empty()) { LOGF(fatal, "cfgFillAcceptanceWeights and cfgAcceptance are mutually exclusive: produce and apply acceptance weights in separate jobs"); } @@ -378,34 +386,34 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis && !eventSeedEstimatorEnabled) { LOGF(fatal, "cfgUserAxis=2 requires an event-seed template configured through cfgNuncSeedsTemplateFile or cfgNuncSeedsTemplate"); } - const bool hasLowerPercentileBoundaries = !cfgEventSeedPercentileLower->empty(); - const bool hasUpperPercentileBoundaries = !cfgEventSeedPercentileUpper->empty(); + const bool hasLowerPercentileBoundaries = !cfgEventClassifierPercentileLower->empty(); + const bool hasUpperPercentileBoundaries = !cfgEventClassifierPercentileUpper->empty(); if (hasLowerPercentileBoundaries != hasUpperPercentileBoundaries) { - LOGF(fatal, "Configure both cfgEventSeedPercentileLower and cfgEventSeedPercentileUpper, or leave both empty"); + LOGF(fatal, "Configure both cfgEventClassifierPercentileLower and cfgEventClassifierPercentileUpper, or leave both empty"); } - eventSeedPercentileAxisEnabled = hasLowerPercentileBoundaries; - if (eventSeedPercentileAxisEnabled) { + eventClassifierPercentileAxisEnabled = hasLowerPercentileBoundaries; + if (eventClassifierPercentileAxisEnabled) { const int nMultiplicityBins = AxisSpec(axisMultiplicity).getNbins(); - if (static_cast(cfgEventSeedPercentileLower->size()) != nMultiplicityBins || - static_cast(cfgEventSeedPercentileUpper->size()) != nMultiplicityBins) { - LOGF(fatal, "Event-seed percentile boundary vectors must each contain exactly %d values, one per axisMultiplicity bin", nMultiplicityBins); + if (static_cast(cfgEventClassifierPercentileLower->size()) != nMultiplicityBins || + static_cast(cfgEventClassifierPercentileUpper->size()) != nMultiplicityBins) { + LOGF(fatal, "Event-classifier percentile boundary vectors must each contain exactly %d values, one per axisMultiplicity bin", nMultiplicityBins); } for (int multBin = 0; multBin < nMultiplicityBins; ++multBin) { - if (!std::isfinite(cfgEventSeedPercentileLower->at(multBin)) || - !std::isfinite(cfgEventSeedPercentileUpper->at(multBin)) || - cfgEventSeedPercentileLower->at(multBin) >= cfgEventSeedPercentileUpper->at(multBin)) { - LOGF(fatal, "Invalid event-seed percentile boundaries in multiplicity bin %d: lower=%g upper=%g", - multBin, cfgEventSeedPercentileLower->at(multBin), cfgEventSeedPercentileUpper->at(multBin)); + if (!std::isfinite(cfgEventClassifierPercentileLower->at(multBin)) || + !std::isfinite(cfgEventClassifierPercentileUpper->at(multBin)) || + cfgEventClassifierPercentileLower->at(multBin) >= cfgEventClassifierPercentileUpper->at(multBin)) { + LOGF(fatal, "Invalid event-classifier percentile boundaries in multiplicity bin %d: lower=%g upper=%g", + multBin, cfgEventClassifierPercentileLower->at(multBin), cfgEventClassifierPercentileUpper->at(multBin)); } } - if (cfgUserAxis == EventSeedAxis) { + if (cfgUserAxis == EventSeedAxis || cfgUserAxis == NMPIAxis) { const std::vector expectedEdges{-1.5, -0.5, 0.5, 1.5, 2.5}; const auto& configuredEdges = AxisSpec(axisUser).binEdges; const float threshold = 1.e-6; if (configuredEdges.size() != expectedEdges.size() || !std::equal(configuredEdges.begin(), configuredEdges.end(), expectedEdges.begin(), [&threshold](double lhs, double rhs) { return std::abs(lhs - rhs) < threshold; })) { - LOGF(fatal, "Percentile-class Seed axis requires axisUser edges {-1.5,-0.5,0.5,1.5,2.5}"); + LOGF(fatal, "Percentile-class event classifier axis requires axisUser edges {-1.5,-0.5,0.5,1.5,2.5}"); } } } @@ -584,7 +592,11 @@ struct TwoParticleCorrelationsMpi { userAxis.emplace_back(axisUser, "m (GeV/c^2)"); userMixingAxis.emplace_back(axisUser, "m (GeV/c^2)"); } else if (cfgUserAxis == EventSeedAxis) { - const char* title = eventSeedPercentileAxisEnabled ? "N_{seed} percentile class" : "N_{seed}"; + const char* title = eventClassifierPercentileAxisEnabled ? "N_{seed} percentile class" : "N_{seed}"; + userAxis.emplace_back(axisUser, title); + userMixingAxis.emplace_back(axisUser, title); + } else if (cfgUserAxis == NMPIAxis) { + const char* title = eventClassifierPercentileAxisEnabled ? "N_{MPI}^{true} percentile class" : "N_{MPI}^{true}"; userAxis.emplace_back(axisUser, title); userMixingAxis.emplace_back(axisUser, title); } @@ -645,7 +657,7 @@ struct TwoParticleCorrelationsMpi { if (!cfgFillAcceptanceWeights && (doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedMultSet)) { auto recoProfiles = std::make_unique(); - addConfigObjectsToObjArray(recoProfiles.get(), mCorrConfigs, eventSeedPercentileAxisEnabled); + addConfigObjectsToObjArray(recoProfiles.get(), mCorrConfigs, cfgUserAxis == EventSeedAxis && eventClassifierPercentileAxisEnabled); fFC.setObject(new FlowContainer("FlowContainer")); fFC->SetXAxis(fPtAxis.get()); fFC->Initialize(recoProfiles.get(), axisMultiplicity, FlowNBootstrap); @@ -659,7 +671,7 @@ struct TwoParticleCorrelationsMpi { } if (!cfgFillAcceptanceWeights && (doprocessMCSameDerived || doprocessSameGenMC)) { auto generatedProfiles = std::make_unique(); - addConfigObjectsToObjArray(generatedProfiles.get(), mCorrConfigs, false); + addConfigObjectsToObjArray(generatedProfiles.get(), mCorrConfigs, cfgUserAxis == NMPIAxis && eventClassifierPercentileAxisEnabled); fFCGen.setObject(new FlowContainer("FlowContainer_gen")); fFCGen->SetXAxis(fPtAxis.get()); fFCGen->Initialize(generatedProfiles.get(), axisMultiplicity, FlowNBootstrap); @@ -692,6 +704,18 @@ struct TwoParticleCorrelationsMpi { const char* getSeedClassSuffix(const SeedClass seedClass) const { + if (cfgUserAxis.value == NMPIAxis) { + switch (seedClass) { + case LowSeedClass: + return "__nMPILow"; + case MiddleSeedClass: + return "__nMPIMiddle"; + case HighSeedClass: + return "__nMPIHigh"; + default: + return ""; + } + } switch (seedClass) { case LowSeedClass: return "__seedLow"; @@ -748,11 +772,9 @@ struct TwoParticleCorrelationsMpi { auto& flowPtContainer = (dt == Gen) ? fFCptGen : fFCpt; const auto fillFlowProfile = [&](const std::string& profileName, const double value, const double weight) { flowContainer->FillProfile(profileName.c_str(), centMult, value, weight, randomNumber); - if constexpr (dt == Reco) { - if (seedClass != InvalidSeedClass) { - const auto classifiedName = getSeedClassProfileName(profileName, seedClass); - flowContainer->FillProfile(classifiedName.c_str(), centMult, value, weight, randomNumber); - } + if (seedClass != InvalidSeedClass) { + const auto classifiedName = getSeedClassProfileName(profileName, seedClass); + flowContainer->FillProfile(classifiedName.c_str(), centMult, value, weight, randomNumber); } }; @@ -1334,17 +1356,35 @@ struct TwoParticleCorrelationsMpi { double getEventSeedUserAxisValue(double multiplicity, double eventSeed) const { - if (!eventSeedPercentileAxisEnabled || eventSeed < 0.0) { + if (!eventClassifierPercentileAxisEnabled || eventSeed < 0.0) { return eventSeed; } const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); - if (multBin < 0 || multBin >= static_cast(cfgEventSeedPercentileLower->size())) { + if (multBin < 0 || multBin >= static_cast(cfgEventClassifierPercentileLower->size())) { + return -1.0; + } + if (eventSeed < cfgEventClassifierPercentileLower->at(multBin)) { + return 0.0; + } + if (eventSeed >= cfgEventClassifierPercentileUpper->at(multBin)) { + return 2.0; + } + return 1.0; + } + + double getTrueNMPIUserAxisValue(double multiplicity, double nMPI) const + { + if (!eventClassifierPercentileAxisEnabled || nMPI < 0.0) { + return nMPI; + } + const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); + if (multBin < 0 || multBin >= static_cast(cfgEventClassifierPercentileLower->size())) { return -1.0; } - if (eventSeed < cfgEventSeedPercentileLower->at(multBin)) { + if (nMPI < cfgEventClassifierPercentileLower->at(multBin)) { return 0.0; } - if (eventSeed >= cfgEventSeedPercentileUpper->at(multBin)) { + if (nMPI >= cfgEventClassifierPercentileUpper->at(multBin)) { return 2.0; } return 1.0; @@ -1352,7 +1392,7 @@ struct TwoParticleCorrelationsMpi { SeedClass getEventSeedClass(const double multiplicity, const double eventSeed) const { - if (!eventSeedPercentileAxisEnabled) { + if (!eventClassifierPercentileAxisEnabled) { return InvalidSeedClass; } const double userAxisValue = getEventSeedUserAxisValue(multiplicity, eventSeed); @@ -1364,6 +1404,20 @@ struct TwoParticleCorrelationsMpi { return static_cast(userAxisValue); } + SeedClass getTrueNMPIClass(const double multiplicity, const double nMPI) const + { + if (!eventClassifierPercentileAxisEnabled) { + return InvalidSeedClass; + } + const double userAxisValue = getTrueNMPIUserAxisValue(multiplicity, nMPI); + if (!std::isfinite(userAxisValue) || + userAxisValue < static_cast(LowSeedClass) || + userAxisValue >= static_cast(SeedClassCount)) { + return InvalidSeedClass; + } + return static_cast(userAxisValue); + } + double getPairAcceptance(double multiplicity, double deltaPhi, double deltaEta, double posZ) const { const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); @@ -1691,7 +1745,7 @@ struct TwoParticleCorrelationsMpi { } template - void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr, std::vector* pendingTriggerFills = nullptr, std::vector* pendingPairFills = nullptr, double eventSeed = -1.0, bool fillEstimatorAcceptanceQA = true, bool fillLoopQA = true) + void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr, std::vector* pendingTriggerFills = nullptr, std::vector* pendingPairFills = nullptr, double userAxisValue = -1.0, bool fillEstimatorAcceptanceQA = true, bool fillLoopQA = true) { if (multiplicity < 0.f) { return; @@ -1773,8 +1827,10 @@ struct TwoParticleCorrelationsMpi { if (pendingTriggerFills) { pendingTriggerFills->push_back({track1.pt(), multiplicity, posZ, triggerWeight}); } else { - target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, eventSeed, triggerWeight); + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, userAxisValue, triggerWeight); } + } else if (cfgUserAxis == NMPIAxis) { + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, userAxisValue, triggerWeight); } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, track1.invMass(), triggerWeight); @@ -1903,8 +1959,10 @@ struct TwoParticleCorrelationsMpi { if (pendingPairFills) { pendingPairFills->push_back({deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, associatedWeight}); } else { - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, eventSeed, associatedWeight); + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, userAxisValue, associatedWeight); } + } else if (cfgUserAxis == NMPIAxis) { + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, userAxisValue, associatedWeight); } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, track1.invMass(), associatedWeight); @@ -2044,6 +2102,9 @@ struct TwoParticleCorrelationsMpi { if (cfgFillAcceptanceWeights) { return; } + if (cfgUserAxis == NMPIAxis && mcCollision.nMPI() < 0) { + return; + } if (!cfgNuncSeedsTemplate.value.empty()) { for (const auto& collision : collisions) { loadCcdbYieldTemplates(collision.timestamp()); @@ -2052,6 +2113,29 @@ struct TwoParticleCorrelationsMpi { } const auto generatedMultiplicity = mcCollision.multiplicity(); + const double trueNMPIUserAxisValue = cfgUserAxis == NMPIAxis + ? getTrueNMPIUserAxisValue(generatedMultiplicity, mcCollision.nMPI()) + : -1.0; + if (cfgUserAxis == NMPIAxis) { + fGFW->Clear(); + fFCptGen->clearVector(); + for (const auto& mcParticle : mcParticles) { + if (!mcParticle.isPhysicalPrimary()) { + continue; + } + const int sign = getParticleSign(mcParticle); + if (sign == 0 || (cfgGeneral.cfgAssociatedCharge != 0 && cfgGeneral.cfgAssociatedCharge * sign < 0)) { + continue; + } + if (!cfgMcTriggerPDGs->empty() && + std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) { + continue; + } + fillPtSums(mcParticle, generatedMultiplicity, mcCollision.posZ()); + fillGFW(mcParticle, generatedMultiplicity, mcCollision.posZ()); + } + fillOutputContainers(generatedMultiplicity, fRndm->Rndm(), getTrueNMPIClass(generatedMultiplicity, mcCollision.nMPI())); + } fillContainerEvent(same, generatedMultiplicity, CorrelationContainer::kCFStepAll); EventSeedEstimate seedEstimate; std::vector pendingTriggerFills; @@ -2062,7 +2146,8 @@ struct TwoParticleCorrelationsMpi { } fillCorrelations(same, mcParticles, mcParticles, generatedMultiplicity, mcCollision.posZ(), 0, 1.0f, &seedEstimate, cfgUserAxis == EventSeedAxis ? &pendingTriggerFills : nullptr, - cfgUserAxis == EventSeedAxis ? &pendingPairFills : nullptr); + cfgUserAxis == EventSeedAxis ? &pendingPairFills : nullptr, + trueNMPIUserAxisValue); finalizeEventSeedEstimate(seedEstimate); if (cfgUserAxis == EventSeedAxis) { flushSeedAxisFills(same, pendingTriggerFills, pendingPairFills, seedEstimate.nuncSeeds()); @@ -2351,6 +2436,9 @@ struct TwoParticleCorrelationsMpi { if (cfgFillAcceptanceWeights) { return; } + if (cfgUserAxis == NMPIAxis && mcCollision.nMPI() < 0) { + return; + } if (cfgVerbosity > 0) { LOGF(info, "processMCSameDerivedT. MC collision: %d, particles1: %d, particles2: %d, collisions: %d", mcCollision.globalIndex(), mcParticles1.size(), mcParticles2.size(), collisions.size()); } @@ -2364,6 +2452,9 @@ struct TwoParticleCorrelationsMpi { multiplicity = collision.multiplicity(); } } + const double trueNMPIUserAxisValue = cfgUserAxis == NMPIAxis + ? getTrueNMPIUserAxisValue(multiplicity, mcCollision.nMPI()) + : -1.0; if (!(doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedCorrected || doprocessSameDerivedMultSet || doprocessSameDerivedMultSetCorrected)) { if constexpr (std::experimental::is_detected::value) { @@ -2389,23 +2480,28 @@ struct TwoParticleCorrelationsMpi { fillPtSums(mcParticle, multiplicity, mcCollision.posZ()); fillGFW(mcParticle, multiplicity, mcCollision.posZ()); } - fillOutputContainers(multiplicity, fRndm->Rndm()); + const SeedClass nMPIClass = cfgUserAxis == NMPIAxis ? getTrueNMPIClass(multiplicity, mcCollision.nMPI()) : InvalidSeedClass; + fillOutputContainers(multiplicity, fRndm->Rndm(), nMPIClass); fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepAll); - fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); + fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f, + nullptr, nullptr, nullptr, trueNMPIUserAxisValue); if (collisions.size() == 0) { return; } fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepVertex); - fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); + fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f, + nullptr, nullptr, nullptr, trueNMPIUserAxisValue); fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepTrackedOnlyPrim); - fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); + fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f, + nullptr, nullptr, nullptr, trueNMPIUserAxisValue); fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepTracked); - fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f); + fillCorrelations(same, mcParticles1, mcParticles2, multiplicity, mcCollision.posZ(), 0, 1.0f, + nullptr, nullptr, nullptr, trueNMPIUserAxisValue); // kCFStepReconstructed and kCFStepCorrected are filled below for every // reconstructed collision associated with this MC collision. @@ -2416,6 +2512,9 @@ struct TwoParticleCorrelationsMpi { 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); + if (cfgUserAxis == NMPIAxis) { + return; + } if (eventSeedEstimatorEnabled && collisions.size() == 0) { registry.fill(HIST("mcValidation/status"), 0.0); }