diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index 8f193a0dabb..39d21e6c32f 100644 --- a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx @@ -75,12 +75,23 @@ struct RadialFlowDecorr { static constexpr int KIntM = 4; // pT-moment order used in the sums (need m up to 3 for c3) static constexpr int KIntK = 4; // weight-power order used in the sums (need k up to 3 for c3) - static constexpr int KNEtaHalfBinsMax = 8; // 0.8 / 0.1: finest (DATA) bins per side + static constexpr int KNEtaHalfBinsMax = 10; // 1.0 / 0.1: finest bins per side (axes span |eta| < KEtaFineMax) static constexpr int KNEtaMax = 2 * KNEtaHalfBinsMax + 1; // + index-0 full-range reference bin - - // Bootstrap: at most KMaxBoot replicas of 8 profiles each, in one registry of at most - // 512 objects. nBoot = min(cfgNBootstrap, KMaxBoot) replicas are filled. + // Observable pT ranges (GeV/c), all filled in the same mean and fluctuation passes with every variation: + // 0 = nominal (0.2, 3], 1 = (0.5, 3] (suffix _pt1), 2 = (0.2, 5] (suffix _pt2). Tracks and MC particles are + // selected in the map range [KPtMin, KPtMax] in every pass; each range is applied by ptRangeMask. + static constexpr int KNPtMax = 3; + static constexpr int KNFlatAxes = 5; // flattening weight maps: (vz, sign, pT, eta, phi); raw maps add the weight-type axis + static constexpr std::size_t KMinEdges = 2; // a binning needs at least two edges + static constexpr std::size_t KNAlignHalves = 2; // FT0 / FV0 alignment: one entry per detector half + static constexpr float KFwdTruthEtaMin = 1.5f; // forward truth signal: charged primaries with |eta| >= this + static constexpr std::array, KNPtMax> KPtRanges = {{{0.2f, 3.f}, {0.5f, 3.f}, {0.2f, 5.f}}}; + + // Bootstrap: at most KMaxBoot replicas; range 0 in histosBoot, the further pT ranges in histosBoot2, each registry + // of at most KMaxRegistryObjects objects. bootstrap.n = min(cfgNBootstrap, KMaxBoot), further capped in init() + // so that the replicas fit the registries. static constexpr int KMaxBoot = 60; + static constexpr int KMaxRegistryObjects = 512; static constexpr float KFloatEpsilon = 1e-6f; static constexpr float KEtaEdgeTolerance = 1e-3f; // slack for cfgCutEta being an integer multiple of the bin width @@ -97,6 +108,7 @@ struct RadialFlowDecorr { static constexpr float KEtaMin = -1.2f; static constexpr float KEtaMax = 1.2f; static constexpr int KNbinsPhi = 64; + static constexpr int KNbinsPhiFlat = 32; // flattening maps // azimuthal harmonics n = 1 ... KNHarm (index ih = n - 1) of the pair maps. C_TA = (C_T + 2 C_A)/3 = c_0 + sum_n beta_n c_n // uses n = 1, 5, 7 (beta_n = 0 for every n divisible by 2 or 3); the others are kept for the region-ratio cross-check static constexpr int KNHarm = 7; @@ -284,7 +296,7 @@ struct RadialFlowDecorr { // DCAxy and DCAz are null checks: with the nominal isPVContributor requirement even // 5 sigma DCAxy or |DCAz| < 0.05 cm removes only ~0.2% of the selected tracks // Eff : efficiency scaled by -10% at all pT (cancels in [pT]: a check) - // EffpTDep : efficiency scaled by -10% at cfgPtMin, linearly to 0 at 10 GeV/c + // EffpTDep : efficiency scaled by -10% at 0.2 GeV/c (KEffVarPtAnchor), linearly to 0 at 10 GeV/c // Flat : DATA: time-split acceptance weights. The Base flattening maps are filled separately for runs // < cfgFlatRunSplit (early) and >= cfgFlatRunSplit (late); each half is weighted with the map of // the other half. cfgFlatRunSplit = 0: the variation is dropped. MC: no flattening weights @@ -340,6 +352,7 @@ struct RadialFlowDecorr { static constexpr float KSystDcaXyScale = 5.f / 7.f; // DCAxy: 5 sigma (nominal 7 sigma) static constexpr float KSystDcaMaxZ = 1.f; // DCAz: |DCAz| < 1 cm (nominal 2 cm) static constexpr float KEffVarPercent = -10.f; // Eff, EffpTDep: efficiency scaled by -10% + static constexpr float KEffVarPtAnchor = 0.2f; // EffpTDep: full variation at 0.2 GeV/c (lowest observable range edge) static constexpr float KEffVarPtZero = 10.f; // EffpTDep: variation falls linearly to 0 at 10 GeV/c static constexpr int KSystITScluster = 4; // NITS: itsNCls >= 4 (nominal 5) static constexpr int KSystTPCcluster = 80; // NTPC: tpcNClsFound >= 80 (nominal 70) @@ -378,9 +391,31 @@ struct RadialFlowDecorr { // Observable eta binning, filled at init(): index 0 = full-range reference bin. std::vector etaLw; std::vector etaUp; - int nEta = 9; // active eta-bin count (set in init()) - int nBoot = 0; // active bootstrap samples (set in init()) - bool doBoot = false; // bootstrap active for this run (base data fluc only) + int nEta = 9; // active eta-bin count (set in init()) + // Observable pT ranges (set in init()). Profiles with an eta-bin index axis carry the range as + // blocks of that axis, index ipt * (nEta + 1) + ieta (etaIdx); the (eta, eta) and (gap, sum) maps of + // range ipt > 0 are separate profiles with the suffix _pt. + struct { + int n = 1; // number of ranges + std::array lo{}; + std::array hi{}; + } ptRange; + int etaIdx(int ipt, int ieta) const { return ipt * (nEta + 1) + ieta; } + uint32_t ptRangeMask(float pt) const + { + uint32_t m = 0; + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + if (pt > ptRange.lo[ipt] && pt <= ptRange.hi[ipt]) { + SETBIT(m, ipt); + } + } + return m; + } + static std::string ptSuffix(int ipt) { return ipt == 0 ? std::string{} : "_pt" + std::to_string(ipt); } + struct { + int n = 0; // active bootstrap samples (set in init()) + bool on = false; // bootstrap active for this run (base data fluc only) + } bootstrap; // Variations filled by each pass (set in init()): maskEff in the MC efficiency pass, // maskFlat in the flattening passes, maskAll in the mean and fluctuation passes. @@ -434,12 +469,8 @@ struct RadialFlowDecorr { } cfgPileup; // Track selection (nominal), on top of isGlobalTrack (Run 3 globalTrack, incl. its 7 sigma pT-dependent DCAxy cut) - Configurable cfgPtMin{"cfgPtMin", 0.2f, "min pT for observables"}; - Configurable cfgPtMax{"cfgPtMax", 5.0f, "max pT for observables"}; - Configurable cfgCutPtLower{"cfgCutPtLower", 0.2f, "Lower pT cut (track selection)"}; - Configurable cfgCutPtUpper{"cfgCutPtUpper", 10.0f, "Higher pT cut (track selection)"}; - Configurable cfgCutEta{"cfgCutEta", 0.8f, "absolute Eta cut"}; - Configurable cfgEtaBinWidth{"cfgEtaBinWidth", 0.1f, "DATA only: width of each narrow observable eta bin (eta units); cfgCutEta must be an integer multiple so eta=0 is a bin edge. MC is pinned to 0.2."}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "|eta| cut of the reco tracks and of the reference eta bin; MC truth particles are kept to |eta| < 1, the edge of all eta axes"}; + Configurable cfgEtaBinWidth{"cfgEtaBinWidth", 0.1f, "Width of each narrow observable eta bin (eta units), MC and DATA; MC bins span |eta| < 1, DATA bins |eta| < cfgCutEta; the width must divide both. Use the same value in the mean and fluctuation passes."}; Configurable cfgCutTracKDcaMaxZ{"cfgCutTracKDcaMaxZ", 2.0f, "Maximum |DCAz| (cm)"}; Configurable cfgITScluster{"cfgITScluster", 5, "Nominal: itsNCls >= value"}; Configurable cfgTPCcluster{"cfgTPCcluster", 70, "Nominal: tpcNClsFound >= value"}; @@ -448,7 +479,7 @@ struct RadialFlowDecorr { // Systematic variations (names stay flat in the JSON) struct : ConfigurableGroup { - Configurable cfgSystMask{"cfgSystMask", static_cast(KAllSystMask), "Bitmask of active variations, bit i = ESystType i: 0 Base (always on), 1 Cent, 2 Chi2TPC, 3 DCAxy, 4 DCAz, 5 Eff, 6 EffpTDep, 7 Flat, 8 ITSLayer, 9 NITS, 10 NTPC, 11 Pileup, 12 Vz, 13 NchBin2, 14 Occupancy, 15 PileupNSigma, 16 VzCentred, 17 VzPos, 18 VzNeg"}; + Configurable cfgSystMask{"cfgSystMask", static_cast(KAllSystMask), "Bitmask of active variations, bit i = ESystType i: 0 Base (always on in the efficiency and flattening passes; may be off in mean/fluctuation split groups, the bootstrap needs it), 1 Cent, 2 Chi2TPC, 3 DCAxy, 4 DCAz, 5 Eff, 6 EffpTDep, 7 Flat, 8 ITSLayer, 9 NITS, 10 NTPC, 11 Pileup, 12 Vz, 13 NchBin2, 14 Occupancy, 15 PileupNSigma, 16 VzCentred, 17 VzPos, 18 VzNeg"}; Configurable cfgSystVzCentre{"cfgSystVzCentre", 0.f, "VzCentred: centre of the |z_vtx - centre| < 8 cm window (cm), e.g. the mean vertex of the sample; 0: VzCentred off"}; Configurable cfgDcaXyP0{"cfgDcaXyP0", 0.0105f, "globalTrack DCAxy cut p0 + p1 / pT^p2 (trackselection dcaSetup=0; dcaSetup=1: 0.004, 0.013, 1)"}; Configurable cfgDcaXyP1{"cfgDcaXyP1", 0.035f, "globalTrack DCAxy cut p1"}; @@ -500,14 +531,20 @@ struct RadialFlowDecorr { AxisSpec nchAxisNchBin2{1, 0., 1., "N_{ch}"}; // NchBin2 (bins of nchWidthNchBin2) AxisSpec centAxisEff{{0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.}, "Centrality (%)"}; // efficiency maps - AxisSpec vzAxis{5, -12.5, 12.5, "Vz"}; - AxisSpec chgAxis{3, -1.5, 1.5}; - AxisSpec pTAxis{{0.0, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 7, 10}, "pT Axis"}; - AxisSpec phiAxis{KNbinsPhi, KPhiMin, TwoPI, "#phi"}; - - // etaFlatAxis is the *fixed* granularity used only for the flattening map, so - // that a flattening map is reusable across the 0.1/0.2 observable binnings. - AxisSpec etaFlatAxis{{-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, "#eta"}; + AxisSpec vzAxis{24, -12., 12., "Vz"}; // flattening maps: 1 cm bins (edges at 0, 8 and 10 match the vertex variations) + AxisSpec chgAxis{2, -2., 2.}; // flattening maps: charge -1, +1 + // One fixed pT axis for every weight map (efficiency, fake, secondary, other-collision, flattening) and the + // pT weight QA. The map passes fill [KPtMin, KPtMax], so the 0-0.15 bin stays empty by construction. + static constexpr std::array KPtMapEdges = {0., 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.5, 3., 3.5, 4., 5., 7., 10.}; + static std::vector effPtEdges() { return {KPtMapEdges.begin(), KPtMapEdges.end()}; } + // flattening maps: a coarser pT axis whose edges are a subset of KPtMapEdges (the phi acceptance varies slowly with pT) + AxisSpec pTAxis{{0., 0.15, 0.2, 0.4, 0.6, 0.8, 1., 3., 5., 7., 10.}, "p_{T} (GeV/c)"}; + AxisSpec phiAxis{KNbinsPhiFlat, KPhiMin, TwoPI, "#phi"}; // flattening maps + + // etaFlatAxis is the *fixed* eta granularity of the flattening maps, the same as the efficiency maps + // (KNbinsEtaFine bins of 0.1 on [-KEtaFineMax, KEtaFineMax]); its edges include those of the 0.1 and + // 0.2 observable binnings, so a map is reusable for both. |eta| > cfgCutEta bins stay empty. + AxisSpec etaFlatAxis{KNbinsEtaFine, -KEtaFineMax, KEtaFineMax, "#eta"}; // etaAxis (physical) and etaBinAxis (integer index) follow the observable // binning and are rebuilt in init(). Placeholders here. @@ -534,6 +571,7 @@ struct RadialFlowDecorr { // 4 per registry (a registry holds at most 512 objects). HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry histosBoot{"HistosBoot", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry histosBoot2{"HistosBoot2", {}, OutputObjHandlingPolicy::AnalysisObject}; // replicas of pT ranges >= 1 HistogramRegistry histosSystA{"HistosSystA", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry histosSystB{"HistosSystB", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry histosSystC{"HistosSystC", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -572,9 +610,9 @@ struct RadialFlowDecorr { TH3* hFake = nullptr; TH3* hEffPtEta = nullptr; TH3* hFakePtEta = nullptr; - THnSparseF* hFlat = nullptr; - THnSparseF* hFlatEarly = nullptr; // Flat variation (data): weights from the early runs, applied to the late runs - THnSparseF* hFlatLate = nullptr; // and from the late runs, applied to the early runs + THnF* hFlat = nullptr; // flattening weights, dense (the builder sets every cell) + THnF* hFlatEarly = nullptr; // Flat variation (data): weights from the early runs, applied to the late runs + THnF* hFlatLate = nullptr; // and from the late runs, applied to the early runs // mean pass, per level: <[pT]> and vs (N_ch, eta bin), vs (N_ch, forward bin) std::array pmean{}; std::array pmeanFwd{}; @@ -637,9 +675,10 @@ struct RadialFlowDecorr { std::array, kNSyst> nchCent{}; // (centrality, N_ch) per variation std::array eff{}; std::array flat{}; - std::array, kNSyst> mean{}; - std::array, kNSyst> fluc{}; - std::array boot{}; + std::array, kNSyst> mean{}; // all pT ranges (eta-bin index axis) + // per pT range: the index-axis profiles are shared with range 0, the maps are separate + std::array, kNLevels>, kNSyst> fluc{}; + std::array, KMaxBoot> boot{}; std::shared_ptr fwdChannel; // (channel: FT0 0-207, FV0 208-255; EChanQA) std::shared_ptr fwdEta; // amplitude-weighted eta per forward bin std::shared_ptr mcTrackOrigin; // (N_ch, pT, ETrackOrigin) of the Base tracks @@ -666,19 +705,9 @@ struct RadialFlowDecorr { std::array, KNHarm>, KNEtaMax> harm{}; void reset() { - for (auto& a : harm) { - for (auto& b : a) { - b.fill(0.); - } - } - for (auto& a : pmwk) { - for (auto& b : a) { - b.fill(0.); - } - } - for (auto& a : wk) { - a.fill(0.); - } + harm = {}; + pmwk = {}; + wk = {}; } void add(int ieta, double pt, double w) { @@ -720,8 +749,9 @@ struct RadialFlowDecorr { } return b; }(); - // recentering offsets: per event and eta bin, q_n = sum w e^{in phi} / W and o_n = sum w pT e^{in phi} / W - static void fillQMean(const std::shared_ptr& h, const PowerSums& s, int ieta, float nCh) + // recentering offsets: per event and eta bin, q_n = sum w e^{in phi} / W and o_n = sum w pT e^{in phi} / W, + // filled at the eta-bin index y = etaIdx(ipt, ieta) + static void fillQMean(const std::shared_ptr& h, const PowerSums& s, int ieta, int y, float nCh) { const double wSum = s.wk[ieta][1]; if (!h || wSum <= KFloatEpsilon) { @@ -729,10 +759,10 @@ struct RadialFlowDecorr { } for (int ih = 0; ih < KNHarm; ++ih) { const auto& q = s.harm[ieta][ih]; - h->Fill(nCh, ieta, KNHarmComp * ih, q[2] / wSum); - h->Fill(nCh, ieta, KNHarmComp * ih + 1, q[3] / wSum); - h->Fill(nCh, ieta, KNHarmComp * ih + 2, q[0] / wSum); - h->Fill(nCh, ieta, KNHarmComp * ih + 3, q[1] / wSum); + h->Fill(nCh, y, KNHarmComp * ih, q[2] / wSum); + h->Fill(nCh, y, KNHarmComp * ih + 1, q[3] / wSum); + h->Fill(nCh, y, KNHarmComp * ih + 2, q[0] / wSum); + h->Fill(nCh, y, KNHarmComp * ih + 3, q[1] / wSum); } } static std::array harmTerms(float phi) @@ -744,11 +774,33 @@ struct RadialFlowDecorr { } return cs; } - std::array, kNSyst> sums{}; + // power sums per (variation, level, pT range) + using PtSums = std::array; + std::array, kNSyst> sums{}; + void resetSums(int v) + { + for (int l = 0; l < kNLevels; ++l) { + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + sums[v][l][ipt].reset(); + } + } + } + // one track (or particle) into the sums of every pT range in ptMask + void addTrack(PtSums& s, uint32_t ptMask, int ieta, double pt, double w, const std::array& cs) + { + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + if (TESTBIT(ptMask, ipt)) { + s[ipt].add(ieta, pt, w); + if (cfgAziTA) { + s[ipt].addHarm(ieta, pt, w, cs); + } + } + } + } // track cut flow of the current event, added to hTrkCount once per event std::array, kNSyst> trkCounts{}; // Flat variation (data): weight map of the other run half for the current event (nullptr: no weights) - const THnSparseF* flatCross = nullptr; + const THnF* flatCross = nullptr; struct FwdSignal { std::array ampl{}; // amplitude (data, MC reco) or number of particles (MC truth) @@ -879,6 +931,7 @@ struct RadialFlowDecorr { countAll(kEvVertexITSTPC); const float centBase = getCentrality(col, cfgEvSel.cfgCentralityChoice); + ev.cent[kSystBase] = centBase; // the track weights use the Base centrality, also when Base itself is off const float multPV = col.multNTracksPV(); ev.occupancy = col.trackOccupancyInTimeRange(); const bool passAddPileup = isPassAddPileup(multPV, nCh, centBase, pupNSigma); @@ -970,7 +1023,7 @@ struct RadialFlowDecorr { return 0; } countAll(kTrkPVContrib); - if (trk.pt() < cfgCutPtLower || trk.pt() > cfgCutPtUpper || std::abs(trk.eta()) > cfgCutEta) { + if (trk.pt() < KPtMin || trk.pt() > KPtMax || std::abs(trk.eta()) > cfgCutEta) { return 0; } countAll(kTrkPtEta); @@ -1038,7 +1091,9 @@ struct RadialFlowDecorr { if (std::abs(pd->Charge()) == 0) { return false; } - if (particle.pt() < cfgCutPtLower || particle.pt() > cfgCutPtUpper || std::abs(particle.eta()) > cfgCutEta) { + // truth particles out to the edge of the eta axes (|eta| < KEtaFineMax = 1), so that tracks migrating + // across |eta| = cfgCutEta between truth and reco stay inside the efficiency maps and observable bins + if (particle.pt() < KPtMin || particle.pt() > KPtMax || std::abs(particle.eta()) >= KEtaFineMax) { return false; } return true; @@ -1079,7 +1134,7 @@ struct RadialFlowDecorr { int etaBinsOf(float eta, std::array& bins) const { int n = 0; - for (int ieta = 0; ieta < nEta && n < 2; ++ieta) { + for (int ieta = 0; ieta < nEta && n < static_cast(bins.size()); ++ieta) { if (eta > etaLw[ieta] && eta <= etaUp[ieta]) { bins[n++] = ieta; } @@ -1151,19 +1206,19 @@ struct RadialFlowDecorr { return 1.f + 0.01f * KEffVarPercent; } if (v == kSystEffPtDep) { - // Linear in pT: full variation at cfgPtMin, none at KEffVarPtZero. - const float shape = std::clamp((KEffVarPtZero - pt) / (KEffVarPtZero - cfgPtMin), 0.f, 1.f); + // Linear in pT: full variation at KEffVarPtAnchor, none at KEffVarPtZero. + const float shape = std::clamp((KEffVarPtZero - pt) / (KEffVarPtZero - KEffVarPtAnchor), 0.f, 1.f); return 1.f + 0.01f * KEffVarPercent * shape; } return 1.f; } - static float getFlatteningWeight(const THnSparseF* h, float vz, float chg, float pt, float eta, float phi) + static float getFlatteningWeight(const THnF* h, float vz, float chg, float pt, float eta, float phi) { if (!h) { return 0.0; } - std::array bins{}; + std::array bins{}; bins[0] = h->GetAxis(0)->FindBin(vz); bins[1] = h->GetAxis(1)->FindBin(chg); bins[2] = h->GetAxis(2)->FindBin(pt); @@ -1279,7 +1334,7 @@ struct RadialFlowDecorr { int fwdSubBin(int det, float eta) const { const auto& e = fwdEdges[det]; - if (e.size() < 2 || eta < e.front() || eta >= e.back()) { + if (e.size() < KMinEdges || eta < e.front() || eta >= e.back()) { return -1; } const int i = static_cast(std::upper_bound(e.begin(), e.end(), eta) - e.begin()) - 1; @@ -1338,7 +1393,7 @@ struct RadialFlowDecorr { auto* alFT0 = ccdb->getForTimeStamp>("FT0/Calib/Align", ts); auto* alFV0 = ccdb->getForTimeStamp>("FV0/Calib/Align", ts); - if (!alFT0 || alFT0->size() < 2 || !alFV0 || alFV0->size() < 2) { + if (!alFT0 || alFT0->size() < KNAlignHalves || !alFV0 || alFV0->size() < KNAlignHalves) { LOGF(fatal, "Could not load FT0/FV0 alignment (FT0/Calib/Align, FV0/Calib/Align) for timestamp %llu", ts); return; } @@ -1448,7 +1503,7 @@ struct RadialFlowDecorr { continue; } const float eta = particle.eta(); - if (std::abs(eta) < 1.5f) { + if (std::abs(eta) < KFwdTruthEtaMin) { continue; } auto* pd = pdg->GetParticle(particle.pdgCode()); @@ -1492,7 +1547,8 @@ struct RadialFlowDecorr { return p->GetBinContent(gb); } - void computeEventObs(const PowerSums& s, int v, int level, const FwdSignal* fwd, float nCh, EventObs& evObs) + // Observables of pT range ipt; the mean-pass references are read at the eta-bin index etaIdx(ipt, ieta). + void computeEventObs(const PowerSums& s, int v, int level, const FwdSignal* fwd, float nCh, EventObs& evObs, int ipt) { const auto& input = inputs[v]; const double nan = std::numeric_limits::quiet_NaN(); @@ -1506,7 +1562,8 @@ struct RadialFlowDecorr { evObs.p1kBar[ieta] = nan; evObs.qTAValid[ieta] = false; // c2 (>=2 tracks), c3 (>=3 tracks) and the pT deviation need a usable reference - const double mmpt = profileValue(input.pmean[level], nCh, ieta + 1, kMeanPt + 1); + const int ybin = etaIdx(ipt, ieta) + 1; + const double mmpt = profileValue(input.pmean[level], nCh, ybin, kMeanPt + 1); if (std::isfinite(mmpt) && mmpt != 0) { const auto corr = calculateC2C3FromSums(s.pmwk[ieta], s.wk[ieta], mmpt); evObs.c2[ieta] = corr.c2; @@ -1523,7 +1580,7 @@ struct RadialFlowDecorr { std::array off{}; if (pq) { for (int c = 0; c < KNHarmComp; ++c) { - const double o = profileValue(pq, nCh, ieta + 1, KNHarmComp * ih + c + 1); + const double o = profileValue(pq, nCh, ybin, KNHarmComp * ih + c + 1); off[c] = std::isfinite(o) ? o : 0.; } } @@ -1536,7 +1593,7 @@ struct RadialFlowDecorr { evObs.qTAValid[ieta] = true; } } - evObs.dN[ieta] = evObs.nW[ieta] - profileValue(input.pmean[level], nCh, ieta + 1, kMeanN + 1); + evObs.dN[ieta] = evObs.nW[ieta] - profileValue(input.pmean[level], nCh, ybin, kMeanN + 1); } evObs.fwdValid.fill(false); if (!fwd || !input.pmeanFwd[level]) { @@ -1556,9 +1613,9 @@ struct RadialFlowDecorr { } } - // Fill the fluctuation profiles of one (variation, level) at x = N_ch, or of one bootstrap - // replica with weight w at x = centrality. - void fillFluc(const FlucHists& h, const EventObs& evObs, double w, float x) + // Fill the fluctuation profiles of one (variation, level, pT range ipt) at x = N_ch, or of one + // bootstrap replica with weight w at x = centrality. The eta-bin index axes are filled at etaIdx(ipt, ieta). + void fillFluc(const FlucHists& h, const EventObs& evObs, double w, float x, int ipt) { auto fill = [w](const std::shared_ptr& p, double a, double b, double c, double val) { if (p) { @@ -1568,15 +1625,16 @@ struct RadialFlowDecorr { // multiplicity, meanpT, C2, C3 vs eta bin for (int ieta = 0; ieta < nEta; ++ieta) { - fill(h.obs, x, ieta, kObsN, evObs.nW[ieta]); + const int y = etaIdx(ipt, ieta); + fill(h.obs, x, y, kObsN, evObs.nW[ieta]); if (std::isfinite(evObs.mean[ieta])) { - fill(h.obs, x, ieta, kObsMeanPt, evObs.mean[ieta]); + fill(h.obs, x, y, kObsMeanPt, evObs.mean[ieta]); } if (std::isfinite(evObs.c2[ieta])) { - fill(h.obs, x, ieta, kObsC2, evObs.c2[ieta]); + fill(h.obs, x, y, kObsC2, evObs.c2[ieta]); } if (std::isfinite(evObs.c3[ieta])) { - fill(h.obs, x, ieta, kObsC3, evObs.c3[ieta]); + fill(h.obs, x, y, kObsC3, evObs.c3[ieta]); } } @@ -1593,13 +1651,13 @@ struct RadialFlowDecorr { float c3SubC = evObs.c2[ietaC] * evObs.p1kBar[ietaA]; // ABB: 2 particles from C, 1 from A if (std::isfinite(c2Sub)) { - fill(h.obs, x, ietaA, kObsC2Sub, c2Sub); + fill(h.obs, x, etaIdx(ipt, ietaA), kObsC2Sub, c2Sub); } if (std::isfinite(c3SubA)) { - fill(h.obs, x, ietaA, kObsC3Sub, c3SubA); + fill(h.obs, x, etaIdx(ipt, ietaA), kObsC3Sub, c3SubA); } if (std::isfinite(c3SubC)) { - fill(h.obs, x, ietaC, kObsC3Sub, c3SubC); + fill(h.obs, x, etaIdx(ipt, ietaC), kObsC3Sub, c3SubC); } } @@ -1679,10 +1737,10 @@ struct RadialFlowDecorr { } for (int ieta = 0; ieta < nEta; ++ieta) { if (std::isfinite(evObs.p1kBar[ieta])) { - fill(h.covFwdPt, x, ieta, k, evObs.p1kBar[ieta] * da); + fill(h.covFwdPt, x, etaIdx(ipt, ieta), k, evObs.p1kBar[ieta] * da); } if (std::isfinite(evObs.dN[ieta])) { - fill(h.covFwdN, x, ieta, k, evObs.dN[ieta] * da); + fill(h.covFwdN, x, etaIdx(ipt, ieta), k, evObs.dN[ieta] * da); } } // tag-tag covariances (all pairs; the diagonal repeats Prof_FwdAmpVar): forward subevents, e.g. A +- C @@ -1876,7 +1934,7 @@ struct RadialFlowDecorr { } for (int c = 0; c < kNEffCat; ++c) { hists.eff[v].maps[c] = reg(v).add((systPrefix(v) + effMapName[c]).c_str(), ";centrality (%);p_{T} (GeV/c);#eta", kTH3F, - {{centAxisEff}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); + {{centAxisEff}, {effPtEdges(), "p_{T} (GeV/c)"}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); } } } @@ -1888,7 +1946,7 @@ struct RadialFlowDecorr { histos.add("etaTruthReco", ";#eta^{MC};#eta^{reco}", kTH2F, {{KNbinsEtaRes, -KEtaFineMax, KEtaFineMax}, {KNbinsEtaRes, -KEtaFineMax, KEtaFineMax}}); histos.add("TruthTracKVz", ";Vz^{MC};Vz^{Reco}", kTH2F, {{KNbinsVz, KVzMin, KVzMax}, {KNbinsVz, KVzMin, KVzMax}}); histos.add("vzResolution", ";Vz^{MC};(Vz^{reco}-Vz^{MC})/Vz^{MC}", kTH2F, {{KNbinsVz, KVzMin, KVzMax}, {100, -0.1, 0.1}}); - hists.mcTrackOrigin = histos.add("MCQA/Hist3D_Nch_pT_origin", ";N_{ch};p_{T} (GeV/c);origin (0 primary, 1 secondary, 2 other MC collision, 3 no label)", kTH3F, {{nchAxisObs}, {KNbinsPtRes, KPtMin, KPtMax}, {kNOrigin, -0.5, static_cast(kNOrigin) - 0.5}}); + hists.mcTrackOrigin = histos.add("MCQA/Hist3D_Nch_pT_origin", ";N_{ch};p_{T} (GeV/c);origin (0 primary, 1 secondary, 2 other MC collision, 3 no label)", kTH3F, {{nchAxisObs}, {effPtEdges(), "p_{T} (GeV/c)"}, {kNOrigin, -0.5, static_cast(kNOrigin) - 0.5}}); } // vertex QA of all MC passes @@ -1935,6 +1993,9 @@ struct RadialFlowDecorr { if (cfgFwdOpt.cfgFwd) { h.fwdMean = r.add((p + "Prof_FwdAmpMean_Nch_fwdbin").c_str(), ";N_{ch};forward bin;#LT a#GT", kTProfile2D, {{nchAxisOf(v)}, {fwdAxis}}); } + dropBinSumw2(h.mean); + dropBinSumw2(h.qMean); + dropBinSumw2(h.fwdMean); } } } @@ -1942,11 +2003,19 @@ struct RadialFlowDecorr { // Fluctuation-pass profiles named prefix + name + suffix, with x axis ax. // withGapSum: also the maps in (Delta eta, Sigma eta). // withHarm: also the c_n (n >= 2) and d_n maps (Base only, not in the bootstrap replicas). - FlucHists declareFluc(HistogramRegistry& r, const std::string& p, const std::string& s, const AxisSpec& ax, bool withGapSum, bool withHarm = false) + // pT range ipt > 0: the eta-bin index profiles are those of range 0 (h0), the maps get the suffix _pt, + // and the forward-only profiles (no pT dependence) are not booked again. + FlucHists declareFluc(HistogramRegistry& r, const std::string& p, const std::string& s, const AxisSpec& ax, bool withGapSum, bool withHarm, int ipt, const FlucHists* h0) { FlucHists h; - auto n = [&](const char* name) { return p + name + s; }; - h.obs = r.add(n("Prof_Obs_Nch_etabin_obs").c_str(), ";N_{ch};#eta-bin;observable (0 N, 1 [p_{T}], 2 c_{2}, 3 c_{3}, 4 C_{2}^{sub} mirror, 5 C_{3}^{sub} mirror)", kTProfile3D, {{ax}, {etaBinAxis}, {kNFlucObs, -0.5, static_cast(kNFlucObs) - 0.5}}); + auto n = [&](const char* name) { return p + name + s + ptSuffix(ipt); }; + if (ipt > 0) { + h.obs = h0->obs; + h.covFwdPt = h0->covFwdPt; + h.covFwdN = h0->covFwdN; + } else { + h.obs = r.add(n("Prof_Obs_Nch_etabin_obs").c_str(), ";N_{ch};#eta-bin;observable (0 N, 1 [p_{T}], 2 c_{2}, 3 c_{3}, 4 C_{2}^{sub} mirror, 5 C_{3}^{sub} mirror)", kTProfile3D, {{ax}, {etaBinAxis}, {kNFlucObs, -0.5, static_cast(kNFlucObs) - 0.5}}); + } h.c2Sub2D = r.add(n("Prof_C2Sub2D_Nch_etaA_etaC").c_str(), ";N_{ch};#eta_{A};#eta_{C}", kTProfile3D, {{ax}, {etaAxis}, {etaAxis}}); h.cNN2D = r.add(n("Prof_CNN2D_Nch_etaA_etaC").c_str(), ";N_{ch};#eta_{A};#eta_{C};#LT #delta N_{A} #delta N_{C}#GT", kTProfile3D, {{ax}, {etaAxis}, {etaAxis}}); h.c3Sub2D = r.add(n("Prof_C3Sub2D_Nch_etaA_etaC").c_str(), ";N_{ch};#eta_{A} (pair);#eta_{C} (single)", kTProfile3D, {{ax}, {etaAxis}, {etaAxis}}); @@ -1962,11 +2031,11 @@ struct RadialFlowDecorr { } } } - if (withGapSum) { + if (withGapSum && ipt == 0) { // nominal range only: an exact re-indexing of the (eta_A, eta_C) maps h.gapSum2D = r.add(n("Prof_GapSum2D_Nch_gap_sum").c_str(), ";N_{ch};#eta_{A} #minus #eta_{C} (Gap);#eta_{A} + #eta_{C} (Sum)", kTProfile3D, {{ax}, {gapAxis}, {sumAxis}}); h.c3GapSum2D = r.add(n("Prof_C3GapSum2D_Nch_gap_sum").c_str(), ";N_{ch};#eta_{pair} #minus #eta_{single};#eta_{pair} + #eta_{single}", kTProfile3D, {{ax}, {gapAxis}, {sumAxis}}); } - if (cfgFwdOpt.cfgFwd) { + if (cfgFwdOpt.cfgFwd && ipt == 0) { h.fwdVar = r.add(n("Prof_FwdAmpVar_Nch_fwdbin").c_str(), ";N_{ch};forward bin;#LT #delta a^{2}#GT", kTProfile2D, {{ax}, {fwdAxis}}); h.covFwdPt = r.add(n("Prof_CovFwdAmpMeanpT_Nch_etabin_fwdbin").c_str(), ";N_{ch};#eta-bin;forward bin;#LT #delta[p_{T}] #delta a#GT", kTProfile3D, {{ax}, {etaBinAxis}, {fwdAxis}}); h.covFwdN = r.add(n("Prof_CovFwdAmpNch_Nch_etabin_fwdbin").c_str(), ";N_{ch};#eta-bin;forward bin;#LT #delta N #delta a#GT", kTProfile3D, {{ax}, {etaBinAxis}, {fwdAxis}}); @@ -1975,6 +2044,29 @@ struct RadialFlowDecorr { return h; } + // Profiles filled with unit weights only (everything but the bootstrap replicas) do not need the per-bin + // sum of squared weights: it equals the bin entries. Dropping it saves a quarter of their memory. + template + static void dropBinSumw2(const std::shared_ptr

& p) + { + if (p) { + p->Sumw2(false); + } + } + static void dropBinSumw2(const FlucHists& h) + { + for (const auto& p : {h.obs, h.c2Sub2D, h.cNN2D, h.c2TA2D, h.c1Sub2D, h.c3Sub2D, h.gapSum2D, h.c3GapSum2D, h.covFwdPt, h.covFwdN, h.covFwdFwd}) { + dropBinSumw2(p); + } + for (const auto& p : h.cn2D) { + dropBinSumw2(p); + } + for (const auto& p : h.dn2D) { + dropBinSumw2(p); + } + dropBinSumw2(h.fwdVar); + } + // Gap/sum maps for every DATA variation; in MC only for Base, to limit the memory. void declareFlucAll(uint32_t mask, bool isMC) { @@ -1985,7 +2077,10 @@ struct RadialFlowDecorr { continue; } for (int l = 0; l < nLev; ++l) { - hists.fluc[v][l] = declareFluc(reg(v), systPrefix(v) + (isMC ? mcDir[l] : std::string{}), "", nchAxisOf(v), !isMC || v == kSystBase, v == kSystBase); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + hists.fluc[v][l][ipt] = declareFluc(reg(v), systPrefix(v) + (isMC ? mcDir[l] : std::string{}), "", nchAxisOf(v), !isMC || v == kSystBase, v == kSystBase, ipt, &hists.fluc[v][l][0]); + dropBinSumw2(hists.fluc[v][l][ipt]); + } } } } @@ -1994,8 +2089,10 @@ struct RadialFlowDecorr { // is the Base centrality in 5% classes, the binning in which the errors are evaluated. void declareBootstrapHists() { - for (int s = 0; s < nBoot; ++s) { - hists.boot[s] = declareFluc(histosBoot, "Bootstrap/", Form("_sample%d", s), centAxisBoot, true); + for (int s = 0; s < bootstrap.n; ++s) { + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + hists.boot[s][ipt] = declareFluc(ipt == 0 ? histosBoot : histosBoot2, "Bootstrap/", Form("_sample%d", s), centAxisBoot, true, false, ipt, &hists.boot[s][0]); + } } } @@ -2021,7 +2118,7 @@ struct RadialFlowDecorr { declareFwdQA(); } declareMean(maskAll, true); - histos.add("Prof_Weights_pT_type", ";p_{T};weight (0 eff, 1 fake, 2 w)", kTProfile2D, {{KNbinsPtRes, KPtMin, KPtMax}, {3, -0.5, 2.5}}); + histos.add("Prof_Weights_pT_type", ";p_{T};weight (0 eff, 1 fake, 2 w)", kTProfile2D, {{effPtEdges(), "p_{T} (GeV/c)"}, {3, -0.5, 2.5}}); } void declareMCFlucHists() @@ -2040,11 +2137,12 @@ struct RadialFlowDecorr { declareFlat(maskFlat, ""); if (TESTBIT(maskAll, kSystFlat) && cfgFlatRunSplit > 0) { auto& h = hists.flat[kSystFlat]; - for (auto* m : {&h.early, &h.late}) { - const char* nm = (m == &h.early) ? "hFlatMapsEarly" : "hFlatMapsLate"; - *m = reg(kSystFlat).add((systPrefix(kSystFlat) + nm).c_str(), ";vz;sign;pt;eta;phi;weight", kTHnSparseF, - {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}, {kNFlatType, -0.5, static_cast(kNFlatType) - 0.5}}); - } + auto book = [&](const char* nm) { + return reg(kSystFlat).add((systPrefix(kSystFlat) + nm).c_str(), ";vz;sign;pt;eta;phi;weight", kTHnSparseF, + {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}, {kNFlatType, -0.5, static_cast(kNFlatType) - 0.5}}); + }; + h.early = book("hFlatMapsEarly"); + h.late = book("hFlatMapsLate"); } } @@ -2056,7 +2154,7 @@ struct RadialFlowDecorr { declareFwdQA(); } declareMean(maskAll, false); - histos.add("Prof_Weights_pt_eta_type", ";p_{T} (GeV/c);#eta;weight (0 eff, 1 fake, 2 w_{#phi})", kTProfile3D, {{KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}, {3, -0.5, 2.5}}); + histos.add("Prof_Weights_pt_eta_type", ";p_{T} (GeV/c);#eta;weight (0 eff, 1 fake, 2 w_{#phi})", kTProfile3D, {{effPtEdges(), "p_{T} (GeV/c)"}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}, {3, -0.5, 2.5}}); } void declareDataFlucHists() @@ -2067,7 +2165,7 @@ struct RadialFlowDecorr { declareFwdQA(); } declareFlucAll(maskAll, false); - if (doBoot) { + if (bootstrap.on) { declareBootstrapHists(); } } @@ -2086,23 +2184,29 @@ struct RadialFlowDecorr { } // Flattening weights w_phi = / N_phi from bin `type` of the last axis of hFlatMaps. - THnSparseF* buildWeightMapFromRaw(const THnSparse* hRaw, int type, const char* mapName) + // Dense weight map (THnF): the builder sets every cell, so a sparse map would hold all of them anyway. + THnF* buildWeightMapFromRaw(const THnSparse* hRaw, int type, const char* mapName) { - if (!hRaw || hRaw->GetNdimensions() != 6) { + if (!hRaw || hRaw->GetNdimensions() != KNFlatAxes + 1) { LOGF(error, "Raw eta-phi map for '%s' is null or not 6D; no flattening will be applied.", mapName); return nullptr; } - std::array nb{}; - for (int i = 0; i < 5; ++i) { + std::array nb{}; + for (int i = 0; i < KNFlatAxes; ++i) { nb[i] = hRaw->GetAxis(i)->GetNbins(); } - auto* hWMap = new THnSparseF(mapName, Form("Flattening Weight Map %s (w_{#phi} = / N_{#phi})", mapName), 5, nb.data()); - for (int i = 0; i < 5; ++i) { + std::array xMin{}, xMax{}; + for (int i = 0; i < KNFlatAxes; ++i) { + xMin[i] = hRaw->GetAxis(i)->GetXmin(); + xMax[i] = hRaw->GetAxis(i)->GetXmax(); + } + auto* hWMap = new THnF(mapName, Form("Flattening Weight Map %s (w_{#phi} = / N_{#phi})", mapName), KNFlatAxes, nb.data(), xMin.data(), xMax.data()); + for (int i = 0; i < KNFlatAxes; ++i) { copyAxis(hWMap->GetAxis(i), hRaw->GetAxis(i)); } - std::array bins{}; - std::array wbins{}; + std::array bins{}; + std::array wbins{}; bins[5] = type + 1; for (int iv = 1; iv <= nb[0]; ++iv) { bins[0] = wbins[0] = iv; @@ -2207,7 +2311,7 @@ struct RadialFlowDecorr { return true; } - THnSparseF* loadFlatMap(int s, const std::string& path, const std::string& dir, int64_t now, const char* name = "hFlatMaps") + THnF* loadFlatMap(int s, const std::string& path, const std::string& dir, int64_t now, const char* name = "hFlatMaps") { auto* lst = ccdb->getForTimeStamp(path, now); if (!lst) { @@ -2242,19 +2346,20 @@ struct RadialFlowDecorr { // The fluctuation pass reads its reference at (N_ch bin, eta bin), so the binning of the // mean pass has to be the same. const auto* p = inputs[v].pmean[l]; - if (p->GetYaxis()->GetNbins() != nEta + 1 || + if (p->GetYaxis()->GetNbins() != ptRange.n * (nEta + 1) || p->GetXaxis()->GetNbins() != nchAxisOf(v).getNbins() || std::abs(p->GetXaxis()->GetBinWidth(1) - nchBinWidthOf(v)) > KFloatEpsilon) { LOGF(fatal, - "Mean pass %s has %d N_ch bins of width %.1f and %d eta bins; this job uses %d bins of " - "width %d and %d eta bins. Rerun the mean pass with the same binning.", + "Mean pass %s has %d N_ch bins of width %.1f and %d (pT range, eta) bins; this job uses %d bins of " + "width %d and %d = %d pT range(s) x %d eta bins. Rerun the mean pass with the same binning " + "(cfgEtaBinWidth).", path.c_str(), p->GetXaxis()->GetNbins(), p->GetXaxis()->GetBinWidth(1), p->GetYaxis()->GetNbins(), - nchAxisOf(v).getNbins(), nchBinWidthOf(v), nEta + 1); + nchAxisOf(v).getNbins(), nchBinWidthOf(v), ptRange.n * (nEta + 1), ptRange.n, nEta + 1); } if (cfgAziTA) { loadProfileFromList(lst, (d + "Prof_QMean_Nch_etabin_comp").c_str(), inputs[v].pqmean[l]); if (!inputs[v].pqmean[l] || inputs[v].pqmean[l]->GetZaxis()->GetNbins() != KNHarmComp * KNHarm || - inputs[v].pqmean[l]->GetYaxis()->GetNbins() != nEta + 1) { + inputs[v].pqmean[l]->GetYaxis()->GetNbins() != ptRange.n * (nEta + 1)) { LOGF(warning, "%s: no recentering offsets Prof_QMean_Nch_etabin_comp (rerun the mean pass with cfgAziTA); flow vectors not recentred", path.c_str()); inputs[v].pqmean[l] = nullptr; } @@ -2300,42 +2405,59 @@ struct RadialFlowDecorr { nchAxisQMean = AxisSpec{KNQMeanNchBins, -KBinOffset, nMax - KBinOffset, "N_{ch}"}; } - // ---- run type: MC is pinned to the reference binning, DATA is configurable ---- - // cfgEtaBinWidth exists only to vary the DATA measurement. MC (efficiency, - // closure) always uses 0.2-wide bins so its mean and fluctuation passes stay - // mutually consistent and reproducible. + // ---- run type ---- + // cfgEtaBinWidth sets the observable eta bins for MC and DATA alike. The efficiency and + // flattening maps have their own fixed eta axes and do not depend on it. The MC mean and + // fluctuation passes must use the same width, since the fluctuation pass reads the mean + // pass output per eta bin. const bool isMcRun = (cfgRunGetEff || cfgRunGetMCFlat || cfgRunMCMean || cfgRunMCFluc); const bool isDataRun = (cfgRunGetDataFlat || cfgRunDataMean || cfgRunDataFluc); if (isMcRun && isDataRun) { LOGF(fatal, - "MC and DATA process switches are both enabled in one job: MC binning is " - "pinned to 0.2 while DATA uses cfgEtaBinWidth. Run them as separate jobs."); + "MC and DATA process switches are both enabled in one job: they read and write " + "different CCDB objects. Run them as separate jobs."); + } + const float effEtaBinWidth = static_cast(cfgEtaBinWidth); + + // ---- observable pT ranges: fixed, all filled in the mean and fluctuation passes ---- + ptRange.n = KNPtMax; + for (int k = 0; k < KNPtMax; ++k) { + ptRange.lo[k] = KPtRanges[k][0]; + ptRange.hi[k] = KPtRanges[k][1]; + } + if (cfgRunGetEff || cfgRunGetMCFlat || cfgRunGetDataFlat) { + LOGF(info, "Efficiency/flattening pass: maps filled over the pT range [%.2f, %.2f] GeV/c", KPtMin, KPtMax); + } + if (cfgRunMCMean || cfgRunMCFluc || cfgRunDataMean || cfgRunDataFluc) { + LOGF(info, "Mean/fluctuation pass: tracks in [%.2f, %.2f] GeV/c, observable pT ranges (%.2f, %.2f], (%.2f, %.2f] (_pt1), (%.2f, %.2f] (_pt2) GeV/c", + KPtMin, KPtMax, ptRange.lo[0], ptRange.hi[0], ptRange.lo[1], ptRange.hi[1], ptRange.lo[2], ptRange.hi[2]); } - const float mcEtaBinWidth = 0.2f; - const float effEtaBinWidth = isMcRun ? mcEtaBinWidth : static_cast(cfgEtaBinWidth); // ---- observable eta binning (width = effEtaBinWidth) -------------------------- { - const float halfEta = cfgCutEta; + // narrow bins: MC over the full axis |eta| < KEtaFineMax (truth particles fill it, reco |eta| < cfgCutEta), + // DATA over |eta| < cfgCutEta (no track beyond it); the reference bin (index 0) is |eta| < cfgCutEta + const float halfEta = isMcRun ? KEtaFineMax : static_cast(cfgCutEta); const float width = effEtaBinWidth; if (width <= 0.f) { LOGF(fatal, "eta bin width must be > 0 (got %.3f)", width); } const int nHalf = static_cast(std::lround(halfEta / width)); - if (nHalf < 1 || std::abs(nHalf * width - halfEta) > KEtaEdgeTolerance) { + const int nHalfCut = static_cast(std::lround(cfgCutEta / width)); + if (nHalf < 1 || std::abs(nHalf * width - halfEta) > KEtaEdgeTolerance || nHalfCut < 1 || + std::abs(nHalfCut * width - cfgCutEta) > KEtaEdgeTolerance || cfgCutEta > halfEta + KEtaEdgeTolerance) { LOGF(fatal, - "cfgCutEta=%.3f is not an integer multiple of the eta bin width=%.3f; " - "eta=0 must be a bin edge for mirror pairing.", - halfEta, width); + "the eta bin width %.3f must divide both cfgCutEta = %.3f and the axis edge %.1f (cfgCutEta <= %.1f); " + "eta = 0 and |eta| = cfgCutEta must be bin edges.", + width, cfgCutEta.value, halfEta, halfEta); } const int nbins = 2 * nHalf; // even -> symmetric about eta = 0 const float lo = -halfEta; - const float hi = halfEta; etaLw.clear(); etaUp.clear(); - etaLw.push_back(lo); // index 0: full-range reference bin - etaUp.push_back(hi); + etaLw.push_back(-cfgCutEta); // index 0: reference bin |eta| < cfgCutEta + etaUp.push_back(cfgCutEta); for (int i = 0; i < nbins; ++i) { etaLw.push_back(lo + i * width); etaUp.push_back(lo + (i + 1) * width); @@ -2352,7 +2474,9 @@ struct RadialFlowDecorr { obsEdges.push_back(etaUp[i]); } etaAxis = AxisSpec{obsEdges, "#eta"}; - etaBinAxis = AxisSpec{nEta + 1, -0.5, static_cast(nEta) + 0.5, "#eta bin Number"}; + // index axis: pT range ipt, eta bin ieta at etaIdx(ipt, ieta) = ipt * (nEta + 1) + ieta + const int nIdx = ptRange.n * (nEta + 1); + etaBinAxis = AxisSpec{nIdx, -0.5, static_cast(nIdx) - 0.5, ptRange.n > 1 ? "p_{T} range #times (n_{#eta} + 1) + #eta bin" : "#eta bin Number"}; // (Delta eta, Sigma eta): fills are integer multiples of w, so 2N-1 bins // of width w centred on them, out to (N-1/2) w. const int nGS = 2 * nbins - 1; @@ -2369,7 +2493,7 @@ struct RadialFlowDecorr { nFwd = kNFwdDet; for (int det = 0; det < kNFwdDet; ++det) { auto& e = fwdEdges[det]; - if (e.size() < 2 || !std::is_sorted(e.begin(), e.end())) { + if (e.size() < KMinEdges || !std::is_sorted(e.begin(), e.end())) { LOGF(fatal, "%s eta edges must be >= 2 increasing values", detName[det]); } fwdOffset[det] = nFwd; @@ -2390,7 +2514,16 @@ struct RadialFlowDecorr { } // ---- variations ---- - const uint32_t active = (static_cast(cfgSyst.cfgSystMask.value) | 1) & KAllSystMask; + // Base is always on in the efficiency and flattening passes (other variations borrow its maps); a + // mean/fluctuation job may leave it out, so that split groups of variations do not all repeat Base. + const bool mapPassMask = cfgRunGetEff || cfgRunGetMCFlat || cfgRunGetDataFlat; + const uint32_t active = (static_cast(cfgSyst.cfgSystMask.value) | (mapPassMask ? 1u : 0u)) & KAllSystMask; + if (active == 0) { + LOGF(fatal, "cfgSystMask = %d selects no variation", cfgSyst.cfgSystMask.value); + } + if (!TESTBIT(active, kSystBase)) { + LOGF(info, "Base is off in this job (cfgSystMask = %d): split group of variations only", cfgSyst.cfgSystMask.value); + } maskAll = maskEff = maskFlat = 0; for (int v = 0; v < kNSyst; ++v) { if (!TESTBIT(active, v)) { @@ -2430,11 +2563,24 @@ struct RadialFlowDecorr { } // bootstrap: Base data fluctuation pass only - doBoot = cfgRunDataFluc && (cfgNBootstrap > 0); - nBoot = std::min(cfgNBootstrap, KMaxBoot); + bootstrap.on = cfgRunDataFluc && (cfgNBootstrap > 0) && TESTBIT(maskAll, kSystBase); // replicas of Base only + bootstrap.n = std::min(cfgNBootstrap, KMaxBoot); if (cfgNBootstrap > KMaxBoot) { LOGF(warning, "cfgNBootstrap=%d capped at %d", cfgNBootstrap.value, KMaxBoot); } + { + // objects per replica, as booked in declareFluc (with gap/sum maps): range 0 in histosBoot, the + // further pT ranges together in histosBoot2 + const int nTA = cfgAziTA ? 2 : 0; + const int nFwdObj = cfgFwdOpt.cfgFwd ? 4 : 0; + const int perReplica0 = 4 + nTA + 2 + nFwdObj; + const int perReplicaExtra = (ptRange.n - 1) * (3 + nTA); // no gap/sum maps for the further ranges + const int nFit = std::min(KMaxRegistryObjects / perReplica0, perReplicaExtra > 0 ? KMaxRegistryObjects / perReplicaExtra : KMaxBoot); + if (bootstrap.on && bootstrap.n > nFit) { + LOGF(warning, "%d bootstrap replicas (%d + %d profiles each) exceed the %d objects of a registry: capped at %d", bootstrap.n, perReplica0, perReplicaExtra, KMaxRegistryObjects, nFit); + bootstrap.n = nFit; + } + } rng.SetSeed(cfgBootstrapSeed); // ---- CCDB ---- @@ -2813,7 +2959,7 @@ struct RadialFlowDecorr { continue; } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + if (pt <= KPtMin || pt > KPtMax) { // map range continue; } float eta = track.eta(), phi = track.phi(); @@ -2875,23 +3021,24 @@ struct RadialFlowDecorr { for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - for (auto& s : sums[v]) { - s.reset(); - } + resetSums(v); } } // --- truth: independent of the track cuts, copied to every selected variation --- - PowerSums& sumsTru = sums[kSystBase][kGen]; + PtSums& sumsTru = sums[kSystBase][kGen]; if (!TESTBIT(ev.mask, kSystBase)) { - sumsTru.reset(); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + sumsTru[ipt].reset(); + } } for (const auto& particle : mcParticles) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) { continue; } float pt = particle.pt(), eta = particle.eta(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } const auto csTru = harmTerms(particle.phi()); @@ -2899,15 +3046,14 @@ struct RadialFlowDecorr { if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; } - sumsTru.add(ieta, pt, 1.); - if (cfgAziTA) { - sumsTru.addHarm(ieta, pt, 1., csTru); - } + addTrack(sumsTru, ptMask, ieta, pt, 1., csTru); } } for (int v = 1; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - sums[v][kGen] = sumsTru; + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + sums[v][kGen][ipt] = sumsTru[ipt]; + } } } @@ -2919,12 +3065,14 @@ struct RadialFlowDecorr { continue; } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } + const bool inPt0 = TESTBIT(ptMask, 0); // QA of range 0 only float eta = track.eta(), phi = track.phi(); float sign = track.sign(); - if (TESTBIT(trackMask, kSystBase)) { + if (TESTBIT(trackMask, kSystBase) && inPt0) { histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); @@ -2938,10 +3086,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(trackMask, v)) { - sums[v][kReco].add(ieta, pt, 1.); - if (cfgAziTA) { - sums[v][kReco].addHarm(ieta, pt, 1., cs); - } + addTrack(sums[v][kReco], ptMask, ieta, pt, 1., cs); } } } @@ -2951,7 +3096,7 @@ struct RadialFlowDecorr { if (weightMask == 0) { continue; } - if (TESTBIT(weightMask, kSystBase)) { + if (TESTBIT(weightMask, kSystBase) && inPt0) { histos.fill(HIST("Prof_Weights_pT_type"), pt, 0, eff); histos.fill(HIST("Prof_Weights_pT_type"), pt, 1, fake); histos.fill(HIST("Prof_Weights_pT_type"), pt, 2, w[kSystBase]); @@ -2965,10 +3110,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(weightMask, v)) { - sums[v][kRecoEffCorr].add(ieta, pt, w[v]); - if (cfgAziTA) { - sums[v][kRecoEffCorr].addHarm(ieta, pt, w[v], cs); - } + addTrack(sums[v][kRecoEffCorr], ptMask, ieta, pt, w[v], cs); } } } @@ -2989,18 +3131,21 @@ struct RadialFlowDecorr { continue; } for (int l = 0; l < kNLevels; ++l) { - const auto& s = sums[v][l]; const auto& h = hists.mean[v][l]; - for (int ieta = 0; ieta < nEta; ++ieta) { - const double n = s.wk[ieta][1]; - h.mean->Fill(nCh, ieta, kMeanN, n); - fillQMean(h.qMean, s, ieta, nCh); - if (n <= KFloatEpsilon) { - continue; - } - const double mpt = s.pmwk[ieta][1][1] / n; - if (std::isfinite(mpt)) { - h.mean->Fill(nCh, ieta, kMeanPt, mpt); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + const auto& s = sums[v][l][ipt]; + for (int ieta = 0; ieta < nEta; ++ieta) { + const int y = etaIdx(ipt, ieta); + const double n = s.wk[ieta][1]; + h.mean->Fill(nCh, y, kMeanN, n); + fillQMean(h.qMean, s, ieta, y, nCh); + if (n <= KFloatEpsilon) { + continue; + } + const double mpt = s.pmwk[ieta][1][1] / n; + if (std::isfinite(mpt)) { + h.mean->Fill(nCh, y, kMeanPt, mpt); + } } } const FwdSignal& fwd = (l == kGen) ? fwdTru : fwdReco; @@ -3054,23 +3199,24 @@ struct RadialFlowDecorr { for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - for (auto& s : sums[v]) { - s.reset(); - } + resetSums(v); } } // --- truth sums: independent of the track cuts, copied to every selected variation --- - PowerSums& sumsTru = sums[kSystBase][kGen]; + PtSums& sumsTru = sums[kSystBase][kGen]; if (!TESTBIT(ev.mask, kSystBase)) { - sumsTru.reset(); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + sumsTru[ipt].reset(); + } } for (const auto& particle : mcParticles) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) { continue; } float pt = particle.pt(), eta = particle.eta(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } const auto csTru = harmTerms(particle.phi()); @@ -3078,15 +3224,14 @@ struct RadialFlowDecorr { if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; } - sumsTru.add(ieta, pt, 1.); - if (cfgAziTA) { - sumsTru.addHarm(ieta, pt, 1., csTru); - } + addTrack(sumsTru, ptMask, ieta, pt, 1., csTru); } } for (int v = 1; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - sums[v][kGen] = sumsTru; + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + sums[v][kGen][ipt] = sumsTru[ipt]; + } } } @@ -3098,12 +3243,14 @@ struct RadialFlowDecorr { continue; } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } + const bool inPt0 = TESTBIT(ptMask, 0); // QA of range 0 only float eta = track.eta(), phi = track.phi(); float sign = track.sign(); - if (TESTBIT(trackMask, kSystBase)) { + if (TESTBIT(trackMask, kSystBase) && inPt0) { histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); @@ -3117,10 +3264,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(trackMask, v)) { - sums[v][kReco].add(ieta, pt, 1.); - if (cfgAziTA) { - sums[v][kReco].addHarm(ieta, pt, 1., cs); - } + addTrack(sums[v][kReco], ptMask, ieta, pt, 1., cs); } } } @@ -3130,7 +3274,7 @@ struct RadialFlowDecorr { if (weightMask == 0) { continue; } - if (TESTBIT(weightMask, kSystBase)) { + if (TESTBIT(weightMask, kSystBase) && inPt0) { fillFlat(kSystBase, vz, sign, pt, eta, phi, w[kSystBase], (1.f - fake) / eff); } @@ -3141,10 +3285,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(weightMask, v)) { - sums[v][kRecoEffCorr].add(ieta, pt, w[v]); - if (cfgAziTA) { - sums[v][kRecoEffCorr].addHarm(ieta, pt, w[v], cs); - } + addTrack(sums[v][kRecoEffCorr], ptMask, ieta, pt, w[v], cs); } } } @@ -3167,8 +3308,10 @@ struct RadialFlowDecorr { } for (int l = 0; l < kNLevels; ++l) { const FwdSignal* fwd = cfgFwdOpt.cfgFwd ? (l == kGen ? &fwdTru : &fwdReco) : nullptr; - computeEventObs(sums[v][l], v, l, fwd, nCh, evObs); - fillFluc(hists.fluc[v][l], evObs, 1., nCh); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + computeEventObs(sums[v][l][ipt], v, l, fwd, nCh, evObs, ipt); + fillFluc(hists.fluc[v][l][ipt], evObs, 1., nCh, ipt); + } } } } @@ -3214,7 +3357,7 @@ struct RadialFlowDecorr { continue; } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) { + if (pt <= KPtMin || pt > KPtMax) { // map range continue; } float eta = track.eta(), phi = track.phi(); @@ -3277,7 +3420,7 @@ struct RadialFlowDecorr { for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - sums[v][0].reset(); + resetSums(v); } } std::array w{}; @@ -3291,12 +3434,14 @@ struct RadialFlowDecorr { if (p < KFloatEpsilon) { continue; } - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } + const bool inPt0 = TESTBIT(ptMask, 0); // QA of range 0 only float eta = track.eta(), phi = track.phi(); float sign = track.sign(); - if (TESTBIT(trackMask, kSystBase)) { + if (TESTBIT(trackMask, kSystBase) && inPt0) { histos.fill(HIST("hP"), p); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); @@ -3308,7 +3453,7 @@ struct RadialFlowDecorr { if (weightMask == 0) { continue; } - if (TESTBIT(weightMask, kSystBase)) { + if (TESTBIT(weightMask, kSystBase) && inPt0) { histos.fill(HIST("Prof_Weights_pt_eta_type"), pt, eta, 0, eff); histos.fill(HIST("Prof_Weights_pt_eta_type"), pt, eta, 1, fake); histos.fill(HIST("Prof_Weights_pt_eta_type"), pt, eta, 2, flatWeight); @@ -3322,10 +3467,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(weightMask, v)) { - sums[v][0].add(ieta, pt, w[v]); - if (cfgAziTA) { - sums[v][0].addHarm(ieta, pt, w[v], cs); - } + addTrack(sums[v][0], ptMask, ieta, pt, w[v], cs); } } } @@ -3344,18 +3486,21 @@ struct RadialFlowDecorr { if (!TESTBIT(ev.mask, v)) { continue; } - const auto& s = sums[v][0]; const auto& h = hists.mean[v][0]; - for (int ieta = 0; ieta < nEta; ++ieta) { - const double n = s.wk[ieta][1]; - h.mean->Fill(nCh, ieta, kMeanN, n); - fillQMean(h.qMean, s, ieta, nCh); - if (n <= KFloatEpsilon) { - continue; - } - const double mpt = s.pmwk[ieta][1][1] / n; - if (std::isfinite(mpt)) { - h.mean->Fill(nCh, ieta, kMeanPt, mpt); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + const auto& s = sums[v][0][ipt]; + for (int ieta = 0; ieta < nEta; ++ieta) { + const int y = etaIdx(ipt, ieta); + const double n = s.wk[ieta][1]; + h.mean->Fill(nCh, y, kMeanN, n); + fillQMean(h.qMean, s, ieta, y, nCh); + if (n <= KFloatEpsilon) { + continue; + } + const double mpt = s.pmwk[ieta][1][1] / n; + if (std::isfinite(mpt)) { + h.mean->Fill(nCh, y, kMeanPt, mpt); + } } } if (cfgFwdOpt.cfgFwd && h.fwdMean) { @@ -3402,7 +3547,7 @@ struct RadialFlowDecorr { for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(ev.mask, v)) { - sums[v][0].reset(); + resetSums(v); } } std::array w{}; @@ -3416,12 +3561,14 @@ struct RadialFlowDecorr { if (p < KFloatEpsilon) { continue; } - if (pt <= cfgPtMin || pt > cfgPtMax) { + const uint32_t ptMask = ptRangeMask(pt); + if (ptMask == 0) { continue; } + const bool inPt0 = TESTBIT(ptMask, 0); // QA of range 0 only float eta = track.eta(), phi = track.phi(); float sign = track.sign(); - if (TESTBIT(trackMask, kSystBase)) { + if (TESTBIT(trackMask, kSystBase) && inPt0) { histos.fill(HIST("hP"), p); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); @@ -3433,7 +3580,7 @@ struct RadialFlowDecorr { if (weightMask == 0) { continue; } - if (TESTBIT(weightMask, kSystBase)) { + if (TESTBIT(weightMask, kSystBase) && inPt0) { fillFlat(kSystBase, vz, sign, pt, eta, phi, w[kSystBase], (1.f - fake) / eff); } @@ -3444,10 +3591,7 @@ struct RadialFlowDecorr { } for (int v = 0; v < kNSyst; ++v) { if (TESTBIT(weightMask, v)) { - sums[v][0].add(ieta, pt, w[v]); - if (cfgAziTA) { - sums[v][0].addHarm(ieta, pt, w[v], cs); - } + addTrack(sums[v][0], ptMask, ieta, pt, w[v], cs); } } } @@ -3463,8 +3607,8 @@ struct RadialFlowDecorr { // --- Poisson bootstrap: one weight per sample per event --- std::array poisW{}; - if (doBoot) { - for (int s = 0; s < nBoot; ++s) { + if (bootstrap.on) { + for (int s = 0; s < bootstrap.n; ++s) { poisW[s] = rng.Poisson(1.0); } } @@ -3475,13 +3619,15 @@ struct RadialFlowDecorr { if (!TESTBIT(ev.mask, v)) { continue; } - computeEventObs(sums[v][0], v, 0, cfgFwdOpt.cfgFwd ? &fwd : nullptr, nCh, evObs); - fillFluc(hists.fluc[v][0], evObs, 1., nCh); - // bootstrap replicas of Base, in the 5% centrality classes - if (v == kSystBase && doBoot) { - for (int s = 0; s < nBoot; ++s) { - if (poisW[s] > 0.) { - fillFluc(hists.boot[s], evObs, poisW[s], ev.cent[kSystBase]); + for (int ipt = 0; ipt < ptRange.n; ++ipt) { + computeEventObs(sums[v][0][ipt], v, 0, cfgFwdOpt.cfgFwd ? &fwd : nullptr, nCh, evObs, ipt); + fillFluc(hists.fluc[v][0][ipt], evObs, 1., nCh, ipt); + // bootstrap replicas of Base, in the 5% centrality classes + if (v == kSystBase && bootstrap.on) { + for (int s = 0; s < bootstrap.n; ++s) { + if (poisW[s] > 0.) { + fillFluc(hists.boot[s][ipt], evObs, poisW[s], ev.cent[kSystBase], ipt); + } } } }