diff --git a/PWGJE/Tasks/jetSpectraEseTask.cxx b/PWGJE/Tasks/jetSpectraEseTask.cxx index 54dbc47901a..faff2f096fc 100644 --- a/PWGJE/Tasks/jetSpectraEseTask.cxx +++ b/PWGJE/Tasks/jetSpectraEseTask.cxx @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -53,6 +54,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; struct JetSpectraEseTask { + Configurable cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"}; Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; Configurable jetR{"jetR", 0.2, "jet resolution parameter"}; Configurable randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"}; @@ -133,6 +135,12 @@ struct JetSpectraEseTask { static constexpr float Acceptance = 0.9f; static constexpr float LowFT0Cut = 1e-8; + Service ccdb; + struct Efficiency { + TH1D* hEff = nullptr; + bool isLoaded = false; + } cfg; + Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax); Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * Scaler) && nabs(aod::jet::eta) < Acceptance - jetR; Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut; @@ -140,10 +148,13 @@ struct JetSpectraEseTask { using ChargedMCDJets = soa::Filtered>; Preslice mcdjetsPerJCollision = o2::aod::jet::collisionId; Preslice tracksPerJCollision = o2::aod::jtrack::collisionId; + Preslice mcdTracksPerJCollision = o2::aod::jtrack::collisionId; + Preslice particlesPerJMcCollision = o2::aod::jmcparticle::mcCollisionId; SliceCache cache; using BinningType = ColumnBinningPolicy; BinningType corrBinning{{binsZVtx, binsCentrality}, true}; + Service pdg; enum class DetID { FT0C, FT0A, @@ -389,6 +400,47 @@ struct JetSpectraEseTask { registry.add("hOccupancy", "Occupancy;Occupancy;entries", {HistType::kTH1F, {{occAxis}}}); registry.add("hPsiOccupancy", "Occupancy;#Psi_{2};entries", {HistType::kTH3F, {{centAxis}, {150, -2.5, 2.5}, {occAxis}}}); } + if (doprocessMCGenTrack) { + LOGF(info, "JetSpectraEseTask::init() - MCGen track"); + registry.add("hTrackPtGen", "", {HistType::kTH1F, {{assocTrackPt}}}); + registry.add("hTrackEtaGen", "", {HistType::kTH1F, {{etaAxis}}}); + registry.add("hTrackPhiGen", "", {HistType::kTH1F, {{phiAxis}}}); + } + if (doprocessMCRecoTrack) { + LOGF(info, "JetSpectraEseTask::init() - MCRec track"); + registry.add("hTrackPtReco", "", {HistType::kTH1F, {{assocTrackPt}}}); + registry.add("hTrackEtaReco", "", {HistType::kTH1F, {{etaAxis}}}); + registry.add("hTrackPhiReco", "", {HistType::kTH1F, {{phiAxis}}}); + } + } + + void loadEfficiency(aod::BCsWithTimestamps::iterator const& bc) + { + uint64_t timestamp = bc.timestamp(); + if (cfg.isLoaded) { + return; + } + if (!cfgEfficiency.value.empty()) { + cfg.hEff = ccdb->getForTimeStamp(cfgEfficiency, timestamp); + if (cfg.hEff == nullptr) { + LOGF(fatal, "Could not load track efficiency from %s", cfgEfficiency.value.c_str()); + } + LOGF(info, "Loaded tracking efficiency from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.hEff); + } + cfg.isLoaded = true; + return; + } + + template + double getEfficiency(TTrack track) + { + double eff{1.0}; + if (cfg.hEff) + eff = cfg.hEff->GetBinContent(cfg.hEff->FindBin(track.pt())); + if (eff == 0) + return -1.; + else + return 1. / eff; } template @@ -463,19 +515,25 @@ struct JetSpectraEseTask { for (const auto& track : tracks) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) continue; + double weff = getEfficiency(track); + if (weff < 0) + continue; auto deta = track.eta() - jet.eta(); auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf); - registry.fill(HIST("thn_jethad_corr_same"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]); - hSameSub[lRndInd]->Fill(centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]); + registry.fill(HIST("thn_jethad_corr_same"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff); + hSameSub[lRndInd]->Fill(centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff); } } for (const auto& track : tracks) { - registry.fill(HIST("trackQA/before/hTrackPt"), centrality, track.pt()); + double weff = getEfficiency(track); + if (weff < 0) + continue; + registry.fill(HIST("trackQA/before/hTrackPt"), centrality, track.pt(), weff); registry.fill(HIST("trackQA/before/hTrackEta"), centrality, track.eta()); registry.fill(HIST("trackQA/before/hTrackPhi"), centrality, track.phi()); if (!jetderiveddatautilities::selectTrack(track, trackSelection)) continue; - registry.fill(HIST("trackQA/after/hTrackPt"), centrality, track.pt()); + registry.fill(HIST("trackQA/after/hTrackPt"), centrality, track.pt(), weff); registry.fill(HIST("trackQA/after/hTrackEta"), centrality, track.eta()); registry.fill(HIST("trackQA/after/hTrackPhi"), centrality, track.phi()); registry.fill(HIST("h3CenttrPhiPsi2"), centrality, RecoDecay::constrainAngle(track.phi() - psi.psi2, -o2::constants::math::PI), qPerc[0]); @@ -489,7 +547,10 @@ struct JetSpectraEseTask { { auto tracksTuple = std::make_tuple(jets, tracks); Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + for (const auto& [c1, jets1, c2, tracks2] : pairData) { + auto bc = c2.template bc_as(); + loadEfficiency(bc); auto c1Tracks = tracks.sliceBy(tracksPerJCollision, c1.globalIndex()); registry.fill(HIST("eventQA/before/hVtxZMixed"), c1.posZ()); registry.fill(HIST("eventQA/before/hVtxZMixed2"), c2.posZ()); @@ -548,23 +609,29 @@ struct JetSpectraEseTask { auto vCorrL = cfgrhoPhi ? corrL(jet) : vCorr; float dPhi{RecoDecay::constrainAngle(jet.phi() - psi.psi2, -o2::constants::math::PI)}; + registry.fill(HIST("hNtrigMixed"), centrality, vCorrL, dPhi, qPerc[0]); for (const auto& track : tracks2) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) continue; + double weff = getEfficiency(track); + if (weff < 0) + continue; auto deta = track.eta() - jet.eta(); auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf); - registry.fill(HIST("hNtrigMixed"), centrality, vCorrL, dPhi, qPerc[0]); - registry.fill(HIST("thn_jethad_corr_mixed"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0]); + registry.fill(HIST("thn_jethad_corr_mixed"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff); } } for (const auto& track : tracks2) { - registry.fill(HIST("trackQA/before/hTrackPtMixed"), centrality, track.pt()); + double weff = getEfficiency(track); + if (weff < 0) + continue; + registry.fill(HIST("trackQA/before/hTrackPtMixed"), centrality, track.pt(), weff); registry.fill(HIST("trackQA/before/hTrackEtaMixed"), centrality, track.eta()); registry.fill(HIST("trackQA/before/hTrackPhiMixed"), centrality, track.phi()); if (!jetderiveddatautilities::selectTrack(track, trackSelection)) continue; - registry.fill(HIST("trackQA/after/hTrackPtMixed"), centrality, track.pt()); + registry.fill(HIST("trackQA/after/hTrackPtMixed"), centrality, track.pt(), weff); registry.fill(HIST("trackQA/after/hTrackEtaMixed"), centrality, track.eta()); registry.fill(HIST("trackQA/after/hTrackPhiMixed"), centrality, track.phi()); } @@ -573,7 +640,7 @@ struct JetSpectraEseTask { void processESEDataCharged(soa::Join::iterator const& collision, soa::Filtered> const& jets, - aod::JetTracks const& tracks) + aod::JetTracks const& tracks, aod::BCsWithTimestamps const&) { registry.fill(HIST("eventQA/hEventCounter"), kFilteredInputEv); registry.fill(HIST("eventQA/before/hVtxZ"), collision.posZ()); @@ -585,13 +652,15 @@ struct JetSpectraEseTask { return; registry.fill(HIST("eventQA/hEventCounter"), kOccupancyCut); + auto bc = collision.bc_as(); + loadEfficiency(bc); jetSpectra(collision, jets, tracks); } PROCESS_SWITCH(JetSpectraEseTask, processESEDataCharged, "process ese collisions", true); void processESEDataChargedMixed(soa::Join const& collisions, soa::Filtered> const& jets, - aod::JetTracks const& tracks) + aod::JetTracks const& tracks, aod::BCsWithTimestamps const&) { jetMixed(collisions, jets, tracks); } @@ -768,6 +837,132 @@ struct JetSpectraEseTask { } PROCESS_SWITCH(JetSpectraEseTask, processMCChargedMatched, "jet MC process: geometrically matched MCP and MCD for response matrix and efficiency", false); + // process for gen/reco tracks for pT efficiency + bool isChargedParticle(int pdgCode) + { + auto pdgParticle = pdg->GetParticle(pdgCode); + return pdgParticle && pdgParticle->Charge() != 0.0; + } + + void processMCGenTrack(soa::Filtered>::iterator const& mcCollision, + soa::SmallGroups const& collisions, + aod::JetTracksMCD const&, + aod::JetParticles const& particles) + { + if (mcCollision.size() < 1) { + return; + } + if (collisions.size() != 1) { + return; + } + if (!(std::abs(mcCollision.posZ()) < vertexZCut)) { + return; + } + + auto collision = collisions.begin(); + if (!(std::abs(collision.posZ()) < vertexZCut)) { + return; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } + if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) { + return; + } + + auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f; + if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) { + return; + } + + auto particlesInCollision = particles.sliceBy(particlesPerJMcCollision, mcCollision.globalIndex()); + for (const auto& particle : particlesInCollision) { + if (!isChargedParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) { + continue; + } + if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) { + continue; + } + + registry.fill(HIST("hTrackPtGen"), particle.pt()); + registry.fill(HIST("hTrackEtaGen"), particle.eta()); + registry.fill(HIST("hTrackPhiGen"), particle.phi()); + } + } + PROCESS_SWITCH(JetSpectraEseTask, processMCGenTrack, "jet MC process: Generated track", false); + void processMCRecoTrack(soa::Filtered>::iterator const& mcCollision, + soa::SmallGroups const& collisions, + aod::JetTracksMCD const& tracks, + aod::JetParticles const&) + { + if (mcCollision.size() < 1) { + return; + } + if (collisions.size() < 1) { + return; + } + if (!(std::abs(mcCollision.posZ()) < vertexZCut)) { + return; + } + + std::vector seenMcParticles; + for (const auto& collision : collisions) { + if (!(std::abs(collision.posZ()) < vertexZCut)) { + continue; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + continue; + } + if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) { + continue; + } + + auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f; + if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) { + continue; + } + + auto tracksInCollision = tracks.sliceBy(mcdTracksPerJCollision, collision.globalIndex()); + for (const auto& track : tracksInCollision) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + if (!track.has_mcParticle()) { + continue; + } + + auto particle = track.mcParticle_as(); + if (!isChargedParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) { + continue; + } + if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) { + continue; + } + if (std::find(seenMcParticles.begin(), seenMcParticles.end(), particle.globalIndex()) != seenMcParticles.end()) { + continue; + } + seenMcParticles.push_back(particle.globalIndex()); + + registry.fill(HIST("hTrackPtReco"), track.pt()); + registry.fill(HIST("hTrackEtaReco"), track.eta()); + registry.fill(HIST("hTrackPhiReco"), track.phi()); + } + } + } + PROCESS_SWITCH(JetSpectraEseTask, processMCRecoTrack, "jet MC process: Reconstructed track", false); + static constexpr float InvalidValue = 999.; // template