Skip to content

[PWGJE] Add new jet substructure observable #15244

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nzardosh merged 1 commit into from
Mar 4, 2026
+151 −52
Merged

[PWGJE] Add new jet substructure observable #15244

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
203 changes: 151 additions & 52 deletions PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@
#include "Framework/ASoA.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/HistogramRegistry.h"
#include <CommonConstants/MathConstants.h>
#include <Framework/AnalysisTask.h>
#include <Framework/ConfigContext.h>
#include <Framework/Configurable.h>
Expand All @@ -42,7 +41,6 @@
#include "TVector3.h"

#include <cmath>
#include <cstdlib>
#include <string>
#include <vector>

Expand All @@ -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",
Expand All @@ -82,59 +86,130 @@ 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<float> 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<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
Configurable<float> jetR{"jetR", 0.4, "jet resolution parameter"};

Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};

// extra angularity knob
Configurable<float> kappa{"kappa", 1.0f, "angularity kappa"};
Configurable<float> alpha{"alpha", 1.0f, "angularity alpha"};

bool doExtraAngularity = false;

std::vector<int> eventSelectionBits;
int trackSelection = -1;

Produces<aod::JetDistanceTable> distJetTable;

template <typename JET, typename TRACKS>
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 <typename TRACKS>
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<float>(std::sqrt(m2)) : 0.f;
}

void init(o2::framework::InitContext&)
{
eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(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);
Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut;

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;
Expand All @@ -146,12 +221,12 @@ struct JetDsSpecSubs {
}
PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "process JE collisions", false);

void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision, soa::Filtered<aod::ChargedJets> const& jets)
void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision,
soa::Filtered<aod::ChargedJets> 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());
Expand All @@ -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<aod::CandidatesDsData>()) {

// 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);
Expand All @@ -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<aod::JetTracks>();

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<float>(axisDistance),
jet.pt(), jet.eta(), jet.phi(),
static_cast<int>(jet.tracks_as<aod::JetTracks>().size()),
static_cast<int>(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);
};

Expand Down
Loading