Skip to content

Commit 9df3864

Browse files
EmilGormEmilGorm
committed
fix pid efficiency logic, counting logic
1 parent 9f1ece1 commit 9df3864

File tree

1 file changed

+107
-67
lines changed

1 file changed

+107
-67
lines changed

‎PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx‎

Lines changed: 107 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -137,11 +137,13 @@ struct FlowGenericFramework {
137137
O2_DEFINE_CONFIGURABLE(cfgUsePID, bool, true, "Enable PID information")
138138
O2_DEFINE_CONFIGURABLE(cfgUseGapMethod, bool, false, "Use gap method in vn-pt calculations")
139139
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
140+
O2_DEFINE_CONFIGURABLE(cfgUsePIDEfficiencies, bool, false, "Use species dependent efficiencies")
140141
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
141142
O2_DEFINE_CONFIGURABLE(cfgPtmin, float, 0.2, "minimum pt (GeV/c)");
142143
O2_DEFINE_CONFIGURABLE(cfgPtmax, float, 10, "maximum pt (GeV/c)");
143144
O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut");
144145
O2_DEFINE_CONFIGURABLE(cfgEtaPtPt, float, 0.4, "eta cut for pt-pt correlations");
146+
O2_DEFINE_CONFIGURABLE(cfgEtaNch, float, 0.4, "eta cut for nch selection");
145147
O2_DEFINE_CONFIGURABLE(cfgUsePIDTotal, bool, false, "use fraction of PID total");
146148
O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)");
147149
struct : ConfigurableGroup {
@@ -223,6 +225,7 @@ struct FlowGenericFramework {
223225

224226
struct Config {
225227
TH1D* mEfficiency = nullptr;
228+
std::vector<TH1D*> mPIDEfficiencies;
226229
std::vector<GFWWeights*> mAcceptance;
227230
bool correctionsLoaded = false;
228231
} cfg;
@@ -845,11 +848,21 @@ struct FlowGenericFramework {
845848
}
846849
}
847850
if (!cfgEfficiency.value.empty()) {
848-
cfg.mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
849-
if (cfg.mEfficiency == nullptr) {
850-
LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str());
851+
if (!cfgUsePIDEfficiencies) {
852+
cfg.mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
853+
if (cfg.mEfficiency == nullptr) {
854+
LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str());
855+
}
856+
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency);
857+
} else {
858+
std::vector<std::string> species = {"ch", "pi", "ka", "pr", "k0", "lambda"};
859+
for (const auto& sp : species) {
860+
cfg.mPIDEfficiencies.push_back(ccdb->getForTimeStamp<TH1D>(cfgEfficiency.value + "/" + sp, timestamp));
861+
if (cfg.mPIDEfficiencies.back() == nullptr)
862+
LOGF(fatal, "Could not load PID efficiency histograms from %s", cfgEfficiency.value + "/" + sp);
863+
LOGF(info, "Loaded PID efficiency histogram from %s", cfgEfficiency.value + "/" + sp);
864+
}
851865
}
852-
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency);
853866
}
854867
cfg.correctionsLoaded = true;
855868
}
@@ -864,11 +877,15 @@ struct FlowGenericFramework {
864877
}
865878

866879
template <typename TTrack>
867-
double getEfficiency(TTrack track)
880+
double getEfficiency(TTrack track, int pidIndex = 0)
868881
{ //-1 ref, 0 ch, 1 pi, 2 ka, 3 pr
869882
double eff = 1.;
870-
if (cfg.mEfficiency)
871-
eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt()));
883+
if (!cfgUsePIDEfficiencies) {
884+
if (cfg.mEfficiency)
885+
eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt()));
886+
} else {
887+
eff = cfg.mPIDEfficiencies[pidIndex]->GetBinContent(cfg.mPIDEfficiencies[pidIndex]->FindBin(track.pt()));
888+
}
872889
if (eff == 0)
873890
return -1.;
874891
else
@@ -1231,7 +1248,7 @@ struct FlowGenericFramework {
12311248

12321249
float total = chtotal;
12331250
if (cfgUsePIDTotal)
1234-
(pidcounter) ? acceptedtracks.pidtotal[pidcounter] : chtotal;
1251+
total = (pidcounter) ? acceptedtracks.pidtotal[pidcounter] : chtotal;
12351252

12361253
if (total == 0.) {
12371254
++pidcounter;
@@ -1336,11 +1353,11 @@ struct FlowGenericFramework {
13361353
int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1;
13371354
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
13381355
if (v0.mK0Short() < cfgPIDCuts.cfgK0SignalMin)
1339-
npt_resonances[0][ptBinIndex] += getEfficiency(v0);
1356+
npt_resonances[0][ptBinIndex]++;
13401357
else if (v0.mK0Short() > cfgPIDCuts.cfgK0SignalMax)
1341-
npt_resonances[2][ptBinIndex] += getEfficiency(v0);
1358+
npt_resonances[2][ptBinIndex]++;
13421359
else
1343-
npt_resonances[1][ptBinIndex] += getEfficiency(v0);
1360+
npt_resonances[1][ptBinIndex]++;
13441361
registry.fill(HIST("K0/hK0s"), 1);
13451362
}
13461363
}
@@ -1351,11 +1368,11 @@ struct FlowGenericFramework {
13511368
int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1;
13521369
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
13531370
if (v0.mLambda() < cfgPIDCuts.cfgLambdaSignalMin)
1354-
npt_resonances[3][ptBinIndex] += getEfficiency(v0);
1371+
npt_resonances[3][ptBinIndex]++;
13551372
else if (v0.mLambda() > cfgPIDCuts.cfgLambdaSignalMax)
1356-
npt_resonances[5][ptBinIndex] += getEfficiency(v0);
1373+
npt_resonances[5][ptBinIndex]++;
13571374
else
1358-
npt_resonances[4][ptBinIndex] += getEfficiency(v0);
1375+
npt_resonances[4][ptBinIndex]++;
13591376
registry.fill(HIST("Lambda/hLambdas"), 1);
13601377
}
13611378
}
@@ -1368,7 +1385,7 @@ struct FlowGenericFramework {
13681385
for (const auto& vec : fractions_resonances) {
13691386
float total = chtotal;
13701387
if (cfgUsePIDTotal)
1371-
(pidcounter) ? std::accumulate(vec.begin(), vec.end(), 0) : chtotal;
1388+
total = (pidcounter) ? std::accumulate(vec.begin(), vec.end(), 0) : chtotal;
13721389

13731390
if (total == 0.) {
13741391
++pidcounter;
@@ -1441,8 +1458,12 @@ struct FlowGenericFramework {
14411458
if (!nchSelected(track))
14421459
return;
14431460

1444-
acceptedTracks.total += getEfficiency(track);
1445-
++acceptedTracks.total_uncorr;
1461+
double weffCh = getEfficiency(track, 0);
1462+
if (std::abs(track.eta()) < cfgEtaNch) {
1463+
if (weffCh > 0)
1464+
acceptedTracks.total += weffCh;
1465+
++acceptedTracks.total_uncorr;
1466+
}
14461467

14471468
if (!trackSelected(track, field))
14481469
return;
@@ -1457,18 +1478,23 @@ struct FlowGenericFramework {
14571478
pidIndex = 3;
14581479
}
14591480

1460-
if (pidIndex)
1461-
acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track);
1462-
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
1463-
1464-
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1465-
acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1466-
if (pidIndex == 1)
1467-
acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1468-
if (pidIndex == 2)
1469-
acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1470-
if (pidIndex == 3)
1471-
acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1481+
if (std::abs(track.eta()) < cfgEtaNch) {
1482+
double weff = getEfficiency(track, pidIndex);
1483+
1484+
if (pidIndex && weff > 0)
1485+
acceptedTracks.pidtotal[pidIndex - 1] += weff;
1486+
1487+
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
1488+
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1489+
if (weffCh > 0)
1490+
acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? weffCh : 1.0;
1491+
if (pidIndex == 1 && weff > 0)
1492+
acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1493+
if (pidIndex == 2 && weff > 0)
1494+
acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1495+
if (pidIndex == 3 && weff > 0)
1496+
acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1497+
}
14721498
}
14731499

14741500
if (cfgFillWeights) {
@@ -1504,21 +1530,24 @@ struct FlowGenericFramework {
15041530
if (std::abs(track.pdgCode()) == kProton)
15051531
pidIndex = 3;
15061532
}
1507-
++acceptedTracks.total;
1508-
++acceptedTracks.total_uncorr;
15091533

1510-
if (pidIndex)
1511-
acceptedTracks.pidtotal[pidIndex - 1] += 1;
1512-
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
1513-
1514-
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1515-
acceptedTracks.nch[ptBinIndex] += 1.0;
1516-
if (pidIndex == 1)
1517-
acceptedTracks.npi[ptBinIndex] += 1.0;
1518-
if (pidIndex == 2)
1519-
acceptedTracks.nka[ptBinIndex] += 1.0;
1520-
if (pidIndex == 3)
1521-
acceptedTracks.npr[ptBinIndex] += 1.0;
1534+
if (std::abs(track.eta()) < cfgEtaNch) {
1535+
++acceptedTracks.total;
1536+
++acceptedTracks.total_uncorr;
1537+
1538+
if (pidIndex)
1539+
acceptedTracks.pidtotal[pidIndex - 1] += 1;
1540+
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
1541+
1542+
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1543+
acceptedTracks.nch[ptBinIndex] += 1.0;
1544+
if (pidIndex == 1)
1545+
acceptedTracks.npi[ptBinIndex] += 1.0;
1546+
if (pidIndex == 2)
1547+
acceptedTracks.nka[ptBinIndex] += 1.0;
1548+
if (pidIndex == 3)
1549+
acceptedTracks.npr[ptBinIndex] += 1.0;
1550+
}
15221551
}
15231552

15241553
fillPtSums<kGen>(track, vtxz, pidIndex);
@@ -1532,32 +1561,39 @@ struct FlowGenericFramework {
15321561
// Select tracks with nominal cuts always
15331562
if (!nchSelected(track))
15341563
return;
1535-
acceptedTracks.total += getEfficiency(track);
1536-
++acceptedTracks.total_uncorr;
1564+
1565+
double weffCh = getEfficiency(track, 0);
1566+
if (std::abs(track.eta()) < cfgEtaNch) {
1567+
if (weffCh > 0)
1568+
acceptedTracks.total += weffCh;
1569+
++acceptedTracks.total_uncorr;
1570+
}
15371571

15381572
if (!trackSelected(track, field))
15391573
return;
15401574
// int pidIndex = 0;
15411575
// if (cfgUsePID) Need PID for v02
15421576
int pidIndex = getNsigmaPID(track);
15431577

1544-
if (pidIndex)
1545-
acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track);
1578+
if (std::abs(track.eta()) < cfgEtaNch) {
1579+
double weff = getEfficiency(track, pidIndex);
1580+
if (pidIndex && weff > 0)
1581+
acceptedTracks.pidtotal[pidIndex - 1] += weff;
15461582

1547-
if (pidIndex)
1548-
acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track);
1549-
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
1583+
int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1;
15501584

1551-
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1552-
acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1553-
if (pidIndex == 1) {
1554-
acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1555-
}
1556-
if (pidIndex == 2) {
1557-
acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1558-
}
1559-
if (pidIndex == 3) {
1560-
acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0;
1585+
if (!(ptBinIndex < 0 || ptBinIndex >= static_cast<int>(o2::analysis::gfw::ptbinning.size()))) {
1586+
if (weffCh > 0)
1587+
acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? weffCh : 1.0;
1588+
if (pidIndex == 1 && weff > 0) {
1589+
acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1590+
}
1591+
if (pidIndex == 2 && weff > 0) {
1592+
acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1593+
}
1594+
if (pidIndex == 3 && weff > 0) {
1595+
acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0;
1596+
}
15611597
}
15621598
}
15631599

@@ -1779,27 +1815,31 @@ struct FlowGenericFramework {
17791815
bool withinPtNch = (track.pt() > ptmins[0] && track.pt() < ptmaxs[0]);
17801816
if (!withinPtPOI && !withinPtRef)
17811817
return;
1818+
double weff = (dt == kGen) ? 1. : getEfficiency(track, pid_index);
1819+
double weffInclusive = (dt == kGen) ? 1. : getEfficiency(track, 0);
1820+
if (weff < 0)
1821+
return;
17821822
double waccRef = (dt == kGen) ? 1. : getAcceptance(track, vtxz, 0);
17831823
double waccPOI = (dt == kGen) ? 1. : withinPtPOI ? getAcceptance(track, vtxz, pid_index + 1)
17841824
: getAcceptance(track, vtxz, 0); //
17851825
if (withinPtRef && withinPtPOI && pid_index)
17861826
waccRef = waccPOI; // if particle is both (then it's overlap), override ref with POI
17871827
if (withinPtRef)
1788-
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccRef, 1);
1828+
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccRef * weffInclusive, 1);
17891829
if (withinPtPOI && pid_index)
1790-
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, (1 << (pid_index + 1)));
1830+
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff, (1 << (pid_index + 1)));
17911831
if (withinPtNch)
1792-
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, 2);
1832+
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff * weffInclusive, 2);
17931833
if (withinPtPOI && withinPtRef && pid_index)
1794-
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, (1 << (pid_index + 5)));
1834+
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff, (1 << (pid_index + 5)));
17951835
if (withinPtNch && withinPtRef)
1796-
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, 32);
1836+
fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff * weffInclusive, 32);
17971837
} else { // Analysing only integrated flow
17981838
bool withinPtRef = (track.pt() > o2::analysis::gfw::ptreflow && track.pt() < o2::analysis::gfw::ptrefup);
17991839
bool withinPtPOI = (track.pt() > o2::analysis::gfw::ptpoilow && track.pt() < o2::analysis::gfw::ptpoiup);
18001840
if (!withinPtPOI && !withinPtRef)
18011841
return;
1802-
double weff = (dt == kGen) ? 1. : getEfficiency(track);
1842+
double weff = (dt == kGen) ? 1. : getEfficiency(track, 0);
18031843
if (weff < 0)
18041844
return;
18051845
if (cfgUseDensityDependentCorrection && withinPtRef && dt != kGen) {
@@ -1886,7 +1926,7 @@ struct FlowGenericFramework {
18861926
}
18871927

18881928
o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ;
1889-
o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz;
1929+
o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz;
18901930

18911931
using GFWCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>>;
18921932
// using GFWTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::pidTOFPi, aod::pidTPCPi, aod::pidTOFKa, aod::pidTPCKa, aod::pidTOFPr, aod::pidTPCPr>>;

0 commit comments

Comments
 (0)