From 25dfd1d8a777ba12e1bcb0b5714128c2a812fe27 Mon Sep 17 00:00:00 2001 From: Monalisa Melo Date: Tue, 3 Mar 2026 16:27:47 -0300 Subject: [PATCH] [PWGJE] Add angularity --- PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx | 203 ++++++++++++++----- 1 file changed, 151 insertions(+), 52 deletions(-) diff --git a/PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx b/PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx index a7810742abe..3e2e12d378a 100644 --- a/PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx +++ b/PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx @@ -30,7 +30,6 @@ #include "Framework/ASoA.h" #include "Framework/AnalysisDataModel.h" #include "Framework/HistogramRegistry.h" -#include #include #include #include @@ -42,7 +41,6 @@ #include "TVector3.h" #include -#include #include #include @@ -67,6 +65,12 @@ DECLARE_SOA_COLUMN(HfY, hfY, float); DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float); DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float); DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float); + +// extra +DECLARE_SOA_COLUMN(JetMass, jetMass, float); +DECLARE_SOA_COLUMN(JetGirth, jetGirth, float); +DECLARE_SOA_COLUMN(JetThrust, jetThrust, float); // lambda_2^1 +DECLARE_SOA_COLUMN(JetLambda11, jetLambda11, float); // lambda_1^1 } // namespace jet_distance DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE", @@ -82,46 +86,117 @@ DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE", jet_distance::HfY, jet_distance::HfMlScore0, jet_distance::HfMlScore1, - jet_distance::HfMlScore2); + jet_distance::HfMlScore2, + jet_distance::JetMass, + jet_distance::JetGirth, + jet_distance::JetThrust, + jet_distance::JetLambda11); } // namespace o2::aod struct JetDsSpecSubs { - HistogramRegistry registry{"registry", - {{"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}}, - {"h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, - {"h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, - {"h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, - {"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, - {"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, - {"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, - {"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{200, 0., 200.}}}}, - {"h_jet_counter", ";# of D_{S} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, - {"h_ds_jet_projection", ";z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}", {HistType::kTH1F, {{1000, 0., 2.}}}}, - {"h_ds_jet_distance_vs_projection", ";#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}", {HistType::kTH2F, {{1000, 0., 1.}, {1000, 0., 2.}}}}, - {"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 1.}}}}, - {"h_ds_jet_pt", ";p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}", {HistType::kTH1F, {{1000, 0., 100.}}}}, - {"h_ds_jet_eta", ";#eta_{T,D_{S} jet};dN/d#eta_{D_{S} jet}", {HistType::kTH1F, {{250, -1., 1.}}}}, - {"h_ds_jet_phi", ";#phi_{T,D_{S} jet};dN/d#phi_{D_{S} jet}", {HistType::kTH1F, {{250, -1., 7.}}}}, - {"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});dN/dm_{D_{S}}", {HistType::kTH1F, {{1000, 0., 6.}}}}, - {"h_ds_eta", ";#eta_{D_{S}} (GeV/c^{2});dN/d#eta_{D_{S}}", {HistType::kTH1F, {{250, -1., 1.}}}}, - {"h_ds_phi", ";#phi_{D_{S}} (GeV/c^{2});dN/d#phi_{D_{S}}", {HistType::kTH1F, {{250, -1., 7.}}}}}}; - Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; + HistogramRegistry registry{ + "registry", + { + {"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}}, + {"h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + {"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + {"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_jet_counter", ";# of D_{S} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, + {"h_ds_jet_projection", ";z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}", {HistType::kTH1F, {{1000, 0., 2.}}}}, + {"h_ds_jet_distance_vs_projection", ";#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}", {HistType::kTH2F, {{1000, 0., 1.}, {1000, 0., 2.}}}}, + {"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 1.}}}}, + {"h_ds_jet_pt", ";p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}", {HistType::kTH1F, {{1000, 0., 100.}}}}, + {"h_ds_jet_eta", ";#eta_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_jet_phi", ";#phi_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + {"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});entries", {HistType::kTH1F, {{1000, 0., 6.}}}}, + {"h_ds_eta", ";#eta_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_phi", ";#phi_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + {"h_ds_jet_mass", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}}, + {"h_ds_jet_lambda11", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"h_ds_jet_lambda12", ";#lambda_{2}^{1} (thrust);entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"h_ds_jet_girth", ";g (#equiv #lambda_{1}^{1}R);entries", {HistType::kTH1F, {{200, 0., 0.5}}}}, + {"h2_dsjet_pt_lambda11", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{1}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}}, + {"h2_dsjet_pt_lambda12", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{2}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}}, + {"h2_dsjet_pt_mass", ";#it{p}_{T,jet} (GeV/#it{c});m_{jet}^{ch} (GeV/#it{c}^{2})", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 50.0}}}}, + {"h2_dsjet_pt_girth", ";#it{p}_{T,jet} (GeV/#it{c});g", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 0.5}}}}, + {"h_ds_jet_lambda_extra", ";#lambda_{#alpha}^{#kappa};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"h2_dsjet_pt_lambda_extra", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{#alpha}^{#kappa}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}}, + }}; + Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; Configurable jetR{"jetR", 0.4, "jet resolution parameter"}; Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; + // extra angularity knob + Configurable kappa{"kappa", 1.0f, "angularity kappa"}; + Configurable alpha{"alpha", 1.0f, "angularity alpha"}; + + bool doExtraAngularity = false; + std::vector eventSelectionBits; int trackSelection = -1; Produces distJetTable; + template + float computeLambda(JET const& jet, TRACKS const& tracks, float a, float k) + { + if (jet.pt() <= 0.f) { + return -1.f; + } + float sum = 0.f; + for (auto const& trk : tracks) { + const float dr = jetutilities::deltaR(jet, trk); + sum += std::pow(trk.pt(), k) * std::pow(dr, a); + } + const float R = jet.r() / 100.f; + const float denom = std::pow(jet.pt(), k) * std::pow(R, a); + if (denom <= 0.f) { + return -1.f; + } + return sum / denom; + } + + template + float computeJetMassFromTracksMass(TRACKS const& tracks) + { + double sumPx = 0.0, sumPy = 0.0, sumPz = 0.0, sumE = 0.0; + + for (auto const& trk : tracks) { + const double pt = trk.pt(); + const double phi = trk.phi(); + const double eta = trk.eta(); + + const double px = pt * std::cos(phi); + const double py = pt * std::sin(phi); + const double pz = pt * std::sinh(eta); + const double p = std::sqrt(px * px + py * py + pz * pz); + + sumPx += px; + sumPy += py; + sumPz += pz; + sumE += p; // massless + } + + const double m2 = sumE * sumE - (sumPx * sumPx + sumPy * sumPy + sumPz * sumPz); + return (m2 > 0.0) ? static_cast(std::sqrt(m2)) : 0.f; + } + void init(o2::framework::InitContext&) { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + + const bool is11 = (std::abs(kappa.value - 1.f) < 1e-6f) && (std::abs(alpha.value - 1.f) < 1e-6f); + const bool is12 = (std::abs(kappa.value - 1.f) < 1e-6f) && (std::abs(alpha.value - 2.f) < 1e-6f); + doExtraAngularity = !(is11 || is12); } Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); @@ -129,12 +204,12 @@ struct JetDsSpecSubs { void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks) { - registry.fill(HIST("h_collisions"), 0.5); if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } registry.fill(HIST("h_collisions"), 1.5); + for (auto const& track : tracks) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { continue; @@ -146,12 +221,12 @@ struct JetDsSpecSubs { } PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "process JE collisions", false); - void processDataCharged(soa::Filtered::iterator const& collision, soa::Filtered const& jets) + void processDataCharged(soa::Filtered::iterator const& collision, + soa::Filtered const& jets) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } - // jets -> charged jet for (auto& jet : jets) { registry.fill(HIST("h_jet_pt"), jet.pt()); registry.fill(HIST("h_jet_eta"), jet.eta()); @@ -165,31 +240,24 @@ struct JetDsSpecSubs { aod::CandidatesDsData const&, aod::JetTracks const&) { - // apply event selection and fill histograms for sanity check registry.fill(HIST("h_collision_counter"), 2.0); - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || + !(std::abs(collision.posZ()) < vertexZCut)) { return; } registry.fill(HIST("h_collision_counter"), 3.0); - // jets -> charged jet with Ds + for (const auto& jet : jets) { - // number of charged jets with Ds registry.fill(HIST("h_jet_counter"), 0.5); - // obtaining jet 3-vector + TVector3 jetVector(jet.px(), jet.py(), jet.pz()); for (const auto& dsCandidate : jet.candidates_as()) { - - // obtaining jet 3-vector TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz()); - // calculating fraction of the jet momentum carried by the Ds along the direction of the jet axis - double zParallel = (jetVector * dsVector) / (jetVector * jetVector); - - // calculating angular distance in eta-phi plane - double axisDistance = jetutilities::deltaR(jet, dsCandidate); + const double zParallel = (jetVector * dsVector) / (jetVector * jetVector); + const double axisDistance = jetutilities::deltaR(jet, dsCandidate); - // filling histograms registry.fill(HIST("h_ds_jet_projection"), zParallel); registry.fill(HIST("h_ds_jet_distance_vs_projection"), axisDistance, zParallel); registry.fill(HIST("h_ds_jet_distance"), axisDistance); @@ -200,25 +268,56 @@ struct JetDsSpecSubs { registry.fill(HIST("h_ds_eta"), dsCandidate.eta()); registry.fill(HIST("h_ds_phi"), dsCandidate.phi()); - // Retrieve ML scores safely - auto scores = dsCandidate.mlScores(); + auto jetTracks = jet.tracks_as(); + + const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f); + const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f); // thrust = λ_2^1 + const float mjet = computeJetMassFromTracksMass(jetTracks); - float s0 = (scores.size() > 0) ? scores[0] : -999.f; - float s1 = (scores.size() > 1) ? scores[1] : -999.f; - float s2 = (scores.size() > 2) ? scores[2] : -999.f; + const float R = jet.r() / 100.f; + const float girth = (lambda11 >= 0.f) ? (lambda11 * R) : -1.f; - distJetTable(axisDistance, + if (lambda11 >= 0.f) { + registry.fill(HIST("h_ds_jet_lambda11"), lambda11); + registry.fill(HIST("h2_dsjet_pt_lambda11"), jet.pt(), lambda11); + } + if (lambda12 >= 0.f) { + registry.fill(HIST("h_ds_jet_lambda12"), lambda12); + registry.fill(HIST("h2_dsjet_pt_lambda12"), jet.pt(), lambda12); + } + registry.fill(HIST("h_ds_jet_mass"), mjet); + registry.fill(HIST("h2_dsjet_pt_mass"), jet.pt(), mjet); + + if (girth >= 0.f) { + registry.fill(HIST("h_ds_jet_girth"), girth); + registry.fill(HIST("h2_dsjet_pt_girth"), jet.pt(), girth); + } + + if (doExtraAngularity) { + const float lambdaExtra = computeLambda(jet, jetTracks, alpha.value, kappa.value); + if (lambdaExtra >= 0.f) { + registry.fill(HIST("h_ds_jet_lambda_extra"), lambdaExtra); + registry.fill(HIST("h2_dsjet_pt_lambda_extra"), jet.pt(), lambdaExtra); + } + } + + auto scores = dsCandidate.mlScores(); + const float s0 = (scores.size() > 0) ? scores[0] : -999.f; + const float s1 = (scores.size() > 1) ? scores[1] : -999.f; + const float s2 = (scores.size() > 2) ? scores[2] : -999.f; + + distJetTable(static_cast(axisDistance), jet.pt(), jet.eta(), jet.phi(), - static_cast(jet.tracks_as().size()), + static_cast(jetTracks.size()), dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(), dsCandidate.m(), dsCandidate.y(), - s0, s1, s2); - break; // get out of candidates' loop after first HF particle is found in jet - } // end of DS candidates loop + s0, s1, s2, + mjet, girth, lambda12, lambda11); - } // end of jets loop - - } // end of process function + break; // only first Ds per jet + } + } + } PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false); };