Skip to content

Commit 5dbe5e9

Browse files
abmodakabmodak
authored
Add track efficiency correction
1 parent bc257c4 commit 5dbe5e9

1 file changed

Lines changed: 167 additions & 8 deletions

File tree

‎PWGCF/TwoParticleCorrelations/Tasks/longrangecorrDerived.cxx‎

Lines changed: 167 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,15 @@
1919
#include "PWGCF/TwoParticleCorrelations/DataModel/LongRangeDerived.h"
2020
#include "PWGUD/Core/SGSelector.h"
2121

22+
#include "Common/CCDB/EventSelectionParams.h"
2223
#include "Common/Core/RecoDecay.h"
24+
#include "Common/Core/TrackSelection.h"
25+
#include "Common/Core/TrackSelectionDefaults.h"
26+
#include "Common/DataModel/EventSelection.h"
27+
#include "Common/DataModel/McCollisionExtra.h"
28+
#include "Common/DataModel/TrackSelectionTables.h"
2329

30+
#include <CCDB/BasicCCDBManager.h>
2431
#include <CommonConstants/MathConstants.h>
2532
#include <Framework/AnalysisDataModel.h>
2633
#include <Framework/AnalysisHelpers.h>
@@ -31,6 +38,7 @@
3138
#include <Framework/HistogramRegistry.h>
3239
#include <Framework/HistogramSpec.h>
3340
#include <Framework/InitContext.h>
41+
#include <Framework/O2DatabasePDGPlugin.h>
3442
#include <Framework/OutputObjHeader.h>
3543
#include <Framework/runDataProcessing.h>
3644

@@ -51,11 +59,16 @@ using namespace o2::aod::fwdtrack;
5159
using namespace o2::aod::evsel;
5260
using namespace o2::constants::math;
5361

62+
auto static constexpr KminCharge = 3.0f;
63+
5464
struct LongrangecorrDerived {
5565

5666
SliceCache cache;
5767
SGSelector sgSelector;
5868
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
69+
TrackSelection myTrackFilter;
70+
Service<o2::framework::O2DatabasePDG> pdg;
71+
Service<o2::ccdb::BasicCCDBManager> ccdb;
5972

6073
struct : ConfigurableGroup {
6174
Configurable<int> cfgNmixedevent{"cfgNmixedevent", 5, "how many events are mixed"};
@@ -66,10 +79,16 @@ struct LongrangecorrDerived {
6679
Configurable<bool> isUseCentEst{"isUseCentEst", false, "Centrality based classification"};
6780
Configurable<int> isUseDataLikeMult{"isUseDataLikeMult", 0, "Data like mult/cent classification"};
6881

82+
Configurable<float> cfgEtaCut{"cfgEtaCut", 0.8f, "Eta range to consider"};
83+
Configurable<float> cfgPtCutMin{"cfgPtCutMin", 0.2f, "minimum accepted track pT"};
84+
Configurable<float> cfgPtCutMax{"cfgPtCutMax", 10.0f, "maximum accepted track pT"};
6985
Configurable<float> cfgTpcMinNclsFound{"cfgTpcMinNclsFound", 50.0f, ""};
7086
Configurable<float> cfgTpcMinNCrossedRows{"cfgTpcMinNCrossedRows", 70.0f, ""};
7187
Configurable<float> cfgTpcMaxChi2PerCluster{"cfgTpcMaxChi2PerCluster", 4.0f, ""};
7288
Configurable<float> cfgTpcMaxDcaZ{"cfgTpcMaxDcaZ", 1.0f, ""};
89+
Configurable<float> cfgLowEffCut{"cfgLowEffCut", 0.001f, "Low efficiency cut"};
90+
Configurable<bool> applyEffCorr{"applyEffCorr", true, "Enable efficiency correction"};
91+
Configurable<std::string> cfgEffccdbPath{"cfgEffccdbPath", "/alice/data/CCDB/Users/a/abmodak/OO/Efficiency", "Browse track eff object from CCDB"};
7392

7493
Configurable<int> cfgMftCluster{"cfgMftCluster", 5, "cut on MFT Cluster"};
7594
Configurable<float> cfgMftDcaxy{"cfgMftDcaxy", 2.0f, "cut on DCA xy for MFT tracks"};
@@ -121,6 +140,10 @@ struct LongrangecorrDerived {
121140
OutputObj<CorrelationContainer> same{"sameEvent"};
122141
OutputObj<CorrelationContainer> mixed{"mixedEvent"};
123142

143+
// corrections
144+
TH3D* hTrkEff = nullptr;
145+
bool fLoadTrkEffCorr = false;
146+
124147
using CollsTable = aod::LRCollisions;
125148
using TrksTable = aod::LRMidTracks;
126149
using MftTrksTable = aod::LRMftTracks;
@@ -211,6 +234,52 @@ struct LongrangecorrDerived {
211234

212235
histos.add("ReassignedMFTtrackAmbDegree", "ReassignedMFTtrackAmbDegree", kTH1D, {cfgAxis.axisMFTAmbDegree});
213236
histos.add("AssignedMFTtrackAmbDegree", "AssignedMFTtrackAmbDegree", kTH1D, {cfgAxis.axisMFTAmbDegree});
237+
238+
if (doprocessTPCtrackEff) {
239+
histos.add("hGenMCdndpt", "hGenMCdndpt", kTH3D, {cfgAxis.axisVtxZ, cfgAxis.axisEtaEfficiency, cfgAxis.axisPtEfficiency});
240+
histos.add("hRecMCdndpt", "hRecMCdndpt", kTH3D, {cfgAxis.axisVtxZ, cfgAxis.axisEtaEfficiency, cfgAxis.axisPtEfficiency});
241+
}
242+
243+
myTrackFilter = getGlobalTrackSelectionRun3ITSMatch(TrackSelection::GlobalTrackRun3ITSMatching::Run3ITSibAny,
244+
TrackSelection::GlobalTrackRun3DCAxyCut::Default);
245+
myTrackFilter.SetPtRange(cfgSel.cfgPtCutMin, cfgSel.cfgPtCutMax);
246+
myTrackFilter.SetEtaRange(-cfgSel.cfgEtaCut, cfgSel.cfgEtaCut);
247+
myTrackFilter.SetMinNCrossedRowsTPC(cfgSel.cfgTpcMinNCrossedRows);
248+
myTrackFilter.SetMinNClustersTPC(cfgSel.cfgTpcMinNclsFound);
249+
myTrackFilter.SetMaxChi2PerClusterTPC(cfgSel.cfgTpcMaxChi2PerCluster);
250+
myTrackFilter.SetMaxDcaZ(cfgSel.cfgTpcMaxDcaZ);
251+
myTrackFilter.print();
252+
}
253+
254+
void loadEffCorrection(uint64_t timestamp)
255+
{
256+
if (fLoadTrkEffCorr) {
257+
return;
258+
}
259+
if (cfgSel.cfgEffccdbPath.value.empty() == false) {
260+
hTrkEff = ccdb->getForTimeStamp<TH3D>(cfgSel.cfgEffccdbPath, timestamp);
261+
if (hTrkEff == nullptr) {
262+
LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgSel.cfgEffccdbPath.value.c_str());
263+
}
264+
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgSel.cfgEffccdbPath.value.c_str(), (void*)hTrkEff);
265+
}
266+
fLoadTrkEffCorr = true;
267+
}
268+
269+
float getTrkEffCorr(float posZ, float eta, float pt)
270+
{
271+
float eff = 1.0f;
272+
if (hTrkEff) {
273+
int zBin = hTrkEff->GetXaxis()->FindBin(posZ);
274+
int etaBin = hTrkEff->GetYaxis()->FindBin(eta);
275+
int ptBin = hTrkEff->GetZaxis()->FindBin(pt);
276+
eff = hTrkEff->GetBinContent(zBin, etaBin, ptBin);
277+
} else {
278+
eff = 1.0f;
279+
}
280+
if (eff < cfgSel.cfgLowEffCut)
281+
eff = 1.0;
282+
return eff;
214283
}
215284

216285
template <typename TTrack>
@@ -336,6 +405,7 @@ struct LongrangecorrDerived {
336405
int fSampleIndex = gRandom->Uniform(0, cfgSel.cfgSampleSize);
337406
for (auto const& triggerTrack : triggers) {
338407
auto trigAmpl = 1.0f;
408+
auto trkeff = 1.0f;
339409
if constexpr (requires { triggerTrack.channelID(); }) {
340410
trigAmpl = triggerTrack.amplitude();
341411
} else {
@@ -352,14 +422,21 @@ struct LongrangecorrDerived {
352422
if (cfgSel.cfgV0Mask != 0 && (cfgSel.cfgV0Mask & (1u << static_cast<uint32_t>(triggerTrack.v0Type()))) == 0u)
353423
continue;
354424
}
425+
426+
if constexpr (requires { triggerTrack.trackType(); }) {
427+
trkeff = getTrkEffCorr(vz, triggerTrack.eta(), triggerTrack.pt());
428+
} else {
429+
trkeff = 1.0;
430+
}
431+
355432
if (!mixing) {
356433
fillTrigTrackQA(triggerTrack);
357434
if constexpr (requires { triggerTrack.channelID(); }) {
358-
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, 1.0, 1.0, eventWeight * trigAmpl);
435+
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, 1.0, 1.0, eventWeight * trigAmpl * trkeff);
359436
} else if constexpr (requires { triggerTrack.invMass(); }) {
360-
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, triggerTrack.pt(), triggerTrack.invMass(), eventWeight * trigAmpl);
437+
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, triggerTrack.pt(), triggerTrack.invMass(), eventWeight * trigAmpl * trkeff);
361438
} else {
362-
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, triggerTrack.pt(), 1.0, eventWeight * trigAmpl);
439+
histos.fill(HIST("Trig_hist"), fSampleIndex, vz, multiplicity, triggerTrack.pt(), 1.0, eventWeight * trigAmpl * trkeff);
363440
}
364441
}
365442
for (auto const& assoTrack : assocs) {
@@ -377,16 +454,16 @@ struct LongrangecorrDerived {
377454
float deltaEta = triggerTrack.eta() - assoTrack.eta();
378455
if (!mixing) {
379456
fillAssocTrackQA(assoTrack);
380-
histos.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta, eventWeight * trigAmpl * assoAmpl);
457+
histos.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta, eventWeight * trigAmpl * assoAmpl * trkeff);
381458
} else {
382-
histos.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta, eventWeight * trigAmpl * assoAmpl);
459+
histos.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta, eventWeight * trigAmpl * assoAmpl * trkeff);
383460
}
384461
if constexpr (requires { triggerTrack.channelID(); }) {
385-
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, 1.0, deltaPhi, deltaEta, 1.0, eventWeight * trigAmpl * assoAmpl);
462+
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, 1.0, deltaPhi, deltaEta, 1.0, eventWeight * trigAmpl * assoAmpl * trkeff);
386463
} else if constexpr (requires { triggerTrack.invMass(); }) {
387-
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, triggerTrack.pt(), deltaPhi, deltaEta, triggerTrack.invMass(), eventWeight * trigAmpl * assoAmpl);
464+
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, triggerTrack.pt(), deltaPhi, deltaEta, triggerTrack.invMass(), eventWeight * trigAmpl * assoAmpl * trkeff);
388465
} else {
389-
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, triggerTrack.pt(), deltaPhi, deltaEta, 1.0, eventWeight * trigAmpl * assoAmpl);
466+
target->getPairHist()->Fill(step, fSampleIndex, vz, multiplicity, triggerTrack.pt(), deltaPhi, deltaEta, 1.0, eventWeight * trigAmpl * assoAmpl * trkeff);
390467
}
391468
} // associated tracks
392469
} // trigger tracks
@@ -782,6 +859,87 @@ struct LongrangecorrDerived {
782859
processMcGenMixed(mccollisions, ft0as, ft0cs);
783860
}
784861

862+
using ColMCTrueTable = soa::Join<aod::McCollisions, aod::McCollsExtra>;
863+
using ColMCRecTable = soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>>;
864+
using TrksMCRecTable = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::McTrackLabels, aod::TrackSelection>>;
865+
Preslice<TrksMCRecTable> perColMidtrack = aod::track::collisionId;
866+
867+
template <typename CheckCol>
868+
bool isEventSelected(CheckCol const& col)
869+
{
870+
if (!col.sel8())
871+
return false;
872+
if (!col.selection_bit(o2::aod::evsel::kNoSameBunchPileup))
873+
return false;
874+
if (!col.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV))
875+
return false;
876+
if (!col.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll))
877+
return false;
878+
if (std::abs(col.posZ()) >= cfgSel.cfgVtxCut)
879+
return false;
880+
return true;
881+
}
882+
883+
template <typename CheckGenPart>
884+
bool isGenPartSelected(CheckGenPart const& particle)
885+
{
886+
if (!particle.isPhysicalPrimary() || !particle.producedByGenerator())
887+
return false;
888+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
889+
if (pdgParticle == nullptr)
890+
return false;
891+
if (std::abs(pdgParticle->Charge()) < KminCharge)
892+
return false;
893+
return true;
894+
}
895+
896+
void processTPCtrackEff(ColMCTrueTable::iterator const& mcCollision, ColMCRecTable const& RecCols,
897+
TrksMCRecTable const& RecTracks, aod::McParticles const& mcparticles)
898+
{
899+
if (std::abs(mcCollision.posZ()) >= cfgSel.cfgVtxCut) {
900+
return;
901+
}
902+
bool atLeastOne = false;
903+
for (const auto& RecCol : RecCols) {
904+
if (!isEventSelected(RecCol))
905+
continue;
906+
if (RecCol.globalIndex() != mcCollision.bestCollisionIndex())
907+
continue;
908+
auto recTracksPart = RecTracks.sliceBy(perColMidtrack, RecCol.globalIndex());
909+
atLeastOne = true;
910+
}
911+
for (const auto& particle : mcparticles) {
912+
if (!isGenPartSelected(particle) ||
913+
std::abs(particle.eta()) > cfgSel.cfgEtaCut ||
914+
particle.pt() < cfgSel.cfgPtCutMin ||
915+
particle.pt() > cfgSel.cfgPtCutMax)
916+
continue;
917+
if (atLeastOne)
918+
histos.fill(HIST("hGenMCdndpt"), mcCollision.posZ(), particle.eta(), particle.pt());
919+
}
920+
for (const auto& RecCol : RecCols) {
921+
if (!isEventSelected(RecCol))
922+
continue;
923+
if (RecCol.globalIndex() != mcCollision.bestCollisionIndex())
924+
continue;
925+
auto recTracksPart = RecTracks.sliceBy(perColMidtrack, RecCol.globalIndex());
926+
for (const auto& track : recTracksPart) {
927+
if (!track.isGlobalTrack())
928+
continue;
929+
if (!myTrackFilter.IsSelected(track))
930+
continue;
931+
if (!track.has_mcParticle())
932+
continue;
933+
auto particle = track.mcParticle();
934+
if (RecCol.mcCollisionId() != particle.mcCollisionId())
935+
continue;
936+
if (particle.isPhysicalPrimary()) {
937+
histos.fill(HIST("hRecMCdndpt"), mcCollision.posZ(), particle.eta(), particle.pt());
938+
}
939+
}
940+
}
941+
}
942+
785943
PROCESS_SWITCH(LongrangecorrDerived, processTpcft0aSE, "same event TPC vs FT0A", false);
786944
PROCESS_SWITCH(LongrangecorrDerived, processTpcft0aME, "mixed event TPC vs FT0A", false);
787945
PROCESS_SWITCH(LongrangecorrDerived, processTpcft0cSE, "same event TPC vs FT0C", false);
@@ -828,6 +986,7 @@ struct LongrangecorrDerived {
828986
PROCESS_SWITCH(LongrangecorrDerived, processMcGenMftft0aME, "mixed MC gen event MFT vs FT0A", false);
829987
PROCESS_SWITCH(LongrangecorrDerived, processMcGenFt0aft0cSE, "same MC gen event FT0A vs FT0C", false);
830988
PROCESS_SWITCH(LongrangecorrDerived, processMcGenFt0aft0cME, "mixed MC gen event FT0A vs FT0C", false);
989+
PROCESS_SWITCH(LongrangecorrDerived, processTPCtrackEff, "process TPC track efficiency", false);
831990
};
832991

833992
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)