Skip to content

Commit bb6032f

Browse files
ArvindiitiArvindiiti
committed
[PWGJE] Update leading-jet selection corrections to include unmatched jets and normalize per selected reconstructed collision
1 parent 335d658 commit bb6032f

1 file changed

Lines changed: 83 additions & 50 deletions

File tree

‎PWGJE/Tasks/nucleiInJets.cxx‎

Lines changed: 83 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -258,6 +258,8 @@ struct nucleiInJets {
258258
using JetMCPartTable = soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>>;
259259
using JetMCDetTable = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>>;
260260

261+
Preslice<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> detectorJetsPerCollision = aod::jet::collisionId;
262+
261263
SliceCache cache;
262264
HistogramRegistry jetHist{"jetHist", {}, OutputObjHandlingPolicy::AnalysisObject};
263265

@@ -817,15 +819,23 @@ struct nucleiInJets {
817819
jetHist.add<TH2>("eff/recmatched/mcCSpectra/gen/perpCone/pt/PtParticleType", "Pt (gen, mcCSpectra, perp cone) vs particletype", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
818820

819821
jetHist.add<TH1>("jetSelCorr/genSel/hLeadingJetPt", "particle-level leading jet selected for jet-selection correction; #it{p}_{T,lead}^{gen} (GeV/#it{c}); Entries", HistType::kTH1F, {{100, 0., 100.}});
820-
jetHist.add<TH1>("jetSelCorr/recoSel/hLeadingJetPtBkgSub", "matched detector-level leading jet selected for jet-selection correction; #it{p}_{T,lead}^{reco,corr} (GeV/#it{c}); Entries", HistType::kTH1F, {{120, -20., 100.}});
822+
jetHist.add<TH1>("jetSelCorr/recoSel/hLeadingJetPtBkgSub", "detector-level leading jet selected for jet-selection correction; #it{p}_{T,lead}^{reco,corr} (GeV/#it{c}); Entries", HistType::kTH1F, {{120, -20., 100.}});
821823
jetHist.add<TH2>("jetSelCorr/recoSel/hMatchedLeadingJetPt", "selected matched leading jet; #it{p}_{T,lead}^{reco,corr} (GeV/#it{c}); #it{p}_{T,lead}^{gen} (GeV/#it{c})", HistType::kTH2F, {{120, -20., 100.}, {100, 0., 100.}});
822824
jetHist.add<TH1>("jetSelCorr/eventNorm/hEventCounts", "event count for leading-jet event-normalization correction; event selection; Entries", HistType::kTH1D, {{2, 0., 2.}});
823825
jetHist.get<TH1>(HIST("jetSelCorr/eventNorm/hEventCounts"))->GetXaxis()->SetBinLabel(1, "gen lead #it{p}_{T} > cut");
824826
jetHist.get<TH1>(HIST("jetSelCorr/eventNorm/hEventCounts"))->GetXaxis()->SetBinLabel(2, "reco lead #it{p}_{T}^{corr} > cut");
827+
jetHist.add<TH1>("jetSelCorr/recoSel/hMatchStatus", "selected detector leading-jet matching; status; reconstructed collisions", HistType::kTH1D, {{3, 0., 3.}});
828+
jetHist.get<TH1>(HIST("jetSelCorr/recoSel/hMatchStatus"))->GetXaxis()->SetBinLabel(1, "no truth match in MC collision");
829+
jetHist.get<TH1>(HIST("jetSelCorr/recoSel/hMatchStatus"))->GetXaxis()->SetBinLabel(2, "truth axis outside acceptance");
830+
jetHist.get<TH1>(HIST("jetSelCorr/recoSel/hMatchStatus"))->GetXaxis()->SetBinLabel(3, "accepted truth axis");
831+
jetHist.add<TH1>("jetSelCorr/eventNorm/hSelectedRecoCollisions", "selected reconstructed collisions per MC collision; selected reconstructed collisions; MC collisions", HistType::kTH1D, {{11, -0.5, 10.5}});
832+
jetHist.add<TH1>("jetSelCorr/matchedTruthAxis/hEventCount", "selected reconstructed collisions with accepted matched truth axis; selection; reconstructed collisions", HistType::kTH1D, {{1, 0., 1.}});
833+
jetHist.add<TH2>("jetSelCorr/matchedTruthAxis/jetCone/pt/PtParticleType", "generated primaries around the matched truth axis; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
834+
jetHist.add<TH2>("jetSelCorr/matchedTruthAxis/perpCone/pt/PtParticleType", "generated primaries in perpendicular cones of the matched truth jet; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
825835
jetHist.add<TH2>("jetSelCorr/genSel/jetCone/pt/PtParticleType", "generated primaries in particle-level leading-jet cone; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
826836
jetHist.add<TH2>("jetSelCorr/genSel/perpCone/pt/PtParticleType", "generated primaries in particle-level perpendicular cone; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
827-
jetHist.add<TH2>("jetSelCorr/recoSel/jetCone/pt/PtParticleType", "generated primaries in matched reco-selected leading-jet cone; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
828-
jetHist.add<TH2>("jetSelCorr/recoSel/perpCone/pt/PtParticleType", "generated primaries in matched reco-selected perpendicular cone; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
837+
jetHist.add<TH2>("jetSelCorr/recoSel/jetCone/pt/PtParticleType", "generated primaries around the selected detector leading-jet axis; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
838+
jetHist.add<TH2>("jetSelCorr/recoSel/perpCone/pt/PtParticleType", "generated primaries in perpendicular cones of the selected detector leading jet; #it{p}_{T}^{gen} (GeV/#it{c}); particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}});
829839

830840
jetHist.add<TH2>("feeddown/antiProton/jetCone/PtOrigin", "reconstructed #bar{p} origin in jet cone; #it{p}_{T}^{rec} (GeV/#it{c}); origin", HistType::kTH2D, {{PtAxis}, {ParticleOriginAxis}});
831841
jetHist.add<TH2>("feeddown/antiProton/jetCone/PtOriginTPC", "reconstructed #bar{p} origin in jet cone, TPC PID; #it{p}_{T}^{rec} (GeV/#it{c}); origin", HistType::kTH2D, {{PtAxis}, {ParticleOriginAxis}});
@@ -2406,7 +2416,7 @@ struct nucleiInJets {
24062416
return true;
24072417
}
24082418

2409-
template <bool RecoSelected, typename ParticlesType>
2419+
template <bool RecoSelected, bool MatchedTruthAxis = false, typename ParticlesType>
24102420
void fillJetSelectionCorrectionParticles(const ParticlesType& mcParticles, double jetEta, double jetPhi)
24112421
{
24122422
const auto perpConePhiJet = getPerpendicuarPhi(jetPhi);
@@ -2430,7 +2440,9 @@ struct nucleiInJets {
24302440
const double delPhi = TVector2::Phi_mpi_pi(jetPhi - mcParticle.phi());
24312441
const double rJet = RecoDecay::sqrtSumOfSquares(delEta, delPhi);
24322442
if (rJet < cfgjetR) {
2433-
if constexpr (RecoSelected) {
2443+
if constexpr (MatchedTruthAxis) {
2444+
jetHist.fill(HIST("jetSelCorr/matchedTruthAxis/jetCone/pt/PtParticleType"), mcParticle.pt(), particleType);
2445+
} else if constexpr (RecoSelected) {
24342446
jetHist.fill(HIST("jetSelCorr/recoSel/jetCone/pt/PtParticleType"), mcParticle.pt(), particleType);
24352447
} else {
24362448
jetHist.fill(HIST("jetSelCorr/genSel/jetCone/pt/PtParticleType"), mcParticle.pt(), particleType);
@@ -2442,7 +2454,9 @@ struct nucleiInJets {
24422454
const double rPerpCone1 = RecoDecay::sqrtSumOfSquares(delEta, delPhiPerpCone1);
24432455
const double rPerpCone2 = RecoDecay::sqrtSumOfSquares(delEta, delPhiPerpCone2);
24442456
if (rPerpCone1 < cfgjetR || rPerpCone2 < cfgjetR) {
2445-
if constexpr (RecoSelected) {
2457+
if constexpr (MatchedTruthAxis) {
2458+
jetHist.fill(HIST("jetSelCorr/matchedTruthAxis/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType);
2459+
} else if constexpr (RecoSelected) {
24462460
jetHist.fill(HIST("jetSelCorr/recoSel/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType);
24472461
} else {
24482462
jetHist.fill(HIST("jetSelCorr/genSel/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType);
@@ -2818,7 +2832,7 @@ struct nucleiInJets {
28182832
int nprocessSimJEEvents = 0;
28192833
void processJetSelectionCorrection(aod::JetMcCollision const& collision,
28202834
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>> const& recocolls,
2821-
JetMCDetTable const&, JetMCPartTable const& mcpjets, aod::JetParticles const& mcParticles)
2835+
JetMCDetTable const& mcdjets, JetMCPartTable const& mcpjets, aod::JetParticles const& mcParticles)
28222836
{
28232837
if (std::abs(collision.posZ()) > cfgMaxZVertex) {
28242838
return;
@@ -2847,61 +2861,80 @@ struct nucleiInJets {
28472861
fillJetSelectionCorrectionParticles<false>(mcParticles, genLeadingJetEta, genLeadingJetPhi);
28482862
}
28492863

2850-
bool hasRecoSelectedLeadingJet = false;
2851-
double recoSelectedLeadingJetPt = -999.;
2852-
double recoSelectedMatchedGenJetPt = -999.;
2853-
double recoSelectedMatchedGenJetEta = -999.;
2854-
double recoSelectedMatchedGenJetPhi = -999.;
2855-
2856-
for (const auto& mcpjet : mcpjets) {
2857-
if (!mcpjet.has_matchedJetGeo()) {
2864+
// Count each reconstructed collision, as in data. The generated selection
2865+
// above is independent: neither branch requires the other to pass the cut.
2866+
int nSelectedRecoCollisions = 0;
2867+
const double jetArea = M_PI * cfgjetR * cfgjetR;
2868+
for (const auto& recocoll : recocolls) {
2869+
if (!isRecoCollisionSelectedForJetSelectionCorrection(recocoll)) {
28582870
continue;
28592871
}
2860-
if (!isConeAxisAccepted(mcpjet.eta())) {
2872+
if (usebkgSubractionMC && !(recocoll.rho() > 0.)) {
28612873
continue;
28622874
}
2863-
for (const auto& mcdjet : mcpjet.template matchedJetGeo_as<JetMCDetTable>()) {
2864-
if (!isConeAxisAccepted(mcdjet.eta())) {
2865-
continue;
2866-
}
2867-
double selectedRecoCollisionRho = -1.;
2868-
bool hasSelectedRecoCollision = false;
2869-
for (const auto& recocoll : recocolls) {
2870-
if (mcdjet.collisionId() != recocoll.globalIndex()) {
2871-
continue;
2872-
}
2873-
if (!isRecoCollisionSelectedForJetSelectionCorrection(recocoll)) {
2874-
continue;
2875-
}
2876-
selectedRecoCollisionRho = recocoll.rho();
2877-
hasSelectedRecoCollision = true;
2878-
break;
2879-
}
2880-
if (!hasSelectedRecoCollision) {
2881-
continue;
2882-
}
28832875

2884-
const double jetArea = M_PI * cfgjetR * cfgjetR;
2885-
const double mcdJetPtBkgSub = usebkgSubractionMC ? mcdjet.pt() - selectedRecoCollisionRho * jetArea : mcdjet.pt();
2886-
if (mcdJetPtBkgSub <= cfgjetPtBkgSubMinMC) {
2876+
const auto jetsInCollision = mcdjets.sliceBy(detectorJetsPerCollision, recocoll.globalIndex());
2877+
auto leadingJet = jetsInCollision.begin();
2878+
bool hasLeadingJet = false;
2879+
double leadingJetPt = cfgjetPtBkgSubMinMC;
2880+
for (auto jet = jetsInCollision.begin(); jet != jetsInCollision.end(); ++jet) {
2881+
if (!isConeAxisAccepted(jet.eta())) {
28872882
continue;
28882883
}
2889-
if (mcdJetPtBkgSub > recoSelectedLeadingJetPt) {
2890-
hasRecoSelectedLeadingJet = true;
2891-
recoSelectedLeadingJetPt = mcdJetPtBkgSub;
2892-
recoSelectedMatchedGenJetPt = mcpjet.pt();
2893-
recoSelectedMatchedGenJetEta = mcpjet.eta();
2894-
recoSelectedMatchedGenJetPhi = mcpjet.phi();
2884+
const double correctedPt = usebkgSubractionMC ? jet.pt() - recocoll.rho() * jetArea : jet.pt();
2885+
if (correctedPt > leadingJetPt) {
2886+
leadingJet = jet;
2887+
leadingJetPt = correctedPt;
2888+
hasLeadingJet = true;
28952889
}
28962890
}
2897-
}
2891+
if (!hasLeadingJet) {
2892+
continue;
2893+
}
28982894

2899-
if (hasRecoSelectedLeadingJet) {
2895+
++nSelectedRecoCollisions;
29002896
jetHist.fill(HIST("jetSelCorr/eventNorm/hEventCounts"), 1.5);
2901-
jetHist.fill(HIST("jetSelCorr/recoSel/hLeadingJetPtBkgSub"), recoSelectedLeadingJetPt);
2902-
jetHist.fill(HIST("jetSelCorr/recoSel/hMatchedLeadingJetPt"), recoSelectedLeadingJetPt, recoSelectedMatchedGenJetPt);
2903-
fillJetSelectionCorrectionParticles<true>(mcParticles, recoSelectedMatchedGenJetEta, recoSelectedMatchedGenJetPhi);
2897+
jetHist.fill(HIST("jetSelCorr/recoSel/hLeadingJetPtBkgSub"), leadingJetPt);
2898+
// All selected detector jets, including unmatched jets, use detector axes
2899+
// Never substitute a matched subleading jet or mix truth and detector axes
2900+
fillJetSelectionCorrectionParticles<true>(mcParticles, leadingJet.eta(), leadingJet.phi());
2901+
2902+
// Preserve a separate matched-truth-axis sample for efficiency comparisons
2903+
// Resolve multiple geometrical matches by the smallest angular distance
2904+
bool hasTruthMatch = false;
2905+
double closestDistance = 1.e9;
2906+
double truthPt = 0.;
2907+
double truthEta = 0.;
2908+
double truthPhi = 0.;
2909+
if (leadingJet.has_matchedJetGeo()) {
2910+
for (const auto& truthJet : leadingJet.template matchedJetGeo_as<JetMCPartTable>()) {
2911+
if (truthJet.mcCollisionId() != collision.globalIndex()) {
2912+
continue;
2913+
}
2914+
const double distance = RecoDecay::sqrtSumOfSquares(leadingJet.eta() - truthJet.eta(), TVector2::Phi_mpi_pi(leadingJet.phi() - truthJet.phi()));
2915+
if (distance < closestDistance) {
2916+
closestDistance = distance;
2917+
hasTruthMatch = true;
2918+
truthPt = truthJet.pt();
2919+
truthEta = truthJet.eta();
2920+
truthPhi = truthJet.phi();
2921+
}
2922+
}
2923+
}
2924+
if (!hasTruthMatch) {
2925+
jetHist.fill(HIST("jetSelCorr/recoSel/hMatchStatus"), 0.5);
2926+
continue;
2927+
}
2928+
jetHist.fill(HIST("jetSelCorr/recoSel/hMatchedLeadingJetPt"), leadingJetPt, truthPt);
2929+
if (!isConeAxisAccepted(truthEta)) {
2930+
jetHist.fill(HIST("jetSelCorr/recoSel/hMatchStatus"), 1.5);
2931+
continue;
2932+
}
2933+
jetHist.fill(HIST("jetSelCorr/recoSel/hMatchStatus"), 2.5);
2934+
jetHist.fill(HIST("jetSelCorr/matchedTruthAxis/hEventCount"), 0.5);
2935+
fillJetSelectionCorrectionParticles<true, true>(mcParticles, truthEta, truthPhi);
29042936
}
2937+
jetHist.fill(HIST("jetSelCorr/eventNorm/hSelectedRecoCollisions"), nSelectedRecoCollisions);
29052938
}
29062939

29072940
void processGenMatched(aod::JetMcCollision const& collision,

0 commit comments

Comments
 (0)