From 6fa8ab725eba2438ece0f4feda759f6d49075036 Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Thu, 24 Sep 2026 20:50:21 +0200 Subject: [PATCH 1/6] [PWGJE] Fix leading-jet selection (Rec) and event normalization --- PWGJE/Tasks/nucleiInJets.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 156020d18b3..2175cf10fd9 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -2495,13 +2495,12 @@ struct nucleiInJets { } // Select the accepted detector-level leading jet before considering matches. - int32_t leadingDetJetId = -1; + int64_t leadingDetJetId = -1; float leadingDetJetPt = -1.f; if (isWithLeadingJet) { for (const auto& mcdjet : mcdjets) { - const float mcdJetPtForResponseSel = usebkgSubractionMC ? mcdjet.pt() - backgroundRho * jetArea : mcdjet.pt(); - if (isConeAxisAccepted(mcdjet.eta()) && mcdJetPtForResponseSel > leadingDetJetPt) { - leadingDetJetPt = mcdJetPtForResponseSel; + if (isConeAxisAccepted(mcdjet.eta()) && mcdjet.pt() > leadingDetJetPt) { + leadingDetJetPt = mcdjet.pt(); leadingDetJetId = mcdjet.globalIndex(); } } @@ -2510,7 +2509,7 @@ struct nucleiInJets { } } - std::vector mcdJetIds{}; + std::vector mcdJetIds{}; std::vector mcdJetPt{}; std::vector mcdJetPhi{}; std::vector mcdJetEta{}; From 617c33b43745705ee08cc0c932f9bd28aa07d94a Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Fri, 25 Sep 2026 11:20:13 +0200 Subject: [PATCH 2/6] Use corrected jet pT for leading-jet selection --- PWGJE/Tasks/nucleiInJets.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 2175cf10fd9..89b90eea3c6 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -2499,8 +2499,9 @@ struct nucleiInJets { float leadingDetJetPt = -1.f; if (isWithLeadingJet) { for (const auto& mcdjet : mcdjets) { - if (isConeAxisAccepted(mcdjet.eta()) && mcdjet.pt() > leadingDetJetPt) { - leadingDetJetPt = mcdjet.pt(); + const float mcdJetPtForResponseSel = usebkgSubractionMC ? mcdjet.pt() - backgroundRho * jetArea : mcdjet.pt(); + if (isConeAxisAccepted(mcdjet.eta()) && mcdJetPtForResponseSel > leadingDetJetPt) { + leadingDetJetPt = mcdJetPtForResponseSel; leadingDetJetId = mcdjet.globalIndex(); } } From 6331fb1eeea32b629f3cff4a1d1b9ba14588480d Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Tue, 29 Sep 2026 19:31:47 +0200 Subject: [PATCH 3/6] [PWGJE] Use int32_t for jet-ID --- PWGJE/Tasks/nucleiInJets.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 89b90eea3c6..156020d18b3 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -2495,7 +2495,7 @@ struct nucleiInJets { } // Select the accepted detector-level leading jet before considering matches. - int64_t leadingDetJetId = -1; + int32_t leadingDetJetId = -1; float leadingDetJetPt = -1.f; if (isWithLeadingJet) { for (const auto& mcdjet : mcdjets) { @@ -2510,7 +2510,7 @@ struct nucleiInJets { } } - std::vector mcdJetIds{}; + std::vector mcdJetIds{}; std::vector mcdJetPt{}; std::vector mcdJetPhi{}; std::vector mcdJetEta{}; From e3c798cb733b2ed7171b72fb594b59ed127bfca5 Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Mon, 19 Oct 2026 09:27:02 +0200 Subject: [PATCH 4/6] [PWGJE] Add configurable MC closure split and antiproton DCA templates --- PWGJE/Tasks/nucleiInJets.cxx | 167 ++++++++++++++++++++++++++++++++--- 1 file changed, 155 insertions(+), 12 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 156020d18b3..806199d2a4d 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -177,6 +177,8 @@ struct nucleiInJets { Configurable useTOFVeto{"useTOFVeto", false, "true: use TOF veto, false: no TOF veto"}; Configurable isRequireHitsInITSLayers{"isRequireHitsInITSLayers", true, "true: at least one hit in the its inner layes"}; Configurable useMcC{"useMcC", true, "use mcC"}; + Configurable mcClosureTrainingFraction{"mcClosureTrainingFraction", 0.5f, "fraction of events assigned to MC correction training; the rest are data-like"}; + Configurable mcClosureSplitSeed{"mcClosureSplitSeed", 0, "TRandom3 seed for MC/data-like split (0: ROOT automatic seed)"}; Configurable usebkgSubractionMC{"usebkgSubractionMC", true, "use rho-area background subtraction for detector-level matched MC jets"}; Configurable cfgjetPtBkgSubMinMC{"cfgjetPtBkgSubMinMC", 10.0f, "minimum detector-level matched MC jet pT after optional background subtraction"}; Configurable useFullConeAcceptance{"useFullConeAcceptance", true, "require jet/perpendicular-cone axes to satisfy |eta_axis| < cfgtrkMaxEta - cfgjetR for jet and efficiency histograms"}; @@ -268,12 +270,17 @@ struct nucleiInJets { Service ccdb; Service pdgDB; TRandom3 randUniform; + void init(o2::framework::InitContext&) { if (doprocessJetTracksData && doprocessJetTracksDataLfPid) { LOGP(fatal, "only one process function should be enabled!!!"); } + if (!(mcClosureTrainingFraction >= 0.f && mcClosureTrainingFraction <= 1.f)) { + LOGP(fatal, "mcClosureTrainingFraction must be between 0 and 1, inclusive"); + } + randUniform.SetSeed(mcClosureSplitSeed); jetHist.print(); // const AxisSpec PtAxis = {100, 0, 10.0}; const AxisSpec PtJetAxis = {100, 0, 100.0}; @@ -642,6 +649,16 @@ struct nucleiInJets { jetHist.add(("tracksInc/dcaxy/rec/" + particle + "/tpcPtVsDcaxy3DPIDTOF").c_str(), "pT vs Dcaxy vs Centrality (PID TOF)", HistType::kTH3F, {{PtAxis}, {100, 0, 100}, dcaxyAxis}); } + // DCAxy templates are needed by both the data and reconstructed-MC processes. + const AxisSpec antiProtonTemplateOriginAxis{4, 0.5, 4.5, "1 primary #bar{p}, 2 material #bar{p}, 3 weak-decay #bar{p}, 4 other/fake"}; + for (const auto& region : {"jetCone", "perpCone"}) { + for (const auto& pid : {"TPC", "TPCTOF", "TPCTOFVeto"}) { + const std::string base = std::string("feeddown/antiProton/template/") + region + "/" + pid + "/"; + jetHist.add(base + "data", "data #bar{p} candidates before DCAxy cut; #it{p}_{T}^{rec}; DCA_{xy}", HistType::kTH2D, {{PtAxis}, {dcaxyAxis}}); + jetHist.add(base + "mc", "MC #bar{p} candidates before DCAxy cut; #it{p}_{T}^{rec}; DCA_{xy}; truth origin", HistType::kTH3D, {{PtAxis}, {dcaxyAxis}, {antiProtonTemplateOriginAxis}}); + } + } + if (isMC) { // inc jetHist.add("recInc/eventStat", "Event statistics (inclusive)", HistType::kTH1F, {{6, 0.f, 6.f}}); @@ -818,6 +835,19 @@ struct nucleiInJets { jetHist.add("eff/recmatched/mcC/gen/perpCone/pt/PtParticleType", "Pt (gen, mcC, perp cone) vs particletype", HistType::kTH2D, {{PtAxis}, {14, -7, 7}}); jetHist.add("eff/recmatched/mcCSpectra/gen/perpCone/pt/PtParticleType", "Pt (gen, mcCSpectra, perp cone) vs particletype", HistType::kTH2D, {{PtAxis}, {14, -7, 7}}); + jetHist.add("eff/recmatched/mcClosure/hSelectedRecoEvents", "matched reconstructed events per closure sample; sample; events", HistType::kTH1D, {{2, 0., 2.}}); + jetHist.get(HIST("eff/recmatched/mcClosure/hSelectedRecoEvents"))->GetXaxis()->SetBinLabel(1, "training"); + jetHist.get(HIST("eff/recmatched/mcClosure/hSelectedRecoEvents"))->GetXaxis()->SetBinLabel(2, "data-like"); + for (const auto& sample : {"training", "dataLike"}) { + for (const auto& region : {"jetCone", "perpCone"}) { + const std::string base = std::string("eff/recmatched/mcClosure/") + sample + "/" + region + "/"; + jetHist.add(base + "hRecAll", "matched primary tracks; reconstructed #it{p}_{T}; true #it{p}_{T}; particle type", HistType::kTH3D, {{PtAxis}, {PtAxis}, {14, -7, 7}}); + jetHist.add(base + "hRecTPCTOFVeto", "matched primary tracks with TPC PID and TOF veto; reconstructed #it{p}_{T}; true #it{p}_{T}; particle type", HistType::kTH3D, {{PtAxis}, {PtAxis}, {14, -7, 7}}); + jetHist.add(base + "hRecTPCTOF", "matched primary tracks with TPC and TOF PID; reconstructed #it{p}_{T}; true #it{p}_{T}; particle type", HistType::kTH3D, {{PtAxis}, {PtAxis}, {14, -7, 7}}); + jetHist.add(base + "hGen", "generated primary particles; true #it{p}_{T}; particle type", HistType::kTH2D, {{PtAxis}, {14, -7, 7}}); + } + } + jetHist.add("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.}}); jetHist.add("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.}}); jetHist.add("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.}}); @@ -953,6 +983,34 @@ struct nucleiInJets { template void fillTrackInfo(const TracksType& trk, const JetType& jets, std::vector& leadingJetPtEtaPhi, float backgroundRho = -1.0f) { + if (isWithLeadingJet && leadingJetPtEtaPhi.size() >= 3 && isConeAxisAccepted(leadingJetPtEtaPhi[1]) && + backgroundRho > 0.f && leadingJetPtEtaPhi[0] - backgroundRho * M_PI * cfgjetR * cfgjetR > cfgjetPtBkgSubMin && + trk.sign() < 0 && trk.pt() <= cMaxPt && isTrackSelectedWithoutDcaxy(trk) && + (!useRapidityCutForPID || isRapiditySelectedForPID(trk.rapidity(o2::constants::physics::MassProton)))) { + const bool tpcPassed = std::abs(trk.tpcNSigmaPr()) < cfgnTPCPIDPr; + if (tpcPassed) { + const bool tofPassed = trk.hasTOF() && std::abs(trk.tofNSigmaPr()) < cfgnTPCPIDPrTOF; + const double deltaEta = leadingJetPtEtaPhi[1] - trk.eta(); + const double deltaPhi = TVector2::Phi_mpi_pi(leadingJetPtEtaPhi[2] - trk.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhi) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPC/data"), trk.pt(), trk.dcaXY()); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); + if (!trk.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); + } + const auto perpPhi = getPerpendicuarPhi(leadingJetPtEtaPhi[2]); + const double deltaPhiPerp1 = TVector2::Phi_mpi_pi(perpPhi[0] - trk.phi()); + const double deltaPhiPerp2 = TVector2::Phi_mpi_pi(perpPhi[1] - trk.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp1) < cfgjetR || RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp2) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPC/data"), trk.pt(), trk.dcaXY()); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); + if (!trk.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); + } + } + } if (!isTrackSelected(trk)) return; if (trk.pt() > cMaxPt) @@ -2483,8 +2541,8 @@ struct nucleiInJets { jetHist.fill(HIST("recmatched/vertexZ"), collision.posZ()); - // Event-wise random splitting for closure test: decide once per event - bool useDataLikeHist = (randUniform.Uniform(0, 1) < 0.5); + // One random choice is shared by the reconstructed and generated spectra of this event. + const bool useTrainingHist = randUniform.Rndm() < mcClosureTrainingFraction; const float backgroundRho = collision.rho(); if (usebkgSubractionMC && backgroundRho <= 0.f) { return; @@ -2580,13 +2638,16 @@ struct nucleiInJets { indexJet = std::distance(mcdJetIds.begin(), leadingMatch); } if (useMcC) { - if (useDataLikeHist) + if (useTrainingHist) jetHist.fill(HIST("recmatched/h2ResponseMatrixLeadingJet"), mcdJetPt.at(indexJet), mcpJetPt.at(indexJet)); else jetHist.fill(HIST("recmatched/mcC/h2ResponseMatrixLeadingJet"), mcdJetPt.at(indexJet), mcpJetPt.at(indexJet)); } else { jetHist.fill(HIST("recmatched/h2ResponseMatrixLeadingJet"), mcdJetPt.at(indexJet), mcpJetPt.at(indexJet)); } + if (useMcC) { + jetHist.fill(HIST("eff/recmatched/mcClosure/hSelectedRecoEvents"), useTrainingHist ? 0.5 : 1.5); + } if (useLeadingJetDetLevelValue) { leadingJetWithPtEtaPhi[0] = mcdJetPt.at(indexJet); leadingJetWithPtEtaPhi[1] = mcdJetEta.at(indexJet); @@ -2597,6 +2658,47 @@ struct nucleiInJets { leadingJetWithPtEtaPhi[2] = mcpJetPhi.at(indexJet); } + // Fill detector-level DCA templates before the DCAxy selection used by the spectra. + const double templateJetPt = mcdJetPt.at(indexJet) - (usebkgSubractionMC ? 0.f : backgroundRho * jetArea); + if (isWithLeadingJet && backgroundRho > 0.f && templateJetPt > cfgjetPtBkgSubMin) { + const double templateJetEta = mcdJetEta.at(indexJet); + const double templateJetPhi = mcdJetPhi.at(indexJet); + const auto perpPhi = getPerpendicuarPhi(templateJetPhi); + for (const auto& track : tracks) { + const auto completeTrack = track.track_as(); + if (completeTrack.sign() >= 0 || completeTrack.pt() > cMaxPt || !isTrackSelectedWithoutDcaxy(completeTrack) || + (useRapidityCutForPID && !isRapiditySelectedForPID(completeTrack.rapidity(o2::constants::physics::MassProton)))) + continue; + if (std::abs(completeTrack.tpcNSigmaPr()) >= cfgnTPCPIDPr) + continue; + int origin = 4; + if (track.has_mcParticle()) { + const auto mcTrack = track.mcParticle_as(); + if (mcTrack.pdgCode() == -2212) + origin = getParticleOriginType(mcTrack); + } + const bool tofPassed = completeTrack.hasTOF() && std::abs(completeTrack.tofNSigmaPr()) < cfgnTPCPIDPrTOF; + const double deltaEta = templateJetEta - completeTrack.eta(); + const double deltaPhi = TVector2::Phi_mpi_pi(templateJetPhi - completeTrack.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhi) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPC/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOF/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + if (!completeTrack.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOFVeto/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + } + const double deltaPhiPerp1 = TVector2::Phi_mpi_pi(perpPhi[0] - completeTrack.phi()); + const double deltaPhiPerp2 = TVector2::Phi_mpi_pi(perpPhi[1] - completeTrack.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp1) < cfgjetR || RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp2) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPC/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOF/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + if (!completeTrack.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOFVeto/mc"), completeTrack.pt(), completeTrack.dcaXY(), origin); + } + } + } + for (const auto& track : tracks) { auto completeTrack = track.track_as(); if (std::fabs(completeTrack.eta()) > cfgtrkMaxEta) @@ -2606,7 +2708,7 @@ struct nucleiInJets { if (!track.has_mcParticle()) continue; auto mcTrack = track.mcParticle_as(); - if (!isRapiditySelectedForPID(mcTrack.y())) + if (useRapidityCutForPID && !isRapiditySelectedForPID(mcTrack.y())) continue; bool isTof(completeTrack.hasTOF()); @@ -2710,7 +2812,7 @@ struct nucleiInJets { const bool isTPCTOFVetoPassed = isTpcPassed && (!isTof || isTofPidPassed); jetHist.fill(HIST("eff/recmatched/pt/PtParticleType"), mcTrack.pt(), jetFlag, particleType); if (useMcC) { - if (useDataLikeHist) + if (useTrainingHist) jetHist.fill(HIST("eff/recmatched/mcC/pt/PtParticleType"), track.pt(), mcTrack.pt(), particleType); else jetHist.fill(HIST("eff/recmatched/mcCSpectra/pt/PtParticleType"), track.pt(), mcTrack.pt(), particleType); @@ -2726,6 +2828,21 @@ struct nucleiInJets { if (jetFlag) { jetHist.fill(HIST("eff/recmatched/jetCone/ptEta/PtEtaParticleType"), mcTrack.pt(), mcTrack.eta(), particleType); + if (useMcC) { + if (useTrainingHist) { + jetHist.fill(HIST("eff/recmatched/mcClosure/training/jetCone/hRecAll"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFVetoPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/jetCone/hRecTPCTOFVeto"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/jetCone/hRecTPCTOF"), track.pt(), mcTrack.pt(), particleType); + } else { + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/jetCone/hRecAll"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFVetoPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/jetCone/hRecTPCTOFVeto"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/jetCone/hRecTPCTOF"), track.pt(), mcTrack.pt(), particleType); + } + } if (isTpcPassed) jetHist.fill(HIST("eff/recmatched/jetCone/ptEta/PtEtaParticleTypeTPC"), mcTrack.pt(), mcTrack.eta(), particleType); if (isTof) @@ -2740,7 +2857,22 @@ struct nucleiInJets { jetHist.fill(HIST("eff/recmatched/perpCone/pt/PtParticleType"), mcTrack.pt(), particleType); jetHist.fill(HIST("eff/recmatched/perpCone/ptEta/PtEtaParticleType"), mcTrack.pt(), mcTrack.eta(), particleType); if (useMcC) { - if (useDataLikeHist) + if (useTrainingHist) { + jetHist.fill(HIST("eff/recmatched/mcClosure/training/perpCone/hRecAll"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFVetoPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/perpCone/hRecTPCTOFVeto"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/perpCone/hRecTPCTOF"), track.pt(), mcTrack.pt(), particleType); + } else { + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/perpCone/hRecAll"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFVetoPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/perpCone/hRecTPCTOFVeto"), track.pt(), mcTrack.pt(), particleType); + if (isTPCTOFPassed) + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/perpCone/hRecTPCTOF"), track.pt(), mcTrack.pt(), particleType); + } + } + if (useMcC) { + if (useTrainingHist) jetHist.fill(HIST("eff/recmatched/perpCone/mcC/pt/PtParticleType"), track.pt(), mcTrack.pt(), particleType); else jetHist.fill(HIST("eff/recmatched/perpCone/mcCSpectra/pt/PtParticleType"), track.pt(), mcTrack.pt(), particleType); @@ -2771,7 +2903,7 @@ struct nucleiInJets { continue; if (std::fabs(mcParticle.eta()) > cfgtrkMaxEta) continue; - if (!isRapiditySelectedForPID(mcParticle.y())) + if (useRapidityCutForPID && !isRapiditySelectedForPID(mcParticle.y())) continue; bool jetFlagMC = false; bool jetFlagPerpConeMC = false; @@ -2806,9 +2938,15 @@ struct nucleiInJets { jetHist.fill(HIST("eff/recmatched/gen/pt/PtParticleType"), mcParticle.pt(), jetFlagMC, particleType); if (jetFlagMC) { jetHist.fill(HIST("eff/recmatched/gen/jetCone/ptEta/PtEtaParticleType"), mcParticle.pt(), mcParticle.eta(), particleType); + if (useMcC) { + if (useTrainingHist) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/jetCone/hGen"), mcParticle.pt(), particleType); + else + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/jetCone/hGen"), mcParticle.pt(), particleType); + } } if (useMcC) { - if (useDataLikeHist) { + if (useTrainingHist) { jetHist.fill(HIST("eff/recmatched/mcC/gen/pt/PtParticleType"), mcParticle.pt(), jetFlagMC, particleType); } else { jetHist.fill(HIST("eff/recmatched/mcCSpectra/gen/pt/PtParticleType"), mcParticle.pt(), jetFlagMC, particleType); @@ -2818,7 +2956,13 @@ struct nucleiInJets { jetHist.fill(HIST("eff/recmatched/gen/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType); jetHist.fill(HIST("eff/recmatched/gen/perpCone/ptEta/PtEtaParticleType"), mcParticle.pt(), mcParticle.eta(), particleType); if (useMcC) { - if (useDataLikeHist) { + if (useTrainingHist) + jetHist.fill(HIST("eff/recmatched/mcClosure/training/perpCone/hGen"), mcParticle.pt(), particleType); + else + jetHist.fill(HIST("eff/recmatched/mcClosure/dataLike/perpCone/hGen"), mcParticle.pt(), particleType); + } + if (useMcC) { + if (useTrainingHist) { jetHist.fill(HIST("eff/recmatched/mcC/gen/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType); } else { jetHist.fill(HIST("eff/recmatched/mcCSpectra/gen/perpCone/pt/PtParticleType"), mcParticle.pt(), particleType); @@ -2957,8 +3101,7 @@ struct nucleiInJets { return; jetHist.fill(HIST("genmatched/vertexZ"), collision.posZ()); - // Event-wise random splitting for closure test: decide once per event - bool useDataLikeHist = (randUniform.Uniform(0, 1) < 0.5); + const bool useTrainingHist = randUniform.Rndm() < mcClosureTrainingFraction; std::vector mcpJetPt{}; std::vector mcpJetPhi{}; @@ -3010,7 +3153,7 @@ struct nucleiInJets { jetHist.fill(HIST("genmatched/hRecMatchedJetEta"), leadingMCPJet.eta(), leadingMCPJet.eta() - mcdjet.eta()); if (useMcC) { - if (useDataLikeHist) { + if (useTrainingHist) { jetHist.fill(HIST("genmatched/hRecMatchedVsGenJetPt"), mcdjet.pt(), leadingMCPJet.pt()); } else { jetHist.fill(HIST("genmatched/mcC/hRecMatchedVsGenJetPt"), mcdjet.pt(), leadingMCPJet.pt()); From 748d79c82974f82ad556cb2ed2f5c9fa289c2a14 Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Tue, 20 Oct 2026 12:57:12 +0200 Subject: [PATCH 5/6] Fix MC closure splitting --- PWGJE/Tasks/nucleiInJets.cxx | 73 +++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 27 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 806199d2a4d..1b690d4bd5c 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -270,6 +270,19 @@ struct nucleiInJets { Service ccdb; Service pdgDB; TRandom3 randUniform; + UInt_t mcClosureSplitSeedEffective = 1; + + bool useTrainingHistForMcCollision(int64_t mcCollisionId) + { + // Use the MC collision ID as part of the seed, so both process functions + // obtain the same decision without storing one entry per collision. + UInt_t collisionSeed = mcClosureSplitSeedEffective ^ static_cast(mcCollisionId); + if (collisionSeed == 0U) { + collisionSeed = 1U; + } + TRandom3 collisionRandom(collisionSeed); + return collisionRandom.Rndm() < mcClosureTrainingFraction; + } void init(o2::framework::InitContext&) { @@ -281,6 +294,10 @@ struct nucleiInJets { LOGP(fatal, "mcClosureTrainingFraction must be between 0 and 1, inclusive"); } randUniform.SetSeed(mcClosureSplitSeed); + mcClosureSplitSeedEffective = randUniform.Integer(0xffffffffU); + if (mcClosureSplitSeedEffective == 0U) { + mcClosureSplitSeedEffective = 1U; + } jetHist.print(); // const AxisSpec PtAxis = {100, 0, 10.0}; const AxisSpec PtJetAxis = {100, 0, 100.0}; @@ -983,31 +1000,33 @@ struct nucleiInJets { template void fillTrackInfo(const TracksType& trk, const JetType& jets, std::vector& leadingJetPtEtaPhi, float backgroundRho = -1.0f) { - if (isWithLeadingJet && leadingJetPtEtaPhi.size() >= 3 && isConeAxisAccepted(leadingJetPtEtaPhi[1]) && - backgroundRho > 0.f && leadingJetPtEtaPhi[0] - backgroundRho * M_PI * cfgjetR * cfgjetR > cfgjetPtBkgSubMin && - trk.sign() < 0 && trk.pt() <= cMaxPt && isTrackSelectedWithoutDcaxy(trk) && - (!useRapidityCutForPID || isRapiditySelectedForPID(trk.rapidity(o2::constants::physics::MassProton)))) { - const bool tpcPassed = std::abs(trk.tpcNSigmaPr()) < cfgnTPCPIDPr; - if (tpcPassed) { - const bool tofPassed = trk.hasTOF() && std::abs(trk.tofNSigmaPr()) < cfgnTPCPIDPrTOF; - const double deltaEta = leadingJetPtEtaPhi[1] - trk.eta(); - const double deltaPhi = TVector2::Phi_mpi_pi(leadingJetPtEtaPhi[2] - trk.phi()); - if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhi) < cfgjetR) { - jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPC/data"), trk.pt(), trk.dcaXY()); - if (tofPassed) - jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); - if (!trk.hasTOF() || tofPassed) - jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); - } - const auto perpPhi = getPerpendicuarPhi(leadingJetPtEtaPhi[2]); - const double deltaPhiPerp1 = TVector2::Phi_mpi_pi(perpPhi[0] - trk.phi()); - const double deltaPhiPerp2 = TVector2::Phi_mpi_pi(perpPhi[1] - trk.phi()); - if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp1) < cfgjetR || RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp2) < cfgjetR) { - jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPC/data"), trk.pt(), trk.dcaXY()); - if (tofPassed) - jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); - if (!trk.hasTOF() || tofPassed) - jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); + if constexpr (!IsMC) { + if (isWithLeadingJet && leadingJetPtEtaPhi.size() >= 3 && isConeAxisAccepted(leadingJetPtEtaPhi[1]) && + backgroundRho > 0.f && leadingJetPtEtaPhi[0] - backgroundRho * M_PI * cfgjetR * cfgjetR > cfgjetPtBkgSubMin && + trk.sign() < 0 && trk.pt() <= cMaxPt && isTrackSelectedWithoutDcaxy(trk) && + (!useRapidityCutForPID || isRapiditySelectedForPID(trk.rapidity(o2::constants::physics::MassProton)))) { + const bool tpcPassed = std::abs(trk.tpcNSigmaPr()) < cfgnTPCPIDPr; + if (tpcPassed) { + const bool tofPassed = trk.hasTOF() && std::abs(trk.tofNSigmaPr()) < cfgnTPCPIDPrTOF; + const double deltaEta = leadingJetPtEtaPhi[1] - trk.eta(); + const double deltaPhi = TVector2::Phi_mpi_pi(leadingJetPtEtaPhi[2] - trk.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhi) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPC/data"), trk.pt(), trk.dcaXY()); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); + if (!trk.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/jetCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); + } + const auto perpPhi = getPerpendicuarPhi(leadingJetPtEtaPhi[2]); + const double deltaPhiPerp1 = TVector2::Phi_mpi_pi(perpPhi[0] - trk.phi()); + const double deltaPhiPerp2 = TVector2::Phi_mpi_pi(perpPhi[1] - trk.phi()); + if (RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp1) < cfgjetR || RecoDecay::sqrtSumOfSquares(deltaEta, deltaPhiPerp2) < cfgjetR) { + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPC/data"), trk.pt(), trk.dcaXY()); + if (tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOF/data"), trk.pt(), trk.dcaXY()); + if (!trk.hasTOF() || tofPassed) + jetHist.fill(HIST("feeddown/antiProton/template/perpCone/TPCTOFVeto/data"), trk.pt(), trk.dcaXY()); + } } } } @@ -2542,7 +2561,7 @@ struct nucleiInJets { jetHist.fill(HIST("recmatched/vertexZ"), collision.posZ()); // One random choice is shared by the reconstructed and generated spectra of this event. - const bool useTrainingHist = randUniform.Rndm() < mcClosureTrainingFraction; + const bool useTrainingHist = useTrainingHistForMcCollision(collision.mcCollision().globalIndex()); const float backgroundRho = collision.rho(); if (usebkgSubractionMC && backgroundRho <= 0.f) { return; @@ -3101,7 +3120,7 @@ struct nucleiInJets { return; jetHist.fill(HIST("genmatched/vertexZ"), collision.posZ()); - const bool useTrainingHist = randUniform.Rndm() < mcClosureTrainingFraction; + const bool useTrainingHist = useTrainingHistForMcCollision(collision.globalIndex()); std::vector mcpJetPt{}; std::vector mcpJetPhi{}; From 18a2281bf30f7db1af00ab762753d4ef0fd92062 Mon Sep 17 00:00:00 2001 From: Arvind Khuntia Date: Tue, 20 Oct 2026 14:05:13 +0200 Subject: [PATCH 6/6] Fix MC closure split and template jet pT --- PWGJE/Tasks/nucleiInJets.cxx | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/PWGJE/Tasks/nucleiInJets.cxx b/PWGJE/Tasks/nucleiInJets.cxx index 1b690d4bd5c..2abfb2f8fa3 100644 --- a/PWGJE/Tasks/nucleiInJets.cxx +++ b/PWGJE/Tasks/nucleiInJets.cxx @@ -48,7 +48,6 @@ #include #include #include -#include #include #include @@ -57,6 +56,7 @@ #include #include #include +#include #include #include @@ -178,7 +178,7 @@ struct nucleiInJets { Configurable isRequireHitsInITSLayers{"isRequireHitsInITSLayers", true, "true: at least one hit in the its inner layes"}; Configurable useMcC{"useMcC", true, "use mcC"}; Configurable mcClosureTrainingFraction{"mcClosureTrainingFraction", 0.5f, "fraction of events assigned to MC correction training; the rest are data-like"}; - Configurable mcClosureSplitSeed{"mcClosureSplitSeed", 0, "TRandom3 seed for MC/data-like split (0: ROOT automatic seed)"}; + Configurable mcClosureSplitSeed{"mcClosureSplitSeed", 1234, "random seed for MC/data-like split (0: automatic seed)"}; Configurable usebkgSubractionMC{"usebkgSubractionMC", true, "use rho-area background subtraction for detector-level matched MC jets"}; Configurable cfgjetPtBkgSubMinMC{"cfgjetPtBkgSubMinMC", 10.0f, "minimum detector-level matched MC jet pT after optional background subtraction"}; Configurable useFullConeAcceptance{"useFullConeAcceptance", true, "require jet/perpendicular-cone axes to satisfy |eta_axis| < cfgtrkMaxEta - cfgjetR for jet and efficiency histograms"}; @@ -269,19 +269,16 @@ struct nucleiInJets { OutputObj zorroSummary{"zorroSummary"}; Service ccdb; Service pdgDB; - TRandom3 randUniform; - UInt_t mcClosureSplitSeedEffective = 1; + uint64_t mcClosureSplitSeedEffective = 1; bool useTrainingHistForMcCollision(int64_t mcCollisionId) { // Use the MC collision ID as part of the seed, so both process functions // obtain the same decision without storing one entry per collision. - UInt_t collisionSeed = mcClosureSplitSeedEffective ^ static_cast(mcCollisionId); - if (collisionSeed == 0U) { - collisionSeed = 1U; - } - TRandom3 collisionRandom(collisionSeed); - return collisionRandom.Rndm() < mcClosureTrainingFraction; + const auto collisionId = static_cast(mcCollisionId); + std::mt19937_64 collisionRandom(mcClosureSplitSeedEffective ^ collisionId); + std::uniform_real_distribution uniform(0.0, 1.0); + return uniform(collisionRandom) < mcClosureTrainingFraction; } void init(o2::framework::InitContext&) @@ -293,10 +290,11 @@ struct nucleiInJets { if (!(mcClosureTrainingFraction >= 0.f && mcClosureTrainingFraction <= 1.f)) { LOGP(fatal, "mcClosureTrainingFraction must be between 0 and 1, inclusive"); } - randUniform.SetSeed(mcClosureSplitSeed); - mcClosureSplitSeedEffective = randUniform.Integer(0xffffffffU); - if (mcClosureSplitSeedEffective == 0U) { - mcClosureSplitSeedEffective = 1U; + if (mcClosureSplitSeed == 0) { + std::random_device randomDevice; + mcClosureSplitSeedEffective = (static_cast(randomDevice()) << 32) ^ randomDevice(); + } else { + mcClosureSplitSeedEffective = static_cast(mcClosureSplitSeed); } jetHist.print(); // const AxisSpec PtAxis = {100, 0, 10.0}; @@ -2678,7 +2676,8 @@ struct nucleiInJets { } // Fill detector-level DCA templates before the DCAxy selection used by the spectra. - const double templateJetPt = mcdJetPt.at(indexJet) - (usebkgSubractionMC ? 0.f : backgroundRho * jetArea); + // mcdJetPt is already background-subtracted when usebkgSubractionMC is enabled. + const double templateJetPt = mcdJetPt.at(indexJet); if (isWithLeadingJet && backgroundRho > 0.f && templateJetPt > cfgjetPtBkgSubMin) { const double templateJetEta = mcdJetEta.at(indexJet); const double templateJetPhi = mcdJetPhi.at(indexJet);