Skip to content

Commit 296df02

Browse files
Luzhiyonggalibuild
Luzhiyongg
and
alibuild
authored
[PWGCF] produce and apply Efficiency correction (#11469)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent ef6dcb2 commit 296df02

File tree

1 file changed

+125
-32
lines changed

1 file changed

+125
-32
lines changed

‎PWGCF/TwoParticleCorrelations/Tasks/diHadronCor.cxx‎

Lines changed: 125 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include "Common/DataModel/Centrality.h"
3737
#include "PWGCF/DataModel/CorrelationsDerived.h"
3838
#include "Common/DataModel/CollisionAssociationTables.h"
39+
#include "Common/DataModel/PIDResponse.h"
3940
#include "PWGCF/Core/CorrelationContainer.h"
4041
#include "PWGCF/Core/PairCuts.h"
4142
#include "PWGCF/GenericFramework/Core/GFWPowerArray.h"
@@ -44,6 +45,7 @@
4445
#include "PWGCF/GenericFramework/Core/GFWWeights.h"
4546
#include "DataFormatsParameters/GRPObject.h"
4647
#include "DataFormatsParameters/GRPMagField.h"
48+
#include <TPDGCode.h>
4749

4850
using namespace o2;
4951
using namespace o2::framework;
@@ -100,34 +102,46 @@ struct DiHadronCor {
100102
O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 2000, "High cut on TPC occupancy")
101103
O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy")
102104
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
105+
O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object")
106+
O2_DEFINE_CONFIGURABLE(cfgVerbosity, bool, false, "Verbose output")
103107

104108
SliceCache cache;
105109
SliceCache cacheNch;
106110

107111
ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"};
108112
ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity axis for histograms"};
109113
ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "centrality axis for histograms"};
110-
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt axis for histograms"};
114+
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 4, 5, 6, 8, 10}, "pt axis for histograms"};
111115
ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"};
112116
ConfigurableAxis axisDeltaEta{"axisDeltaEta", {48, -2.4, 2.4}, "delta eta axis for histograms"};
113-
ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"};
114-
ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt associated axis for histograms"};
117+
ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 4, 5, 6, 8, 10}, "pt trigger axis for histograms"};
118+
ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 4, 5, 6, 8, 10}, "pt associated axis for histograms"};
115119
ConfigurableAxis axisVtxMix{"axisVtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"};
116120
ConfigurableAxis axisMultMix{"axisMultMix", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for mixed event histograms"};
117121
ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"};
118122

119123
ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"};
120124
ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"};
121-
ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, "pt axis for efficiency histograms"};
125+
ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 4, 5, 6, 8, 10}, "pt axis for efficiency histograms"};
122126

123127
// make the filters and cuts.
124128
Filter collisionFilter = (nabs(aod::collision::posZ) < cfgCutVtxZ) && (aod::evsel::sel8) == true && (aod::cent::centFT0C > cfgCentFT0CMin) && (aod::cent::centFT0C < cfgCentFT0CMax);
125129
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
130+
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSel, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>; // aod::CentFT0Cs
131+
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::McTrackLabels>>;
132+
133+
// Filter for MCParticle
134+
Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax);
135+
using MyMcParticles = soa::Filtered<aod::McParticles>;
136+
137+
// Filter for MCcollisions
138+
Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVtxZ;
139+
using MyMcCollisions = soa::Filtered<aod::McCollisions>;
140+
141+
Preslice<aod::Tracks> perCollision = aod::track::collisionId;
126142

127143
// Corrections
128-
TH1D* mEfficiency = nullptr;
129-
GFWWeights* mAcceptance = nullptr;
130-
TObjArray* mAcceptanceList = nullptr;
144+
TH3D* mEfficiency = nullptr;
131145
bool correctionsLoaded = false;
132146

133147
// Define the outputs
@@ -150,6 +164,9 @@ struct DiHadronCor {
150164
MixedEvent = 3
151165
};
152166

167+
// persistent caches
168+
std::vector<float> efficiencyAssociatedCache;
169+
153170
// Additional Event selection cuts - Copy from flowGenericFramework.cxx
154171
TF1* fMultPVCutLow = nullptr;
155172
TF1* fMultPVCutHigh = nullptr;
@@ -159,9 +176,6 @@ struct DiHadronCor {
159176
TF1* fT0AV0AMean = nullptr;
160177
TF1* fT0AV0ASigma = nullptr;
161178

162-
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSel, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>; // aod::CentFT0Cs
163-
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA>>;
164-
165179
void init(InitContext&)
166180
{
167181
const AxisSpec axisPhi{72, 0.0, constants::math::TwoPI, "#varphi"};
@@ -224,6 +238,7 @@ struct DiHadronCor {
224238
registry.add("deltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}});
225239
registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}});
226240
registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}});
241+
registry.add("EtaCorrected", "EtaCorrected", {HistType::kTH1D, {axisEta}});
227242
registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}});
228243
registry.add("pTCorrected", "pTCorrected", {HistType::kTH1D, {axisPtTrigger}});
229244
registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}});
@@ -237,16 +252,20 @@ struct DiHadronCor {
237252

238253
registry.add("eventcount", "bin", {HistType::kTH1F, {{4, 0, 4, "bin"}}}); // histogram to see how many events are in the same and mixed event
239254

255+
if (!cfgEfficiency.value.empty()) {
256+
efficiencyAssociatedCache.reserve(512);
257+
}
258+
240259
std::vector<AxisSpec> corrAxis = {{axisSample, "Sample"},
241260
{axisVertex, "z-vtx (cm)"},
242261
{axisPtTrigger, "p_{T} (GeV/c)"},
243262
{axisPtAssoc, "p_{T} (GeV/c)"},
244263
{axisDeltaPhi, "#Delta#varphi (rad)"},
245264
{axisDeltaEta, "#Delta#eta"}};
246265
std::vector<AxisSpec> effAxis = {
247-
{axisVertexEfficiency, "z-vtx (cm)"},
248-
{axisPtEfficiency, "p_{T} (GeV/c)"},
249266
{axisEtaEfficiency, "#eta"},
267+
{axisPtEfficiency, "p_{T} (GeV/c)"},
268+
{axisVertexEfficiency, "z-vtx (cm)"},
250269
};
251270
std::vector<AxisSpec> userAxis;
252271

@@ -298,12 +317,18 @@ struct DiHadronCor {
298317
return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu));
299318
}
300319

301-
void loadCorrections(uint64_t timestamp)
320+
void loadEfficiency(uint64_t timestamp)
302321
{
303-
if (correctionsLoaded)
322+
if (correctionsLoaded) {
304323
return;
324+
}
305325
if (cfgEfficiency.value.empty() == false) {
306-
mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
326+
if (cfgLocalEfficiency > 0) {
327+
TFile* fEfficiencyTrigger = TFile::Open(cfgEfficiency.value.c_str(), "READ");
328+
mEfficiency = reinterpret_cast<TH3D*>(fEfficiencyTrigger->Get("ccdb_object"));
329+
} else {
330+
mEfficiency = ccdb->getForTimeStamp<TH3D>(cfgEfficiency, timestamp);
331+
}
307332
if (mEfficiency == nullptr) {
308333
LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiency.value.c_str());
309334
}
@@ -312,13 +337,17 @@ struct DiHadronCor {
312337
correctionsLoaded = true;
313338
}
314339

315-
bool setCurrentParticleWeights(float& weight_nue, float pt)
340+
bool getEfficiencyCorrection(float& weight_nue, float eta, float pt, float posZ)
316341
{
317342
float eff = 1.;
318-
if (mEfficiency)
319-
eff = mEfficiency->GetBinContent(mEfficiency->FindBin(pt));
320-
else
343+
if (mEfficiency) {
344+
int etaBin = mEfficiency->GetXaxis()->FindBin(eta);
345+
int ptBin = mEfficiency->GetYaxis()->FindBin(pt);
346+
int zBin = mEfficiency->GetZaxis()->FindBin(posZ);
347+
eff = mEfficiency->GetBinContent(etaBin, ptBin, zBin);
348+
} else {
321349
eff = 1.0;
350+
}
322351
if (eff == 0)
323352
return false;
324353
weight_nue = 1. / eff;
@@ -336,13 +365,15 @@ struct DiHadronCor {
336365
registry.fill(HIST("zVtx"), collision.posZ());
337366

338367
float weff1 = 1;
368+
float vtxz = collision.posZ();
339369
for (auto const& track1 : tracks) {
340370
if (!trackSelected(track1))
341371
continue;
342-
if (!setCurrentParticleWeights(weff1, track1.pt()))
372+
if (!getEfficiencyCorrection(weff1, track1.eta(), track1.pt(), vtxz))
343373
continue;
344374
registry.fill(HIST("Phi"), RecoDecay::constrainAngle(track1.phi(), 0.0));
345375
registry.fill(HIST("Eta"), track1.eta());
376+
registry.fill(HIST("EtaCorrected"), track1.eta(), weff1);
346377
registry.fill(HIST("pT"), track1.pt());
347378
registry.fill(HIST("pTCorrected"), track1.pt(), weff1);
348379
}
@@ -376,6 +407,16 @@ struct DiHadronCor {
376407
template <CorrelationContainer::CFStep step, typename TTracks, typename TTracksAssoc>
377408
void fillCorrelations(TTracks tracks1, TTracksAssoc tracks2, float posZ, int system, int magneticField, float cent) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms
378409
{
410+
// Cache efficiency for particles (too many FindBin lookups)
411+
if (mEfficiency) {
412+
efficiencyAssociatedCache.clear();
413+
efficiencyAssociatedCache.reserve(tracks2.size());
414+
for (const auto& track2 : tracks2) {
415+
float weff = 1.;
416+
getEfficiencyCorrection(weff, track2.eta(), track2.pt(), posZ);
417+
efficiencyAssociatedCache.push_back(weff);
418+
}
419+
}
379420

380421
if (system == SameEvent) {
381422
registry.fill(HIST("Centrality_used"), cent);
@@ -384,25 +425,26 @@ struct DiHadronCor {
384425

385426
int fSampleIndex = gRandom->Uniform(0, cfgSampleSize);
386427

387-
float weff1 = 1;
388-
float weff2 = 1;
428+
float triggerWeight = 1.0f;
429+
float associatedWeight = 1.0f;
389430
// loop over all tracks
390431
for (auto const& track1 : tracks1) {
391432

392433
if (!trackSelected(track1))
393434
continue;
394-
if (!setCurrentParticleWeights(weff1, track1.pt()))
435+
if (!getEfficiencyCorrection(triggerWeight, track1.eta(), track1.pt(), posZ))
395436
continue;
396437
if (system == SameEvent) {
397-
registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt());
438+
registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), triggerWeight);
398439
}
399440

400441
for (auto const& track2 : tracks2) {
401442

402443
if (!trackSelected(track2))
403444
continue;
404-
if (!setCurrentParticleWeights(weff2, track2.pt()))
405-
continue;
445+
if (mEfficiency) {
446+
associatedWeight = efficiencyAssociatedCache[track2.filteredIndex()];
447+
}
406448

407449
if (track1.pt() <= track2.pt())
408450
continue; // skip if the trigger pt is less than the associate pt
@@ -435,12 +477,12 @@ struct DiHadronCor {
435477
// fill the right sparse and histograms
436478
if (system == SameEvent) {
437479

438-
same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
439-
registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta);
480+
same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta, triggerWeight * associatedWeight);
481+
registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta, triggerWeight * associatedWeight);
440482
} else if (system == MixedEvent) {
441483

442-
mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
443-
registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta);
484+
mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta, triggerWeight * associatedWeight);
485+
registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta, triggerWeight * associatedWeight);
444486
}
445487
}
446488
}
@@ -528,7 +570,7 @@ struct DiHadronCor {
528570

529571
registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin
530572

531-
loadCorrections(bc.timestamp());
573+
loadEfficiency(bc.timestamp());
532574
fillYield(collision, tracks);
533575

534576
if (cfgSelCollByNch && (tracks.size() < cfgCutMultMin || tracks.size() >= cfgCutMultMax)) {
@@ -561,7 +603,7 @@ struct DiHadronCor {
561603
for (auto const& [collision1, tracks1, collision2, tracks2] : pair) {
562604
registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin
563605
auto bc = collision1.bc_as<aod::BCsWithTimestamps>();
564-
loadCorrections(bc.timestamp());
606+
loadEfficiency(bc.timestamp());
565607

566608
if (cfgSelCollByNch && (tracks1.size() < cfgCutMultMin || tracks1.size() >= cfgCutMultMax))
567609
continue;
@@ -587,6 +629,57 @@ struct DiHadronCor {
587629
}
588630

589631
PROCESS_SWITCH(DiHadronCor, processMixed, "Process mixed events", true);
632+
633+
int getSpecies(int pdgCode)
634+
{
635+
switch (pdgCode) {
636+
case 211: // pion
637+
case -211:
638+
return 0;
639+
case 321: // Kaon
640+
case -321:
641+
return 1;
642+
case 2212: // proton
643+
case -2212:
644+
return 2;
645+
default: // NOTE. The efficiency histogram is hardcoded to contain 4 species. Anything special will have the last slot.
646+
return 3;
647+
}
648+
}
649+
650+
void processMCEfficiency(MyMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions>> const& collisions, MyMcParticles const& mcParticles, AodTracks const& tracks)
651+
{
652+
if (cfgSelCollByNch && (tracks.size() < cfgCutMultMin || tracks.size() >= cfgCutMultMax)) {
653+
return;
654+
}
655+
// Primaries
656+
for (const auto& mcParticle : mcParticles) {
657+
if (mcParticle.isPhysicalPrimary()) {
658+
same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), 0., mcCollision.posZ());
659+
}
660+
}
661+
for (const auto& collision : collisions) {
662+
auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex());
663+
if (cfgVerbosity) {
664+
LOGF(info, " Reconstructed collision at vtx-z = %f", collision.posZ());
665+
LOGF(info, " which has %d tracks", groupedTracks.size());
666+
}
667+
668+
for (const auto& track : groupedTracks) {
669+
if (track.has_mcParticle()) {
670+
const auto& mcParticle = track.mcParticle();
671+
if (mcParticle.isPhysicalPrimary()) {
672+
same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), 0., mcCollision.posZ());
673+
}
674+
same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoAll, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), 0., mcCollision.posZ());
675+
} else {
676+
// fake track
677+
same->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, track.eta(), track.pt(), 0, 0., mcCollision.posZ());
678+
}
679+
}
680+
}
681+
}
682+
PROCESS_SWITCH(DiHadronCor, processMCEfficiency, "MC: Extract efficiencies", false);
590683
};
591684

592685
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)