diff --git a/PWGCF/GenericFramework/Core/CMakeLists.txt b/PWGCF/GenericFramework/Core/CMakeLists.txt old mode 100755 new mode 100644 diff --git a/PWGCF/GenericFramework/Core/FlowPtContainer.cxx b/PWGCF/GenericFramework/Core/FlowPtContainer.cxx index d2d5a1740f6..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 @@ -31,16 +32,154 @@ #include #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; + +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])) == 0) { + 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 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 += tokenIsPt ? "0" : token; + groupIsPt.push_back(tokenIsPt); + } + 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}) != 0) { + name += "_flow"; + } + const int ptOffset = hasFlow ? 1 : 0; + for (int i = 0; i < ptOrder; ++i) { + if ((activeFactors & (std::uint64_t{1} << (ptOffset + i))) != 0) { + 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 +193,9 @@ FlowPtContainer::FlowPtContainer() : fCMTermList(nullptr), fillCounter(0), fEventWeight(EventWeight::UnityWeight), fUseCentralMoments(true), - fUseGap(false) + fUseGap(false), + fPtCentralMomentList(nullptr), + fPtCumulantList(nullptr) { } FlowPtContainer::~FlowPtContainer() @@ -66,6 +207,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 +223,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 +241,349 @@ 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 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; + } + 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 (size_t i = 0; i < numberOfRegions; ++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()); + 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 (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; + } + 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({.config = std::move(config), + .isPt = std::move(*isPt), + .factors = std::move(factors), + .ptTerms = std::move(ptTerms), + .ptOrder = ptOrder, + .hasFlow = hasFlow, + .allFactors = 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({.region = getPtRegion(request, 0), + .ptIndex = request.config.ptInd[group], + .pTDif = request.config.pTDif, + .profile = 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{ + .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. + 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}) == 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); + } + } + } + profileConfig.Head = makeFullTopologyName(request.config.Head, activeFactors, request.allFactors, request.hasFlow, request.ptOrder); + parsed.profiles.push_back({.config = std::move(profileConfig), + .isPt = request.isPt, + .factors = request.factors, + .activeFactors = activeFactors, + .usesFullTopology = 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({.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{ + .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({.region = regions.GetNames()[region], + .order = static_cast(request.factors[group].size()), + .meanProfile = 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({.maxOrder = order, + .values = std::vector((order + 1) * (order + 1), 0.)}); + } + mpar = maxPtOrder; + 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) { @@ -618,6 +1104,94 @@ 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]; + 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) { + 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 +1263,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)) != 0 ? 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 +2373,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 +2402,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 +2433,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 +2813,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 +2855,7 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fCorrList, tCorr); } + merged = true; } if (tCov) { if (!fCovList) { @@ -1932,20 +2863,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 +2871,7 @@ Long64_t FlowPtContainer::Merge(TCollection* collist) // o2-linter: disable=name } else { mergeBSLists(fSubList, tSub); } + merged = true; } if (tSubCM) { if (!fSubCMList) { @@ -1960,6 +2879,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..6bd33833f7e 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(); + [[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.}; 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> + ;