Skip to content

Commit 4554a2e

Browse files
monamelopmonamelop
authored
[PWGJE] Add new jet substructure observable (#15244)
1 parent 4cbfcbd commit 4554a2e

File tree

1 file changed

+151
-52
lines changed

1 file changed

+151
-52
lines changed

‎PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx‎

Lines changed: 151 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,6 @@
3030
#include "Framework/ASoA.h"
3131
#include "Framework/AnalysisDataModel.h"
3232
#include "Framework/HistogramRegistry.h"
33-
#include <CommonConstants/MathConstants.h>
3433
#include <Framework/AnalysisTask.h>
3534
#include <Framework/ConfigContext.h>
3635
#include <Framework/Configurable.h>
@@ -42,7 +41,6 @@
4241
#include "TVector3.h"
4342

4443
#include <cmath>
45-
#include <cstdlib>
4644
#include <string>
4745
#include <vector>
4846

@@ -67,6 +65,12 @@ DECLARE_SOA_COLUMN(HfY, hfY, float);
6765
DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float);
6866
DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float);
6967
DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float);
68+
69+
// extra
70+
DECLARE_SOA_COLUMN(JetMass, jetMass, float);
71+
DECLARE_SOA_COLUMN(JetGirth, jetGirth, float);
72+
DECLARE_SOA_COLUMN(JetThrust, jetThrust, float); // lambda_2^1
73+
DECLARE_SOA_COLUMN(JetLambda11, jetLambda11, float); // lambda_1^1
7074
} // namespace jet_distance
7175

7276
DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE",
@@ -82,59 +86,130 @@ DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE",
8286
jet_distance::HfY,
8387
jet_distance::HfMlScore0,
8488
jet_distance::HfMlScore1,
85-
jet_distance::HfMlScore2);
89+
jet_distance::HfMlScore2,
90+
jet_distance::JetMass,
91+
jet_distance::JetGirth,
92+
jet_distance::JetThrust,
93+
jet_distance::JetLambda11);
8694
} // namespace o2::aod
8795

8896
struct JetDsSpecSubs {
89-
HistogramRegistry registry{"registry",
90-
{{"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}},
91-
{"h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
92-
{"h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
93-
{"h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
94-
{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
95-
{"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
96-
{"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
97-
{"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{200, 0., 200.}}}},
98-
{"h_jet_counter", ";# of D_{S} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}},
99-
{"h_ds_jet_projection", ";z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}", {HistType::kTH1F, {{1000, 0., 2.}}}},
100-
{"h_ds_jet_distance_vs_projection", ";#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}", {HistType::kTH2F, {{1000, 0., 1.}, {1000, 0., 2.}}}},
101-
{"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 1.}}}},
102-
{"h_ds_jet_pt", ";p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}", {HistType::kTH1F, {{1000, 0., 100.}}}},
103-
{"h_ds_jet_eta", ";#eta_{T,D_{S} jet};dN/d#eta_{D_{S} jet}", {HistType::kTH1F, {{250, -1., 1.}}}},
104-
{"h_ds_jet_phi", ";#phi_{T,D_{S} jet};dN/d#phi_{D_{S} jet}", {HistType::kTH1F, {{250, -1., 7.}}}},
105-
{"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});dN/dm_{D_{S}}", {HistType::kTH1F, {{1000, 0., 6.}}}},
106-
{"h_ds_eta", ";#eta_{D_{S}} (GeV/c^{2});dN/d#eta_{D_{S}}", {HistType::kTH1F, {{250, -1., 1.}}}},
107-
{"h_ds_phi", ";#phi_{D_{S}} (GeV/c^{2});dN/d#phi_{D_{S}}", {HistType::kTH1F, {{250, -1., 7.}}}}}};
108-
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
97+
HistogramRegistry registry{
98+
"registry",
99+
{
100+
{"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}},
101+
{"h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
102+
{"h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
103+
{"h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
104+
{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
105+
{"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
106+
{"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
107+
{"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{200, 0., 200.}}}},
108+
{"h_jet_counter", ";# of D_{S} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}},
109+
{"h_ds_jet_projection", ";z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}", {HistType::kTH1F, {{1000, 0., 2.}}}},
110+
{"h_ds_jet_distance_vs_projection", ";#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}", {HistType::kTH2F, {{1000, 0., 1.}, {1000, 0., 2.}}}},
111+
{"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 1.}}}},
112+
{"h_ds_jet_pt", ";p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}", {HistType::kTH1F, {{1000, 0., 100.}}}},
113+
{"h_ds_jet_eta", ";#eta_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
114+
{"h_ds_jet_phi", ";#phi_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
115+
{"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});entries", {HistType::kTH1F, {{1000, 0., 6.}}}},
116+
{"h_ds_eta", ";#eta_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
117+
{"h_ds_phi", ";#phi_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
118+
{"h_ds_jet_mass", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}},
119+
{"h_ds_jet_lambda11", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
120+
{"h_ds_jet_lambda12", ";#lambda_{2}^{1} (thrust);entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
121+
{"h_ds_jet_girth", ";g (#equiv #lambda_{1}^{1}R);entries", {HistType::kTH1F, {{200, 0., 0.5}}}},
122+
{"h2_dsjet_pt_lambda11", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{1}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}},
123+
{"h2_dsjet_pt_lambda12", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{2}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}},
124+
{"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}}}},
125+
{"h2_dsjet_pt_girth", ";#it{p}_{T,jet} (GeV/#it{c});g", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 0.5}}}},
126+
{"h_ds_jet_lambda_extra", ";#lambda_{#alpha}^{#kappa};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
127+
{"h2_dsjet_pt_lambda_extra", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{#alpha}^{#kappa}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}},
128+
}};
109129

130+
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
110131
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
111132
Configurable<float> jetR{"jetR", 0.4, "jet resolution parameter"};
112133

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

137+
// extra angularity knob
138+
Configurable<float> kappa{"kappa", 1.0f, "angularity kappa"};
139+
Configurable<float> alpha{"alpha", 1.0f, "angularity alpha"};
140+
141+
bool doExtraAngularity = false;
142+
116143
std::vector<int> eventSelectionBits;
117144
int trackSelection = -1;
118145

119146
Produces<aod::JetDistanceTable> distJetTable;
120147

148+
template <typename JET, typename TRACKS>
149+
float computeLambda(JET const& jet, TRACKS const& tracks, float a, float k)
150+
{
151+
if (jet.pt() <= 0.f) {
152+
return -1.f;
153+
}
154+
float sum = 0.f;
155+
for (auto const& trk : tracks) {
156+
const float dr = jetutilities::deltaR(jet, trk);
157+
sum += std::pow(trk.pt(), k) * std::pow(dr, a);
158+
}
159+
const float R = jet.r() / 100.f;
160+
const float denom = std::pow(jet.pt(), k) * std::pow(R, a);
161+
if (denom <= 0.f) {
162+
return -1.f;
163+
}
164+
return sum / denom;
165+
}
166+
167+
template <typename TRACKS>
168+
float computeJetMassFromTracksMass(TRACKS const& tracks)
169+
{
170+
double sumPx = 0.0, sumPy = 0.0, sumPz = 0.0, sumE = 0.0;
171+
172+
for (auto const& trk : tracks) {
173+
const double pt = trk.pt();
174+
const double phi = trk.phi();
175+
const double eta = trk.eta();
176+
177+
const double px = pt * std::cos(phi);
178+
const double py = pt * std::sin(phi);
179+
const double pz = pt * std::sinh(eta);
180+
const double p = std::sqrt(px * px + py * py + pz * pz);
181+
182+
sumPx += px;
183+
sumPy += py;
184+
sumPz += pz;
185+
sumE += p; // massless
186+
}
187+
188+
const double m2 = sumE * sumE - (sumPx * sumPx + sumPy * sumPy + sumPz * sumPz);
189+
return (m2 > 0.0) ? static_cast<float>(std::sqrt(m2)) : 0.f;
190+
}
191+
121192
void init(o2::framework::InitContext&)
122193
{
123194
eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
124195
trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
196+
197+
const bool is11 = (std::abs(kappa.value - 1.f) < 1e-6f) && (std::abs(alpha.value - 1.f) < 1e-6f);
198+
const bool is12 = (std::abs(kappa.value - 1.f) < 1e-6f) && (std::abs(alpha.value - 2.f) < 1e-6f);
199+
doExtraAngularity = !(is11 || is12);
125200
}
126201

127202
Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
128203
Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut;
129204

130205
void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks)
131206
{
132-
133207
registry.fill(HIST("h_collisions"), 0.5);
134208
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
135209
return;
136210
}
137211
registry.fill(HIST("h_collisions"), 1.5);
212+
138213
for (auto const& track : tracks) {
139214
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
140215
continue;
@@ -146,12 +221,12 @@ struct JetDsSpecSubs {
146221
}
147222
PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "process JE collisions", false);
148223

149-
void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision, soa::Filtered<aod::ChargedJets> const& jets)
224+
void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision,
225+
soa::Filtered<aod::ChargedJets> const& jets)
150226
{
151227
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
152228
return;
153229
}
154-
// jets -> charged jet
155230
for (auto& jet : jets) {
156231
registry.fill(HIST("h_jet_pt"), jet.pt());
157232
registry.fill(HIST("h_jet_eta"), jet.eta());
@@ -165,31 +240,24 @@ struct JetDsSpecSubs {
165240
aod::CandidatesDsData const&,
166241
aod::JetTracks const&)
167242
{
168-
// apply event selection and fill histograms for sanity check
169243
registry.fill(HIST("h_collision_counter"), 2.0);
170-
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) {
244+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) ||
245+
!(std::abs(collision.posZ()) < vertexZCut)) {
171246
return;
172247
}
173248
registry.fill(HIST("h_collision_counter"), 3.0);
174-
// jets -> charged jet with Ds
249+
175250
for (const auto& jet : jets) {
176-
// number of charged jets with Ds
177251
registry.fill(HIST("h_jet_counter"), 0.5);
178-
// obtaining jet 3-vector
252+
179253
TVector3 jetVector(jet.px(), jet.py(), jet.pz());
180254

181255
for (const auto& dsCandidate : jet.candidates_as<aod::CandidatesDsData>()) {
182-
183-
// obtaining jet 3-vector
184256
TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz());
185257

186-
// calculating fraction of the jet momentum carried by the Ds along the direction of the jet axis
187-
double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
188-
189-
// calculating angular distance in eta-phi plane
190-
double axisDistance = jetutilities::deltaR(jet, dsCandidate);
258+
const double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
259+
const double axisDistance = jetutilities::deltaR(jet, dsCandidate);
191260

192-
// filling histograms
193261
registry.fill(HIST("h_ds_jet_projection"), zParallel);
194262
registry.fill(HIST("h_ds_jet_distance_vs_projection"), axisDistance, zParallel);
195263
registry.fill(HIST("h_ds_jet_distance"), axisDistance);
@@ -200,25 +268,56 @@ struct JetDsSpecSubs {
200268
registry.fill(HIST("h_ds_eta"), dsCandidate.eta());
201269
registry.fill(HIST("h_ds_phi"), dsCandidate.phi());
202270

203-
// Retrieve ML scores safely
204-
auto scores = dsCandidate.mlScores();
271+
auto jetTracks = jet.tracks_as<aod::JetTracks>();
272+
273+
const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f);
274+
const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f); // thrust = λ_2^1
275+
const float mjet = computeJetMassFromTracksMass(jetTracks);
205276

206-
float s0 = (scores.size() > 0) ? scores[0] : -999.f;
207-
float s1 = (scores.size() > 1) ? scores[1] : -999.f;
208-
float s2 = (scores.size() > 2) ? scores[2] : -999.f;
277+
const float R = jet.r() / 100.f;
278+
const float girth = (lambda11 >= 0.f) ? (lambda11 * R) : -1.f;
209279

210-
distJetTable(axisDistance,
280+
if (lambda11 >= 0.f) {
281+
registry.fill(HIST("h_ds_jet_lambda11"), lambda11);
282+
registry.fill(HIST("h2_dsjet_pt_lambda11"), jet.pt(), lambda11);
283+
}
284+
if (lambda12 >= 0.f) {
285+
registry.fill(HIST("h_ds_jet_lambda12"), lambda12);
286+
registry.fill(HIST("h2_dsjet_pt_lambda12"), jet.pt(), lambda12);
287+
}
288+
registry.fill(HIST("h_ds_jet_mass"), mjet);
289+
registry.fill(HIST("h2_dsjet_pt_mass"), jet.pt(), mjet);
290+
291+
if (girth >= 0.f) {
292+
registry.fill(HIST("h_ds_jet_girth"), girth);
293+
registry.fill(HIST("h2_dsjet_pt_girth"), jet.pt(), girth);
294+
}
295+
296+
if (doExtraAngularity) {
297+
const float lambdaExtra = computeLambda(jet, jetTracks, alpha.value, kappa.value);
298+
if (lambdaExtra >= 0.f) {
299+
registry.fill(HIST("h_ds_jet_lambda_extra"), lambdaExtra);
300+
registry.fill(HIST("h2_dsjet_pt_lambda_extra"), jet.pt(), lambdaExtra);
301+
}
302+
}
303+
304+
auto scores = dsCandidate.mlScores();
305+
const float s0 = (scores.size() > 0) ? scores[0] : -999.f;
306+
const float s1 = (scores.size() > 1) ? scores[1] : -999.f;
307+
const float s2 = (scores.size() > 2) ? scores[2] : -999.f;
308+
309+
distJetTable(static_cast<float>(axisDistance),
211310
jet.pt(), jet.eta(), jet.phi(),
212-
static_cast<int>(jet.tracks_as<aod::JetTracks>().size()),
311+
static_cast<int>(jetTracks.size()),
213312
dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(),
214313
dsCandidate.m(), dsCandidate.y(),
215-
s0, s1, s2);
216-
break; // get out of candidates' loop after first HF particle is found in jet
217-
} // end of DS candidates loop
314+
s0, s1, s2,
315+
mjet, girth, lambda12, lambda11);
218316

219-
} // end of jets loop
220-
221-
} // end of process function
317+
break; // only first Ds per jet
318+
}
319+
}
320+
}
222321
PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false);
223322
};
224323

0 commit comments

Comments
 (0)