Skip to content

Commit 6093b67

Browse files
somabhattasomabhatta
committed
[PWGCF] Add configurable pT-dependent material budget efficiency systematic for Radial Decorrelation code
1 parent 2e196f9 commit 6093b67

1 file changed

Lines changed: 46 additions & 23 deletions

File tree

‎PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx‎

Lines changed: 46 additions & 23 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"};
@@ -491,15 +493,36 @@ struct RadialFlowDecorr {
491493
return static_cast<float>(val);
492494
}
493495

494-
// Reject missing or invalid efficiency/fake corrections consistently in all passes.
496+
float efficiencyVariationFactor(float pt) const
497+
{
498+
if (cfgSystType.value != kSystEff) {
499+
return 1.f;
500+
}
501+
// Empirical shape: full variation at cfgPtMin, zero at cfgPtMax.
502+
const float shape = cfgEffVarPtDep.value ? (cfgPtMax.value - pt) / (cfgPtMax.value - cfgPtMin.value) : 1.f;
503+
return 1.f + 0.01f * cfgEffVarPercent.value * shape;
504+
}
505+
495506
bool getValidEffFake(float mult, float pt, float eta, bool useEff,
496507
float& eff, float& fake) const
497508
{
498509
eff = getEfficiency(mult, pt, eta, 0, useEff);
499510
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;
511+
512+
// Validate the nominal corrections before applying any variation.
513+
if (!std::isfinite(eff) || !std::isfinite(fake) ||
514+
eff <= KFloatEpsilon || eff > 1.f ||
515+
fake < 0.f || fake >= 1.f) {
516+
return false;
517+
}
518+
519+
// Apply the variation once, in DATA mean/fluctuation jobs.
520+
if (useEff && (cfgRunDataMean.value || cfgRunDataFluc.value)) {
521+
eff *= efficiencyVariationFactor(pt);
522+
}
523+
524+
// The varied efficiency is an effective correction denominator.
525+
return std::isfinite(eff) && eff > KFloatEpsilon;
503526
}
504527

505528
float getFlatteningWeight(float vz, float chg, float pt, float eta, float phi, bool useFlat) const
@@ -1487,7 +1510,7 @@ struct RadialFlowDecorr {
14871510

14881511
histos.fill(HIST("MCGen/Prof_Cent_Nchrec"), cent, sumWiTruth[0]);
14891512
histos.fill(HIST("MCGen/Prof_Mult_Nchrec"), multPV, sumWiTruth[0]);
1490-
if (sumWiTruth[0] > 1.0f) {
1513+
if (sumWiTruth[0] > KFloatEpsilon) {
14911514
histos.fill(HIST("MCGen/Prof_Cent_MeanpT"), cent, sumWiptiTruth[0] / sumWiTruth[0]);
14921515
histos.fill(HIST("MCGen/Prof_Mult_MeanpT"), multPV, sumWiptiTruth[0] / sumWiTruth[0]);
14931516
}
@@ -1554,26 +1577,26 @@ struct RadialFlowDecorr {
15541577
float nRecoAB = sumWiReco[ietaA] + sumWiReco[ietaC];
15551578
float nCorrAB = sumWiRecoEffCorr[ietaA] + sumWiRecoEffCorr[ietaC];
15561579

1557-
if (nTruAB > 0) {
1580+
if (sumWiTruth[ietaA] > KFloatEpsilon && sumWiTruth[ietaC] > KFloatEpsilon) {
15581581
histos.fill(HIST("Prof2D_MeanpTSub_Tru"), cent, ietaA, ietaC, (sumWiptiTruth[ietaA] + sumWiptiTruth[ietaC]) / nTruAB);
15591582
}
1560-
if (nRecoAB > 0) {
1583+
if (sumWiReco[ietaA] > KFloatEpsilon && sumWiReco[ietaC] > KFloatEpsilon) {
15611584
histos.fill(HIST("Prof2D_MeanpTSub_Reco"), cent, ietaA, ietaC, (sumWiptiReco[ietaA] + sumWiptiReco[ietaC]) / nRecoAB);
15621585
}
1563-
if (nCorrAB > 0) {
1586+
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon && sumWiRecoEffCorr[ietaC] > KFloatEpsilon) {
15641587
histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr"), cent, ietaA, ietaC, (sumWiptiRecoEffCorr[ietaA] + sumWiptiRecoEffCorr[ietaC]) / nCorrAB);
15651588
}
15661589
}
15671590

1568-
if (sumWiTruth[ietaA] > 0) {
1591+
if (sumWiTruth[ietaA] > KFloatEpsilon) {
15691592
histos.fill(HIST("pmeanTru_nch_etabin"), multPV, ietaA, sumWiptiTruth[ietaA] / sumWiTruth[ietaA]);
15701593
histos.fill(HIST("pmeanMultTru_nch_etabin"), multPV, ietaA, sumWiTruth[ietaA]);
15711594
}
1572-
if (sumWiReco[ietaA] > 0) {
1595+
if (sumWiReco[ietaA] > KFloatEpsilon) {
15731596
histos.fill(HIST("pmeanReco_nch_etabin"), multPV, ietaA, sumWiptiReco[ietaA] / sumWiReco[ietaA]);
15741597
histos.fill(HIST("pmeanMultReco_nch_etabin"), multPV, ietaA, sumWiReco[ietaA]);
15751598
}
1576-
if (sumWiRecoEffCorr[ietaA] > 0) {
1599+
if (sumWiRecoEffCorr[ietaA] > KFloatEpsilon) {
15771600
histos.fill(HIST("pmeanRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiptiRecoEffCorr[ietaA] / sumWiRecoEffCorr[ietaA]);
15781601
histos.fill(HIST("pmeanMultRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiRecoEffCorr[ietaA]);
15791602
}
@@ -1706,7 +1729,7 @@ struct RadialFlowDecorr {
17061729
float mmptRecoEffCor = state.pmeanRecoEffcorrNchEtabinStep2->GetBinContent(ibx, iby);
17071730

17081731
// truth
1709-
meanTru[ieta] = (sumWkTru[ieta][1] >= 1.0) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1732+
meanTru[ieta] = (sumWkTru[ieta][1] > KFloatEpsilon) ? (sumPmwkTru[ieta][1][1] / sumWkTru[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17101733
c2Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
17111734
c3Tru[ieta] = std::numeric_limits<double>::quiet_NaN();
17121735
p1kBarTru[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1719,7 +1742,7 @@ struct RadialFlowDecorr {
17191742
}
17201743
}
17211744
// reco
1722-
meanReco[ieta] = (sumWkReco[ieta][1] >= 1.0) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1745+
meanReco[ieta] = (sumWkReco[ieta][1] > KFloatEpsilon) ? (sumPmwkReco[ieta][1][1] / sumWkReco[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17231746
c2Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
17241747
c3Reco[ieta] = std::numeric_limits<double>::quiet_NaN();
17251748
p1kBarReco[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1732,7 +1755,7 @@ struct RadialFlowDecorr {
17321755
}
17331756
}
17341757
// reco, efficiency-corrected
1735-
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] >= 1.0) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
1758+
meanRecoEffCor[ieta] = (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) ? (sumPmwkRecoEffCor[ieta][1][1] / sumWkRecoEffCor[ieta][1]) : std::numeric_limits<double>::quiet_NaN();
17361759
c2RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
17371760
c3RecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
17381761
p1kBarRecoEffCor[ieta] = std::numeric_limits<double>::quiet_NaN();
@@ -1755,15 +1778,15 @@ struct RadialFlowDecorr {
17551778
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkRecoEffCor[ieta][1]);
17561779
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkRecoEffCor[ieta][1]);
17571780

1758-
if (sumWkTru[ieta][1] > 1.0f) {
1781+
if (sumWkTru[ieta][1] > KFloatEpsilon) {
17591782
histos.fill(HIST("MCGen/Prof_Cent_NEta_MeanpT"), cent, ieta, meanTru[ieta]);
17601783
histos.fill(HIST("MCGen/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanTru[ieta]);
17611784
}
1762-
if (sumWkReco[ieta][1] > 1.0f) {
1785+
if (sumWkReco[ieta][1] > KFloatEpsilon) {
17631786
histos.fill(HIST("MCReco/Prof_Cent_NEta_MeanpT"), cent, ieta, meanReco[ieta]);
17641787
histos.fill(HIST("MCReco/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanReco[ieta]);
17651788
}
1766-
if (sumWkRecoEffCor[ieta][1] > 1.0f) {
1789+
if (sumWkRecoEffCor[ieta][1] > KFloatEpsilon) {
17671790
histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_MeanpT"), cent, ieta, meanRecoEffCor[ieta]);
17681791
histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanRecoEffCor[ieta]);
17691792
}
@@ -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();

0 commit comments

Comments
 (0)