Skip to content

[PWGCF] Add configurable pT-dependent material budget efficiency systematic for Radial Decorrelation Analysis #18198

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fsii merged 2 commits into from
Oct 7, 2026
Merged

[PWGCF] Add configurable pT-dependent material budget efficiency systematic for Radial Decorrelation Analysis #18198

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6093b67
[PWGCF] Add configurable pT-dependent material budget efficiency syst…
somabhatta Oct 6, 2026
da7d395
[PWGCF] fixing sumgap bins
somabhatta Oct 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 60 additions & 36 deletions PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,6 +187,8 @@ struct RadialFlowDecorr {

Configurable<bool> cfgFlat{"cfgFlat", false, "Whether to use flattening weights"};
Configurable<bool> cfgEff{"cfgEff", false, "Whether to use Efficiency weights"};
Configurable<float> cfgEffVarPercent{"cfgEffVarPercent", 0.f, "Signed relative efficiency variation (%): +10 or -10"};
Configurable<bool> cfgEffVarPtDep{"cfgEffVarPtDep", false, "false: constant variation; true: linear decrease with pT"};
Configurable<bool> cfgZDC{"cfgZDC", false, "Whether to use ZDC for pileup histograms"};

Configurable<std::string> cfgCCDBurl{"cfgCCDBurl", "https://alice-ccdb.cern.ch", "ccdb url"};
Expand All @@ -213,13 +215,8 @@ struct RadialFlowDecorr {
// binning and are rebuilt in init(). Placeholders here.
AxisSpec etaAxis{9, -0.9, 0.9, "#eta"};
AxisSpec etaBinAxis{10, -0.5, 9.5, "#eta bin Number"};

AxisSpec gapAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
"Gap"};
AxisSpec sumAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
"Sum"};
AxisSpec gapAxis{15, -1.5, 1.5, "#Delta#eta"};
AxisSpec sumAxis{15, -1.5, 1.5, "#Sigma#eta"};

// --- process switches --------------------------------------------------------
Configurable<bool> cfgRunGetEff{"cfgRunGetEff", false, "Run MC pass to build efficiency/fake maps"};
Expand Down Expand Up @@ -491,15 +488,36 @@ struct RadialFlowDecorr {
return static_cast<float>(val);
}

// Reject missing or invalid efficiency/fake corrections consistently in all passes.
float efficiencyVariationFactor(float pt) const
{
if (cfgSystType.value != kSystEff) {
return 1.f;
}
// Empirical shape: full variation at cfgPtMin, zero at cfgPtMax.
const float shape = cfgEffVarPtDep.value ? (cfgPtMax.value - pt) / (cfgPtMax.value - cfgPtMin.value) : 1.f;
return 1.f + 0.01f * cfgEffVarPercent.value * shape;
}

bool getValidEffFake(float mult, float pt, float eta, bool useEff,
float& eff, float& fake) const
{
eff = getEfficiency(mult, pt, eta, 0, useEff);
fake = getEfficiency(mult, pt, eta, 1, useEff);
return std::isfinite(eff) && std::isfinite(fake) &&
eff > KFloatEpsilon && eff <= 1.f &&
fake >= 0.f && fake < 1.f;

// Validate the nominal corrections before applying any variation.
if (!std::isfinite(eff) || !std::isfinite(fake) ||
eff <= KFloatEpsilon || eff > 1.f ||
fake < 0.f || fake >= 1.f) {
return false;
}

// Apply the variation once, in DATA mean/fluctuation jobs.
if (useEff && (cfgRunDataMean.value || cfgRunDataFluc.value)) {
eff *= efficiencyVariationFactor(pt);
}

// The varied efficiency is an effective correction denominator.
return std::isfinite(eff) && eff > KFloatEpsilon;
}

float getFlatteningWeight(float vz, float chg, float pt, float eta, float phi, bool useFlat) const
Expand Down Expand Up @@ -1073,8 +1091,12 @@ struct RadialFlowDecorr {
}
etaAxis = AxisSpec{obsEdges, "#eta"};
etaBinAxis = AxisSpec{nEta + 1, -0.5, static_cast<double>(nEta) + 0.5, "#eta bin Number"};
LOGF(info, "Observable eta binning (%s): %d bins of width %.2f over |eta|<%.2f (+ reference), nEta=%d",
isMcRun ? "MC pinned" : "DATA", nbins, width, halfEta, nEta);
// (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;
const double gsMax = (nbins - 0.5) * width;
gapAxis = AxisSpec{nGS, -gsMax, gsMax, "#Delta#eta"};
sumAxis = AxisSpec{nGS, -gsMax, gsMax, "#Sigma#eta"};
}

// bootstrap active only for the base data fluctuation pass
Expand Down Expand Up @@ -1487,7 +1509,7 @@ struct RadialFlowDecorr {

histos.fill(HIST("MCGen/Prof_Cent_Nchrec"), cent, sumWiTruth[0]);
histos.fill(HIST("MCGen/Prof_Mult_Nchrec"), multPV, sumWiTruth[0]);
if (sumWiTruth[0] > 1.0f) {
if (sumWiTruth[0] > KFloatEpsilon) {
histos.fill(HIST("MCGen/Prof_Cent_MeanpT"), cent, sumWiptiTruth[0] / sumWiTruth[0]);
histos.fill(HIST("MCGen/Prof_Mult_MeanpT"), multPV, sumWiptiTruth[0] / sumWiTruth[0]);
}
Expand Down Expand Up @@ -1554,26 +1576,26 @@ struct RadialFlowDecorr {
float nRecoAB = sumWiReco[ietaA] + sumWiReco[ietaC];
float nCorrAB = sumWiRecoEffCorr[ietaA] + sumWiRecoEffCorr[ietaC];

if (nTruAB > 0) {
if (sumWiTruth[ietaA] > KFloatEpsilon && sumWiTruth[ietaC] > KFloatEpsilon) {
histos.fill(HIST("Prof2D_MeanpTSub_Tru"), cent, ietaA, ietaC, (sumWiptiTruth[ietaA] + sumWiptiTruth[ietaC]) / nTruAB);
}
if (nRecoAB > 0) {
if (sumWiReco[ietaA] > KFloatEpsilon && sumWiReco[ietaC] > KFloatEpsilon) {
histos.fill(HIST("Prof2D_MeanpTSub_Reco"), cent, ietaA, ietaC, (sumWiptiReco[ietaA] + sumWiptiReco[ietaC]) / nRecoAB);
}
if (nCorrAB > 0) {
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon && sumWiRecoEffCorr[ietaC] > KFloatEpsilon) {
histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr"), cent, ietaA, ietaC, (sumWiptiRecoEffCorr[ietaA] + sumWiptiRecoEffCorr[ietaC]) / nCorrAB);
}
}

if (sumWiTruth[ietaA] > 0) {
if (sumWiTruth[ietaA] > KFloatEpsilon) {
histos.fill(HIST("pmeanTru_nch_etabin"), multPV, ietaA, sumWiptiTruth[ietaA] / sumWiTruth[ietaA]);
histos.fill(HIST("pmeanMultTru_nch_etabin"), multPV, ietaA, sumWiTruth[ietaA]);
}
if (sumWiReco[ietaA] > 0) {
if (sumWiReco[ietaA] > KFloatEpsilon) {
histos.fill(HIST("pmeanReco_nch_etabin"), multPV, ietaA, sumWiptiReco[ietaA] / sumWiReco[ietaA]);
histos.fill(HIST("pmeanMultReco_nch_etabin"), multPV, ietaA, sumWiReco[ietaA]);
}
if (sumWiRecoEffCorr[ietaA] > 0) {
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon) {
histos.fill(HIST("pmeanRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiptiRecoEffCorr[ietaA] / sumWiRecoEffCorr[ietaA]);
histos.fill(HIST("pmeanMultRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiRecoEffCorr[ietaA]);
}
Expand Down Expand Up @@ -1706,7 +1728,7 @@ struct RadialFlowDecorr {
float mmptRecoEffCor = state.pmeanRecoEffcorrNchEtabinStep2->GetBinContent(ibx, iby);

// truth
meanTru[ieta] = (sumWkTru[ieta][1] >= 1.0) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
meanTru[ieta] = (sumWkTru[ieta][1] > KFloatEpsilon) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
c2Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
c3Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
p1kBarTru[ieta] = std::numeric_limits<double>::quiet_NaN();
Expand All @@ -1719,7 +1741,7 @@ struct RadialFlowDecorr {
}
}
// reco
meanReco[ieta] = (sumWkReco[ieta][1] >= 1.0) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
meanReco[ieta] = (sumWkReco[ieta][1] > KFloatEpsilon) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
c2Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
c3Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
p1kBarReco[ieta] = std::numeric_limits<double>::quiet_NaN();
Expand All @@ -1732,7 +1754,7 @@ struct RadialFlowDecorr {
}
}
// reco, efficiency-corrected
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] >= 1.0) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
c2RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
c3RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
p1kBarRecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
Expand All @@ -1755,15 +1777,15 @@ struct RadialFlowDecorr {
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkRecoEffCor[ieta][1]);
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkRecoEffCor[ieta][1]);

if (sumWkTru[ieta][1] > 1.0f) {
if (sumWkTru[ieta][1] > KFloatEpsilon) {
histos.fill(HIST("MCGen/Prof_Cent_NEta_MeanpT"), cent, ieta, meanTru[ieta]);
histos.fill(HIST("MCGen/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanTru[ieta]);
}
if (sumWkReco[ieta][1] > 1.0f) {
if (sumWkReco[ieta][1] > KFloatEpsilon) {
histos.fill(HIST("MCReco/Prof_Cent_NEta_MeanpT"), cent, ieta, meanReco[ieta]);
histos.fill(HIST("MCReco/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanReco[ieta]);
}
if (sumWkRecoEffCor[ieta][1] > 1.0f) {
if (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) {
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_MeanpT"), cent, ieta, meanRecoEffCor[ieta]);
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanRecoEffCor[ieta]);
}
Expand Down Expand Up @@ -1867,8 +1889,9 @@ struct RadialFlowDecorr {
for (int ietaC = 1; ietaC < nEta; ++ietaC) {
float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f;
float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f;
float gap = etaValA - etaValB;
float sum = (etaValA + etaValB);
const float w = etaUp[ietaA] - etaLw[ietaA];
const float gap = (ietaA - ietaC) * w;
const float sum = 2.f * etaLw[1] + (ietaA + ietaC - 1) * w; // = cA + cC on the lattice

float c2SubTru = (ietaA == ietaC) ? static_cast<float>(c2Tru[ietaA]) : p1kBarTru[ietaA] * p1kBarTru[ietaC];
float c2SubReco = (ietaA == ietaC) ? static_cast<float>(c2Reco[ietaA]) : p1kBarReco[ietaA] * p1kBarReco[ietaC];
Expand Down Expand Up @@ -2059,7 +2082,7 @@ struct RadialFlowDecorr {
}
}

if (sumWi[0] >= 1.0f) {
if (sumWi[0] > KFloatEpsilon) {
histos.fill(HIST("Prof_Cent_Nchrec"), cent, sumWi[0]);
histos.fill(HIST("Prof_Mult_Nchrec"), coll.multNTracksPV(), sumWi[0]);
histos.fill(HIST("Prof_Cent_MeanpT"), cent, sumWipti[0] / sumWi[0]);
Expand All @@ -2068,7 +2091,7 @@ struct RadialFlowDecorr {

for (int ietaA = 0; ietaA < nEta; ++ietaA) {
for (int ietaC = 0; ietaC < nEta; ++ietaC) {
if ((sumWi[ietaA] < 1.0f) || (sumWi[ietaC] < 1.0f)) {
if ((sumWi[ietaA] < KFloatEpsilon) || (sumWi[ietaC] < KFloatEpsilon)) {
continue;
}
double wCorrAB = sumWi[ietaA] + sumWi[ietaC];
Expand All @@ -2077,7 +2100,7 @@ struct RadialFlowDecorr {
histos.fill(HIST("Prof2D_MeanpTSub"), cent, ietaA, ietaC, mptsub);
}
}
if (sumWi[ietaA] >= 1.0f) {
if (sumWi[ietaA] > KFloatEpsilon) {
double mpt = sumWipti[ietaA] / sumWi[ietaA];
if (std::isfinite(mpt)) {
histos.fill(HIST("pmean_nch_etabin"), coll.multNTracksPV(), ietaA, mpt);
Expand Down Expand Up @@ -2177,7 +2200,9 @@ struct RadialFlowDecorr {
if (!getValidEffFake(coll.multNTracksPV(), pt, eta, cfgEff, eff, fake)) {
continue;
}

float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat);

float w = flatWeight * (1.0f - fake) / eff;
if (!std::isfinite(w) || w <= 0.f) {
continue;
Expand All @@ -2203,9 +2228,7 @@ struct RadialFlowDecorr {
float mmpt = state.pmeanNchEtabinStep2->GetBinContent(ibx, iby);

meanMult[ieta] = sumwk[ieta][1];

// mean pT of the bin: valid only with >=1 (weighted) track
mean[ieta] = (sumwk[ieta][1] >= 1.0) ? (sumpmwk[ieta][1][1] / sumwk[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
mean[ieta] = (sumwk[ieta][1] > KFloatEpsilon) ? (sumpmwk[ieta][1][1] / sumwk[ieta][1]) : std::numeric_limits<double>::quiet_NaN();

// c2 (>=2 tracks), c3 (>=3 tracks) and the pT deviation need a usable reference
c2[ieta] = std::numeric_limits<double>::quiet_NaN();
Expand Down Expand Up @@ -2285,8 +2308,9 @@ struct RadialFlowDecorr {
for (int ietaC = 1; ietaC < nEta; ++ietaC) {
float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f;
float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f;
float gap = etaValA - etaValB;
float sum = (etaValA + etaValB);
const float w = etaUp[ietaA] - etaLw[ietaA];
const float gap = (ietaA - ietaC) * w;
const float sum = 2.f * etaLw[1] + (ietaA + ietaC - 1) * w; // = cA + cC on the lattice

float c2Sub = (ietaA == ietaC) ? static_cast<float>(c2[ietaA]) : p1kBar[ietaA] * p1kBar[ietaC];
// C3 sub-event map, ORDERED: cell (A,C) = c2[A] p1kBar[C] = <dpT dpT>_A <dpT>_C,
Expand Down
Loading