Skip to content

Commit f0078b6

Browse files
joachimckhjoachimckh
authored
[PWGJE] tracking efficiency (#15520)
1 parent 58d02bd commit f0078b6

File tree

1 file changed

+205
-10
lines changed

1 file changed

+205
-10
lines changed

‎PWGJE/Tasks/jetSpectraEseTask.cxx‎

Lines changed: 205 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@
3232
#include <Framework/HistogramSpec.h>
3333
#include <Framework/InitContext.h>
3434
#include <Framework/Logger.h>
35+
#include <Framework/O2DatabasePDGPlugin.h>
3536
#include <Framework/OutputObjHeader.h>
3637
#include <Framework/runDataProcessing.h>
3738

@@ -53,6 +54,7 @@ using namespace o2::framework;
5354
using namespace o2::framework::expressions;
5455

5556
struct JetSpectraEseTask {
57+
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
5658
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
5759
Configurable<float> jetR{"jetR", 0.2, "jet resolution parameter"};
5860
Configurable<float> randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"};
@@ -133,17 +135,26 @@ struct JetSpectraEseTask {
133135
static constexpr float Acceptance = 0.9f;
134136
static constexpr float LowFT0Cut = 1e-8;
135137

138+
Service<ccdb::BasicCCDBManager> ccdb;
139+
struct Efficiency {
140+
TH1D* hEff = nullptr;
141+
bool isLoaded = false;
142+
} cfg;
143+
136144
Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax);
137145
Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * Scaler) && nabs(aod::jet::eta) < Acceptance - jetR;
138146
Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut;
139147
Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;
140148
using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets, aod::ChargedMCDetectorLevelJetEventWeights>>;
141149
Preslice<ChargedMCDJets> mcdjetsPerJCollision = o2::aod::jet::collisionId;
142150
Preslice<aod::JetTracks> tracksPerJCollision = o2::aod::jtrack::collisionId;
151+
Preslice<aod::JetTracksMCD> mcdTracksPerJCollision = o2::aod::jtrack::collisionId;
152+
Preslice<aod::JetParticles> particlesPerJMcCollision = o2::aod::jmcparticle::mcCollisionId;
143153

144154
SliceCache cache;
145155
using BinningType = ColumnBinningPolicy<aod::jcollision::PosZ, aod::jcollision::CentFT0C>;
146156
BinningType corrBinning{{binsZVtx, binsCentrality}, true};
157+
Service<o2::framework::O2DatabasePDG> pdg;
147158

148159
enum class DetID { FT0C,
149160
FT0A,
@@ -389,6 +400,47 @@ struct JetSpectraEseTask {
389400
registry.add("hOccupancy", "Occupancy;Occupancy;entries", {HistType::kTH1F, {{occAxis}}});
390401
registry.add("hPsiOccupancy", "Occupancy;#Psi_{2};entries", {HistType::kTH3F, {{centAxis}, {150, -2.5, 2.5}, {occAxis}}});
391402
}
403+
if (doprocessMCGenTrack) {
404+
LOGF(info, "JetSpectraEseTask::init() - MCGen track");
405+
registry.add("hTrackPtGen", "", {HistType::kTH1F, {{assocTrackPt}}});
406+
registry.add("hTrackEtaGen", "", {HistType::kTH1F, {{etaAxis}}});
407+
registry.add("hTrackPhiGen", "", {HistType::kTH1F, {{phiAxis}}});
408+
}
409+
if (doprocessMCRecoTrack) {
410+
LOGF(info, "JetSpectraEseTask::init() - MCRec track");
411+
registry.add("hTrackPtReco", "", {HistType::kTH1F, {{assocTrackPt}}});
412+
registry.add("hTrackEtaReco", "", {HistType::kTH1F, {{etaAxis}}});
413+
registry.add("hTrackPhiReco", "", {HistType::kTH1F, {{phiAxis}}});
414+
}
415+
}
416+
417+
void loadEfficiency(aod::BCsWithTimestamps::iterator const& bc)
418+
{
419+
uint64_t timestamp = bc.timestamp();
420+
if (cfg.isLoaded) {
421+
return;
422+
}
423+
if (!cfgEfficiency.value.empty()) {
424+
cfg.hEff = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
425+
if (cfg.hEff == nullptr) {
426+
LOGF(fatal, "Could not load track efficiency from %s", cfgEfficiency.value.c_str());
427+
}
428+
LOGF(info, "Loaded tracking efficiency from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.hEff);
429+
}
430+
cfg.isLoaded = true;
431+
return;
432+
}
433+
434+
template <typename TTrack>
435+
double getEfficiency(TTrack track)
436+
{
437+
double eff{1.0};
438+
if (cfg.hEff)
439+
eff = cfg.hEff->GetBinContent(cfg.hEff->FindBin(track.pt()));
440+
if (eff == 0)
441+
return -1.;
442+
else
443+
return 1. / eff;
392444
}
393445

394446
template <typename TCollision, typename TJets>
@@ -463,19 +515,25 @@ struct JetSpectraEseTask {
463515
for (const auto& track : tracks) {
464516
if (!jetderiveddatautilities::selectTrack(track, trackSelection))
465517
continue;
518+
double weff = getEfficiency(track);
519+
if (weff < 0)
520+
continue;
466521
auto deta = track.eta() - jet.eta();
467522
auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf);
468-
registry.fill(HIST("thn_jethad_corr_same"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]);
469-
hSameSub[lRndInd]->Fill(centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]);
523+
registry.fill(HIST("thn_jethad_corr_same"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
524+
hSameSub[lRndInd]->Fill(centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
470525
}
471526
}
472527
for (const auto& track : tracks) {
473-
registry.fill(HIST("trackQA/before/hTrackPt"), centrality, track.pt());
528+
double weff = getEfficiency(track);
529+
if (weff < 0)
530+
continue;
531+
registry.fill(HIST("trackQA/before/hTrackPt"), centrality, track.pt(), weff);
474532
registry.fill(HIST("trackQA/before/hTrackEta"), centrality, track.eta());
475533
registry.fill(HIST("trackQA/before/hTrackPhi"), centrality, track.phi());
476534
if (!jetderiveddatautilities::selectTrack(track, trackSelection))
477535
continue;
478-
registry.fill(HIST("trackQA/after/hTrackPt"), centrality, track.pt());
536+
registry.fill(HIST("trackQA/after/hTrackPt"), centrality, track.pt(), weff);
479537
registry.fill(HIST("trackQA/after/hTrackEta"), centrality, track.eta());
480538
registry.fill(HIST("trackQA/after/hTrackPhi"), centrality, track.phi());
481539
registry.fill(HIST("h3CenttrPhiPsi2"), centrality, RecoDecay::constrainAngle(track.phi() - psi.psi2, -o2::constants::math::PI), qPerc[0]);
@@ -489,7 +547,10 @@ struct JetSpectraEseTask {
489547
{
490548
auto tracksTuple = std::make_tuple(jets, tracks);
491549
Pair<TCollisions, TJets, TTracks, BinningType> pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache};
550+
492551
for (const auto& [c1, jets1, c2, tracks2] : pairData) {
552+
auto bc = c2.template bc_as<aod::BCsWithTimestamps>();
553+
loadEfficiency(bc);
493554
auto c1Tracks = tracks.sliceBy(tracksPerJCollision, c1.globalIndex());
494555
registry.fill(HIST("eventQA/before/hVtxZMixed"), c1.posZ());
495556
registry.fill(HIST("eventQA/before/hVtxZMixed2"), c2.posZ());
@@ -548,23 +609,29 @@ struct JetSpectraEseTask {
548609
auto vCorrL = cfgrhoPhi ? corrL(jet) : vCorr;
549610
float dPhi{RecoDecay::constrainAngle(jet.phi() - psi.psi2, -o2::constants::math::PI)};
550611

612+
registry.fill(HIST("hNtrigMixed"), centrality, vCorrL, dPhi, qPerc[0]);
551613
for (const auto& track : tracks2) {
552614
if (!jetderiveddatautilities::selectTrack(track, trackSelection))
553615
continue;
616+
double weff = getEfficiency(track);
617+
if (weff < 0)
618+
continue;
554619
auto deta = track.eta() - jet.eta();
555620
auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf);
556-
registry.fill(HIST("hNtrigMixed"), centrality, vCorrL, dPhi, qPerc[0]);
557-
registry.fill(HIST("thn_jethad_corr_mixed"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]);
621+
registry.fill(HIST("thn_jethad_corr_mixed"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
558622
}
559623
}
560624

561625
for (const auto& track : tracks2) {
562-
registry.fill(HIST("trackQA/before/hTrackPtMixed"), centrality, track.pt());
626+
double weff = getEfficiency(track);
627+
if (weff < 0)
628+
continue;
629+
registry.fill(HIST("trackQA/before/hTrackPtMixed"), centrality, track.pt(), weff);
563630
registry.fill(HIST("trackQA/before/hTrackEtaMixed"), centrality, track.eta());
564631
registry.fill(HIST("trackQA/before/hTrackPhiMixed"), centrality, track.phi());
565632
if (!jetderiveddatautilities::selectTrack(track, trackSelection))
566633
continue;
567-
registry.fill(HIST("trackQA/after/hTrackPtMixed"), centrality, track.pt());
634+
registry.fill(HIST("trackQA/after/hTrackPtMixed"), centrality, track.pt(), weff);
568635
registry.fill(HIST("trackQA/after/hTrackEtaMixed"), centrality, track.eta());
569636
registry.fill(HIST("trackQA/after/hTrackPhiMixed"), centrality, track.phi());
570637
}
@@ -573,7 +640,7 @@ struct JetSpectraEseTask {
573640

574641
void processESEDataCharged(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
575642
soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& jets,
576-
aod::JetTracks const& tracks)
643+
aod::JetTracks const& tracks, aod::BCsWithTimestamps const&)
577644
{
578645
registry.fill(HIST("eventQA/hEventCounter"), kFilteredInputEv);
579646
registry.fill(HIST("eventQA/before/hVtxZ"), collision.posZ());
@@ -585,13 +652,15 @@ struct JetSpectraEseTask {
585652
return;
586653
registry.fill(HIST("eventQA/hEventCounter"), kOccupancyCut);
587654

655+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
656+
loadEfficiency(bc);
588657
jetSpectra(collision, jets, tracks);
589658
}
590659
PROCESS_SWITCH(JetSpectraEseTask, processESEDataCharged, "process ese collisions", true);
591660

592661
void processESEDataChargedMixed(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs> const& collisions,
593662
soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& jets,
594-
aod::JetTracks const& tracks)
663+
aod::JetTracks const& tracks, aod::BCsWithTimestamps const&)
595664
{
596665
jetMixed(collisions, jets, tracks);
597666
}
@@ -768,6 +837,132 @@ struct JetSpectraEseTask {
768837
}
769838
PROCESS_SWITCH(JetSpectraEseTask, processMCChargedMatched, "jet MC process: geometrically matched MCP and MCD for response matrix and efficiency", false);
770839

840+
// process for gen/reco tracks for pT efficiency
841+
bool isChargedParticle(int pdgCode)
842+
{
843+
auto pdgParticle = pdg->GetParticle(pdgCode);
844+
return pdgParticle && pdgParticle->Charge() != 0.0;
845+
}
846+
847+
void processMCGenTrack(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
848+
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
849+
aod::JetTracksMCD const&,
850+
aod::JetParticles const& particles)
851+
{
852+
if (mcCollision.size() < 1) {
853+
return;
854+
}
855+
if (collisions.size() != 1) {
856+
return;
857+
}
858+
if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
859+
return;
860+
}
861+
862+
auto collision = collisions.begin();
863+
if (!(std::abs(collision.posZ()) < vertexZCut)) {
864+
return;
865+
}
866+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
867+
return;
868+
}
869+
if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) {
870+
return;
871+
}
872+
873+
auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f;
874+
if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) {
875+
return;
876+
}
877+
878+
auto particlesInCollision = particles.sliceBy(particlesPerJMcCollision, mcCollision.globalIndex());
879+
for (const auto& particle : particlesInCollision) {
880+
if (!isChargedParticle(particle.pdgCode())) {
881+
continue;
882+
}
883+
if (!particle.isPhysicalPrimary()) {
884+
continue;
885+
}
886+
if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) {
887+
continue;
888+
}
889+
if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) {
890+
continue;
891+
}
892+
893+
registry.fill(HIST("hTrackPtGen"), particle.pt());
894+
registry.fill(HIST("hTrackEtaGen"), particle.eta());
895+
registry.fill(HIST("hTrackPhiGen"), particle.phi());
896+
}
897+
}
898+
PROCESS_SWITCH(JetSpectraEseTask, processMCGenTrack, "jet MC process: Generated track", false);
899+
void processMCRecoTrack(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
900+
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
901+
aod::JetTracksMCD const& tracks,
902+
aod::JetParticles const&)
903+
{
904+
if (mcCollision.size() < 1) {
905+
return;
906+
}
907+
if (collisions.size() < 1) {
908+
return;
909+
}
910+
if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
911+
return;
912+
}
913+
914+
std::vector<int> seenMcParticles;
915+
for (const auto& collision : collisions) {
916+
if (!(std::abs(collision.posZ()) < vertexZCut)) {
917+
continue;
918+
}
919+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
920+
continue;
921+
}
922+
if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) {
923+
continue;
924+
}
925+
926+
auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f;
927+
if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) {
928+
continue;
929+
}
930+
931+
auto tracksInCollision = tracks.sliceBy(mcdTracksPerJCollision, collision.globalIndex());
932+
for (const auto& track : tracksInCollision) {
933+
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
934+
continue;
935+
}
936+
if (!track.has_mcParticle()) {
937+
continue;
938+
}
939+
940+
auto particle = track.mcParticle_as<aod::JetParticles>();
941+
if (!isChargedParticle(particle.pdgCode())) {
942+
continue;
943+
}
944+
if (!particle.isPhysicalPrimary()) {
945+
continue;
946+
}
947+
if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) {
948+
continue;
949+
}
950+
if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) {
951+
continue;
952+
}
953+
if (std::find(seenMcParticles.begin(), seenMcParticles.end(), particle.globalIndex()) != seenMcParticles.end()) {
954+
continue;
955+
}
956+
seenMcParticles.push_back(particle.globalIndex());
957+
958+
registry.fill(HIST("hTrackPtReco"), track.pt());
959+
registry.fill(HIST("hTrackEtaReco"), track.eta());
960+
registry.fill(HIST("hTrackPhiReco"), track.phi());
961+
}
962+
}
963+
}
964+
PROCESS_SWITCH(JetSpectraEseTask, processMCRecoTrack, "jet MC process: Reconstructed track", false);
965+
771966
static constexpr float InvalidValue = 999.;
772967

773968
// template <bool FillAllPsi, bool FillHist, typename EPCol>

0 commit comments

Comments
 (0)