Skip to content

Commit 4322868

Browse files
monamelopmonamelop
committed
[PWGJE] Fix O2 linter: task name and file name
1 parent 95672e7 commit 4322868

File tree

1 file changed

+345
-0
lines changed

1 file changed

+345
-0
lines changed

‎PWGJE/Tasks/jetDsSpecSubs.cxx‎

Lines changed: 345 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,345 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// Jet substructure and spectrum task for D_s mesons
13+
//
14+
// This task is used to reconstruct and analyse jets containing charged D_s
15+
// mesons
16+
//
17+
/// \author Monalisa Melo <monalisa.melo@cern.ch>, Universidade de São Paulo
18+
//
19+
20+
#include "PWGHF/Core/DecayChannels.h"
21+
#include "PWGJE/Core/JetDerivedDataUtilities.h"
22+
#include "PWGJE/Core/JetUtilities.h"
23+
#include "PWGJE/DataModel/Jet.h"
24+
#include "PWGJE/DataModel/JetReducedData.h"
25+
#include "PWGJE/DataModel/JetSubtraction.h"
26+
27+
#include "Common/Core/RecoDecay.h"
28+
#include "Common/DataModel/TrackSelectionTables.h"
29+
30+
#include "Framework/ASoA.h"
31+
#include "Framework/AnalysisDataModel.h"
32+
#include "Framework/HistogramRegistry.h"
33+
#include <Framework/AnalysisTask.h>
34+
#include <Framework/ConfigContext.h>
35+
#include <Framework/Configurable.h>
36+
#include <Framework/DataProcessorSpec.h>
37+
#include <Framework/HistogramSpec.h>
38+
#include <Framework/InitContext.h>
39+
#include <Framework/runDataProcessing.h>
40+
41+
#include "TVector3.h"
42+
43+
#include <cmath>
44+
#include <string>
45+
#include <vector>
46+
47+
using namespace o2;
48+
using namespace o2::framework;
49+
using namespace o2::framework::expressions;
50+
51+
namespace o2::aod
52+
{
53+
namespace jet_distance
54+
{
55+
DECLARE_SOA_COLUMN(JetHfDist, jetHfDist, float);
56+
DECLARE_SOA_COLUMN(JetPt, jetPt, float);
57+
DECLARE_SOA_COLUMN(JetEta, jetEta, float);
58+
DECLARE_SOA_COLUMN(JetPhi, jetPhi, float);
59+
DECLARE_SOA_COLUMN(JetNConst, jetNConst, int);
60+
DECLARE_SOA_COLUMN(HfPt, hfPt, float);
61+
DECLARE_SOA_COLUMN(HfEta, hfEta, float);
62+
DECLARE_SOA_COLUMN(HfPhi, hfPhi, float);
63+
DECLARE_SOA_COLUMN(HfMass, hfMass, float);
64+
DECLARE_SOA_COLUMN(HfY, hfY, float);
65+
DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float);
66+
DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float);
67+
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
74+
} // namespace jet_distance
75+
76+
DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE",
77+
jet_distance::JetHfDist,
78+
jet_distance::JetPt,
79+
jet_distance::JetEta,
80+
jet_distance::JetPhi,
81+
jet_distance::JetNConst,
82+
jet_distance::HfPt,
83+
jet_distance::HfEta,
84+
jet_distance::HfPhi,
85+
jet_distance::HfMass,
86+
jet_distance::HfY,
87+
jet_distance::HfMlScore0,
88+
jet_distance::HfMlScore1,
89+
jet_distance::HfMlScore2,
90+
jet_distance::JetMass,
91+
jet_distance::JetGirth,
92+
jet_distance::JetThrust,
93+
jet_distance::JetLambda11);
94+
} // namespace o2::aod
95+
96+
struct JetDsSpecSubs {
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", ";type;counts", {HistType::kTH1F, {{2, 0., 2.}}}},
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+
{"hSparse_ds", ";m_{D_{S}};p_{T,D_{S}};z^{D_{S},jet}_{||};#DeltaR_{D_{S},jet}", {HistType::kTHnSparseF, {{60, 1.7, 2.1}, {60, 0., 100.}, {60, 0., 2.}, {60, 0., 1.0}}}},
116+
{"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});entries", {HistType::kTH1F, {{1000, 0., 6.}}}},
117+
{"h_ds_eta", ";#eta_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
118+
{"h_ds_phi", ";#phi_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
119+
{"h_ds_jet_mass", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}},
120+
{"h_ds_jet_lambda11", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
121+
{"h_ds_jet_lambda12", ";#lambda_{2}^{1} (thrust);entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
122+
{"h_ds_jet_girth", ";g (#equiv #lambda_{1}^{1}R);entries", {HistType::kTH1F, {{200, 0., 0.5}}}},
123+
{"h2_dsjet_pt_lambda11", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{1}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}},
124+
{"h2_dsjet_pt_lambda12", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{2}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}},
125+
{"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}}}},
126+
{"h2_dsjet_pt_girth", ";#it{p}_{T,jet} (GeV/#it{c});g", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 0.5}}}},
127+
}};
128+
129+
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
130+
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
131+
Configurable<float> jetR{"jetR", 0.4, "jet resolution parameter"};
132+
133+
Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
134+
Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};
135+
136+
std::vector<int> eventSelectionBits;
137+
int trackSelection = -1;
138+
139+
Produces<aod::JetDistanceTable> distJetTable;
140+
141+
template <typename JET, typename TRACKS>
142+
float computeLambda(JET const& jet, TRACKS const& tracks, float a, float k)
143+
{
144+
if (jet.pt() <= 0.f) {
145+
return -1.f;
146+
}
147+
float sum = 0.f;
148+
for (auto const& trk : tracks) {
149+
const float dr = jetutilities::deltaR(jet, trk);
150+
sum += std::pow(trk.pt(), k) * std::pow(dr, a);
151+
}
152+
const float R = jet.r() / 100.f;
153+
const float denom = std::pow(jet.pt(), k) * std::pow(R, a);
154+
if (denom <= 0.f) {
155+
return -1.f;
156+
}
157+
return sum / denom;
158+
}
159+
160+
template <typename TRACKS>
161+
float computeJetMassFromTracksMass(TRACKS const& tracks)
162+
{
163+
double sumPx = 0.0, sumPy = 0.0, sumPz = 0.0, sumE = 0.0;
164+
165+
for (auto const& trk : tracks) {
166+
const double pt = trk.pt();
167+
const double phi = trk.phi();
168+
const double eta = trk.eta();
169+
170+
const double px = pt * std::cos(phi);
171+
const double py = pt * std::sin(phi);
172+
const double pz = pt * std::sinh(eta);
173+
const double p = std::sqrt(px * px + py * py + pz * pz);
174+
175+
sumPx += px;
176+
sumPy += py;
177+
sumPz += pz;
178+
sumE += p; // massless
179+
}
180+
181+
const double m2 = sumE * sumE - (sumPx * sumPx + sumPy * sumPy + sumPz * sumPz);
182+
return (m2 > 0.0) ? static_cast<float>(std::sqrt(m2)) : 0.f;
183+
}
184+
185+
void init(o2::framework::InitContext&)
186+
{
187+
eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
188+
trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
189+
190+
auto h = registry.get<TH1>(HIST("h_jet_counter"));
191+
h->GetXaxis()->SetBinLabel(1, "All jets");
192+
h->GetXaxis()->SetBinLabel(2, "Ds-tagged jets");
193+
}
194+
195+
Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
196+
Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut;
197+
198+
void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks)
199+
{
200+
registry.fill(HIST("h_collisions"), 0.5);
201+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
202+
return;
203+
}
204+
205+
registry.fill(HIST("h_collisions"), 1.5);
206+
207+
for (auto const& track : tracks) {
208+
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
209+
continue;
210+
}
211+
registry.fill(HIST("h_track_pt"), track.pt());
212+
registry.fill(HIST("h_track_eta"), track.eta());
213+
registry.fill(HIST("h_track_phi"), track.phi());
214+
}
215+
}
216+
PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "process JE collisions", false);
217+
218+
void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision,
219+
soa::Filtered<aod::ChargedJets> const& jets)
220+
{
221+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
222+
return;
223+
}
224+
for (const auto& jet : jets) {
225+
registry.fill(HIST("h_jet_pt"), jet.pt());
226+
registry.fill(HIST("h_jet_eta"), jet.eta());
227+
registry.fill(HIST("h_jet_phi"), jet.phi());
228+
}
229+
}
230+
PROCESS_SWITCH(JetDsSpecSubs, processDataCharged, "charged jets in data", false);
231+
232+
void processDataChargedSubstructure(aod::JetCollision const& collision,
233+
soa::Join<aod::DsChargedJets, aod::DsChargedJetConstituents> const& jets,
234+
aod::CandidatesDsData const&,
235+
aod::JetTracks const&)
236+
{
237+
registry.fill(HIST("h_collision_counter"), 2.0);
238+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) ||
239+
!(std::abs(collision.posZ()) < vertexZCut)) {
240+
return;
241+
}
242+
registry.fill(HIST("h_collision_counter"), 3.0);
243+
244+
for (const auto& jet : jets) {
245+
246+
registry.fill(HIST("h_jet_counter"), 0.5);
247+
248+
bool hasDs = false;
249+
250+
TVector3 jetVector(jet.px(), jet.py(), jet.pz());
251+
252+
// Compute jet-level quantities once (independent of Ds candidates)
253+
auto jetTracks = jet.tracks_as<aod::JetTracks>();
254+
255+
const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f);
256+
const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f);
257+
const float mjet = computeJetMassFromTracksMass(jetTracks);
258+
259+
const float R = jet.r() / 100.f;
260+
const float girth = (lambda11 >= 0.f) ? (lambda11 * R) : -1.f;
261+
262+
// Loop over Ds candidates (particle level)
263+
for (const auto& dsCandidate : jet.candidates_as<aod::CandidatesDsData>()) {
264+
265+
hasDs = true;
266+
267+
TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz());
268+
269+
// zParallel defined as longitudinal momentum fraction along the jet axis
270+
const double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
271+
const double axisDistance = jetutilities::deltaR(jet, dsCandidate);
272+
273+
// --- Ds-level observables ---
274+
registry.fill(HIST("h_ds_jet_projection"), zParallel);
275+
registry.fill(HIST("h_ds_jet_distance_vs_projection"), axisDistance, zParallel);
276+
registry.fill(HIST("h_ds_jet_distance"), axisDistance);
277+
278+
registry.fill(HIST("h_ds_mass"), dsCandidate.m());
279+
registry.fill(HIST("h_ds_eta"), dsCandidate.eta());
280+
registry.fill(HIST("h_ds_phi"), dsCandidate.phi());
281+
282+
const float mass = dsCandidate.m();
283+
const float pt = dsCandidate.pt();
284+
const float z = zParallel;
285+
const float dR = axisDistance;
286+
287+
// Main THnSparse: invariant mass, pT, z, and ΔR
288+
registry.fill(HIST("hSparse_ds"), mass, pt, z, dR);
289+
290+
// --- output table ---
291+
auto scores = dsCandidate.mlScores();
292+
constexpr int kScore0 = 0;
293+
constexpr int kScore1 = 1;
294+
constexpr int kScore2 = 2;
295+
296+
const float s0 = (scores.size() > kScore0) ? scores[kScore0] : -999.f;
297+
const float s1 = (scores.size() > kScore1) ? scores[kScore1] : -999.f;
298+
const float s2 = (scores.size() > kScore2) ? scores[kScore2] : -999.f;
299+
300+
distJetTable(static_cast<float>(axisDistance),
301+
jet.pt(), jet.eta(), jet.phi(),
302+
static_cast<int>(jetTracks.size()),
303+
dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(),
304+
dsCandidate.m(), dsCandidate.y(),
305+
s0, s1, s2,
306+
mjet, girth, lambda12, lambda11);
307+
}
308+
309+
// Jet-level quantities (filled once per jet containing at least one Ds)
310+
if (hasDs) {
311+
312+
registry.fill(HIST("h_jet_counter"), 1.5);
313+
314+
// Jet properties
315+
registry.fill(HIST("h_ds_jet_pt"), jet.pt());
316+
registry.fill(HIST("h_ds_jet_eta"), jet.eta());
317+
registry.fill(HIST("h_ds_jet_phi"), jet.phi());
318+
319+
// Jet substructure observables
320+
if (lambda11 >= 0.f) {
321+
registry.fill(HIST("h_ds_jet_lambda11"), lambda11);
322+
registry.fill(HIST("h2_dsjet_pt_lambda11"), jet.pt(), lambda11);
323+
}
324+
325+
if (lambda12 >= 0.f) {
326+
registry.fill(HIST("h_ds_jet_lambda12"), lambda12);
327+
registry.fill(HIST("h2_dsjet_pt_lambda12"), jet.pt(), lambda12);
328+
}
329+
330+
registry.fill(HIST("h_ds_jet_mass"), mjet);
331+
registry.fill(HIST("h2_dsjet_pt_mass"), jet.pt(), mjet);
332+
333+
if (girth >= 0.f) {
334+
registry.fill(HIST("h_ds_jet_girth"), girth);
335+
registry.fill(HIST("h2_dsjet_pt_girth"), jet.pt(), girth);
336+
}
337+
}
338+
}
339+
}
340+
PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false);
341+
};
342+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
343+
{
344+
return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc)};
345+
}

0 commit comments

Comments
 (0)