[PWGJE] Add new plots for Ds jets

PWGJE/Tasks/CMakeLists.txt

@@ -1,4 +1,4 @@
 # Copyright 2019-2026 CERN and copyright holders of ALICE O2.
 # See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 # All rights not expressly granted are reserved.
 #
@@ -28,19 +28,19 @@
 )
 endif()
 
 o2physics_add_dpl_workflow(emc-cellmonitor
                     SOURCES emcCellMonitor.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(emc-clustermonitor
                     SOURCES emcClusterMonitor.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(emc-eventselection-qa
                     SOURCES emcEventSelectionQA.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(emc-vertexselection-qa
                     SOURCES emcVertexSelectionQA.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
@@ -48,11 +48,11 @@
                     SOURCES emcalGammaGammaBcWise.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(emc-pi0-energyscale-calib
                     SOURCES emcalPi0EnergyScaleCalib.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(emc-tmmonitor
                     SOURCES emcTmMonitor.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
@@ -60,7 +60,7 @@
                     SOURCES mcGeneratorStudies.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 o2physics_add_dpl_workflow(photon-isolation-qa
                     SOURCES photonIsolationQA.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
@@ -395,8 +395,8 @@
                         SOURCES substructureDebug.cxx
                         PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
                         COMPONENT_NAME Analysis)
-    o2physics_add_dpl_workflow(jet-ds-spectrum-subs
-                        SOURCES jetDsSpectrumAndSubstructure.cxx
+    o2physics_add_dpl_workflow(jet-ds-spec-subs
+                        SOURCES jetDsSpecSubs.cxx
                         PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
                         COMPONENT_NAME Analysis)
     o2physics_add_dpl_workflow(jet-d0-ang-substructure

PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx → PWGJE/Tasks/jetDsSpecSubs.cxx

@@ -14,26 +14,31 @@
 // This task is used to reconstruct and analyse jets containing charged D_s
 // mesons
 //
-/// \author Monalisa Melo <monalisa.melo@cern.ch>
+/// \author Monalisa Melo <monalisa.melo@cern.ch>, Universidade de São Paulo
 //
 
+#include "PWGHF/Core/DecayChannels.h"
 #include "PWGJE/Core/JetDerivedDataUtilities.h"
 #include "PWGJE/Core/JetUtilities.h"
 #include "PWGJE/DataModel/Jet.h"
 #include "PWGJE/DataModel/JetReducedData.h"
+#include "PWGJE/DataModel/JetSubtraction.h"
 
-#include <Framework/ASoA.h>
-#include <Framework/AnalysisDataModel.h>
-#include <Framework/AnalysisHelpers.h>
+#include "Common/Core/RecoDecay.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+#include "Framework/ASoA.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/HistogramRegistry.h"
 #include <Framework/AnalysisTask.h>
 #include <Framework/ConfigContext.h>
 #include <Framework/Configurable.h>
-#include <Framework/HistogramRegistry.h>
+#include <Framework/DataProcessorSpec.h>
 #include <Framework/HistogramSpec.h>
 #include <Framework/InitContext.h>
 #include <Framework/runDataProcessing.h>
 
-#include <TVector3.h>
+#include "TVector3.h"
 
 #include <cmath>
 #include <string>
@@ -100,13 +105,14 @@ struct JetDsSpecSubs {
       {"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_jet_counter", ";type;counts", {HistType::kTH1F, {{2, 0., 2.}}}},
       {"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.}}}},
+      {"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}}}},
       {"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.}}}},
@@ -118,8 +124,6 @@ struct JetDsSpecSubs {
       {"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"};
@@ -129,12 +133,6 @@ struct JetDsSpecSubs {
   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;
 
@@ -189,9 +187,9 @@ struct JetDsSpecSubs {
     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);
+    auto h = registry.get<TH1>(HIST("h_jet_counter"));
+    h->GetXaxis()->SetBinLabel(1, "All jets");
+    h->GetXaxis()->SetBinLabel(2, "Ds-tagged jets");
   }
 
   Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
@@ -203,6 +201,7 @@ struct JetDsSpecSubs {
     if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
       return;
     }
+
     registry.fill(HIST("h_collisions"), 1.5);
 
     for (auto const& track : tracks) {
@@ -222,7 +221,7 @@ struct JetDsSpecSubs {
     if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
       return;
     }
-    for (auto& jet : jets) {
+    for (const auto& jet : jets) {
       registry.fill(HIST("h_jet_pt"), jet.pt());
       registry.fill(HIST("h_jet_eta"), jet.eta());
       registry.fill(HIST("h_jet_phi"), jet.phi());
@@ -243,77 +242,104 @@ struct JetDsSpecSubs {
     registry.fill(HIST("h_collision_counter"), 3.0);
 
     for (const auto& jet : jets) {
+
       registry.fill(HIST("h_jet_counter"), 0.5);
 
+      bool hasDs = false;
+
       TVector3 jetVector(jet.px(), jet.py(), jet.pz());
 
+      // Compute jet-level quantities once (independent of Ds candidates)
+      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);
+      const float mjet = computeJetMassFromTracksMass(jetTracks);
+
+      const float R = jet.r() / 100.f;
+      const float girth = (lambda11 >= 0.f) ? (lambda11 * R) : -1.f;
+
+      // Loop over Ds candidates (particle level)
       for (const auto& dsCandidate : jet.candidates_as<aod::CandidatesDsData>()) {
+
+        hasDs = true;
+
         TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz());
 
+        // zParallel defined as longitudinal momentum fraction along the jet axis
         const double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
         const double axisDistance = jetutilities::deltaR(jet, dsCandidate);
 
+        // --- Ds-level observables ---
         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);
-        registry.fill(HIST("h_ds_jet_pt"), jet.pt());
-        registry.fill(HIST("h_ds_jet_eta"), jet.eta());
-        registry.fill(HIST("h_ds_jet_phi"), jet.phi());
+
         registry.fill(HIST("h_ds_mass"), dsCandidate.m());
         registry.fill(HIST("h_ds_eta"), dsCandidate.eta());
         registry.fill(HIST("h_ds_phi"), dsCandidate.phi());
 
-        auto jetTracks = jet.tracks_as<aod::JetTracks>();
+        const float mass = dsCandidate.m();
+        const float pt = dsCandidate.pt();
+        const float z = zParallel;
+        const float dR = axisDistance;
+
+        // Main THnSparse: invariant mass, pT, z, and ΔR
+        registry.fill(HIST("hSparse_ds"), mass, pt, z, dR);
+
+        // --- output table ---
+        auto scores = dsCandidate.mlScores();
+        constexpr int kScore0 = 0;
+        constexpr int kScore1 = 1;
+        constexpr int kScore2 = 2;
+
+        const float s0 = (scores.size() > kScore0) ? scores[kScore0] : -999.f;
+        const float s1 = (scores.size() > kScore1) ? scores[kScore1] : -999.f;
+        const float s2 = (scores.size() > kScore2) ? scores[kScore2] : -999.f;
+
+        distJetTable(static_cast<float>(axisDistance),
+                     jet.pt(), jet.eta(), jet.phi(),
+                     static_cast<int>(jetTracks.size()),
+                     dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(),
+                     dsCandidate.m(), dsCandidate.y(),
+                     s0, s1, s2,
+                     mjet, girth, lambda12, lambda11);
+      }
+
+      // Jet-level quantities (filled once per jet containing at least one Ds)
+      if (hasDs) {
 
-        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);
+        registry.fill(HIST("h_jet_counter"), 1.5);
 
-        const float R = jet.r() / 100.f;
-        const float girth = (lambda11 >= 0.f) ? (lambda11 * R) : -1.f;
+        // Jet properties
+        registry.fill(HIST("h_ds_jet_pt"), jet.pt());
+        registry.fill(HIST("h_ds_jet_eta"), jet.eta());
+        registry.fill(HIST("h_ds_jet_phi"), jet.phi());
 
+        // Jet substructure observables
         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>(jetTracks.size()),
-                     dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(),
-                     dsCandidate.m(), dsCandidate.y(),
-                     s0, s1, s2,
-                     mjet, girth, lambda12, lambda11);
-
-        break; // only first Ds per jet
       }
     }
   }
   PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false);
 };
-
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc, TaskName{"jet-ds-spectrum-subs"})}; }
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc)};
+}