From 7a556a51e369462dbeb7e65c052b50bb1c97f5af Mon Sep 17 00:00:00 2001 From: Zoltan Varga Date: Fri, 27 Mar 2026 04:15:12 +0100 Subject: [PATCH 1/4] Add Lund plane task --- PWGJE/Tasks/CMakeLists.txt | 4 + PWGJE/Tasks/jetLundPlane.cxx | 672 +++++++++++++++++++++++++++++++++++ 2 files changed, 676 insertions(+) create mode 100644 PWGJE/Tasks/jetLundPlane.cxx diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 55ec23b2c44..0275a7a66ed 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -317,6 +317,10 @@ if(FastJet_FOUND) SOURCES jetLundReclustering.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(jet-lund-plane + SOURCES jetLundPlane.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(hf-fragmentation-function SOURCES hfFragmentationFunction.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore diff --git a/PWGJE/Tasks/jetLundPlane.cxx b/PWGJE/Tasks/jetLundPlane.cxx new file mode 100644 index 00000000000..80be0959a97 --- /dev/null +++ b/PWGJE/Tasks/jetLundPlane.cxx @@ -0,0 +1,672 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file jetLundPlane.cxx +/// \brief Task for jet Lund plane. Creates histograms for offline unfolding (including QA histos), and optionally tables. +/// \author Zoltan Varga + +#include +#include +#include +#include +#include + +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" + +#include "PWGJE/Core/JetFinder.h" +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" + +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + + +// Mini-AOD (tables) +namespace o2::aod +{ +DECLARE_SOA_COLUMN(MiniCollId, miniCollId, uint64_t); // collision global index +DECLARE_SOA_COLUMN(MiniJetId, miniJetId, uint64_t); // jet global index (within its table) +DECLARE_SOA_COLUMN(Level, level, uint8_t); // 0=reco(det), 1=truth(part) +DECLARE_SOA_COLUMN(JetRint, jetRint, int32_t); // jet.r() as stored (int R*100) +DECLARE_SOA_COLUMN(JetPt, jetPt, float); +DECLARE_SOA_COLUMN(JetEta, jetEta, float); +DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); + +// Per-splitting observables (primary branch) +DECLARE_SOA_COLUMN(SplitId, splitId, uint16_t); +DECLARE_SOA_COLUMN(DeltaR, deltaR, float); +DECLARE_SOA_COLUMN(PtSoft, ptSoft, float); +DECLARE_SOA_COLUMN(PtHard, ptHard, float); +DECLARE_SOA_COLUMN(SoftEta, softEta, float); +DECLARE_SOA_COLUMN(SoftPhi, softPhi, float); + +// Jet-jet matching (MC) +DECLARE_SOA_COLUMN(DetJetId, detJetId, uint64_t); +DECLARE_SOA_COLUMN(PartJetId, partJetId, uint64_t); +DECLARE_SOA_COLUMN(MatchDR, matchDR, float); +DECLARE_SOA_COLUMN(MatchRelPt, matchRelPt, float); + +DECLARE_SOA_TABLE(MiniJets, "AOD", "MINIJET", + MiniCollId, MiniJetId, Level, JetRint, JetPt, JetEta, JetPhi); + +DECLARE_SOA_TABLE(MiniSplittings, "AOD", "MINISPL", + MiniCollId, MiniJetId, Level, SplitId, DeltaR, PtSoft, PtHard, SoftEta, SoftPhi, JetPt); + +DECLARE_SOA_TABLE(MiniJetMatches, "AOD", "MINIMCH", + MiniCollId, DetJetId, PartJetId, MatchDR, MatchRelPt); +} // namespace o2::aod + + +namespace +{ +constexpr float kTiny = 1e-12f; + +template +float jetRfromTable(const JetT& jet, float fallbackR) +{ + if constexpr (requires { jet.r(); }) { + return static_cast(jet.r()) * 0.01f; // r stored as int(R*100) + } + return fallbackR; +} + +struct SplittingObs { + float deltaR{}; + float ptSoft{}; + float ptHard{}; + float kT{}; + float lnRoverDR{}; + float lnkt{}; + float z{}; + float softEta{}; + float softPhi{}; +}; + +struct SplitMatchPair { + size_t recoIdx{}; + size_t truthIdx{}; + float dR{}; +}; + +inline float deltaPhi(float phi1, float phi2) +{ + return std::remainder(phi1 - phi2, 2.f * static_cast(M_PI)); +} + +inline float splittingMatchDistance(const SplittingObs& a, const SplittingObs& b) +{ + return std::hypot(a.softEta - b.softEta, deltaPhi(a.softPhi, b.softPhi)); +} + +std::vector buildUniqueSplittingMatches(const std::vector& recoSpl, + const std::vector& truthSpl, + float maxDR) +{ + std::vector matches; + matches.reserve(std::min(recoSpl.size(), truthSpl.size())); + + if (recoSpl.empty() || truthSpl.empty()) { + return matches; + } + + // Mutual-nearest-neighbor matching in the soft-prong (eta,phi) coordinates: + // 1) for each truth splitting, find the closest reco splitting within maxDR + // 2) for that reco splitting, find the closest truth splitting within maxDR + // 3) keep the pair only if the reco splitting points back to the original truth splitting + for (size_t truthIdx = 0; truthIdx < truthSpl.size(); ++truthIdx) { + float bestRecoDR = maxDR; + size_t bestRecoIdx = recoSpl.size(); + + for (size_t recoIdx = 0; recoIdx < recoSpl.size(); ++recoIdx) { + const float dR = splittingMatchDistance(recoSpl[recoIdx], truthSpl[truthIdx]); + if (dR < bestRecoDR) { + bestRecoDR = dR; + bestRecoIdx = recoIdx; + } + } + + if (bestRecoIdx == recoSpl.size()) { + continue; + } + + float bestTruthDR = maxDR; + size_t reverseTruthIdx = truthSpl.size(); + for (size_t candTruthIdx = 0; candTruthIdx < truthSpl.size(); ++candTruthIdx) { + const float dR = splittingMatchDistance(recoSpl[bestRecoIdx], truthSpl[candTruthIdx]); + if (dR < bestTruthDR) { + bestTruthDR = dR; + reverseTruthIdx = candTruthIdx; + } + } + + if (reverseTruthIdx == truthIdx) { + matches.push_back({bestRecoIdx, truthIdx, bestRecoDR}); + } + } + + return matches; +} + +template +std::vector buildFastJetInputs(ConstituentRangeT&& constituents, float trackPtMin) +{ + std::vector fjInputs; + fjInputs.reserve(64); + for (auto const& c : constituents) { + if (c.pt() < trackPtMin) { + continue; + } + const float mPi = 0.13957f; // GeV/c^2 + const float e = std::sqrt(c.p() * c.p() + mPi * mPi); + fjInputs.emplace_back(c.px(), c.py(), c.pz(), e); + } + return fjInputs; +} + +std::vector primaryDeclusteringSplittings(fastjet::PseudoJet jetCA, float jetR) +{ + std::vector out; + out.reserve(32); + + fastjet::PseudoJet p1, p2; + while (jetCA.has_parents(p1, p2)) { + if (p2.perp() > p1.perp()) { + std::swap(p1, p2); + } + const float pt1 = p1.perp(); + const float pt2 = p2.perp(); + const float pt = pt1 + pt2; + if (pt <= kTiny) { + break; + } + const float z = pt2 / pt; + const float dR = p1.delta_R(p2); + const float kT = pt2 * std::sin(dR); + + if (dR > kTiny && kT > kTiny && z > kTiny) { + SplittingObs s; + s.deltaR = dR; + s.ptSoft = pt2; + s.ptHard = pt1; + s.kT = kT; + s.lnRoverDR = std::log(jetR / dR); + s.lnkt = std::log(kT); + s.z = z; + s.softEta = p2.eta(); + s.softPhi = p2.phi_std(); + out.push_back(s); + } + jetCA = p1; // follow hard branch (primary Lund) + } + return out; +} + +} // namespace + +struct JetLundPlaneUnfolding +{ + // Config + Configurable vertexZCut{"vertexZCut", 10.f, "|z_vtx| cut"}; + Configurable jetPtMin{"jetPtMin", 20.f, "min reco jet pT"}; + Configurable jetEtaMin{"jetEtaMin", -0.5f, "min jet eta"}; + Configurable jetEtaMax{"jetEtaMax", 0.5f, "max jet eta"}; + Configurable jetR{"jetR", 0.4f, "jet radius (must match derived tables)"}; + Configurable trackPtMin{"trackPtMin", 0.15f, "min constituent pT"}; + + Configurable nBinsJetPt{"nBinsJetPt", 200, "jet pT bins"}; + Configurable jetPtMax{"jetPtMax", 200.f, "jet pT max"}; + + Configurable nBinsLund{"nBinsLund", 120, "lund bins (lnR/DR and lnkt)"}; + Configurable lnRoverDRMin{"lnRoverDRMin", -1.f, "min ln(R/DR)"}; + Configurable lnRoverDRMax{"lnRoverDRMax", 6.f, "max ln(R/DR)"}; + Configurable lnKtMin{"lnKtMin", -6.f, "min ln(kT)"}; + Configurable lnKtMax{"lnKtMax", 6.f, "max ln(kT)"}; + + // switches (runtime) + Configurable writeMiniAOD{"writeMiniAOD", true, "write mini-AOD tables (jets, splittings, matching)"}; + Configurable doData{"doData", true, "run reco/data processing"}; + Configurable doMC{"doMC", false, "run MC processing (det+part+response)"}; + + // matching knobs (applied on top of matching relations) + Configurable matchMaxDR{"matchMaxDR", 0.2f, "max ΔR between det and part jet"}; + Configurable splitMatchMaxDR{"splitMatchMaxDR", 0.1f, "max ΔR between reco and truth splittings (soft prong)"}; + Configurable matchUseRelPt{"matchUseRelPt", true, "apply relative pT compatibility cut"}; + Configurable matchMaxRelPtDiff{"matchMaxRelPtDiff", 0.5f, "max |pTdet-pTpart|/pTpart"}; + + // Registry + HistogramRegistry registry{"registry"}; + + // Mini-AOD outputs (optional) + Produces outMiniJets; + Produces outMiniSplittings; + Produces outMiniJetMatches; + + // FastJet reclustering setup (C/A) + JetFinder reclusterer; + std::vector jetReclustered; + + void init(InitContext const&) + { + const int nbPt = nBinsJetPt.value; + const float ptMax = jetPtMax.value; + const int nbL = nBinsLund.value; + + registry.add("hEventCount", "Event counter;step;counts", HistType::kTH1F, {{10, 0.f, 10.f}}); + registry.add("hJetCountSummary", "Jet count summary;category;counts", HistType::kTH1F, {{6, 0.5f, 6.5f}}); + + // Jet spectra / matching QA + registry.add("hJetPtRecoAll", "Reco jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + registry.add("hJetPtRecoMatched", "Reco matched jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + registry.add("hJetPtRecoFake", "Reco unmatched (fake) jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + + registry.add("hJetPtTruthAll", "Truth jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + registry.add("hJetPtTruthMatched", "Truth matched jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + registry.add("hJetPtTruthMiss", "Truth unmatched (miss) jets;p_{T}^{jet} (GeV/c);counts", HistType::kTH1F, {{nbPt, 0.f, ptMax}}); + + // Jet pT response (1D unfolding) + registry.add("hJetPtResponse", "Jet pT response;p_{T}^{det};p_{T}^{part}", HistType::kTH2F, + {{nbPt, 0.f, ptMax}, {nbPt, 0.f, ptMax}}); + + // Lund reco / truth (3D) + registry.add("hLundReco3D", "Primary Lund (reco);ln(R/#DeltaR);ln(k_{T});p_{T}^{jet}", + HistType::kTH3F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, + {nbL, lnKtMin.value, lnKtMax.value}, + {nbPt, 0.f, ptMax}}); + + registry.add("hLundTruth3D", "Primary Lund (truth);ln(R/#DeltaR);ln(k_{T});p_{T}^{jet}", + HistType::kTH3F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, + {nbL, lnKtMin.value, lnKtMax.value}, + {nbPt, 0.f, ptMax}}); + + // 6D response for 3D unfolding: + // (det lnR/DR, det lnkt, det jetpt, part lnR/DR, part lnkt, part jetpt) + registry.add("hLundResponse6D", + "Lund response 6D;" + "ln(R/#DeltaR)_{det};ln(k_{T})_{det};p_{T}^{det};" + "ln(R/#DeltaR)_{part};ln(k_{T})_{part};p_{T}^{part}", + HistType::kTHnSparseF, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, + {nbL, lnKtMin.value, lnKtMax.value}, + {nbPt, 0.f, ptMax}, + {nbL, lnRoverDRMin.value, lnRoverDRMax.value}, + {nbL, lnKtMin.value, lnKtMax.value}, + {nbPt, 0.f, ptMax}}); + + // Early QA histograms for declustering and matching + registry.add("hNSplittingsReco", "Number of primary splittings (reco);N_{split};counts", HistType::kTH1F, {{40, -0.5f, 39.5f}}); + registry.add("hNSplittingsTruth", "Number of primary splittings (truth);N_{split};counts", HistType::kTH1F, {{40, -0.5f, 39.5f}}); + registry.add("hNSplittingsVsJetPtReco", "Reco primary splittings;#it{p}_{T}^{jet} (GeV/#it{c});N_{split}", HistType::kTH2F, {{nbPt, 0.f, ptMax}, {40, -0.5f, 39.5f}}); + registry.add("hNSplittingsVsJetPtTruth", "Truth primary splittings;#it{p}_{T}^{jet} (GeV/#it{c});N_{split}", HistType::kTH2F, {{nbPt, 0.f, ptMax}, {40, -0.5f, 39.5f}}); + + registry.add("hDeltaRReco", "Reco splitting #DeltaR;#DeltaR;counts", HistType::kTH1F, {{120, 0.f, 1.2f}}); + registry.add("hDeltaRTruth", "Truth splitting #DeltaR;#DeltaR;counts", HistType::kTH1F, {{120, 0.f, 1.2f}}); + registry.add("hLnRoverDRReco", "Reco splitting ln(R/#DeltaR);ln(R/#DeltaR);counts", HistType::kTH1F, {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}}); + registry.add("hLnRoverDRTruth", "Truth splitting ln(R/#DeltaR);ln(R/#DeltaR);counts", HistType::kTH1F, {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}}); + registry.add("hLnKtReco", "Reco splitting ln(k_{T});ln(k_{T});counts", HistType::kTH1F, {{nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hLnKtTruth", "Truth splitting ln(k_{T});ln(k_{T});counts", HistType::kTH1F, {{nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hKtReco", "Reco splitting k_{T};k_{T};counts", HistType::kTH1F, {{200, 0.f, 20.f}}); + registry.add("hKtTruth", "Truth splitting k_{T};k_{T};counts", HistType::kTH1F, {{200, 0.f, 20.f}}); + + registry.add("hJetMatchDR", "Matched jet #DeltaR;#DeltaR(det,part);counts", HistType::kTH1F, {{100, 0.f, 1.f}}); + registry.add("hJetMatchRelPt", "Matched jet relative #it{p}_{T} difference;|#it{p}_{T}^{det}-#it{p}_{T}^{part}|/#it{p}_{T}^{part};counts", HistType::kTH1F, {{100, 0.f, 2.f}}); + registry.add("hJetPtResidual", "Jet #it{p}_{T} residual;(#it{p}_{T}^{det}-#it{p}_{T}^{part})/#it{p}_{T}^{part};counts", HistType::kTH1F, {{160, -2.f, 2.f}}); + registry.add("hLnRoverDRResidual", "Splitting ln(R/#DeltaR) residual;ln(R/#DeltaR)_{det}-ln(R/#DeltaR)_{part};counts", HistType::kTH1F, {{160, -4.f, 4.f}}); + registry.add("hLnKtResidual", "Splitting ln(k_{T}) residual;ln(k_{T})_{det}-ln(k_{T})_{part};counts", HistType::kTH1F, {{160, -6.f, 6.f}}); + registry.add("hSplitMatchDR", "Matched splitting #DeltaR (soft prong);#DeltaR_{match};counts", HistType::kTH1F, {{120, 0.f, 1.2f}}); + + registry.add("hSplitTruthDen", "Truth splitting denominator;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hSplitTruthMatched", "Truth matched splitting numerator;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hSplitTruthMiss", "Truth missed splittings;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + + registry.add("hSplitRecoDen", "Reco splitting denominator;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hSplitRecoMatched", "Reco matched splitting numerator;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + registry.add("hSplitRecoFake", "Reco fake splittings;ln(R/#DeltaR);ln(k_{T})", HistType::kTH2F, + {{nbL, lnRoverDRMin.value, lnRoverDRMax.value}, {nbL, lnKtMin.value, lnKtMax.value}}); + + // reclusterer config + reclusterer.isReclustering = true; + reclusterer.algorithm = fastjet::cambridge_algorithm; + reclusterer.jetR = 1.0; // recluster radius for declustering tree + reclusterer.recombScheme = fastjet::E_scheme; + reclusterer.strategy = fastjet::Best; + reclusterer.areaType = fastjet::active_area; + } + + // Filters for reco jets (data or MC reco) + Filter collisionFilter = (nabs(aod::jcollision::posZ) < vertexZCut.node()); + Filter jetFilter = (aod::jet::pt > jetPtMin.node()) && + (aod::jet::eta > jetEtaMin.node()) && + (aod::jet::eta < jetEtaMax.node()) && + (aod::jet::r == nround(jetR.node() * 100.f)); + + // Type aliases + using RecoJets = soa::Join; + using DetJetsMatched = soa::Join; + using PartJetsMatched = soa::Join; + + template + bool passJetFiducial(JetT const& jet, int rWanted) const + { + return (jet.r() == rWanted) && (jet.pt() > jetPtMin.value) && (jet.eta() > jetEtaMin.value) && (jet.eta() < jetEtaMax.value); + } + + template + uint64_t getCollisionGlobalIndex(JetT const& jet, aod::JetCollisions const& collisions) const + { + if constexpr (requires { jet.collisionId(); }) { + return collisions.iteratorAt(jet.collisionId()).globalIndex(); + } else if constexpr (requires { jet.mcCollisionId(); }) { + return static_cast(jet.mcCollisionId()); + } else { + return jet.globalIndex(); + } + } + + void fillJetCountSummary(int bin) + { + registry.fill(HIST("hJetCountSummary"), static_cast(bin)); + } + + // Helpers to fill Lund from a jet row + template + std::vector getPrimarySplittings(JetRowT const& jet, ConstituentTableT const&) + { + auto fjInputs = buildFastJetInputs(jet.template tracks_as(), trackPtMin.value); + if (fjInputs.size() < 2) { + return {}; + } + + jetReclustered.clear(); + fastjet::ClusterSequenceArea cs(reclusterer.findJets(fjInputs, jetReclustered)); + jetReclustered = sorted_by_pt(jetReclustered); + if (jetReclustered.empty()) { + return {}; + } + + const float rjet = jetRfromTable(jet, jetR.value); + return primaryDeclusteringSplittings(jetReclustered[0], rjet); + } + + void fillPrimaryLund3DFromSplittings(std::vector const& spl, float jetPt, bool isTruth, float weight = 1.f) + { + for (auto const& s : spl) { + if (isTruth) { + registry.fill(HIST("hLundTruth3D"), s.lnRoverDR, s.lnkt, jetPt, weight); + } else { + registry.fill(HIST("hLundReco3D"), s.lnRoverDR, s.lnkt, jetPt, weight); + } + } + } + + template + void fillPrimaryLund3D(JetRowT const& jet, ConstituentTableT const& constituents, bool isTruth, float weight = 1.f) + { + fillPrimaryLund3DFromSplittings(getPrimarySplittings(jet, constituents), jet.pt(), isTruth, weight); + } + + void fillSplittingQAHists(std::vector const& splittings, bool isTruth, float jetPt) + { + if (isTruth) { + registry.fill(HIST("hNSplittingsTruth"), static_cast(splittings.size())); + registry.fill(HIST("hNSplittingsVsJetPtTruth"), jetPt, static_cast(splittings.size())); + } else { + registry.fill(HIST("hNSplittingsReco"), static_cast(splittings.size())); + registry.fill(HIST("hNSplittingsVsJetPtReco"), jetPt, static_cast(splittings.size())); + } + + for (auto const& s : splittings) { + if (isTruth) { + registry.fill(HIST("hDeltaRTruth"), s.deltaR); + registry.fill(HIST("hLnRoverDRTruth"), s.lnRoverDR); + registry.fill(HIST("hLnKtTruth"), s.lnkt); + registry.fill(HIST("hKtTruth"), s.kT); + } else { + registry.fill(HIST("hDeltaRReco"), s.deltaR); + registry.fill(HIST("hLnRoverDRReco"), s.lnRoverDR); + registry.fill(HIST("hLnKtReco"), s.lnkt); + registry.fill(HIST("hKtReco"), s.kT); + } + } + } + + std::vector fillSplittingCorrectionHists(std::vector const& recoSpl, std::vector const& truthSpl) + { + for (auto const& s : truthSpl) { + registry.fill(HIST("hSplitTruthDen"), s.lnRoverDR, s.lnkt); + } + for (auto const& s : recoSpl) { + registry.fill(HIST("hSplitRecoDen"), s.lnRoverDR, s.lnkt); + } + + const auto matches = buildUniqueSplittingMatches(recoSpl, truthSpl, splitMatchMaxDR.value); + std::vector recoUsed(recoSpl.size(), false); + std::vector truthUsed(truthSpl.size(), false); + + for (const auto& m : matches) { + recoUsed[m.recoIdx] = true; + truthUsed[m.truthIdx] = true; + registry.fill(HIST("hSplitRecoMatched"), recoSpl[m.recoIdx].lnRoverDR, recoSpl[m.recoIdx].lnkt); + registry.fill(HIST("hSplitTruthMatched"), truthSpl[m.truthIdx].lnRoverDR, truthSpl[m.truthIdx].lnkt); + registry.fill(HIST("hSplitMatchDR"), m.dR); + } + + for (size_t i = 0; i < truthSpl.size(); ++i) { + if (!truthUsed[i]) { + registry.fill(HIST("hSplitTruthMiss"), truthSpl[i].lnRoverDR, truthSpl[i].lnkt); + } + } + for (size_t i = 0; i < recoSpl.size(); ++i) { + if (!recoUsed[i]) { + registry.fill(HIST("hSplitRecoFake"), recoSpl[i].lnRoverDR, recoSpl[i].lnkt); + } + } + return matches; + } + + // DATA / RECO PROCESSING + void processData(soa::Filtered::iterator const& collision, + soa::Filtered const& jets, + aod::JetTracks const& tracks) + { + if (!doData.value) { + return; + } + registry.fill(HIST("hEventCount"), 0.5); + for (auto const& jet : jets) { + registry.fill(HIST("hJetPtRecoAll"), jet.pt()); + auto spl = getPrimarySplittings(jet, tracks); + fillPrimaryLund3DFromSplittings(spl, jet.pt(), /*isTruth*/ false); + fillSplittingQAHists(spl, /*isTruth*/ false, jet.pt()); + + if (writeMiniAOD.value) { + const uint64_t collId = collision.globalIndex(); + const uint64_t jetId = jet.globalIndex(); + outMiniJets(collId, jetId, /*level*/ (uint8_t)0, jet.r(), jet.pt(), jet.eta(), jet.phi()); + uint16_t sid = 0; + for (auto const& s : spl) { + outMiniSplittings(collId, jetId, /*level*/ (uint8_t)0, sid++, s.deltaR, s.ptSoft, s.ptHard, s.softEta, s.softPhi, jet.pt()); + } + } + } + } + PROCESS_SWITCH(JetLundPlaneUnfolding, processData, "Reco/data Lund + jet spectra", true); + + // MC PROCESSING (det + part + response) + void processMC(DetJetsMatched const& detJets, + PartJetsMatched const& partJets, + aod::JetCollisions const& collisions, + aod::JetTracks const& tracks, + aod::JetParticles const& particles) + { + if (!doMC.value) { + return; + } + registry.fill(HIST("hEventCount"), 1.5); + + const int rWanted = static_cast(std::lround(jetR.value * 100.f)); + std::unordered_map truthMatchedById; + auto h6 = registry.get(HIST("hLundResponse6D")); + + // Loop over detector-level jets and select the best truth match from the already-produced matching relation + for (auto const& detJet : detJets) { + if (!passJetFiducial(detJet, rWanted)) { + continue; + } + + registry.fill(HIST("hJetPtRecoAll"), detJet.pt()); + fillJetCountSummary(1); + + const uint64_t detCollId = getCollisionGlobalIndex(detJet, collisions); + const uint64_t detId = detJet.globalIndex(); + auto detSpl = getPrimarySplittings(detJet, tracks); + fillPrimaryLund3DFromSplittings(detSpl, detJet.pt(), /*isTruth*/ false); + fillSplittingQAHists(detSpl, /*isTruth*/ false, detJet.pt()); + + if (writeMiniAOD.value) { + outMiniJets(detCollId, detId, /*level*/ (uint8_t)0, detJet.r(), detJet.pt(), detJet.eta(), detJet.phi()); + uint16_t sidTmp = 0; + for (auto const& s : detSpl) { + outMiniSplittings(detCollId, detId, /*level*/ (uint8_t)0, sidTmp++, s.deltaR, s.ptSoft, s.ptHard, s.softEta, s.softPhi, detJet.pt()); + } + } + + if (!detJet.has_matchedJetGeo()) { + registry.fill(HIST("hJetPtRecoFake"), detJet.pt()); + fillJetCountSummary(3); + continue; + } + + float bestDR = 1e9f; + float bestRelPt = -1.f; + bool foundMatch = false; + + uint64_t bestPartId = 0; + float bestPartPt = 0.f; + std::vector bestPartSpl; + + for (auto const& candPartJet : detJet.template matchedJetGeo_as()) { + if (!passJetFiducial(candPartJet, rWanted)) { + continue; + } + + const float dR = std::hypot(detJet.eta() - candPartJet.eta(), + std::remainder(detJet.phi() - candPartJet.phi(), 2.f * (float)M_PI)); + if (dR > matchMaxDR.value) { + continue; + } + + const float rel = std::abs(detJet.pt() - candPartJet.pt()) / std::max(candPartJet.pt(), kTiny); + if (matchUseRelPt.value && rel > matchMaxRelPtDiff.value) { + continue; + } + + if (dR < bestDR) { + bestDR = dR; + bestRelPt = rel; + bestPartId = candPartJet.globalIndex(); + bestPartPt = candPartJet.pt(); + bestPartSpl = getPrimarySplittings(candPartJet, particles); + foundMatch = true; + } + } + + if (!foundMatch) { + registry.fill(HIST("hJetPtRecoFake"), detJet.pt()); + fillJetCountSummary(3); + continue; + } + + registry.fill(HIST("hJetPtRecoMatched"), detJet.pt()); + fillJetCountSummary(2); + registry.fill(HIST("hJetPtTruthMatched"), bestPartPt); + fillJetCountSummary(5); + registry.fill(HIST("hJetPtResponse"), detJet.pt(), bestPartPt); + registry.fill(HIST("hJetMatchDR"), bestDR); + registry.fill(HIST("hJetMatchRelPt"), bestRelPt); + registry.fill(HIST("hJetPtResidual"), (detJet.pt() - bestPartPt) / std::max(bestPartPt, kTiny)); + truthMatchedById[bestPartId] = true; + + if (writeMiniAOD.value) { + outMiniJetMatches(detCollId, detId, bestPartId, bestDR, bestRelPt); + } + + const auto splitMatches = fillSplittingCorrectionHists(detSpl, bestPartSpl); + + // Fill the 6D response and residual QA using the unique geometric splitting matches + for (const auto& m : splitMatches) { + const auto& detS = detSpl[m.recoIdx]; + const auto& partS = bestPartSpl[m.truthIdx]; + double x[6] = {detS.lnRoverDR, detS.lnkt, detJet.pt(), + partS.lnRoverDR, partS.lnkt, bestPartPt}; + h6->Fill(x); + registry.fill(HIST("hLnRoverDRResidual"), detS.lnRoverDR - partS.lnRoverDR); + registry.fill(HIST("hLnKtResidual"), detS.lnkt - partS.lnkt); + } + } + + // Loop over all truth jets to fill the inclusive truth spectra and identify misses + for (auto const& partJet : partJets) { + if (!passJetFiducial(partJet, rWanted)) { + continue; + } + + registry.fill(HIST("hJetPtTruthAll"), partJet.pt()); + fillJetCountSummary(4); + auto spl = getPrimarySplittings(partJet, particles); + fillPrimaryLund3DFromSplittings(spl, partJet.pt(), /*isTruth*/ true); + fillSplittingQAHists(spl, /*isTruth*/ true, partJet.pt()); + + const uint64_t collId = getCollisionGlobalIndex(partJet, collisions); + const uint64_t jetId = partJet.globalIndex(); + + if (writeMiniAOD.value) { + outMiniJets(collId, jetId, /*level*/ (uint8_t)1, partJet.r(), partJet.pt(), partJet.eta(), partJet.phi()); + uint16_t sid = 0; + for (auto const& s : spl) { + outMiniSplittings(collId, jetId, /*level*/ (uint8_t)1, sid++, s.deltaR, s.ptSoft, s.ptHard, s.softEta, s.softPhi, partJet.pt()); + } + } + + if (!truthMatchedById[jetId]) { + registry.fill(HIST("hJetPtTruthMiss"), partJet.pt()); + fillJetCountSummary(6); + } + } + } + PROCESS_SWITCH(JetLundPlaneUnfolding, processMC, "MC det+part + responses", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"jet-lund-plane"})}; +} \ No newline at end of file From 45b23eb8f32e20aaaf89c8bb58ac2c61c1acd570 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 27 Mar 2026 07:24:18 +0000 Subject: [PATCH 2/4] Please consider the following formatting changes --- PWGJE/Tasks/jetLundPlane.cxx | 47 +++++++++++++++++------------------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/PWGJE/Tasks/jetLundPlane.cxx b/PWGJE/Tasks/jetLundPlane.cxx index 80be0959a97..d6b6eb8d7a8 100644 --- a/PWGJE/Tasks/jetLundPlane.cxx +++ b/PWGJE/Tasks/jetLundPlane.cxx @@ -13,39 +13,38 @@ /// \brief Task for jet Lund plane. Creates histograms for offline unfolding (including QA histos), and optionally tables. /// \author Zoltan Varga -#include -#include -#include -#include -#include +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/Core/JetFinder.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" -#include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" -#include "PWGJE/Core/JetFinder.h" -#include "PWGJE/Core/FastJetUtilities.h" -#include "PWGJE/DataModel/Jet.h" -#include "PWGJE/DataModel/JetReducedData.h" - #include -#include #include +#include + +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; - // Mini-AOD (tables) namespace o2::aod { -DECLARE_SOA_COLUMN(MiniCollId, miniCollId, uint64_t); // collision global index -DECLARE_SOA_COLUMN(MiniJetId, miniJetId, uint64_t); // jet global index (within its table) -DECLARE_SOA_COLUMN(Level, level, uint8_t); // 0=reco(det), 1=truth(part) -DECLARE_SOA_COLUMN(JetRint, jetRint, int32_t); // jet.r() as stored (int R*100) +DECLARE_SOA_COLUMN(MiniCollId, miniCollId, uint64_t); // collision global index +DECLARE_SOA_COLUMN(MiniJetId, miniJetId, uint64_t); // jet global index (within its table) +DECLARE_SOA_COLUMN(Level, level, uint8_t); // 0=reco(det), 1=truth(part) +DECLARE_SOA_COLUMN(JetRint, jetRint, int32_t); // jet.r() as stored (int R*100) DECLARE_SOA_COLUMN(JetPt, jetPt, float); DECLARE_SOA_COLUMN(JetEta, jetEta, float); DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); @@ -74,7 +73,6 @@ DECLARE_SOA_TABLE(MiniJetMatches, "AOD", "MINIMCH", MiniCollId, DetJetId, PartJetId, MatchDR, MatchRelPt); } // namespace o2::aod - namespace { constexpr float kTiny = 1e-12f; @@ -221,8 +219,7 @@ std::vector primaryDeclusteringSplittings(fastjet::PseudoJet jetCA } // namespace -struct JetLundPlaneUnfolding -{ +struct JetLundPlaneUnfolding { // Config Configurable vertexZCut{"vertexZCut", 10.f, "|z_vtx| cut"}; Configurable jetPtMin{"jetPtMin", 20.f, "min reco jet pT"}; @@ -367,11 +364,11 @@ struct JetLundPlaneUnfolding // Type aliases using RecoJets = soa::Join; using DetJetsMatched = soa::Join; + aod::ChargedMCDetectorLevelJetConstituents, + aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>; using PartJetsMatched = soa::Join; + aod::ChargedMCParticleLevelJetConstituents, + aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>; template bool passJetFiducial(JetT const& jet, int rWanted) const @@ -669,4 +666,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc, TaskName{"jet-lund-plane"})}; -} \ No newline at end of file +} From 84994069d95514072893139446d0fa94fab2a3c7 Mon Sep 17 00:00:00 2001 From: Zoltan Varga <91464454+zovarga@users.noreply.github.com> Date: Fri, 27 Mar 2026 08:50:16 +0100 Subject: [PATCH 3/4] Update jetLundPlane.cxx --- PWGJE/Tasks/jetLundPlane.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGJE/Tasks/jetLundPlane.cxx b/PWGJE/Tasks/jetLundPlane.cxx index d6b6eb8d7a8..e55127a3aad 100644 --- a/PWGJE/Tasks/jetLundPlane.cxx +++ b/PWGJE/Tasks/jetLundPlane.cxx @@ -33,6 +33,7 @@ #include #include #include +#include using namespace o2; using namespace o2::framework; @@ -576,7 +577,7 @@ struct JetLundPlaneUnfolding { } const float dR = std::hypot(detJet.eta() - candPartJet.eta(), - std::remainder(detJet.phi() - candPartJet.phi(), 2.f * (float)M_PI)); + std::remainder(detJet.phi() - candPartJet.phi(), 2.f * static_cast(M_PI))); if (dR > matchMaxDR.value) { continue; } From 053aa434e3cc0dff8b960ad1accc94c3df928605 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 27 Mar 2026 07:50:55 +0000 Subject: [PATCH 4/4] Please consider the following formatting changes --- PWGJE/Tasks/jetLundPlane.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGJE/Tasks/jetLundPlane.cxx b/PWGJE/Tasks/jetLundPlane.cxx index e55127a3aad..9ed4acfb398 100644 --- a/PWGJE/Tasks/jetLundPlane.cxx +++ b/PWGJE/Tasks/jetLundPlane.cxx @@ -32,8 +32,8 @@ #include #include #include -#include #include +#include using namespace o2; using namespace o2::framework;