Skip to content

Commit c05dd43

Browse files
somabhattasomabhatta
authored
[PWGCF] Add configurable pT-dependent material budget efficiency systematic for Radial Decorrelation Analysis (#18198)
1 parent 4aa762b commit c05dd43

1 file changed

Lines changed: 60 additions & 36 deletions

File tree

‎PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx‎

Lines changed: 60 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -187,6 +187,8 @@ struct RadialFlowDecorr {
187187

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

192194
Configurable<std::string> cfgCCDBurl{"cfgCCDBurl", "https://alice-ccdb.cern.ch", "ccdb url"};
@@ -213,13 +215,8 @@ struct RadialFlowDecorr {
213215
// binning and are rebuilt in init(). Placeholders here.
214216
AxisSpec etaAxis{9, -0.9, 0.9, "#eta"};
215217
AxisSpec etaBinAxis{10, -0.5, 9.5, "#eta bin Number"};
216-
217-
AxisSpec gapAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
218-
0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
219-
"Gap"};
220-
AxisSpec sumAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
221-
0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
222-
"Sum"};
218+
AxisSpec gapAxis{15, -1.5, 1.5, "#Delta#eta"};
219+
AxisSpec sumAxis{15, -1.5, 1.5, "#Sigma#eta"};
223220

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

494-
// Reject missing or invalid efficiency/fake corrections consistently in all passes.
491+
float efficiencyVariationFactor(float pt) const
492+
{
493+
if (cfgSystType.value != kSystEff) {
494+
return 1.f;
495+
}
496+
// Empirical shape: full variation at cfgPtMin, zero at cfgPtMax.
497+
const float shape = cfgEffVarPtDep.value ? (cfgPtMax.value - pt) / (cfgPtMax.value - cfgPtMin.value) : 1.f;
498+
return 1.f + 0.01f * cfgEffVarPercent.value * shape;
499+
}
500+
495501
bool getValidEffFake(float mult, float pt, float eta, bool useEff,
496502
float& eff, float& fake) const
497503
{
498504
eff = getEfficiency(mult, pt, eta, 0, useEff);
499505
fake = getEfficiency(mult, pt, eta, 1, useEff);
500-
return std::isfinite(eff) && std::isfinite(fake) &&
501-
eff > KFloatEpsilon && eff <= 1.f &&
502-
fake >= 0.f && fake < 1.f;
506+
507+
// Validate the nominal corrections before applying any variation.
508+
if (!std::isfinite(eff) || !std::isfinite(fake) ||
509+
eff <= KFloatEpsilon || eff > 1.f ||
510+
fake < 0.f || fake >= 1.f) {
511+
return false;
512+
}
513+
514+
// Apply the variation once, in DATA mean/fluctuation jobs.
515+
if (useEff && (cfgRunDataMean.value || cfgRunDataFluc.value)) {
516+
eff *= efficiencyVariationFactor(pt);
517+
}
518+
519+
// The varied efficiency is an effective correction denominator.
520+
return std::isfinite(eff) && eff > KFloatEpsilon;
503521
}
504522

505523
float getFlatteningWeight(float vz, float chg, float pt, float eta, float phi, bool useFlat) const
@@ -1073,8 +1091,12 @@ struct RadialFlowDecorr {
10731091
}
10741092
etaAxis = AxisSpec{obsEdges, "#eta"};
10751093
etaBinAxis = AxisSpec{nEta + 1, -0.5, static_cast<double>(nEta) + 0.5, "#eta bin Number"};
1076-
LOGF(info, "Observable eta binning (%s): %d bins of width %.2f over |eta|<%.2f (+ reference), nEta=%d",
1077-
isMcRun ? "MC pinned" : "DATA", nbins, width, halfEta, nEta);
1094+
// (Delta eta, Sigma eta): fills are integer multiples of w, so 2N-1 bins
1095+
// of width w centred on them, out to (N-1/2) w.
1096+
const int nGS = 2 * nbins - 1;
1097+
const double gsMax = (nbins - 0.5) * width;
1098+
gapAxis = AxisSpec{nGS, -gsMax, gsMax, "#Delta#eta"};
1099+
sumAxis = AxisSpec{nGS, -gsMax, gsMax, "#Sigma#eta"};
10781100
}
10791101

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

14881510
histos.fill(HIST("MCGen/Prof_Cent_Nchrec"), cent, sumWiTruth[0]);
14891511
histos.fill(HIST("MCGen/Prof_Mult_Nchrec"), multPV, sumWiTruth[0]);
1490-
if (sumWiTruth[0] > 1.0f) {
1512+
if (sumWiTruth[0] > KFloatEpsilon) {
14911513
histos.fill(HIST("MCGen/Prof_Cent_MeanpT"), cent, sumWiptiTruth[0] / sumWiTruth[0]);
14921514
histos.fill(HIST("MCGen/Prof_Mult_MeanpT"), multPV, sumWiptiTruth[0] / sumWiTruth[0]);
14931515
}
@@ -1554,26 +1576,26 @@ struct RadialFlowDecorr {
15541576
float nRecoAB = sumWiReco[ietaA] + sumWiReco[ietaC];
15551577
float nCorrAB = sumWiRecoEffCorr[ietaA] + sumWiRecoEffCorr[ietaC];
15561578

1557-
if (nTruAB > 0) {
1579+
if (sumWiTruth[ietaA] > KFloatEpsilon && sumWiTruth[ietaC] > KFloatEpsilon) {
15581580
histos.fill(HIST("Prof2D_MeanpTSub_Tru"), cent, ietaA, ietaC, (sumWiptiTruth[ietaA] + sumWiptiTruth[ietaC]) / nTruAB);
15591581
}
1560-
if (nRecoAB > 0) {
1582+
if (sumWiReco[ietaA] > KFloatEpsilon && sumWiReco[ietaC] > KFloatEpsilon) {
15611583
histos.fill(HIST("Prof2D_MeanpTSub_Reco"), cent, ietaA, ietaC, (sumWiptiReco[ietaA] + sumWiptiReco[ietaC]) / nRecoAB);
15621584
}
1563-
if (nCorrAB > 0) {
1585+
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon && sumWiRecoEffCorr[ietaC] > KFloatEpsilon) {
15641586
histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr"), cent, ietaA, ietaC, (sumWiptiRecoEffCorr[ietaA] + sumWiptiRecoEffCorr[ietaC]) / nCorrAB);
15651587
}
15661588
}
15671589

1568-
if (sumWiTruth[ietaA] > 0) {
1590+
if (sumWiTruth[ietaA] > KFloatEpsilon) {
15691591
histos.fill(HIST("pmeanTru_nch_etabin"), multPV, ietaA, sumWiptiTruth[ietaA] / sumWiTruth[ietaA]);
15701592
histos.fill(HIST("pmeanMultTru_nch_etabin"), multPV, ietaA, sumWiTruth[ietaA]);
15711593
}
1572-
if (sumWiReco[ietaA] > 0) {
1594+
if (sumWiReco[ietaA] > KFloatEpsilon) {
15731595
histos.fill(HIST("pmeanReco_nch_etabin"), multPV, ietaA, sumWiptiReco[ietaA] / sumWiReco[ietaA]);
15741596
histos.fill(HIST("pmeanMultReco_nch_etabin"), multPV, ietaA, sumWiReco[ietaA]);
15751597
}
1576-
if (sumWiRecoEffCorr[ietaA] > 0) {
1598+
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon) {
15771599
histos.fill(HIST("pmeanRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiptiRecoEffCorr[ietaA] / sumWiRecoEffCorr[ietaA]);
15781600
histos.fill(HIST("pmeanMultRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiRecoEffCorr[ietaA]);
15791601
}
@@ -1706,7 +1728,7 @@ struct RadialFlowDecorr {
17061728
float mmptRecoEffCor = state.pmeanRecoEffcorrNchEtabinStep2->GetBinContent(ibx, iby);
17071729

17081730
// truth
1709-
meanTru[ieta] = (sumWkTru[ieta][1] >= 1.0) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1731+
meanTru[ieta] = (sumWkTru[ieta][1] > KFloatEpsilon) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17101732
c2Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
17111733
c3Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
17121734
p1kBarTru[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1719,7 +1741,7 @@ struct RadialFlowDecorr {
17191741
}
17201742
}
17211743
// reco
1722-
meanReco[ieta] = (sumWkReco[ieta][1] >= 1.0) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1744+
meanReco[ieta] = (sumWkReco[ieta][1] > KFloatEpsilon) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17231745
c2Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
17241746
c3Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
17251747
p1kBarReco[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1732,7 +1754,7 @@ struct RadialFlowDecorr {
17321754
}
17331755
}
17341756
// reco, efficiency-corrected
1735-
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] >= 1.0) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1757+
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17361758
c2RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
17371759
c3RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
17381760
p1kBarRecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1755,15 +1777,15 @@ struct RadialFlowDecorr {
17551777
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkRecoEffCor[ieta][1]);
17561778
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkRecoEffCor[ieta][1]);
17571779

1758-
if (sumWkTru[ieta][1] > 1.0f) {
1780+
if (sumWkTru[ieta][1] > KFloatEpsilon) {
17591781
histos.fill(HIST("MCGen/Prof_Cent_NEta_MeanpT"), cent, ieta, meanTru[ieta]);
17601782
histos.fill(HIST("MCGen/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanTru[ieta]);
17611783
}
1762-
if (sumWkReco[ieta][1] > 1.0f) {
1784+
if (sumWkReco[ieta][1] > KFloatEpsilon) {
17631785
histos.fill(HIST("MCReco/Prof_Cent_NEta_MeanpT"), cent, ieta, meanReco[ieta]);
17641786
histos.fill(HIST("MCReco/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanReco[ieta]);
17651787
}
1766-
if (sumWkRecoEffCor[ieta][1] > 1.0f) {
1788+
if (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) {
17671789
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_MeanpT"), cent, ieta, meanRecoEffCor[ieta]);
17681790
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanRecoEffCor[ieta]);
17691791
}
@@ -1867,8 +1889,9 @@ struct RadialFlowDecorr {
18671889
for (int ietaC = 1; ietaC < nEta; ++ietaC) {
18681890
float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f;
18691891
float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f;
1870-
float gap = etaValA - etaValB;
1871-
float sum = (etaValA + etaValB);
1892+
const float w = etaUp[ietaA] - etaLw[ietaA];
1893+
const float gap = (ietaA - ietaC) * w;
1894+
const float sum = 2.f * etaLw[1] + (ietaA + ietaC - 1) * w; // = cA + cC on the lattice
18721895

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

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

20692092
for (int ietaA = 0; ietaA < nEta; ++ietaA) {
20702093
for (int ietaC = 0; ietaC < nEta; ++ietaC) {
2071-
if ((sumWi[ietaA] < 1.0f) || (sumWi[ietaC] < 1.0f)) {
2094+
if ((sumWi[ietaA] < KFloatEpsilon) || (sumWi[ietaC] < KFloatEpsilon)) {
20722095
continue;
20732096
}
20742097
double wCorrAB = sumWi[ietaA] + sumWi[ietaC];
@@ -2077,7 +2100,7 @@ struct RadialFlowDecorr {
20772100
histos.fill(HIST("Prof2D_MeanpTSub"), cent, ietaA, ietaC, mptsub);
20782101
}
20792102
}
2080-
if (sumWi[ietaA] >= 1.0f) {
2103+
if (sumWi[ietaA] > KFloatEpsilon) {
20812104
double mpt = sumWipti[ietaA] / sumWi[ietaA];
20822105
if (std::isfinite(mpt)) {
20832106
histos.fill(HIST("pmean_nch_etabin"), coll.multNTracksPV(), ietaA, mpt);
@@ -2177,7 +2200,9 @@ struct RadialFlowDecorr {
21772200
if (!getValidEffFake(coll.multNTracksPV(), pt, eta, cfgEff, eff, fake)) {
21782201
continue;
21792202
}
2203+
21802204
float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat);
2205+
21812206
float w = flatWeight * (1.0f - fake) / eff;
21822207
if (!std::isfinite(w) || w <= 0.f) {
21832208
continue;
@@ -2203,9 +2228,7 @@ struct RadialFlowDecorr {
22032228
float mmpt = state.pmeanNchEtabinStep2->GetBinContent(ibx, iby);
22042229

22052230
meanMult[ieta] = sumwk[ieta][1];
2206-
2207-
// mean pT of the bin: valid only with >=1 (weighted) track
2208-
mean[ieta] = (sumwk[ieta][1] >= 1.0) ? (sumpmwk[ieta][1][1] / sumwk[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
2231+
mean[ieta] = (sumwk[ieta][1] > KFloatEpsilon) ? (sumpmwk[ieta][1][1] / sumwk[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
22092232

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

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

0 commit comments

Comments
 (0)