1414// This task is used to reconstruct and analyse jets containing charged D_s
1515// mesons
1616//
17- // / \author Monalisa Melo <monalisa.melo@cern.ch>
17+ // / \author Monalisa Melo <monalisa.melo@cern.ch>, Universidade de São Paulo
1818//
1919
20+ #include " PWGHF/Core/DecayChannels.h"
2021#include " PWGJE/Core/JetDerivedDataUtilities.h"
2122#include " PWGJE/Core/JetUtilities.h"
2223#include " PWGJE/DataModel/Jet.h"
2324#include " PWGJE/DataModel/JetReducedData.h"
25+ #include " PWGJE/DataModel/JetSubtraction.h"
2426
25- #include < Framework/ASoA.h>
26- #include < Framework/AnalysisDataModel.h>
27- #include < Framework/AnalysisHelpers.h>
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"
2833#include < Framework/AnalysisTask.h>
2934#include < Framework/ConfigContext.h>
3035#include < Framework/Configurable.h>
31- #include < Framework/HistogramRegistry .h>
36+ #include < Framework/DataProcessorSpec .h>
3237#include < Framework/HistogramSpec.h>
3338#include < Framework/InitContext.h>
3439#include < Framework/runDataProcessing.h>
3540
36- #include < TVector3.h>
41+ #include " TVector3.h"
3742
3843#include < cmath>
3944#include < string>
@@ -100,13 +105,14 @@ struct JetDsSpecSubs {
100105 {" h_jet_eta" , " jet #eta;#eta_{jet};entries" , {HistType::kTH1F , {{100 , -1.0 , 1.0 }}}},
101106 {" h_jet_phi" , " jet #phi;#phi_{jet};entries" , {HistType::kTH1F , {{80 , -1.0 , 7 .}}}},
102107 {" h_collision_counter" , " # of collisions;" , {HistType::kTH1F , {{200 , 0 ., 200 .}}}},
103- {" h_jet_counter" , " ;# of D_{S} jets; " , {HistType::kTH1F , {{6 , 0 ., 3.0 }}}},
108+ {" h_jet_counter" , " ;type;counts " , {HistType::kTH1F , {{2 , 0 ., 2 . }}}},
104109 {" h_ds_jet_projection" , " ;z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}" , {HistType::kTH1F , {{1000 , 0 ., 2 .}}}},
105110 {" h_ds_jet_distance_vs_projection" , " ;#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}" , {HistType::kTH2F , {{1000 , 0 ., 1 .}, {1000 , 0 ., 2 .}}}},
106111 {" h_ds_jet_distance" , " ;#DeltaR_{D_{S},jet};dN/d(#DeltaR)" , {HistType::kTH1F , {{1000 , 0 ., 1 .}}}},
107112 {" h_ds_jet_pt" , " ;p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}" , {HistType::kTH1F , {{1000 , 0 ., 100 .}}}},
108113 {" h_ds_jet_eta" , " ;#eta_{D_{S} jet};entries" , {HistType::kTH1F , {{250 , -1 ., 1 .}}}},
109114 {" 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 }}}},
110116 {" h_ds_mass" , " ;m_{D_{S}} (GeV/c^{2});entries" , {HistType::kTH1F , {{1000 , 0 ., 6 .}}}},
111117 {" h_ds_eta" , " ;#eta_{D_{S}};entries" , {HistType::kTH1F , {{250 , -1 ., 1 .}}}},
112118 {" h_ds_phi" , " ;#phi_{D_{S}};entries" , {HistType::kTH1F , {{250 , -1 ., 7 .}}}},
@@ -118,8 +124,6 @@ struct JetDsSpecSubs {
118124 {" h2_dsjet_pt_lambda12" , " ;#it{p}_{T,jet} (GeV/#it{c});#lambda_{2}^{1}" , {HistType::kTH2F , {{100 , 0 ., 100 .}, {200 , 0 ., 1.0 }}}},
119125 {" 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 }}}},
120126 {" h2_dsjet_pt_girth" , " ;#it{p}_{T,jet} (GeV/#it{c});g" , {HistType::kTH2F , {{100 , 0 ., 100 .}, {200 , 0 ., 0.5 }}}},
121- {" h_ds_jet_lambda_extra" , " ;#lambda_{#alpha}^{#kappa};entries" , {HistType::kTH1F , {{200 , 0 ., 1.0 }}}},
122- {" h2_dsjet_pt_lambda_extra" , " ;#it{p}_{T,jet} (GeV/#it{c});#lambda_{#alpha}^{#kappa}" , {HistType::kTH2F , {{100 , 0 ., 100 .}, {200 , 0 ., 1.0 }}}},
123127 }};
124128
125129 Configurable<float > vertexZCut{" vertexZCut" , 10 .0f , " Accepted z-vertex range" };
@@ -129,12 +133,6 @@ struct JetDsSpecSubs {
129133 Configurable<std::string> eventSelections{" eventSelections" , " sel8" , " choose event selection" };
130134 Configurable<std::string> trackSelections{" trackSelections" , " globalTracks" , " set track selections" };
131135
132- // extra angularity knob
133- Configurable<float > kappa{" kappa" , 1 .0f , " angularity kappa" };
134- Configurable<float > alpha{" alpha" , 1 .0f , " angularity alpha" };
135-
136- bool doExtraAngularity = false ;
137-
138136 std::vector<int > eventSelectionBits;
139137 int trackSelection = -1 ;
140138
@@ -189,9 +187,9 @@ struct JetDsSpecSubs {
189187 eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits (static_cast <std::string>(eventSelections));
190188 trackSelection = jetderiveddatautilities::initialiseTrackSelection (static_cast <std::string>(trackSelections));
191189
192- const bool is11 = ( std::abs (kappa. value - 1 . f ) < 1e- 6f ) && ( std::abs (alpha. value - 1 . f ) < 1e- 6f );
193- const bool is12 = ( std::abs (kappa. value - 1 . f ) < 1e- 6f ) && ( std::abs (alpha. value - 2 . f ) < 1e- 6f );
194- doExtraAngularity = !(is11 || is12 );
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 " );
195193 }
196194
197195 Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100 .0f );
@@ -203,6 +201,7 @@ struct JetDsSpecSubs {
203201 if (!jetderiveddatautilities::selectCollision (collision, eventSelectionBits)) {
204202 return ;
205203 }
204+
206205 registry.fill (HIST (" h_collisions" ), 1.5 );
207206
208207 for (auto const & track : tracks) {
@@ -222,7 +221,7 @@ struct JetDsSpecSubs {
222221 if (!jetderiveddatautilities::selectCollision (collision, eventSelectionBits)) {
223222 return ;
224223 }
225- for (auto & jet : jets) {
224+ for (const auto & jet : jets) {
226225 registry.fill (HIST (" h_jet_pt" ), jet.pt ());
227226 registry.fill (HIST (" h_jet_eta" ), jet.eta ());
228227 registry.fill (HIST (" h_jet_phi" ), jet.phi ());
@@ -243,77 +242,104 @@ struct JetDsSpecSubs {
243242 registry.fill (HIST (" h_collision_counter" ), 3.0 );
244243
245244 for (const auto & jet : jets) {
245+
246246 registry.fill (HIST (" h_jet_counter" ), 0.5 );
247247
248+ bool hasDs = false ;
249+
248250 TVector3 jetVector (jet.px (), jet.py (), jet.pz ());
249251
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)
250263 for (const auto & dsCandidate : jet.candidates_as <aod::CandidatesDsData>()) {
264+
265+ hasDs = true ;
266+
251267 TVector3 dsVector (dsCandidate.px (), dsCandidate.py (), dsCandidate.pz ());
252268
269+ // zParallel defined as longitudinal momentum fraction along the jet axis
253270 const double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
254271 const double axisDistance = jetutilities::deltaR (jet, dsCandidate);
255272
273+ // --- Ds-level observables ---
256274 registry.fill (HIST (" h_ds_jet_projection" ), zParallel);
257275 registry.fill (HIST (" h_ds_jet_distance_vs_projection" ), axisDistance, zParallel);
258276 registry.fill (HIST (" h_ds_jet_distance" ), axisDistance);
259- registry.fill (HIST (" h_ds_jet_pt" ), jet.pt ());
260- registry.fill (HIST (" h_ds_jet_eta" ), jet.eta ());
261- registry.fill (HIST (" h_ds_jet_phi" ), jet.phi ());
277+
262278 registry.fill (HIST (" h_ds_mass" ), dsCandidate.m ());
263279 registry.fill (HIST (" h_ds_eta" ), dsCandidate.eta ());
264280 registry.fill (HIST (" h_ds_phi" ), dsCandidate.phi ());
265281
266- auto jetTracks = jet.tracks_as <aod::JetTracks>();
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) {
267311
268- const float lambda11 = computeLambda (jet, jetTracks, 1 .f , 1 .f );
269- const float lambda12 = computeLambda (jet, jetTracks, 2 .f , 1 .f ); // thrust = λ_2^1
270- const float mjet = computeJetMassFromTracksMass (jetTracks);
312+ registry.fill (HIST (" h_jet_counter" ), 1.5 );
271313
272- const float R = jet.r () / 100 .f ;
273- const float girth = (lambda11 >= 0 .f ) ? (lambda11 * R) : -1 .f ;
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 ());
274318
319+ // Jet substructure observables
275320 if (lambda11 >= 0 .f ) {
276321 registry.fill (HIST (" h_ds_jet_lambda11" ), lambda11);
277322 registry.fill (HIST (" h2_dsjet_pt_lambda11" ), jet.pt (), lambda11);
278323 }
324+
279325 if (lambda12 >= 0 .f ) {
280326 registry.fill (HIST (" h_ds_jet_lambda12" ), lambda12);
281327 registry.fill (HIST (" h2_dsjet_pt_lambda12" ), jet.pt (), lambda12);
282328 }
329+
283330 registry.fill (HIST (" h_ds_jet_mass" ), mjet);
284331 registry.fill (HIST (" h2_dsjet_pt_mass" ), jet.pt (), mjet);
285332
286333 if (girth >= 0 .f ) {
287334 registry.fill (HIST (" h_ds_jet_girth" ), girth);
288335 registry.fill (HIST (" h2_dsjet_pt_girth" ), jet.pt (), girth);
289336 }
290-
291- if (doExtraAngularity) {
292- const float lambdaExtra = computeLambda (jet, jetTracks, alpha.value , kappa.value );
293- if (lambdaExtra >= 0 .f ) {
294- registry.fill (HIST (" h_ds_jet_lambda_extra" ), lambdaExtra);
295- registry.fill (HIST (" h2_dsjet_pt_lambda_extra" ), jet.pt (), lambdaExtra);
296- }
297- }
298-
299- auto scores = dsCandidate.mlScores ();
300- const float s0 = (scores.size () > 0 ) ? scores[0 ] : -999 .f ;
301- const float s1 = (scores.size () > 1 ) ? scores[1 ] : -999 .f ;
302- const float s2 = (scores.size () > 2 ) ? scores[2 ] : -999 .f ;
303-
304- distJetTable (static_cast <float >(axisDistance),
305- jet.pt (), jet.eta (), jet.phi (),
306- static_cast <int >(jetTracks.size ()),
307- dsCandidate.pt (), dsCandidate.eta (), dsCandidate.phi (),
308- dsCandidate.m (), dsCandidate.y (),
309- s0, s1, s2,
310- mjet, girth, lambda12, lambda11);
311-
312- break ; // only first Ds per jet
313337 }
314338 }
315339 }
316340 PROCESS_SWITCH (JetDsSpecSubs, processDataChargedSubstructure, " charged HF jet substructure" , false );
317341};
318-
319- WorkflowSpec defineDataProcessing (ConfigContext const & cfgc) { return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc, TaskName{" jet-ds-spectrum-subs" })}; }
342+ WorkflowSpec defineDataProcessing (ConfigContext const & cfgc)
343+ {
344+ return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc)};
345+ }
0 commit comments