From a678aada0e7053b5b77af7fb1d90f5c58a8b6362 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Thu, 22 Oct 2026 23:27:03 +0200 Subject: [PATCH 1/4] Add region-aware initialisation mimicing GFW flow initialisation --- PWGCF/GenericFramework/Core/CMakeLists.txt | 58 +- .../GenericFramework/Core/FlowPtContainer.cxx | 935 +++++++++++++++++- PWGCF/GenericFramework/Core/FlowPtContainer.h | 60 +- .../Core/GenericFrameworkLinkDef.h | 4 + 4 files changed, 1003 insertions(+), 54 deletions(-) mode change 100755 => 100644 PWGCF/GenericFramework/Core/CMakeLists.txt diff --git a/PWGCF/GenericFramework/Core/CMakeLists.txt b/PWGCF/GenericFramework/Core/CMakeLists.txt old mode 100755 new mode 100644 index c82c4a5c751..fe3b5ab1c0e --- a/PWGCF/GenericFramework/Core/CMakeLists.txt +++ b/PWGCF/GenericFramework/Core/CMakeLists.txt @@ -1,35 +1,35 @@ -# Copyright 2019-2020 CERN and copyright holders of ALICE O2. -# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -# All rights not expressly granted are reserved. +#Copyright 2019 - 2020 CERN and copyright holders of ALICE O2. +#See https: // alice-o2.web.cern.ch/copyright for details of the copyright holders. +#All rights not expressly granted are reserved. # -# This software is distributed under the terms of the GNU General Public -# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +#This software is distributed under the terms of the GNU General Public +#License v3(GPL Version 3), copied verbatim in the file "COPYING". # -# In applying this license CERN does not waive the privileges and immunities -# granted to it by virtue of its status as an Intergovernmental Organization +#In applying this license CERN does not waive the privileges and immunities +#granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. o2physics_add_library(GFWCore - SOURCES GFWPowerArray.cxx - GFWCumulant.cxx - GFW.cxx - ProfileSubset.cxx - FlowContainer.cxx - GFWWeights.cxx - GFWWeightsList.cxx - FlowPtContainer.cxx - BootstrapProfile.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) + SOURCES GFWPowerArray.cxx + GFWCumulant.cxx + GFW.cxx + ProfileSubset.cxx + FlowContainer.cxx + GFWWeights.cxx + GFWWeightsList.cxx + FlowPtContainer.cxx + BootstrapProfile.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) -o2physics_target_root_dictionary(GFWCore - HEADERS GFWPowerArray.h - GFWCumulant.h - GFW.h - ProfileSubset.h - FlowContainer.h - GFWWeights.h - GFWWeightsList.h - GFWConfig.h - FlowPtContainer.h - BootstrapProfile.h - LINKDEF GenericFrameworkLinkDef.h) + o2physics_target_root_dictionary(GFWCore + HEADERS GFWPowerArray.h + GFWCumulant.h + GFW.h + ProfileSubset.h + FlowContainer.h + GFWWeights.h + GFWWeightsList.h + GFWConfig.h + FlowPtContainer.h + BootstrapProfile.h + LINKDEF GenericFrameworkLinkDef.h) diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx index d2d5a1740f6..040dc9ff0c6 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx @@ -31,16 +31,154 @@ #include #include +#include +#include #include #include #include #include +#include +#include +#include +#include +#include +#include #include +#include +#include +#include #include #include using namespace o2::analysis::genericframework; using namespace o2::analysis::genericframework::eventweight; +using o2::analysis::genericframework::flowptcontainer::PtCorrConfig; +using o2::analysis::genericframework::flowptcontainer::PtMomentConfig; +using o2::analysis::genericframework::flowptcontainer::PtMomentGroup; + +namespace +{ +constexpr std::size_t TwoParticleOrder = 2; +constexpr int ParityDivisor = 2; + +bool parseInteger(std::string_view token, int& value) +{ + if (!token.empty() && token.front() == '+') { + token.remove_prefix(1); + } + const auto result = std::from_chars(token.data(), token.data() + token.size(), value); + return !token.empty() && result.ec == std::errc{} && result.ptr == token.data() + token.size(); +} + +std::optional>> replacePtTerms(std::string& expression, const std::string& head) +{ + std::vector> isPt; + size_t cursor = 0; + while (true) { + auto openBrace = expression.find('{', cursor); + if (openBrace == std::string::npos) { + break; + } + if (openBrace > 0 && !std::isspace(static_cast(expression[openBrace - 1]))) { + expression.insert(openBrace++, 1, ' '); + } + const auto closeBrace = expression.find('}', openBrace + 1); + if (closeBrace == std::string::npos || expression.find('{', openBrace + 1) < closeBrace) { + LOGF(error, "Invalid pT correlation %s: unmatched brace", head.c_str()); + return std::nullopt; + } + + std::string terms = expression.substr(openBrace + 1, closeBrace - openBrace - 1); + std::replace(terms.begin(), terms.end(), ',', ' '); + std::replace(terms.begin(), terms.end(), ';', ' '); + std::istringstream stream(terms); + std::string replacement; + std::string token; + std::vector groupIsPt; + while (stream >> token) { + int harmonic = 0; + const bool isPt = token == "pt"; + if (!isPt && !parseInteger(token, harmonic)) { + LOGF(error, "Invalid pT correlation %s: unknown observable %s", head.c_str(), token.c_str()); + return std::nullopt; + } + if (!replacement.empty()) { + replacement += ' '; + } + replacement += isPt ? "0" : token; + groupIsPt.push_back(isPt); + } + if (groupIsPt.empty()) { + LOGF(error, "Invalid pT correlation %s: empty observable group", head.c_str()); + return std::nullopt; + } + + expression.replace(openBrace + 1, closeBrace - openBrace - 1, replacement); + cursor = openBrace + replacement.size() + 2; + isPt.push_back(std::move(groupIsPt)); + } + + if (isPt.empty()) { + LOGF(error, "Invalid pT correlation %s: no observable groups", head.c_str()); + return std::nullopt; + } + return isPt; +} + +std::string makeFullTopologyName(const std::string& head, std::uint64_t activeFactors, std::uint64_t allFactors, bool hasFlow, int ptOrder) +{ + if (activeFactors == allFactors) { + return head; + } + + std::string name = head + "_full"; + if (hasFlow && (activeFactors & std::uint64_t{1})) { + name += "_flow"; + } + const int ptOffset = hasFlow ? 1 : 0; + for (int i = 0; i < ptOrder; ++i) { + if (activeFactors & (std::uint64_t{1} << (ptOffset + i))) { + name += "_pt" + std::to_string(i + 1); + } + } + return name; +} + +double binomialCoefficient(int n, int k) +{ + k = std::min(k, n - k); + double result = 1.; + for (int i = 1; i <= k; ++i) { + result *= static_cast(n - k + i) / i; + } + return result; +} + +std::string makePostprocessName(std::string_view observable, const char* head, int ind) +{ + return std::string{head} + "_" + std::string{observable} + "_" + std::to_string(ind); +} + +bool sameMomentConfigs(const std::vector& first, const std::vector& second) +{ + if (first.size() != second.size()) { + return false; + } + for (size_t config = 0; config < first.size(); ++config) { + if (first[config].head != second[config].head || first[config].groups.size() != second[config].groups.size()) { + return false; + } + for (size_t group = 0; group < first[config].groups.size(); ++group) { + const auto& lhs = first[config].groups[group]; + const auto& rhs = second[config].groups[group]; + if (lhs.region != rhs.region || lhs.order != rhs.order || lhs.meanProfile != rhs.meanProfile) { + return false; + } + } + } + return true; +} +} // namespace FlowPtContainer::FlowPtContainer() : fCMTermList(nullptr), fCorrList(nullptr), @@ -54,7 +192,9 @@ FlowPtContainer::FlowPtContainer() : fCMTermList(nullptr), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), - fUseGap(false) + fUseGap(false), + fPtCentralMomentList(nullptr), + fPtCumulantList(nullptr) { } FlowPtContainer::~FlowPtContainer() @@ -66,6 +206,8 @@ FlowPtContainer::~FlowPtContainer() delete fSubCMList; delete fCumulantList; delete fCentralMomentList; + delete fPtCentralMomentList; + delete fPtCumulantList; }; FlowPtContainer::FlowPtContainer(const char* name) : TNamed(name, name), fCMTermList(nullptr), @@ -80,7 +222,9 @@ FlowPtContainer::FlowPtContainer(const char* name) : TNamed(name, name), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), - fUseGap(false) + fUseGap(false), + fPtCentralMomentList(nullptr), + fPtCumulantList(nullptr) { } FlowPtContainer::FlowPtContainer(const char* name, const char* title) : TNamed(name, title), @@ -96,8 +240,317 @@ FlowPtContainer::FlowPtContainer(const char* name, const char* title) : TNamed(n fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), - fUseGap(false) + fUseGap(false), + fPtCentralMomentList(nullptr), + fPtCumulantList(nullptr) +{ +} +void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWRegions& regions, const GFWCorrConfigs& configs, const int& nsub) { + const int numberOfRegions = regions.GetSize(); + const int numberOfConfigs = configs.GetSize(); + if (numberOfRegions < 0 || numberOfConfigs < 0) { + LOGF(error, "Cannot initialize %s: GFW configuration vectors have different sizes", GetName()); + return; + } + if (nsub < 0) { + LOGF(error, "Cannot initialize %s: number of subsamples cannot be negative", GetName()); + return; + } + + std::vector multiBins = axis.binEdges; + int numberOfBins = axis.nBins.value_or(0); + const bool isVariableBinning = numberOfBins <= 0; + if (isVariableBinning) { + numberOfBins = multiBins.size() - 1; + } + const size_t expectedEdges = isVariableBinning ? static_cast(numberOfBins + 1) : 2u; + if (numberOfBins <= 0 || multiBins.size() != expectedEdges) { + LOGF(error, "Cannot initialize %s: invalid multiplicity axis", GetName()); + return; + } + + GFW parser; + for (int i = 0; i < numberOfRegions; ++i) { + parser.AddRegion(regions.GetNames()[i], regions.GetEtaMin()[i], regions.GetEtaMax()[i], regions.GetpTDifs()[i] ? 2 : 1, regions.GetBitmasks()[i]); + } + if (parser.fRegions.size() != static_cast(numberOfRegions)) { + LOGF(error, "Cannot initialize %s: invalid GFW region", GetName()); + return; + } + + struct ParsedPtRequest { + GFW::CorrConfig config; + std::vector> isPt; + std::vector> factors; + std::vector> ptTerms; + int ptOrder; + bool hasFlow; + std::uint64_t allFactors; + }; + + std::vector requests; + requests.reserve(numberOfConfigs); + std::unordered_set correlationHeads; + std::vector regionMaxOrders(numberOfRegions, 0); + int maxPtOrder = 0; + for (int i = 0; i < numberOfConfigs; ++i) { + if (configs.GetpTDifs()[i] != 0) { + LOGF(error, "Invalid pT correlation %s: differential configurations are not supported", configs.GetHeads()[i].c_str()); + return; + } + std::string expression = configs.GetCorrs()[i]; + auto isPt = replacePtTerms(expression, configs.GetHeads()[i]); + if (!isPt) { + return; + } + auto config = parser.GetCorrelatorConfig(expression, configs.GetHeads()[i], configs.GetpTDifs()[i] != 0); + if (config.Regs.size() != isPt->size() || config.Hars.size() != isPt->size() || config.Overlap.size() != isPt->size() || config.ptInd.size() != isPt->size()) { + LOGF(error, "Invalid pT correlation %s", config.Head.c_str()); + return; + } + if (std::any_of(config.ptInd.begin(), config.ptInd.end(), [](int ptIndex) { return ptIndex >= 0; })) { + LOGF(error, "Invalid pT correlation %s: explicit pT-bin indices are not supported", config.Head.c_str()); + return; + } + int ptOrder = 0; + bool hasFlow = false; + std::vector> ptTerms; + std::vector> factors(isPt->size()); + for (size_t group = 0; group < isPt->size(); ++group) { + if (config.Regs[group].empty() || config.Hars[group].size() != isPt->at(group).size()) { + LOGF(error, "Invalid pT correlation %s", config.Head.c_str()); + return; + } + const bool groupHasPt = std::any_of(isPt->at(group).begin(), isPt->at(group).end(), [](bool value) { return value; }); + const bool groupHasFlow = std::any_of(isPt->at(group).begin(), isPt->at(group).end(), [](bool value) { return !value; }); + if (groupHasPt && groupHasFlow) { + LOGF(error, "Invalid pT correlation %s: mixed flow and pt terms in one group are not supported", config.Head.c_str()); + return; + } + if (groupHasPt && (config.Regs[group].size() != 1 || config.Overlap[group] >= 0)) { + LOGF(error, "Invalid pT correlation %s: pt groups must contain exactly one region without an explicit overlap", config.Head.c_str()); + return; + } + factors[group].resize(isPt->at(group).size()); + for (size_t term = 0; term < isPt->at(group).size(); ++term) { + if (isPt->at(group)[term]) { + ptTerms.emplace_back(group, term); + ++ptOrder; + } else { + hasFlow = true; + } + } + } + if (ptOrder == 0) { + LOGF(error, "Invalid pT correlation %s: configuration does not contain pt", config.Head.c_str()); + return; + } + if (config.Head.empty() || !correlationHeads.insert(config.Head).second) { + LOGF(error, "Cannot initialize %s: correlation heads must be non-empty and unique", GetName()); + return; + } + + const int factorCount = ptOrder + (hasFlow ? 1 : 0); + if (factorCount >= std::numeric_limits::digits) { + LOGF(error, "Invalid pT correlation %s: too many observable factors", config.Head.c_str()); + return; + } + const int ptOffset = hasFlow ? 1 : 0; + for (int ptIndex = 0; ptIndex < ptOrder; ++ptIndex) { + const auto [group, term] = ptTerms[ptIndex]; + factors[group][term] = ptOffset + ptIndex; + } + if (hasFlow) { + for (size_t group = 0; group < isPt->size(); ++group) { + if (!isPt->at(group).front()) { + std::fill(factors[group].begin(), factors[group].end(), 0); + } + } + } + const std::uint64_t allFactors = (std::uint64_t{1} << factorCount) - 1; + for (size_t group = 0; group < isPt->size(); ++group) { + if (!isPt->at(group).front()) { + continue; + } + const int order = isPt->at(group).size(); + for (const int& region : config.Regs[group]) { + regionMaxOrders[region] = std::max(regionMaxOrders[region], order); + } + if (config.Overlap[group] >= 0) { + regionMaxOrders[config.Overlap[group]] = std::max(regionMaxOrders[config.Overlap[group]], order); + } + } + + maxPtOrder = std::max(maxPtOrder, ptOrder); + requests.push_back({std::move(config), std::move(*isPt), std::move(factors), std::move(ptTerms), ptOrder, hasFlow, allFactors}); + } + + struct NaturalMean { + int region; + int ptIndex; + bool pTDif; + std::string profile; + }; + std::vector naturalMeans; + auto getPtRegion = [](const ParsedPtRequest& request, int ptIndex) { + const auto group = request.ptTerms[ptIndex].first; + return request.config.Regs[group].front(); + }; + auto findNaturalMean = [&](const ParsedPtRequest& request, int ptIndex) { + const auto group = request.ptTerms[ptIndex].first; + const int region = getPtRegion(request, ptIndex); + return std::find_if(naturalMeans.begin(), naturalMeans.end(), [&](const auto& mean) { + return mean.region == region && mean.ptIndex == request.config.ptInd[group] && mean.pTDif == request.config.pTDif; + }); + }; + + // Prefer explicitly requested one-particle observables as the shared native means, + // irrespective of their order in the configuration list. + for (const auto& request : requests) { + if (request.hasFlow || request.ptOrder != 1 || findNaturalMean(request, 0) != naturalMeans.end()) { + continue; + } + const auto group = request.ptTerms.front().first; + naturalMeans.push_back({getPtRegion(request, 0), request.config.ptInd[group], request.config.pTDif, request.config.Head}); + } + + std::vector parsedConfigs; + parsedConfigs.reserve(requests.size()); + std::vector momentConfigs; + momentConfigs.reserve(requests.size()); + for (const auto& request : requests) { + const int ptOffset = request.hasFlow ? 1 : 0; + PtCorrConfig parsed{request.config.Head, request.allFactors, request.hasFlow, {}, std::vector(request.ptOrder)}; + + // Factors in one single-region pT group are exchangeable. Their full-topology + // terms therefore depend only on how many factors are active. + std::vector> factorClasses; + if (request.hasFlow) { + factorClasses.push_back({0}); + } + for (size_t group = 0; group < request.isPt.size(); ++group) { + if (!request.isPt[group].front()) { + continue; + } + factorClasses.push_back(request.factors[group]); + } + + std::vector activeMasks{0}; + for (const auto& factors : factorClasses) { + const size_t previousSize = activeMasks.size(); + for (size_t i = 0; i < previousSize; ++i) { + std::uint64_t addedFactors = 0; + for (const int& factor : factors) { + addedFactors |= std::uint64_t{1} << factor; + activeMasks.push_back(activeMasks[i] | addedFactors); + } + } + } + std::sort(activeMasks.begin(), activeMasks.end()); + parsed.profiles.reserve(activeMasks.size() + request.ptOrder); + + // Inactive factors stay in the tuple as harmonic-zero multiplicity insertions, + // so every subset retains the event weight of the full correlation. + auto addFullTopologyProfile = [&](std::uint64_t activeFactors) { + auto profileConfig = request.config; + if (request.hasFlow && !(activeFactors & std::uint64_t{1})) { + for (size_t group = 0; group < request.isPt.size(); ++group) { + if (!request.isPt[group].front()) { + std::fill(profileConfig.Hars[group].begin(), profileConfig.Hars[group].end(), 0); + } + } + } + profileConfig.Head = makeFullTopologyName(request.config.Head, activeFactors, request.allFactors, request.hasFlow, request.ptOrder); + parsed.profiles.push_back({std::move(profileConfig), request.isPt, request.factors, activeFactors, true}); + }; + + addFullTopologyProfile(request.allFactors); + for (const auto& activeFactors : activeMasks) { + if (activeFactors != 0 && activeFactors != request.allFactors) { + addFullTopologyProfile(activeFactors); + } + } + + // Each distinct pT source has one naturally weighted mean. Explicit one-particle + // requests provide that mean directly; otherwise the first consumer books it. + for (int ptIndex = 0; ptIndex < request.ptOrder; ++ptIndex) { + const std::uint64_t factor = std::uint64_t{1} << (ptOffset + ptIndex); + if (factor == request.allFactors) { + parsed.meanProfiles[ptIndex] = request.config.Head; + continue; + } + auto existingMean = findNaturalMean(request, ptIndex); + if (existingMean != naturalMeans.end()) { + parsed.meanProfiles[ptIndex] = existingMean->profile; + continue; + } + + const auto group = request.ptTerms[ptIndex].first; + GFW::CorrConfig meanConfig; + meanConfig.Regs = {{getPtRegion(request, ptIndex)}}; + meanConfig.Hars = {{0}}; + meanConfig.Overlap = {-1}; + meanConfig.ptInd = {request.config.ptInd[group]}; + meanConfig.pTDif = request.config.pTDif; + meanConfig.Head = request.config.Head + "_mean_pt" + std::to_string(ptIndex + 1); + parsed.meanProfiles[ptIndex] = meanConfig.Head; + naturalMeans.push_back({meanConfig.Regs.front().front(), meanConfig.ptInd.front(), meanConfig.pTDif, meanConfig.Head}); + parsed.profiles.push_back({std::move(meanConfig), {{true}}, {{ptOffset + ptIndex}}, factor, false}); + } + + if (!request.hasFlow) { + PtMomentConfig momentConfig{request.config.Head, {}}; + momentConfig.groups.reserve(request.isPt.size()); + for (size_t group = 0; group < request.isPt.size(); ++group) { + const int firstFactor = request.factors[group].front(); + const int region = request.config.Regs[group].front(); + momentConfig.groups.push_back({regions.GetNames()[region], static_cast(request.factors[group].size()), parsed.meanProfiles[firstFactor]}); + } + momentConfigs.push_back(std::move(momentConfig)); + } + + parsedConfigs.push_back(std::move(parsed)); + } + + auto* profiles = new TList(); + profiles->SetOwner(kTRUE); + std::unordered_set profileNames; + for (const auto& parsed : parsedConfigs) { + for (const auto& profile : parsed.profiles) { + const auto& head = profile.config.Head; + if (!profileNames.insert(head).second) { + LOGF(error, "Cannot initialize %s: generated profile name %s is not unique", GetName(), head.c_str()); + delete profiles; + return; + } + if (isVariableBinning) { + profiles->Add(new BootstrapProfile(head.c_str(), head.c_str(), numberOfBins, multiBins.data())); + } else { + profiles->Add(new BootstrapProfile(head.c_str(), head.c_str(), numberOfBins, multiBins.front(), multiBins.back())); + } + } + } + if (nsub > 0) { + for (int i = 0; i < profiles->GetEntries(); ++i) { + dynamic_cast(profiles->At(i))->InitializeSubsamples(nsub); + } + } + + delete fCorrList; + fCorrList = profiles; + fPtRegions = regions; + fPtCorrConfigs = std::move(parsedConfigs); + fPtMomentConfigs = std::move(momentConfigs); + clearPtPostprocessLists(); + fPtRegionMoments.clear(); + fPtRegionMoments.reserve(regionMaxOrders.size()); + for (const int& order : regionMaxOrders) { + fPtRegionMoments.push_back({order, std::vector((order + 1) * (order + 1), 0.)}); + } + mpar = maxPtOrder; + nSubevents = numberOfRegions; + LOGF(info, "Container %s initialized with %i pT correlations, %i profiles and %i subsamples", GetName(), numberOfConfigs, profiles->GetEntries(), nsub); } void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const int& maxOrder, const GFWCorrConfigs& configs, const int& nsub) { @@ -618,6 +1071,93 @@ void FlowPtContainer::fill(const double& w, const double& pt) sumP[i] += std::pow(w, i % (mpar + 1)) * std::pow(pt, i / (mpar + 1)); } } +void FlowPtContainer::fill(const double& eta, const double& w, const double& pt, int mask) +{ + int maxOrder = 0; + for (const auto& moments : fPtRegionMoments) { + maxOrder = std::max(maxOrder, moments.maxOrder); + } + if (maxOrder == 0) { + return; + } + + std::vector weightPowers(maxOrder + 1, 1.); + std::vector ptPowers(maxOrder + 1, 1.); + for (int power = 1; power <= maxOrder; ++power) { + weightPowers[power] = weightPowers[power - 1] * w; + ptPowers[power] = ptPowers[power - 1] * pt; + } + + for (size_t region = 0; region < fPtRegionMoments.size(); ++region) { + auto& moments = fPtRegionMoments[region]; + if (moments.maxOrder == 0 || !(fPtRegions.GetEtaMin()[region] < eta && eta < fPtRegions.GetEtaMax()[region]) || + !(fPtRegions.GetBitmasks()[region] & mask)) { + continue; + } + for (int weightPower = 1; weightPower <= moments.maxOrder; ++weightPower) { + for (int ptPower = 0; ptPower <= weightPower; ++ptPower) { + moments.values[weightPower * (moments.maxOrder + 1) + ptPower] += weightPowers[weightPower] * ptPowers[ptPower]; + } + } + } +} +void FlowPtContainer::clearPtRegionMoments() +{ + for (auto& moments : fPtRegionMoments) { // o2-linter: disable=const-ref-in-for-loop (contents are reset) + std::fill(moments.values.begin(), moments.values.end(), 0.); + } +} +void FlowPtContainer::clearPtPostprocessLists() +{ + delete fPtCentralMomentList; + fPtCentralMomentList = nullptr; + delete fPtCumulantList; + fPtCumulantList = nullptr; +} +double FlowPtContainer::getPtMoment(int region, int weightPower, int ptPower) const +{ + if (region < 0 || region >= static_cast(fPtRegionMoments.size())) { + return 0.; + } + const auto& moments = fPtRegionMoments[region]; + if (weightPower < 1 || weightPower > moments.maxOrder || ptPower < 0 || ptPower > weightPower) { + return 0.; + } + return moments.values[weightPower * (moments.maxOrder + 1) + ptPower]; +} +double FlowPtContainer::recursivePtCorrelation(int poi, int ref, int overlap, std::vector& weightPowers, std::vector& ptPowers) const +{ + if (weightPowers.empty() || weightPowers.size() != ptPowers.size()) { + return 0.; + } + if (weightPowers.front() != 1 && overlap >= 0) { + poi = overlap; + } + if (weightPowers.size() == 1) { + return getPtMoment(poi, weightPowers.front(), ptPowers.front()); + } + if (weightPowers.size() == TwoParticleOrder) { + const double product = getPtMoment(poi, weightPowers[0], ptPowers[0]) * getPtMoment(ref, weightPowers[1], ptPowers[1]); + const double selfCorrelation = overlap >= 0 ? getPtMoment(overlap, weightPowers[0] + weightPowers[1], ptPowers[0] + ptPowers[1]) : 0.; + return product - selfCorrelation; + } + + const int lastWeightPower = weightPowers.back(); + const int lastPtPower = ptPowers.back(); + weightPowers.pop_back(); + ptPowers.pop_back(); + double result = recursivePtCorrelation(poi, ref, overlap, weightPowers, ptPowers) * getPtMoment(ref, lastWeightPower, lastPtPower); + for (int term = static_cast(weightPowers.size()) - 1; term >= 0; --term) { + weightPowers[term] += lastWeightPower; + ptPowers[term] += lastPtPower; + result -= recursivePtCorrelation(poi, ref, overlap, weightPowers, ptPowers); + weightPowers[term] -= lastWeightPower; + ptPowers[term] -= lastPtPower; + } + weightPowers.push_back(lastWeightPower); + ptPowers.push_back(lastPtPower); + return result; +} void FlowPtContainer::fillSub(const double& w, const double& pt, int subIndex) { for (size_t i = 0; i < insub[subIndex].size(); ++i) { @@ -689,6 +1229,57 @@ void FlowPtContainer::fillPtProfiles(const double& centmult, const double& rn) } } } +void FlowPtContainer::fillPtProfiles(GFW& gfw, const double& centmult, const double& rn) +{ + if (!fCorrList) { + return; + } + clearPtPostprocessLists(); + for (const auto& correlation : fPtCorrConfigs) { + for (const auto& profileConfig : correlation.profiles) { + std::complex numerator{1., 0.}; + double denominator = 1.; + for (size_t group = 0; group < profileConfig.config.Regs.size(); ++group) { + if (profileConfig.isPt[group].front()) { + const int poi = profileConfig.config.Regs[group].front(); + const int ref = profileConfig.config.Regs[group].size() > 1 ? profileConfig.config.Regs[group][1] : poi; + int overlap = profileConfig.config.Overlap[group]; + if (overlap < 0 && poi == ref) { + overlap = ref; + } + std::vector weightPowers(profileConfig.isPt[group].size(), 1); + std::vector ptPowers; + ptPowers.reserve(profileConfig.isPt[group].size()); + for (const int& factor : profileConfig.factors[group]) { + ptPowers.push_back((profileConfig.activeFactors & (std::uint64_t{1} << factor)) ? 1 : 0); + } + auto denominatorPtPowers = std::vector(ptPowers.size(), 0); + denominator *= recursivePtCorrelation(poi, ref, overlap, weightPowers, denominatorPtPowers); + numerator *= recursivePtCorrelation(poi, ref, overlap, weightPowers, ptPowers); + } else { + GFW::CorrConfig groupConfig; + groupConfig.Regs = {profileConfig.config.Regs[group]}; + groupConfig.Hars = {profileConfig.config.Hars[group]}; + groupConfig.Overlap = {profileConfig.config.Overlap[group]}; + groupConfig.ptInd = {profileConfig.config.ptInd[group]}; + groupConfig.pTDif = profileConfig.config.pTDif; + numerator *= gfw.Calculate(groupConfig, 0, false); + denominator *= gfw.Calculate(groupConfig, 0, true).real(); + } + } + if (!std::isfinite(numerator.real()) || !std::isfinite(numerator.imag()) || !std::isfinite(denominator) || denominator <= 0.) { + continue; + } + auto* profile = dynamic_cast(fCorrList->FindObject(profileConfig.config.Head.c_str())); + if (!profile) { + LOGF(error, "pT profile %s has not been booked", profileConfig.config.Head.c_str()); + continue; + } + const double eventWeight = fEventWeight == EventWeight::UnityWeight ? 1. : denominator; + profile->FillProfile(centmult, numerator.real() / denominator, eventWeight, rn); + } + } +} bool FlowPtContainer::addPtProfile(const char* name, int observableOrder) { if (!fCorrList || !name || (name[0] == 0) || observableOrder < 1 || observableOrder > mpar) { @@ -1748,6 +2339,7 @@ double FlowPtContainer::orderedAddition(std::vector vec) } void FlowPtContainer::rebinMulti(int nbins) { + clearPtPostprocessLists(); if (fCMTermList) { for (int i = 0; i < fCMTermList->GetEntries(); i++) { dynamic_cast(fCMTermList->At(i))->RebinMulti(nbins); @@ -1776,6 +2368,7 @@ void FlowPtContainer::rebinMulti(int nbins) } void FlowPtContainer::rebinMulti(int nbins, double* binedges) { + clearPtPostprocessLists(); if (fCMTermList) { for (int i = 0; i < fCMTermList->GetEntries(); i++) { dynamic_cast(fCMTermList->At(i))->RebinMulti(nbins, binedges); @@ -1806,6 +2399,295 @@ TH1* FlowPtContainer::getCorrHist(int ind, int m) { return dynamic_cast(fCorrList->FindObject(Form("mpt%i", m)))->getHist(ind); } +TH1* FlowPtContainer::getPtCentralMomentHist(const char* head, int ind) +{ + if (!head || head[0] == 0 || !fCorrList) { + return nullptr; + } + + const std::string resultName = makePostprocessName("ptcentral", head, ind); + if (fPtCentralMomentList) { + if (auto* result = dynamic_cast(fPtCentralMomentList->FindObject(resultName.c_str()))) { + return result; + } + } + + const auto config = std::find_if(fPtMomentConfigs.begin(), fPtMomentConfigs.end(), [&](const auto& candidate) { return candidate.head == head; }); + if (config == fPtMomentConfigs.end()) { + LOGF(error, "Cannot calculate pT central moment %s: no pure-pT postprocessing configuration was stored", head); + return nullptr; + } + + int order = 0; + for (const auto& group : config->groups) { + order += group.order; + } + if (order < 1 || order >= std::numeric_limits::digits) { + LOGF(error, "Cannot calculate pT central moment %s: invalid order %i", head, order); + return nullptr; + } + + auto getProfileHist = [&](const std::string& name) -> std::unique_ptr { + auto* profile = dynamic_cast(fCorrList->FindObject(name.c_str())); + if (!profile) { + LOGF(error, "Cannot calculate pT central moment %s: profile %s is missing", head, name.c_str()); + return nullptr; + } + std::unique_ptr histogram{profile->getHist(ind)}; + if (histogram) { + histogram->SetDirectory(nullptr); + } + return histogram; + }; + + auto shape = getProfileHist(config->head); + if (!shape) { + return nullptr; + } + auto* result = dynamic_cast(shape->Clone(resultName.c_str())); + result->SetDirectory(nullptr); + result->Reset(); + + const std::uint64_t allFactors = (std::uint64_t{1} << order) - 1; + std::vector activeCounts(config->groups.size(), 0); + bool finished = false; + while (!finished) { + int activeOrder = 0; + int factorOffset = 0; + double coefficient = 1.; + std::uint64_t activeFactors = 0; + std::map meanPowers; + for (size_t groupIndex = 0; groupIndex < config->groups.size(); ++groupIndex) { + const auto& group = config->groups[groupIndex]; + const int active = activeCounts[groupIndex]; + activeOrder += active; + coefficient *= binomialCoefficient(group.order, active); + for (int factor = 0; factor < active; ++factor) { + activeFactors |= std::uint64_t{1} << (factorOffset + factor); + } + if (active < group.order) { + meanPowers[group.meanProfile] += group.order - active; + } + factorOffset += group.order; + } + + std::unique_ptr term; + if (activeOrder == 0) { + term.reset(dynamic_cast(shape->Clone())); + term->SetDirectory(nullptr); + term->Reset(); + for (int bin = 1; bin <= term->GetNbinsX(); ++bin) { + if (shape->GetBinContent(bin) != 0. || shape->GetBinError(bin) != 0.) { + term->SetBinContent(bin, 1.); + } + } + } else { + term = getProfileHist(makeFullTopologyName(config->head, activeFactors, allFactors, false, order)); + if (!term) { + delete result; + return nullptr; + } + term->SetDirectory(nullptr); + } + + for (const auto& [meanName, power] : meanPowers) { + auto mean = getProfileHist(meanName); + if (!mean) { + delete result; + return nullptr; + } + mean->SetDirectory(nullptr); + if (power > 1) { + for (int bin = 0; bin <= mean->GetNbinsX() + 1; ++bin) { + const double value = mean->GetBinContent(bin); + const double error = mean->GetBinError(bin); + mean->SetBinContent(bin, std::pow(value, power)); + mean->SetBinError(bin, std::abs(power * std::pow(value, power - 1) * error)); + } + } + term->Multiply(mean.get()); + } + + if ((order - activeOrder) % ParityDivisor != 0) { + coefficient = -coefficient; + } + result->Add(term.get(), coefficient); + + size_t groupIndex = 0; + for (; groupIndex < activeCounts.size(); ++groupIndex) { + if (++activeCounts[groupIndex] <= config->groups[groupIndex].order) { + break; + } + activeCounts[groupIndex] = 0; + } + finished = groupIndex == activeCounts.size(); + } + + if (!fPtCentralMomentList) { + fPtCentralMomentList = new TList(); + fPtCentralMomentList->SetOwner(kTRUE); + } + fPtCentralMomentList->Add(result); + return result; +} +TH1* FlowPtContainer::getPtCumulantHist(const char* head, int ind) +{ + if (!head || head[0] == 0 || !fCorrList) { + return nullptr; + } + + const std::string resultName = makePostprocessName("ptcumulant", head, ind); + if (fPtCumulantList) { + if (auto* result = dynamic_cast(fPtCumulantList->FindObject(resultName.c_str()))) { + return result; + } + } + + const auto config = std::find_if(fPtMomentConfigs.begin(), fPtMomentConfigs.end(), [&](const auto& candidate) { return candidate.head == head; }); + if (config == fPtMomentConfigs.end()) { + LOGF(error, "Cannot calculate pT cumulant %s: no pure-pT postprocessing configuration was stored", head); + return nullptr; + } + + std::vector factorGroups; + for (size_t groupIndex = 0; groupIndex < config->groups.size(); ++groupIndex) { + factorGroups.insert(factorGroups.end(), config->groups[groupIndex].order, groupIndex); + } + const int order = factorGroups.size(); + if (order < 1) { + return nullptr; + } + + if (order == 1) { + auto* profile = dynamic_cast(fCorrList->FindObject(config->groups.front().meanProfile.c_str())); + if (!profile) { + LOGF(error, "Cannot calculate pT cumulant %s: profile %s is missing", head, config->groups.front().meanProfile.c_str()); + return nullptr; + } + auto* result = profile->getHist(ind); + if (!result) { + return nullptr; + } + result->SetName(resultName.c_str()); + result->SetDirectory(nullptr); + if (!fPtCumulantList) { + fPtCumulantList = new TList(); + fPtCumulantList->SetOwner(kTRUE); + } + fPtCumulantList->Add(result); + return result; + } + + auto* centralMoment = getPtCentralMomentHist(head, ind); + if (!centralMoment) { + return nullptr; + } + auto* result = dynamic_cast(centralMoment->Clone(resultName.c_str())); + result->SetDirectory(nullptr); + result->Reset(); + + using Topology = std::vector>; + auto topology = [](const PtMomentConfig& candidate, const std::vector* counts = nullptr) { + Topology value; + for (size_t groupIndex = 0; groupIndex < candidate.groups.size(); ++groupIndex) { + const int count = counts ? counts->at(groupIndex) : candidate.groups[groupIndex].order; + if (count > 0) { + value.emplace_back(candidate.groups[groupIndex].region, count); + } + } + std::sort(value.begin(), value.end()); + return value; + }; + auto findMomentConfig = [&](const std::vector& counts) { + const auto required = topology(*config, &counts); + return std::find_if(fPtMomentConfigs.begin(), fPtMomentConfigs.end(), [&](const auto& candidate) { return topology(candidate) == required; }); + }; + auto describeTopology = [&](const std::vector& counts) { + std::string description; + for (const auto& [region, count] : topology(*config, &counts)) { + if (!description.empty()) { + description += " "; + } + description += region + "{" + std::to_string(count) + "pt}"; + } + return description; + }; + + bool failed = false; + std::vector blockForFactor(order, 0); + std::function addPartitions = [&](int factor, int highestBlock) { + if (failed) { + return; + } + if (factor < order) { + for (int block = 0; block <= highestBlock + 1; ++block) { + blockForFactor[factor] = block; + addPartitions(factor + 1, std::max(highestBlock, block)); + } + return; + } + + const int numberOfBlocks = highestBlock + 1; + std::vector> blockCounts(numberOfBlocks, std::vector(config->groups.size(), 0)); + std::vector blockSizes(numberOfBlocks, 0); + for (int factorIndex = 0; factorIndex < order; ++factorIndex) { + const int block = blockForFactor[factorIndex]; + ++blockCounts[block][factorGroups[factorIndex]]; + ++blockSizes[block]; + } + if (std::any_of(blockSizes.begin(), blockSizes.end(), [](int size) { return size == 1; })) { + return; + } + + std::unique_ptr term; + for (int block = 0; block < numberOfBlocks; ++block) { + const PtMomentConfig* blockConfig = &*config; + if (numberOfBlocks > 1) { + const auto found = findMomentConfig(blockCounts[block]); + if (found == fPtMomentConfigs.end()) { + LOGF(error, "Cannot calculate pT cumulant %s: missing lower-order correlation %s", head, describeTopology(blockCounts[block]).c_str()); + failed = true; + return; + } + blockConfig = &*found; + } + auto* blockMoment = getPtCentralMomentHist(blockConfig->head.c_str(), ind); + if (!blockMoment) { + failed = true; + return; + } + if (!term) { + term.reset(dynamic_cast(blockMoment->Clone())); + term->SetDirectory(nullptr); + } else { + term->Multiply(blockMoment); + } + } + + double coefficient = 1.; + for (int i = 2; i < numberOfBlocks; ++i) { + coefficient *= i; + } + if (numberOfBlocks % ParityDivisor == 0) { + coefficient = -coefficient; + } + result->Add(term.get(), coefficient); + }; + + blockForFactor.front() = 0; + addPartitions(1, 0); + if (failed) { + delete result; + return nullptr; + } + + if (!fPtCumulantList) { + fPtCumulantList = new TList(); + fPtCumulantList->SetOwner(kTRUE); + } + fPtCumulantList->Add(result); + return result; +} TH1* FlowPtContainer::getCentralMomentHist(int ind, int m) { if (!fCentralMomentList) { @@ -1897,27 +2779,41 @@ void FlowPtContainer::calculateCumulantHists(std::vector inh, int ind) } Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name/function-variable (Keep name consistent with ROOT streamer merging) { - if (!fCorrList || !fCMTermList) { - return 0; - } + delete fCumulantList; + fCumulantList = nullptr; + delete fCentralMomentList; + fCentralMomentList = nullptr; + clearPtPostprocessLists(); + Long64_t nmerged = 0; TIter allPTC(collist); FlowPtContainer* lPTC = nullptr; while ((lPTC = dynamic_cast(allPTC()))) { + if (!lPTC->fPtMomentConfigs.empty()) { + if (fPtMomentConfigs.empty()) { + fPtMomentConfigs = lPTC->fPtMomentConfigs; + } else if (!sameMomentConfigs(fPtMomentConfigs, lPTC->fPtMomentConfigs)) { + LOGF(error, "Cannot merge FlowPtContainer %s: pT postprocessing configurations differ", GetName()); + continue; + } + } else if (!fPtMomentConfigs.empty()) { + LOGF(error, "Cannot merge FlowPtContainer %s: incoming pT postprocessing configuration is missing", GetName()); + continue; + } + TList* tCMTerm = lPTC->fCMTermList; TList* tCorr = lPTC->fCorrList; TList* tCov = lPTC->fCovList; TList* tSub = lPTC->fSubList; TList* tSubCM = lPTC->fSubCMList; - TList* tCum = lPTC->fCumulantList; - TList* tCM = lPTC->fCentralMomentList; + bool merged = false; if (tCMTerm) { if (!fCMTermList) { fCMTermList = dynamic_cast(tCMTerm->Clone()); } else { mergeBSLists(fCMTermList, tCMTerm); } - nmerged++; + merged = true; } if (tCorr) { if (!fCorrList) { @@ -1925,6 +2821,7 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fCorrList, tCorr); } + merged = true; } if (tCov) { if (!fCovList) { @@ -1932,20 +2829,7 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fCovList, tCov); } - } - if (tCum) { - if (!fCumulantList) { - fCumulantList = dynamic_cast(tCum->Clone()); - } else { - mergeBSLists(fCumulantList, tCum); - } - } - if (tCM) { - if (!fCentralMomentList) { - fCentralMomentList = dynamic_cast(tCM->Clone()); - } else { - mergeBSLists(fCentralMomentList, tCM); - } + merged = true; } if (tSub) { if (!fSubList) { @@ -1953,6 +2837,7 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fSubList, tSub); } + merged = true; } if (tSubCM) { if (!fSubCMList) { @@ -1960,6 +2845,10 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fSubCMList, tSubCM); } + merged = true; + } + if (merged) { + ++nmerged; } } return nmerged; diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.h b/PWGCF/GenericFramework/Core/FlowPtContainer.h index 904153b7887..147babc095a 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.h +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.h @@ -16,6 +16,7 @@ #ifndef PWGCF_GENERICFRAMEWORK_CORE_FLOWPTCONTAINER_H_ #define PWGCF_GENERICFRAMEWORK_CORE_FLOWPTCONTAINER_H_ +#include "GFW.h" #include "GFWConfig.h" #include @@ -31,6 +32,7 @@ #include #include #include +#include #include #include @@ -42,6 +44,43 @@ enum EventWeight { }; }; +namespace o2::analysis::genericframework::flowptcontainer +{ +struct PtProfileConfig { + GFW::CorrConfig config; + std::vector> isPt; + std::vector> factors; + std::uint64_t activeFactors; // Flow, when present, is bit 0; pT terms follow in expression order. + bool usesFullTopology; +}; + +struct PtCorrConfig { + std::string head; + std::uint64_t allFactors; + bool hasFlow; + std::vector profiles; + std::vector meanProfiles; +}; + +struct PtRegionMoments { + int maxOrder; + std::vector values; +}; + +struct PtMomentGroup { + std::string region; + int order; + std::string meanProfile; + ClassDefNV(PtMomentGroup, 1); +}; + +struct PtMomentConfig { + std::string head; + std::vector groups; + ClassDefNV(PtMomentConfig, 1); +}; +} // namespace o2::analysis::genericframework::flowptcontainer + class FlowPtContainer : public TNamed { public: @@ -50,6 +89,7 @@ class FlowPtContainer : public TNamed explicit FlowPtContainer(const char* name); ~FlowPtContainer() override; FlowPtContainer(const char* name, const char* title); + void initialise(const o2::framework::AxisSpec& axis, const o2::analysis::genericframework::GFWRegions& regions, const o2::analysis::genericframework::GFWCorrConfigs& configs, const int& nsub = 10); void initialise(const o2::framework::AxisSpec& axis, const int& maxOrder, const o2::analysis::genericframework::GFWCorrConfigs& configs, const int& nsub = 10); void initialise(int nbinsx, double* xbins, const int& maxOrder, const o2::analysis::genericframework::GFWCorrConfigs& configs, const int& nsub = 10); void initialise(int nbinsx, double xlow, double xhigh, const int& maxOrder, const o2::analysis::genericframework::GFWCorrConfigs& configs, const int& nsub = 10); @@ -58,6 +98,7 @@ class FlowPtContainer : public TNamed void initialiseSubevent(int nbinsx, double* xbins, const int& maxOrder, const int& nsubev = 2, const int& nsub = 10); void initialiseSubevent(int nbinsx, double xlow, double xhigh, const int& maxOrder, const int& nsubev = 2, const int& nsub = 10); void fill(const double& w, const double& pt); + void fill(const double& eta, const double& w, const double& pt, int mask); void fillSub(const double& w, const double& pt, int subIndex); void fillSub1(const double& w, const double& pt) { fillSub(w, pt, 0); } void fillSub2(const double& w, const double& pt) { fillSub(w, pt, nSubevents - 1); } @@ -68,6 +109,7 @@ class FlowPtContainer : public TNamed void calculateSubeventCorrelations(); void calculateCMTerms(); void fillPtProfiles(const double& centmult, const double& rn); + void fillPtProfiles(GFW& gfw, const double& centmult, const double& rn); // Book and fill a separate pT observable with a weight calculated by the task. // The task can supply a separate event weight for the same observable. bool addPtProfile(const char* name, int observableOrder); @@ -133,6 +175,9 @@ class FlowPtContainer : public TNamed TH1* getCentralMomentHist(int ind, int m); TH1* getCumulantHist(int ind, int m); TH1* getCorrHist(int ind, int m); + // Postprocess a pure-pT configuration by head. ind = -1 selects the main profile. + TH1* getPtCentralMomentHist(const char* head, int ind = -1); + TH1* getPtCumulantHist(const char* head, int ind = -1); int getMpar() { return mpar; } Long64_t Merge(TCollection* collist); double orderedAddition(std::vector vec); @@ -157,6 +202,7 @@ class FlowPtContainer : public TNamed arr.resize(3 * 3 * 5 * 5, {0.0, 0.0}); warr.clear(); warr.resize(3 * 3 * 5 * 5, 0.0); + clearPtRegionMoments(); }; TList* fCMTermList; @@ -229,10 +275,20 @@ class FlowPtContainer : public TNamed private: static const int centralMomentMaxOrder = 4; - std::vector> subevents; + o2::analysis::genericframework::GFWRegions fPtRegions; //! + std::vector fPtCorrConfigs; //! + std::vector fPtRegionMoments; //! + std::vector fPtMomentConfigs; + TList* fPtCentralMomentList; //! + TList* fPtCumulantList; //! + std::vector> subevents; //! + void clearPtPostprocessLists(); + void clearPtRegionMoments(); + double getPtMoment(int region, int weightPower, int ptPower) const; + double recursivePtCorrelation(int poi, int ref, int overlap, std::vector& weightPowers, std::vector& ptPowers) const; void getSubevents(int k, int n, std::vector& current, std::vector>& subevents); static constexpr std::array FactorialArray = {1., 1., 2., 6., 24., 120., 720., 5040., 40320.}; static constexpr std::array SignArray = {1, -1, 1, -1, 1, -1, 1, -1, 1}; - ClassDefOverride(FlowPtContainer, 2); + ClassDefOverride(FlowPtContainer, 3); }; #endif // PWGCF_GENERICFRAMEWORK_CORE_FLOWPTCONTAINER_H_ diff --git a/PWGCF/GenericFramework/Core/GenericFrameworkLinkDef.h b/PWGCF/GenericFramework/Core/GenericFrameworkLinkDef.h index e3907fec165..a9a3da8c05e 100644 --- a/PWGCF/GenericFramework/Core/GenericFrameworkLinkDef.h +++ b/PWGCF/GenericFramework/Core/GenericFrameworkLinkDef.h @@ -27,3 +27,7 @@ #pragma link C++ class o2::analysis::genericframework::GFWBinningCuts + ; #pragma link C++ class o2::analysis::genericframework::GFWRegions + ; #pragma link C++ class o2::analysis::genericframework::GFWCorrConfigs + ; +#pragma link C++ class o2::analysis::genericframework::flowptcontainer::PtMomentGroup + ; +#pragma link C++ class o2::analysis::genericframework::flowptcontainer::PtMomentConfig + ; +#pragma link C++ class std::vector < o2::analysis::genericframework::flowptcontainer::PtMomentGroup> + ; +#pragma link C++ class std::vector < o2::analysis::genericframework::flowptcontainer::PtMomentConfig> + ; From 3ce06e68f33768823c3684959948dd81462c748b Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Thu, 22 Oct 2026 23:43:56 +0200 Subject: [PATCH 2/4] revert(cmake): restore GFWCore formatting --- PWGCF/GenericFramework/Core/CMakeLists.txt | 58 +++++++++++----------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/PWGCF/GenericFramework/Core/CMakeLists.txt b/PWGCF/GenericFramework/Core/CMakeLists.txt index fe3b5ab1c0e..c82c4a5c751 100644 --- a/PWGCF/GenericFramework/Core/CMakeLists.txt +++ b/PWGCF/GenericFramework/Core/CMakeLists.txt @@ -1,35 +1,35 @@ -#Copyright 2019 - 2020 CERN and copyright holders of ALICE O2. -#See https: // alice-o2.web.cern.ch/copyright for details of the copyright holders. -#All rights not expressly granted are reserved. +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. # -#This software is distributed under the terms of the GNU General Public -#License v3(GPL Version 3), copied verbatim in the file "COPYING". +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". # -#In applying this license CERN does not waive the privileges and immunities -#granted to it by virtue of its status as an Intergovernmental Organization +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. o2physics_add_library(GFWCore - SOURCES GFWPowerArray.cxx - GFWCumulant.cxx - GFW.cxx - ProfileSubset.cxx - FlowContainer.cxx - GFWWeights.cxx - GFWWeightsList.cxx - FlowPtContainer.cxx - BootstrapProfile.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) + SOURCES GFWPowerArray.cxx + GFWCumulant.cxx + GFW.cxx + ProfileSubset.cxx + FlowContainer.cxx + GFWWeights.cxx + GFWWeightsList.cxx + FlowPtContainer.cxx + BootstrapProfile.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) - o2physics_target_root_dictionary(GFWCore - HEADERS GFWPowerArray.h - GFWCumulant.h - GFW.h - ProfileSubset.h - FlowContainer.h - GFWWeights.h - GFWWeightsList.h - GFWConfig.h - FlowPtContainer.h - BootstrapProfile.h - LINKDEF GenericFrameworkLinkDef.h) +o2physics_target_root_dictionary(GFWCore + HEADERS GFWPowerArray.h + GFWCumulant.h + GFW.h + ProfileSubset.h + FlowContainer.h + GFWWeights.h + GFWWeightsList.h + GFWConfig.h + FlowPtContainer.h + BootstrapProfile.h + LINKDEF GenericFrameworkLinkDef.h) From aa46acdc6c0dbd72b0f10516e85b98b32556e471 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Thu, 22 Oct 2026 23:49:00 +0200 Subject: [PATCH 3/4] fix linter --- .../GenericFramework/Core/FlowPtContainer.cxx | 38 ++++++++++--------- PWGCF/GenericFramework/Core/FlowPtContainer.h | 12 +++--- 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx index 040dc9ff0c6..885efe557c4 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx @@ -97,16 +97,16 @@ std::optional>> replacePtTerms(std::string& expres std::vector groupIsPt; while (stream >> token) { int harmonic = 0; - const bool isPt = token == "pt"; - if (!isPt && !parseInteger(token, harmonic)) { + const bool tokenIsPt = token == "pt"; + if (!tokenIsPt && !parseInteger(token, harmonic)) { LOGF(error, "Invalid pT correlation %s: unknown observable %s", head.c_str(), token.c_str()); return std::nullopt; } if (!replacement.empty()) { replacement += ' '; } - replacement += isPt ? "0" : token; - groupIsPt.push_back(isPt); + replacement += tokenIsPt ? "0" : token; + groupIsPt.push_back(tokenIsPt); } if (groupIsPt.empty()) { LOGF(error, "Invalid pT correlation %s: empty observable group", head.c_str()); @@ -247,9 +247,11 @@ FlowPtContainer::FlowPtContainer(const char* name, const char* title) : TNamed(n } void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWRegions& regions, const GFWCorrConfigs& configs, const int& nsub) { - const int numberOfRegions = regions.GetSize(); - const int numberOfConfigs = configs.GetSize(); - if (numberOfRegions < 0 || numberOfConfigs < 0) { + const auto numberOfRegions = regions.GetNames().size(); + const auto numberOfConfigs = configs.GetCorrs().size(); + if (regions.GetEtaMin().size() != numberOfRegions || regions.GetEtaMax().size() != numberOfRegions || + regions.GetpTDifs().size() != numberOfRegions || regions.GetBitmasks().size() != numberOfRegions || + configs.GetHeads().size() != numberOfConfigs || configs.GetpTDifs().size() != numberOfConfigs) { LOGF(error, "Cannot initialize %s: GFW configuration vectors have different sizes", GetName()); return; } @@ -271,10 +273,10 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR } GFW parser; - for (int i = 0; i < numberOfRegions; ++i) { + for (size_t i = 0; i < numberOfRegions; ++i) { parser.AddRegion(regions.GetNames()[i], regions.GetEtaMin()[i], regions.GetEtaMax()[i], regions.GetpTDifs()[i] ? 2 : 1, regions.GetBitmasks()[i]); } - if (parser.fRegions.size() != static_cast(numberOfRegions)) { + if (parser.fRegions.size() != numberOfRegions) { LOGF(error, "Cannot initialize %s: invalid GFW region", GetName()); return; } @@ -284,9 +286,9 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR std::vector> isPt; std::vector> factors; std::vector> ptTerms; - int ptOrder; - bool hasFlow; - std::uint64_t allFactors; + int ptOrder{}; + bool hasFlow{}; + std::uint64_t allFactors{}; }; std::vector requests; @@ -294,7 +296,7 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR std::unordered_set correlationHeads; std::vector regionMaxOrders(numberOfRegions, 0); int maxPtOrder = 0; - for (int i = 0; i < numberOfConfigs; ++i) { + for (size_t i = 0; i < numberOfConfigs; ++i) { if (configs.GetpTDifs()[i] != 0) { LOGF(error, "Invalid pT correlation %s: differential configurations are not supported", configs.GetHeads()[i].c_str()); return; @@ -387,9 +389,9 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR } struct NaturalMean { - int region; - int ptIndex; - bool pTDif; + int region{}; + int ptIndex{}; + bool pTDif{}; std::string profile; }; std::vector naturalMeans; @@ -549,8 +551,8 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR fPtRegionMoments.push_back({order, std::vector((order + 1) * (order + 1), 0.)}); } mpar = maxPtOrder; - nSubevents = numberOfRegions; - LOGF(info, "Container %s initialized with %i pT correlations, %i profiles and %i subsamples", GetName(), numberOfConfigs, profiles->GetEntries(), nsub); + nSubevents = static_cast(numberOfRegions); + LOGF(info, "Container %s initialized with %zu pT correlations, %i profiles and %i subsamples", GetName(), numberOfConfigs, profiles->GetEntries(), nsub); } void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const int& maxOrder, const GFWCorrConfigs& configs, const int& nsub) { diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.h b/PWGCF/GenericFramework/Core/FlowPtContainer.h index 147babc095a..9a079f8dae4 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.h +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.h @@ -50,26 +50,26 @@ struct PtProfileConfig { GFW::CorrConfig config; std::vector> isPt; std::vector> factors; - std::uint64_t activeFactors; // Flow, when present, is bit 0; pT terms follow in expression order. - bool usesFullTopology; + std::uint64_t activeFactors{}; // Flow, when present, is bit 0; pT terms follow in expression order. + bool usesFullTopology{}; }; struct PtCorrConfig { std::string head; - std::uint64_t allFactors; - bool hasFlow; + std::uint64_t allFactors{}; + bool hasFlow{}; std::vector profiles; std::vector meanProfiles; }; struct PtRegionMoments { - int maxOrder; + int maxOrder{}; std::vector values; }; struct PtMomentGroup { std::string region; - int order; + int order{}; std::string meanProfile; ClassDefNV(PtMomentGroup, 1); }; From 5f29d36d84ba760188e25474afacf1d0867d9427 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Fri, 23 Oct 2026 00:10:06 +0200 Subject: [PATCH 4/4] codecheck fix --- .../GenericFramework/Core/FlowPtContainer.cxx | 68 ++++++++++++++----- PWGCF/GenericFramework/Core/FlowPtContainer.h | 2 +- 2 files changed, 51 insertions(+), 19 deletions(-) diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx index 885efe557c4..d76b8a6bb39 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx @@ -16,6 +16,7 @@ #include "FlowPtContainer.h" #include "PWGCF/GenericFramework/Core/BootstrapProfile.h" +#include "PWGCF/GenericFramework/Core/GFW.h" #include "PWGCF/GenericFramework/Core/GFWConfig.h" #include @@ -45,6 +46,7 @@ #include #include #include +#include #include #include #include @@ -54,7 +56,6 @@ using namespace o2::analysis::genericframework; using namespace o2::analysis::genericframework::eventweight; using o2::analysis::genericframework::flowptcontainer::PtCorrConfig; using o2::analysis::genericframework::flowptcontainer::PtMomentConfig; -using o2::analysis::genericframework::flowptcontainer::PtMomentGroup; namespace { @@ -79,7 +80,7 @@ std::optional>> replacePtTerms(std::string& expres if (openBrace == std::string::npos) { break; } - if (openBrace > 0 && !std::isspace(static_cast(expression[openBrace - 1]))) { + if (openBrace > 0 && std::isspace(static_cast(expression[openBrace - 1])) == 0) { expression.insert(openBrace++, 1, ' '); } const auto closeBrace = expression.find('}', openBrace + 1); @@ -132,12 +133,12 @@ std::string makeFullTopologyName(const std::string& head, std::uint64_t activeFa } std::string name = head + "_full"; - if (hasFlow && (activeFactors & std::uint64_t{1})) { + if (hasFlow && (activeFactors & std::uint64_t{1}) != 0) { name += "_flow"; } const int ptOffset = hasFlow ? 1 : 0; for (int i = 0; i < ptOrder; ++i) { - if (activeFactors & (std::uint64_t{1} << (ptOffset + i))) { + if ((activeFactors & (std::uint64_t{1} << (ptOffset + i))) != 0) { name += "_pt" + std::to_string(i + 1); } } @@ -274,7 +275,7 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR GFW parser; for (size_t i = 0; i < numberOfRegions; ++i) { - parser.AddRegion(regions.GetNames()[i], regions.GetEtaMin()[i], regions.GetEtaMax()[i], regions.GetpTDifs()[i] ? 2 : 1, regions.GetBitmasks()[i]); + parser.AddRegion(regions.GetNames()[i], regions.GetEtaMin()[i], regions.GetEtaMax()[i], regions.GetpTDifs()[i] != 0 ? 2 : 1, regions.GetBitmasks()[i]); } if (parser.fRegions.size() != numberOfRegions) { LOGF(error, "Cannot initialize %s: invalid GFW region", GetName()); @@ -385,7 +386,13 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR } maxPtOrder = std::max(maxPtOrder, ptOrder); - requests.push_back({std::move(config), std::move(*isPt), std::move(factors), std::move(ptTerms), ptOrder, hasFlow, allFactors}); + requests.push_back({.config = std::move(config), + .isPt = std::move(*isPt), + .factors = std::move(factors), + .ptTerms = std::move(ptTerms), + .ptOrder = ptOrder, + .hasFlow = hasFlow, + .allFactors = allFactors}); } struct NaturalMean { @@ -414,7 +421,10 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR continue; } const auto group = request.ptTerms.front().first; - naturalMeans.push_back({getPtRegion(request, 0), request.config.ptInd[group], request.config.pTDif, request.config.Head}); + naturalMeans.push_back({.region = getPtRegion(request, 0), + .ptIndex = request.config.ptInd[group], + .pTDif = request.config.pTDif, + .profile = request.config.Head}); } std::vector parsedConfigs; @@ -423,7 +433,12 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR momentConfigs.reserve(requests.size()); for (const auto& request : requests) { const int ptOffset = request.hasFlow ? 1 : 0; - PtCorrConfig parsed{request.config.Head, request.allFactors, request.hasFlow, {}, std::vector(request.ptOrder)}; + PtCorrConfig parsed{ + .head = request.config.Head, + .allFactors = request.allFactors, + .hasFlow = request.hasFlow, + .profiles = {}, + .meanProfiles = std::vector(request.ptOrder)}; // Factors in one single-region pT group are exchangeable. Their full-topology // terms therefore depend only on how many factors are active. @@ -456,7 +471,7 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR // so every subset retains the event weight of the full correlation. auto addFullTopologyProfile = [&](std::uint64_t activeFactors) { auto profileConfig = request.config; - if (request.hasFlow && !(activeFactors & std::uint64_t{1})) { + if (request.hasFlow && (activeFactors & std::uint64_t{1}) == 0) { for (size_t group = 0; group < request.isPt.size(); ++group) { if (!request.isPt[group].front()) { std::fill(profileConfig.Hars[group].begin(), profileConfig.Hars[group].end(), 0); @@ -464,7 +479,11 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR } } profileConfig.Head = makeFullTopologyName(request.config.Head, activeFactors, request.allFactors, request.hasFlow, request.ptOrder); - parsed.profiles.push_back({std::move(profileConfig), request.isPt, request.factors, activeFactors, true}); + parsed.profiles.push_back({.config = std::move(profileConfig), + .isPt = request.isPt, + .factors = request.factors, + .activeFactors = activeFactors, + .usesFullTopology = true}); }; addFullTopologyProfile(request.allFactors); @@ -497,17 +516,28 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR meanConfig.pTDif = request.config.pTDif; meanConfig.Head = request.config.Head + "_mean_pt" + std::to_string(ptIndex + 1); parsed.meanProfiles[ptIndex] = meanConfig.Head; - naturalMeans.push_back({meanConfig.Regs.front().front(), meanConfig.ptInd.front(), meanConfig.pTDif, meanConfig.Head}); - parsed.profiles.push_back({std::move(meanConfig), {{true}}, {{ptOffset + ptIndex}}, factor, false}); + naturalMeans.push_back({.region = meanConfig.Regs.front().front(), + .ptIndex = meanConfig.ptInd.front(), + .pTDif = meanConfig.pTDif, + .profile = meanConfig.Head}); + parsed.profiles.push_back({.config = std::move(meanConfig), + .isPt = {{true}}, + .factors = {{ptOffset + ptIndex}}, + .activeFactors = factor, + .usesFullTopology = false}); } if (!request.hasFlow) { - PtMomentConfig momentConfig{request.config.Head, {}}; + PtMomentConfig momentConfig{ + .head = request.config.Head, + .groups = {}}; momentConfig.groups.reserve(request.isPt.size()); for (size_t group = 0; group < request.isPt.size(); ++group) { const int firstFactor = request.factors[group].front(); const int region = request.config.Regs[group].front(); - momentConfig.groups.push_back({regions.GetNames()[region], static_cast(request.factors[group].size()), parsed.meanProfiles[firstFactor]}); + momentConfig.groups.push_back({.region = regions.GetNames()[region], + .order = static_cast(request.factors[group].size()), + .meanProfile = parsed.meanProfiles[firstFactor]}); } momentConfigs.push_back(std::move(momentConfig)); } @@ -548,7 +578,8 @@ void FlowPtContainer::initialise(const o2::framework::AxisSpec& axis, const GFWR fPtRegionMoments.clear(); fPtRegionMoments.reserve(regionMaxOrders.size()); for (const int& order : regionMaxOrders) { - fPtRegionMoments.push_back({order, std::vector((order + 1) * (order + 1), 0.)}); + fPtRegionMoments.push_back({.maxOrder = order, + .values = std::vector((order + 1) * (order + 1), 0.)}); } mpar = maxPtOrder; nSubevents = static_cast(numberOfRegions); @@ -1092,8 +1123,9 @@ void FlowPtContainer::fill(const double& eta, const double& w, const double& pt, for (size_t region = 0; region < fPtRegionMoments.size(); ++region) { auto& moments = fPtRegionMoments[region]; - if (moments.maxOrder == 0 || !(fPtRegions.GetEtaMin()[region] < eta && eta < fPtRegions.GetEtaMax()[region]) || - !(fPtRegions.GetBitmasks()[region] & mask)) { + const bool inEtaRange = fPtRegions.GetEtaMin()[region] < eta && eta < fPtRegions.GetEtaMax()[region]; + const bool maskMatches = (fPtRegions.GetBitmasks()[region] & mask) != 0; + if (moments.maxOrder == 0 || !inEtaRange || !maskMatches) { continue; } for (int weightPower = 1; weightPower <= moments.maxOrder; ++weightPower) { @@ -1253,7 +1285,7 @@ void FlowPtContainer::fillPtProfiles(GFW& gfw, const double& centmult, const dou std::vector ptPowers; ptPowers.reserve(profileConfig.isPt[group].size()); for (const int& factor : profileConfig.factors[group]) { - ptPowers.push_back((profileConfig.activeFactors & (std::uint64_t{1} << factor)) ? 1 : 0); + ptPowers.push_back((profileConfig.activeFactors & (std::uint64_t{1} << factor)) != 0 ? 1 : 0); } auto denominatorPtPowers = std::vector(ptPowers.size(), 0); denominator *= recursivePtCorrelation(poi, ref, overlap, weightPowers, denominatorPtPowers); diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.h b/PWGCF/GenericFramework/Core/FlowPtContainer.h index 9a079f8dae4..6bd33833f7e 100644 --- a/PWGCF/GenericFramework/Core/FlowPtContainer.h +++ b/PWGCF/GenericFramework/Core/FlowPtContainer.h @@ -284,7 +284,7 @@ class FlowPtContainer : public TNamed std::vector> subevents; //! void clearPtPostprocessLists(); void clearPtRegionMoments(); - double getPtMoment(int region, int weightPower, int ptPower) const; + [[nodiscard]] double getPtMoment(int region, int weightPower, int ptPower) const; double recursivePtCorrelation(int poi, int ref, int overlap, std::vector& weightPowers, std::vector& ptPowers) const; void getSubevents(int k, int n, std::vector& current, std::vector>& subevents); static constexpr std::array FactorialArray = {1., 1., 2., 6., 24., 120., 720., 5040., 40320.};