From 92931f48e2f394113ac699b6e456ac93bb0fb4c3 Mon Sep 17 00:00:00 2001 From: Preeti Dhankher Date: Mon, 9 Mar 2026 13:26:54 +0100 Subject: [PATCH 1/5] Add new D0 substructure task --- PWGJE/Tasks/CMakeLists.txt | 4 + PWGJE/Tasks/jetD0Substructure.cxx | 266 ++++++++++++++++++++++++++++++ 2 files changed, 270 insertions(+) create mode 100644 PWGJE/Tasks/jetD0Substructure.cxx diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index f05c735ef21..d78d84fbafb 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -389,6 +389,10 @@ if(FastJet_FOUND) SOURCES jetDsSpectrumAndSubstructure.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(jet-d0-substructure + SOURCES jetD0Substructure.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(bjet-cent-mult SOURCES bjetCentMult.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore diff --git a/PWGJE/Tasks/jetD0Substructure.cxx b/PWGJE/Tasks/jetD0Substructure.cxx new file mode 100644 index 00000000000..18d2dbd6328 --- /dev/null +++ b/PWGJE/Tasks/jetD0Substructure.cxx @@ -0,0 +1,266 @@ +// Copyright 2019-2020 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. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// \file jetD0Substructure.cxx +// +// \brief Analysis task for the reconstruction and study of charged jets +// containing D_0 mesons in pp collisions. +// \inherited from D0 fragmentation and Ds +// \P. Dhankher + + +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetSubtraction.h" +#include "PWGJE/DataModel/JetReducedData.h" + +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGHF/Core/DecayChannels.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" +#include +#include "Framework/HistogramRegistry.h" +#include +#include +#include +#include +#include +#include +#include + + +#include "TVector3.h" +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +// Definition of a custom AOD table to store jet–D0 quantities +namespace o2::aod +{ +namespace jet_obj +{ +// Jet-related quantities +DECLARE_SOA_COLUMN(JetHfDist, jetHfDist, float); +DECLARE_SOA_COLUMN(JetPt, jetPt, float); +DECLARE_SOA_COLUMN(JetEta, jetEta, float); +DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); +DECLARE_SOA_COLUMN(JetNConst, jetNConst, int); +DECLARE_SOA_COLUMN(JetAng, jetAng, float); +// D0 candidate quantities +DECLARE_SOA_COLUMN(HfPt, hfPt, float); +DECLARE_SOA_COLUMN(HfEta, hfEta, float); +DECLARE_SOA_COLUMN(HfPhi, hfPhi, float); +DECLARE_SOA_COLUMN(HfMass, hfMass, float); +DECLARE_SOA_COLUMN(HfY, hfY, float); +// ML scores +DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float); +DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float); +DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float); +} // namespace jet_obj +// AOD table definition +DECLARE_SOA_TABLE(JetObjTable, "AOD", "JETOBJTABLE", + jet_obj::JetHfDist, + jet_obj::JetPt, + jet_obj::JetEta, + jet_obj::JetPhi, + jet_obj::JetNConst, + jet_obj::JetAng, + jet_obj::HfPt, + jet_obj::HfEta, + jet_obj::HfPhi, + jet_obj::HfMass, + jet_obj::HfY, + jet_obj::HfMlScore0, + jet_obj::HfMlScore1, + jet_obj::HfMlScore2); +} +struct JetD0Substructure { + /** + * Histogram registry + * + * Contains: + * - Event and track histograms + * - Jet kinematic distributions + * - D0–jet substructure observables + */ + 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, {{2, 0., 2.}}}}, + {"h_jet_counter", ";# of D^{0} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, + {"h_d0_jet_projection", ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_jet_distance_vs_projection", ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}}, + {"h_d0_jet_distance", ";#DeltaR_{D^{0},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_jet_pt", ";p_{T,D^{0} jet};dN/dp_{T,D^{0} jet}", {HistType::kTH1F, {{200, 0., 10.}}}}, + {"h_d0_jet_eta", ";#eta_{T,D^{0} jet};dN/d#eta_{D^{0} jet}", {HistType::kTH1F, {{250, -5., 5.}}}}, + {"h_d0_jet_phi", ";#phi_{T,D^{0} jet};dN/d#phi_{D^{0} jet}", {HistType::kTH1F, {{250, -10., 10.}}}}, + {"h_d0_mass", ";m_{D^{0}} (GeV/c^{2});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_eta", ";#eta_{D^{0}} (GeV/c^{2});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, + {"h_d0_phi", ";#phi_{D^{0}} (GeV/c^{2});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, + {"h_d0_ang", ";#lambda_{#kappa}^{#alpha};counts", {HistType::kTH1F, {{100, 0., 1.}}}}} + }; + + // Configurables + Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; + + Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; + Configurable jetR{"jetR", 0.4, "jet resolution parameter"}; + + Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; + Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; + Configurable kappa{"kappa", 1.0, "angularity kappa"}; //to do: configurable from json + Configurable alpha{"alpha", 1.0, "angularity alpha"}; + + std::vector eventSelectionBits; + int trackSelection = -1; + + // Output table producer + Produces ObjJetTable; + + float angularity; + float leadingConstituentPt; + + void init(o2::framework::InitContext&) + { + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + } + + Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); + Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; + + template + void jetCalculateAngularity(T const& jet, U const& /*tracks*/) + { + angularity = 0.0; + leadingConstituentPt = 0.0; + for (auto& constituent : jet.template tracks_as()) { + if (constituent.pt() >= leadingConstituentPt) { + leadingConstituentPt = constituent.pt(); + } + angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent), alpha); + } + angularity /= (jet.pt() * (jet.r() / 100.f)); + } + // Collision-level and for Tracks QA + 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); + + // Loop over tracks and apply track selection + for (auto const& track : tracks) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + registry.fill(HIST("h_track_pt"), track.pt()); + registry.fill(HIST("h_track_eta"), track.eta()); + registry.fill(HIST("h_track_phi"), track.phi()); + } + } + PROCESS_SWITCH(JetD0Substructure, processCollisions, "process JE collisions", false); + + // Charged jet processing (no HF requirement) + void processDataCharged(soa::Filtered::iterator const& collision, soa::Filtered 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()); + registry.fill(HIST("h_jet_phi"), jet.phi()); + } + } + PROCESS_SWITCH(JetD0Substructure, processDataCharged, "charged jets in data", false); + + void processDataChargedSubstructure(aod::JetCollision const& collision, + soa::Join const& jets, + aod::CandidatesD0Data const&, aod::JetTracks const& tracks) + { + // apply event selection and fill histograms for sanity check + registry.fill(HIST("h_collision_counter"), 0.5); + + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { + return; + } + registry.fill(HIST("h_collision_counter"), 1.5); + + // Loop over jets containing D0 candidates + for (const auto& jet : jets) { + //number of charged jets with D0 + registry.fill(HIST("h_jet_counter"), 0.5); + // obtaining jet 3-vector + TVector3 jetVector(jet.px(), jet.py(), jet.pz()); + + //Loop over D0 candidates associated to the jet + for (const auto& d0Candidate : jet.candidates_as()) { + // obtaining jet 3-vector + TVector3 d0Vector(d0Candidate.px(), d0Candidate.py(), d0Candidate.pz()); + + // calculating fraction of the jet momentum carried by the D0 along the direction of the jet axis + double zParallel = (jetVector * d0Vector) / (jetVector * jetVector); + + // calculating angular distance in eta-phi plane + double axisDistance = jetutilities::deltaR(jet, d0Candidate); + + jetCalculateAngularity(jet,tracks); + + // filling histograms + registry.fill(HIST("h_d0_jet_projection"), zParallel); + registry.fill(HIST("h_d0_jet_distance_vs_projection"), axisDistance, zParallel); + registry.fill(HIST("h_d0_jet_distance"), axisDistance); + registry.fill(HIST("h_d0_jet_pt"), jet.pt()); + registry.fill(HIST("h_d0_jet_eta"), jet.eta()); + registry.fill(HIST("h_d0_jet_phi"), jet.phi()); + registry.fill(HIST("h_d0_mass"), d0Candidate.m()); + registry.fill(HIST("h_d0_eta"), d0Candidate.eta()); + registry.fill(HIST("h_d0_phi"), d0Candidate.phi()); + registry.fill(HIST("h_d0_ang"), angularity); // add more axis + + + // filling table + ObjJetTable(axisDistance, + jet.pt(), jet.eta(), jet.phi(), jet.tracks_as().size(), angularity, + d0Candidate.pt(), d0Candidate.eta(), d0Candidate.phi(), d0Candidate.m(), d0Candidate.y(), d0Candidate.mlScores()[0], d0Candidate.mlScores()[1], d0Candidate.mlScores()[2]); + + break; // get out of candidates' loop after first HF particle is found in jet + } // end of D0 candidates loop + + } // end of jets loop + + } // end of process function + PROCESS_SWITCH(JetD0Substructure, processDataChargedSubstructure, "charged HF jet substructure", false); + + + +}; +// Workflow definition +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-substructure"})}; } From 89724d5c45ac0fe3b54c11befc8747200afd190f Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 9 Mar 2026 12:36:48 +0000 Subject: [PATCH 2/5] Please consider the following formatting changes --- PWGJE/Tasks/jetD0Substructure.cxx | 323 +++++++++++++++--------------- 1 file changed, 158 insertions(+), 165 deletions(-) diff --git a/PWGJE/Tasks/jetD0Substructure.cxx b/PWGJE/Tasks/jetD0Substructure.cxx index 18d2dbd6328..b1309d2e64f 100644 --- a/PWGJE/Tasks/jetD0Substructure.cxx +++ b/PWGJE/Tasks/jetD0Substructure.cxx @@ -16,32 +16,30 @@ // \inherited from D0 fragmentation and Ds // \P. Dhankher - -#include "PWGJE/DataModel/Jet.h" -#include "PWGJE/DataModel/JetSubtraction.h" -#include "PWGJE/DataModel/JetReducedData.h" - +#include "PWGHF/Core/DecayChannels.h" #include "PWGJE/Core/JetDerivedDataUtilities.h" #include "PWGJE/Core/JetUtilities.h" -#include "PWGHF/Core/DecayChannels.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" +#include "PWGJE/DataModel/JetSubtraction.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/TrackSelectionTables.h" #include "Framework/ASoA.h" #include "Framework/AnalysisDataModel.h" -#include #include "Framework/HistogramRegistry.h" #include +#include +#include #include +#include #include #include -#include -#include #include - #include "TVector3.h" + #include #include #include @@ -76,23 +74,23 @@ DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float); } // namespace jet_obj // AOD table definition DECLARE_SOA_TABLE(JetObjTable, "AOD", "JETOBJTABLE", - jet_obj::JetHfDist, - jet_obj::JetPt, - jet_obj::JetEta, - jet_obj::JetPhi, - jet_obj::JetNConst, - jet_obj::JetAng, - jet_obj::HfPt, - jet_obj::HfEta, - jet_obj::HfPhi, - jet_obj::HfMass, - jet_obj::HfY, - jet_obj::HfMlScore0, - jet_obj::HfMlScore1, - jet_obj::HfMlScore2); -} + jet_obj::JetHfDist, + jet_obj::JetPt, + jet_obj::JetEta, + jet_obj::JetPhi, + jet_obj::JetNConst, + jet_obj::JetAng, + jet_obj::HfPt, + jet_obj::HfEta, + jet_obj::HfPhi, + jet_obj::HfMass, + jet_obj::HfY, + jet_obj::HfMlScore0, + jet_obj::HfMlScore1, + jet_obj::HfMlScore2); +} // namespace o2::aod struct JetD0Substructure { - /** + /** * Histogram registry * * Contains: @@ -100,167 +98,162 @@ struct JetD0Substructure { * - Jet kinematic distributions * - D0–jet substructure observables */ - 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, {{2, 0., 2.}}}}, - {"h_jet_counter", ";# of D^{0} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, - {"h_d0_jet_projection", ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_jet_distance_vs_projection", ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}}, - {"h_d0_jet_distance", ";#DeltaR_{D^{0},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_jet_pt", ";p_{T,D^{0} jet};dN/dp_{T,D^{0} jet}", {HistType::kTH1F, {{200, 0., 10.}}}}, - {"h_d0_jet_eta", ";#eta_{T,D^{0} jet};dN/d#eta_{D^{0} jet}", {HistType::kTH1F, {{250, -5., 5.}}}}, - {"h_d0_jet_phi", ";#phi_{T,D^{0} jet};dN/d#phi_{D^{0} jet}", {HistType::kTH1F, {{250, -10., 10.}}}}, - {"h_d0_mass", ";m_{D^{0}} (GeV/c^{2});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_eta", ";#eta_{D^{0}} (GeV/c^{2});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, - {"h_d0_phi", ";#phi_{D^{0}} (GeV/c^{2});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, - {"h_d0_ang", ";#lambda_{#kappa}^{#alpha};counts", {HistType::kTH1F, {{100, 0., 1.}}}}} - }; - - // Configurables - Configurable 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, {{2, 0., 2.}}}}, + {"h_jet_counter", ";# of D^{0} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, + {"h_d0_jet_projection", ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_jet_distance_vs_projection", ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}}, + {"h_d0_jet_distance", ";#DeltaR_{D^{0},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_jet_pt", ";p_{T,D^{0} jet};dN/dp_{T,D^{0} jet}", {HistType::kTH1F, {{200, 0., 10.}}}}, + {"h_d0_jet_eta", ";#eta_{T,D^{0} jet};dN/d#eta_{D^{0} jet}", {HistType::kTH1F, {{250, -5., 5.}}}}, + {"h_d0_jet_phi", ";#phi_{T,D^{0} jet};dN/d#phi_{D^{0} jet}", {HistType::kTH1F, {{250, -10., 10.}}}}, + {"h_d0_mass", ";m_{D^{0}} (GeV/c^{2});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_eta", ";#eta_{D^{0}} (GeV/c^{2});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, + {"h_d0_phi", ";#phi_{D^{0}} (GeV/c^{2});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, + {"h_d0_ang", ";#lambda_{#kappa}^{#alpha};counts", {HistType::kTH1F, {{100, 0., 1.}}}}}}; - Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; - Configurable jetR{"jetR", 0.4, "jet resolution parameter"}; + // Configurables + Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; - Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; - Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; - Configurable kappa{"kappa", 1.0, "angularity kappa"}; //to do: configurable from json - Configurable alpha{"alpha", 1.0, "angularity alpha"}; - - std::vector eventSelectionBits; - int trackSelection = -1; + Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; + Configurable jetR{"jetR", 0.4, "jet resolution parameter"}; - // Output table producer - Produces ObjJetTable; - - float angularity; - float leadingConstituentPt; - - void init(o2::framework::InitContext&) - { - eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); - trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); - } - - Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); - Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; + Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; + Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; + Configurable kappa{"kappa", 1.0, "angularity kappa"}; // to do: configurable from json + Configurable alpha{"alpha", 1.0, "angularity alpha"}; + + std::vector eventSelectionBits; + int trackSelection = -1; + + // Output table producer + Produces ObjJetTable; - template - void jetCalculateAngularity(T const& jet, U const& /*tracks*/) - { - angularity = 0.0; - leadingConstituentPt = 0.0; - for (auto& constituent : jet.template tracks_as()) { - if (constituent.pt() >= leadingConstituentPt) { - leadingConstituentPt = constituent.pt(); - } - angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent), alpha); - } - angularity /= (jet.pt() * (jet.r() / 100.f)); + float angularity; + float leadingConstituentPt; + + void init(o2::framework::InitContext&) + { + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + } + + Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); + Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; + + template + void jetCalculateAngularity(T const& jet, U const& /*tracks*/) + { + angularity = 0.0; + leadingConstituentPt = 0.0; + for (auto& constituent : jet.template tracks_as()) { + if (constituent.pt() >= leadingConstituentPt) { + leadingConstituentPt = constituent.pt(); + } + angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent), alpha); } - // Collision-level and for Tracks QA - void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks) - { + angularity /= (jet.pt() * (jet.r() / 100.f)); + } + // Collision-level and for Tracks QA + 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); + registry.fill(HIST("h_collisions"), 0.5); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } + registry.fill(HIST("h_collisions"), 1.5); - // Loop over tracks and apply track selection - for (auto const& track : tracks) { - if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; - } - registry.fill(HIST("h_track_pt"), track.pt()); - registry.fill(HIST("h_track_eta"), track.eta()); - registry.fill(HIST("h_track_phi"), track.phi()); - } + // Loop over tracks and apply track selection + for (auto const& track : tracks) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + registry.fill(HIST("h_track_pt"), track.pt()); + registry.fill(HIST("h_track_eta"), track.eta()); + registry.fill(HIST("h_track_phi"), track.phi()); } - PROCESS_SWITCH(JetD0Substructure, processCollisions, "process JE collisions", false); + } + PROCESS_SWITCH(JetD0Substructure, processCollisions, "process JE collisions", false); - // Charged jet processing (no HF requirement) - void processDataCharged(soa::Filtered::iterator const& collision, soa::Filtered 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()); - registry.fill(HIST("h_jet_phi"), jet.phi()); - } + // Charged jet processing (no HF requirement) + void processDataCharged(soa::Filtered::iterator const& collision, soa::Filtered 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()); + registry.fill(HIST("h_jet_phi"), jet.phi()); } - PROCESS_SWITCH(JetD0Substructure, processDataCharged, "charged jets in data", false); + } + PROCESS_SWITCH(JetD0Substructure, processDataCharged, "charged jets in data", false); - void processDataChargedSubstructure(aod::JetCollision const& collision, - soa::Join const& jets, - aod::CandidatesD0Data const&, aod::JetTracks const& tracks) - { - // apply event selection and fill histograms for sanity check - registry.fill(HIST("h_collision_counter"), 0.5); + void processDataChargedSubstructure(aod::JetCollision const& collision, + soa::Join const& jets, + aod::CandidatesD0Data const&, aod::JetTracks const& tracks) + { + // apply event selection and fill histograms for sanity check + registry.fill(HIST("h_collision_counter"), 0.5); - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { - return; - } - registry.fill(HIST("h_collision_counter"), 1.5); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { + return; + } + registry.fill(HIST("h_collision_counter"), 1.5); - // Loop over jets containing D0 candidates - for (const auto& jet : jets) { - //number of charged jets with D0 - registry.fill(HIST("h_jet_counter"), 0.5); - // obtaining jet 3-vector - TVector3 jetVector(jet.px(), jet.py(), jet.pz()); + // Loop over jets containing D0 candidates + for (const auto& jet : jets) { + // number of charged jets with D0 + registry.fill(HIST("h_jet_counter"), 0.5); + // obtaining jet 3-vector + TVector3 jetVector(jet.px(), jet.py(), jet.pz()); - //Loop over D0 candidates associated to the jet - for (const auto& d0Candidate : jet.candidates_as()) { - // obtaining jet 3-vector - TVector3 d0Vector(d0Candidate.px(), d0Candidate.py(), d0Candidate.pz()); + // Loop over D0 candidates associated to the jet + for (const auto& d0Candidate : jet.candidates_as()) { + // obtaining jet 3-vector + TVector3 d0Vector(d0Candidate.px(), d0Candidate.py(), d0Candidate.pz()); - // calculating fraction of the jet momentum carried by the D0 along the direction of the jet axis - double zParallel = (jetVector * d0Vector) / (jetVector * jetVector); + // calculating fraction of the jet momentum carried by the D0 along the direction of the jet axis + double zParallel = (jetVector * d0Vector) / (jetVector * jetVector); - // calculating angular distance in eta-phi plane - double axisDistance = jetutilities::deltaR(jet, d0Candidate); - - jetCalculateAngularity(jet,tracks); - - // filling histograms - registry.fill(HIST("h_d0_jet_projection"), zParallel); - registry.fill(HIST("h_d0_jet_distance_vs_projection"), axisDistance, zParallel); - registry.fill(HIST("h_d0_jet_distance"), axisDistance); - registry.fill(HIST("h_d0_jet_pt"), jet.pt()); - registry.fill(HIST("h_d0_jet_eta"), jet.eta()); - registry.fill(HIST("h_d0_jet_phi"), jet.phi()); - registry.fill(HIST("h_d0_mass"), d0Candidate.m()); - registry.fill(HIST("h_d0_eta"), d0Candidate.eta()); - registry.fill(HIST("h_d0_phi"), d0Candidate.phi()); - registry.fill(HIST("h_d0_ang"), angularity); // add more axis - + // calculating angular distance in eta-phi plane + double axisDistance = jetutilities::deltaR(jet, d0Candidate); - // filling table - ObjJetTable(axisDistance, - jet.pt(), jet.eta(), jet.phi(), jet.tracks_as().size(), angularity, - d0Candidate.pt(), d0Candidate.eta(), d0Candidate.phi(), d0Candidate.m(), d0Candidate.y(), d0Candidate.mlScores()[0], d0Candidate.mlScores()[1], d0Candidate.mlScores()[2]); - - break; // get out of candidates' loop after first HF particle is found in jet - } // end of D0 candidates loop + jetCalculateAngularity(jet, tracks); - } // end of jets loop + // filling histograms + registry.fill(HIST("h_d0_jet_projection"), zParallel); + registry.fill(HIST("h_d0_jet_distance_vs_projection"), axisDistance, zParallel); + registry.fill(HIST("h_d0_jet_distance"), axisDistance); + registry.fill(HIST("h_d0_jet_pt"), jet.pt()); + registry.fill(HIST("h_d0_jet_eta"), jet.eta()); + registry.fill(HIST("h_d0_jet_phi"), jet.phi()); + registry.fill(HIST("h_d0_mass"), d0Candidate.m()); + registry.fill(HIST("h_d0_eta"), d0Candidate.eta()); + registry.fill(HIST("h_d0_phi"), d0Candidate.phi()); + registry.fill(HIST("h_d0_ang"), angularity); // add more axis - } // end of process function - PROCESS_SWITCH(JetD0Substructure, processDataChargedSubstructure, "charged HF jet substructure", false); + // filling table + ObjJetTable(axisDistance, + jet.pt(), jet.eta(), jet.phi(), jet.tracks_as().size(), angularity, + d0Candidate.pt(), d0Candidate.eta(), d0Candidate.phi(), d0Candidate.m(), d0Candidate.y(), d0Candidate.mlScores()[0], d0Candidate.mlScores()[1], d0Candidate.mlScores()[2]); + break; // get out of candidates' loop after first HF particle is found in jet + } // end of D0 candidates loop + } // end of jets loop + } // end of process function + PROCESS_SWITCH(JetD0Substructure, processDataChargedSubstructure, "charged HF jet substructure", false); }; // Workflow definition WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-substructure"})}; } From 64f6f89e0691a1298711aceb6a75b92bf6cd5a6a Mon Sep 17 00:00:00 2001 From: Preeti Dhankher Date: Tue, 10 Mar 2026 16:45:36 +0100 Subject: [PATCH 3/5] Renaming the task to avoid confusion with other task with similar name --- PWGJE/Tasks/CMakeLists.txt | 4 +- ...structure.cxx => jetD0AngSubstructure.cxx} | 72 +++---------------- 2 files changed, 13 insertions(+), 63 deletions(-) rename PWGJE/Tasks/{jetD0Substructure.cxx => jetD0AngSubstructure.cxx} (71%) diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index d78d84fbafb..e338093a33c 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -389,8 +389,8 @@ if(FastJet_FOUND) SOURCES jetDsSpectrumAndSubstructure.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(jet-d0-substructure - SOURCES jetD0Substructure.cxx + o2physics_add_dpl_workflow(jet-d0-ang-substructure + SOURCES jetD0AngSubstructure.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(bjet-cent-mult diff --git a/PWGJE/Tasks/jetD0Substructure.cxx b/PWGJE/Tasks/jetD0AngSubstructure.cxx similarity index 71% rename from PWGJE/Tasks/jetD0Substructure.cxx rename to PWGJE/Tasks/jetD0AngSubstructure.cxx index b1309d2e64f..b8a9fc3fea0 100644 --- a/PWGJE/Tasks/jetD0Substructure.cxx +++ b/PWGJE/Tasks/jetD0AngSubstructure.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -// \file jetD0Substructure.cxx +// \file JetD0AngSubstructure.cxx // // \brief Analysis task for the reconstruction and study of charged jets // containing D_0 mesons in pp collisions. @@ -89,7 +89,7 @@ DECLARE_SOA_TABLE(JetObjTable, "AOD", "JETOBJTABLE", jet_obj::HfMlScore1, jet_obj::HfMlScore2); } // namespace o2::aod -struct JetD0Substructure { +struct JetD0AngSubstructure { /** * Histogram registry * @@ -99,14 +99,7 @@ struct JetD0Substructure { * - D0–jet substructure observables */ 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, {{2, 0., 2.}}}}, + {{"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{2, 0., 2.}}}}, {"h_jet_counter", ";# of D^{0} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, {"h_d0_jet_projection", ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}}, {"h_d0_jet_distance_vs_projection", ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}}, @@ -114,9 +107,9 @@ struct JetD0Substructure { {"h_d0_jet_pt", ";p_{T,D^{0} jet};dN/dp_{T,D^{0} jet}", {HistType::kTH1F, {{200, 0., 10.}}}}, {"h_d0_jet_eta", ";#eta_{T,D^{0} jet};dN/d#eta_{D^{0} jet}", {HistType::kTH1F, {{250, -5., 5.}}}}, {"h_d0_jet_phi", ";#phi_{T,D^{0} jet};dN/d#phi_{D^{0} jet}", {HistType::kTH1F, {{250, -10., 10.}}}}, - {"h_d0_mass", ";m_{D^{0}} (GeV/c^{2});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_eta", ";#eta_{D^{0}} (GeV/c^{2});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, - {"h_d0_phi", ";#phi_{D^{0}} (GeV/c^{2});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, + {"h_d0_mass", ";m_{D^{0}});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, + {"h_d0_eta", ";#eta_{D^{0}});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, + {"h_d0_phi", ";#phi_{D^{0}});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, {"h_d0_ang", ";#lambda_{#kappa}^{#alpha};counts", {HistType::kTH1F, {{100, 0., 1.}}}}}}; // Configurables @@ -137,7 +130,6 @@ struct JetD0Substructure { Produces ObjJetTable; float angularity; - float leadingConstituentPt; void init(o2::framework::InitContext&) { @@ -145,58 +137,15 @@ struct JetD0Substructure { trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); } - Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); - Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; - template void jetCalculateAngularity(T const& jet, U const& /*tracks*/) { angularity = 0.0; - leadingConstituentPt = 0.0; for (auto& constituent : jet.template tracks_as()) { - if (constituent.pt() >= leadingConstituentPt) { - leadingConstituentPt = constituent.pt(); - } - angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent), alpha); - } - angularity /= (jet.pt() * (jet.r() / 100.f)); - } - // Collision-level and for Tracks QA - 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); - - // Loop over tracks and apply track selection - for (auto const& track : tracks) { - if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; - } - registry.fill(HIST("h_track_pt"), track.pt()); - registry.fill(HIST("h_track_eta"), track.eta()); - registry.fill(HIST("h_track_phi"), track.phi()); + angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent)/(jet.r() / 100.f), alpha); } + angularity /= std::pow(jet.pt(), kappa); } - PROCESS_SWITCH(JetD0Substructure, processCollisions, "process JE collisions", false); - - // Charged jet processing (no HF requirement) - void processDataCharged(soa::Filtered::iterator const& collision, soa::Filtered 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()); - registry.fill(HIST("h_jet_phi"), jet.phi()); - } - } - PROCESS_SWITCH(JetD0Substructure, processDataCharged, "charged jets in data", false); void processDataChargedSubstructure(aod::JetCollision const& collision, soa::Join const& jets, @@ -253,7 +202,8 @@ struct JetD0Substructure { } // end of jets loop } // end of process function - PROCESS_SWITCH(JetD0Substructure, processDataChargedSubstructure, "charged HF jet substructure", false); + PROCESS_SWITCH(JetD0AngSubstructure, processDataChargedSubstructure, "charged HF jet substructure", false); }; // Workflow definition -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-substructure"})}; } +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-ang-substructure"})}; } + From a0bd868f5b567c5551fffff19aeb74e27ec04fea Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 10 Mar 2026 15:47:23 +0000 Subject: [PATCH 4/5] Please consider the following formatting changes --- PWGJE/Tasks/jetD0AngSubstructure.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGJE/Tasks/jetD0AngSubstructure.cxx b/PWGJE/Tasks/jetD0AngSubstructure.cxx index b8a9fc3fea0..05c1923dd96 100644 --- a/PWGJE/Tasks/jetD0AngSubstructure.cxx +++ b/PWGJE/Tasks/jetD0AngSubstructure.cxx @@ -142,7 +142,7 @@ struct JetD0AngSubstructure { { angularity = 0.0; for (auto& constituent : jet.template tracks_as()) { - angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent)/(jet.r() / 100.f), alpha); + angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent) / (jet.r() / 100.f), alpha); } angularity /= std::pow(jet.pt(), kappa); } @@ -206,4 +206,3 @@ struct JetD0AngSubstructure { }; // Workflow definition WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-ang-substructure"})}; } - From 37349f2816eb0bcfcbee04acf3e33331534b9280 Mon Sep 17 00:00:00 2001 From: Lieuwe Huisman <153300723+Lolle2001@users.noreply.github.com> Date: Tue, 21 Apr 2026 14:55:22 +0200 Subject: [PATCH 5/5] [PWGJE] Updated jetD0AngSubstructure.cxx - Added task to analyse D0 candidates from experimental data. Histograms for different observables such as angularity etc. are made for quick checking purposes. - Added task to match detector level and particle level Monte Carlo particles. - Added an additional task to calculate the Monte Carlo efficiency. - Using templated analyse functions to easily define processfunctions for other particles than D0 if desired. --- PWGJE/Tasks/jetD0AngSubstructure.cxx | 878 ++++++++++++++++++++++++--- 1 file changed, 781 insertions(+), 97 deletions(-) diff --git a/PWGJE/Tasks/jetD0AngSubstructure.cxx b/PWGJE/Tasks/jetD0AngSubstructure.cxx index 9650d7ea4bb..c4d4058e3db 100644 --- a/PWGJE/Tasks/jetD0AngSubstructure.cxx +++ b/PWGJE/Tasks/jetD0AngSubstructure.cxx @@ -9,35 +9,44 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -// \file JetD0AngSubstructure.cxx +/// \file JetD0AngSubstructure.cxx // -// \brief Analysis task for the reconstruction and study of charged jets +/// \brief Analysis task for the reconstruction and study of charged jets // containing D_0 mesons in pp collisions. -// \inherited from D0 fragmentation and Ds -// \P. Dhankher - -#include "PWGJE/Core/JetDerivedDataUtilities.h" -#include "PWGJE/Core/JetUtilities.h" -#include "PWGJE/DataModel/Jet.h" - +// inherited from D0 fragmentation and Ds +/// \author P. Dhankher +/// \author L.J. Huisman +// #include #include -#include #include #include #include +#include #include #include #include #include - +// +#include +// +#include +#include +#include +#include +#include +// +#include +// +#include +// #include - +// #include #include #include #include - +// using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -45,65 +54,276 @@ using namespace o2::framework::expressions; // Definition of a custom AOD table to store jet–D0 quantities namespace o2::aod { -namespace jet_obj -{ +/* // Jet-related quantities -DECLARE_SOA_COLUMN(JetHfDist, jetHfDist, float); -DECLARE_SOA_COLUMN(JetPt, jetPt, float); -DECLARE_SOA_COLUMN(JetEta, jetEta, float); -DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); -DECLARE_SOA_COLUMN(JetNConst, jetNConst, int); -DECLARE_SOA_COLUMN(JetAng, jetAng, float); +*/ +DECLARE_SOA_COLUMN(expJetHfDist, ExpJetHfDist, float); +DECLARE_SOA_COLUMN(expJetPt, ExpJetPt, float); +DECLARE_SOA_COLUMN(expJetEta, ExpJetEta, float); +DECLARE_SOA_COLUMN(expJetPhi, ExpJetPhi, float); +DECLARE_SOA_COLUMN(expJetNConst, ExpJetNConst, int); +DECLARE_SOA_COLUMN(expJetAng, ExpJetAng, float); // D0 candidate quantities -DECLARE_SOA_COLUMN(HfPt, hfPt, float); -DECLARE_SOA_COLUMN(HfEta, hfEta, float); -DECLARE_SOA_COLUMN(HfPhi, hfPhi, float); -DECLARE_SOA_COLUMN(HfMass, hfMass, float); -DECLARE_SOA_COLUMN(HfY, hfY, float); +DECLARE_SOA_COLUMN(expHfPt, ExpHfPt, float); +DECLARE_SOA_COLUMN(expHfEta, ExpHfEta, float); +DECLARE_SOA_COLUMN(expHfPhi, ExpHfPhi, float); +DECLARE_SOA_COLUMN(expHfMass, ExpHfMass, float); +DECLARE_SOA_COLUMN(expHfY, ExpHfY, float); +// ML scores +DECLARE_SOA_COLUMN(expHfMlScore0, ExpHfMlScore0, float); +DECLARE_SOA_COLUMN(expHfMlScore1, ExpHfMlScore1, float); +DECLARE_SOA_COLUMN(expHfMlScore2, ExpHfMlScore2, float); + +/* +// MCP quantities (Particle Level) +*/ +// Jets +DECLARE_SOA_COLUMN(mcpJetHfDist, McpJetHfDist, float); +DECLARE_SOA_COLUMN(mcpJetPt, McpJetPt, float); +DECLARE_SOA_COLUMN(mcpJetEta, McpJetEta, float); +DECLARE_SOA_COLUMN(mcpJetPhi, McpJetPhi, float); +DECLARE_SOA_COLUMN(mcpJetNConst, McpJetNConst, float); +DECLARE_SOA_COLUMN(mcpJetAng, McpJetAng, float); +// D0 candidates (Heavy Flavour) +DECLARE_SOA_COLUMN(mcpHfPt, McpHfPt, float); +DECLARE_SOA_COLUMN(mcpHfEta, McpHfEta, float); +DECLARE_SOA_COLUMN(mcpHfPhi, McpHfPhi, float); +// DECLARE_SOA_COLUMN(HfMass, mcpHfMass, float); +DECLARE_SOA_COLUMN(mcpHfY, McpHfY, float); +DECLARE_SOA_COLUMN(mcpHfPrompt, McpHfPrompt, bool); +DECLARE_SOA_COLUMN(mcpHfMatch, McpHfMatch, bool); + +/* +// MCD quantities (Detector Level) +*/ +// Jets +DECLARE_SOA_COLUMN(mcdJetHfDist, McdJetHfDist, float); +DECLARE_SOA_COLUMN(mcdJetPt, McdJetPt, float); +DECLARE_SOA_COLUMN(mcdJetEta, McdJetEta, float); +DECLARE_SOA_COLUMN(mcdJetPhi, McdJetPhi, float); +DECLARE_SOA_COLUMN(mcdJetNConst, McdJetNConst, float); +DECLARE_SOA_COLUMN(mcdJetAng, McdJetAng, float); +// D0 candidates (Heavy Flavour) +DECLARE_SOA_COLUMN(mcdHfPt, McdHfPt, float); +DECLARE_SOA_COLUMN(mcdHfEta, McdHfEta, float); +DECLARE_SOA_COLUMN(mcdHfPhi, McdHfPhi, float); +DECLARE_SOA_COLUMN(mcdHfMass, McdHfMass, float); +DECLARE_SOA_COLUMN(mcdHfY, McdHfY, float); +DECLARE_SOA_COLUMN(mcdHfPrompt, McdHfPrompt, bool); +DECLARE_SOA_COLUMN(mcdHfMatch, McdHfMatch, bool); +// Other +DECLARE_SOA_COLUMN(mcdHfMatchedFrom, McdHfMatchedFrom, int); +DECLARE_SOA_COLUMN(mcdHfSelectedAs, McdHfSelectedAs, int); // ML scores -DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float); -DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float); -DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float); -} // namespace jet_obj +DECLARE_SOA_COLUMN(mcdHfMlScore0, McdHfMlScore0, float); +DECLARE_SOA_COLUMN(mcdHfMlScore1, McdHfMlScore1, float); +DECLARE_SOA_COLUMN(mcdHfMlScore2, McdHfMlScore2, float); + +/* // AOD table definition -DECLARE_SOA_TABLE(JetObjTable, "AOD", "JETOBJTABLE", - jet_obj::JetHfDist, - jet_obj::JetPt, - jet_obj::JetEta, - jet_obj::JetPhi, - jet_obj::JetNConst, - jet_obj::JetAng, - jet_obj::HfPt, - jet_obj::HfEta, - jet_obj::HfPhi, - jet_obj::HfMass, - jet_obj::HfY, - jet_obj::HfMlScore0, - jet_obj::HfMlScore1, - jet_obj::HfMlScore2); +*/ +DECLARE_SOA_TABLE(EXPJetObjTable, "AOD", "EXPJETOBJTABLE", + expJetHfDist, + expJetPt, + expJetEta, + expJetPhi, + expJetNConst, + expJetAng, + expHfPt, + expHfEta, + expHfPhi, + expHfMass, + expHfY, + expHfMlScore0, + expHfMlScore1, + expHfMlScore2); + +DECLARE_SOA_TABLE(MCPJetObjTable, "AOD", "MCPJETOBJTABLE", + mcpJetHfDist, + mcpJetPt, + mcpJetEta, + mcpJetPhi, + mcpJetNConst, + mcpJetAng, + mcpHfPt, + mcpHfEta, + mcpHfPhi, + mcpHfY, + mcpHfPrompt, + mcpHfMatch); + +DECLARE_SOA_TABLE(MCDJetObjTable, "AOD", "MCDJETOBJTABLE", + mcdJetHfDist, + mcdJetPt, + mcdJetEta, + mcdJetPhi, + mcdJetNConst, + mcdJetAng, + mcdHfPt, + mcdHfEta, + mcdHfPhi, + mcdHfMass, + mcdHfY, + mcdHfPrompt, + mcdHfMatch, + mcdHfMlScore0, + mcdHfMlScore1, + mcdHfMlScore2, + mcdHfMatchedFrom, + mcdHfSelectedAs); + +DECLARE_SOA_TABLE(MatchJetDistanceTable, "AOD", "MATCHTABLE", + mcpJetHfDist, + mcpJetPt, + mcpJetEta, + mcpJetPhi, + mcpJetNConst, + mcpJetAng, + mcpHfPt, + mcpHfEta, + mcpHfPhi, + mcpHfY, + mcpHfPrompt, + mcdJetHfDist, + mcdJetPt, + mcdJetEta, + mcdJetPhi, + mcdJetNConst, + mcdJetAng, + mcdHfPt, + mcdHfEta, + mcdHfPhi, + mcdHfMass, + mcdHfY, + mcdHfPrompt, + mcdHfMlScore0, + mcdHfMlScore1, + mcdHfMlScore2, + mcdHfMatchedFrom, + mcdHfSelectedAs); + } // namespace o2::aod + +// Helps to avoid typos in histogram names when using the "HIST" macro +/* +hf(l): Heavy flavour +ex(p): Experimental data (detector level) +mcd: Monte Carlo data (detector level) +mcp: Monte Carlo data (particle level) +*/ +namespace histnames +{ +#define HNAME(name) constexpr const char* name = #name; +/* +// Experimental Data (analyseDataChargedSubstructure) +*/ +HNAME(ex_col); // Collision Counter +HNAME(ex_jet); // Jet Counter +HNAME(ex_jet_pt); // Jet pT +HNAME(ex_jet_eta); // Jet eta +HNAME(ex_jet_phi); // Jet phi +HNAME(ex_jet_ang); // Jet angularity +HNAME(ex_jet_proj); // Projection of HF candidate momentum on jet axis +HNAME(ex_jet_dist); // Jet–HF candidate angular distance +HNAME(ex_jet_dist_proj); // Jet–HF candidate distance vs projection +HNAME(ex_hfl_pt); // HF candidate pT +HNAME(ex_hfl_mass); // HF candidate mass +HNAME(ex_hfl_eta); // HF candidate eta +HNAME(ex_hfl_phi); // HF candidate phi +/* +// Monte Carlo Data Efficiency (analyseMonteCarloEfficiency) +*/ +HNAME(mc_eff_col); // Collision Counter +HNAME(mc_eff_jet); // Jet Counter +HNAME(mc_eff_det_jet_pt); // Detector level jet pT +HNAME(mc_eff_det_jet_eta); // Detector level jet eta +HNAME(mc_eff_det_jet_phi); // Detector level jet phi +HNAME(mc_eff_det_jet_ang); // Detector level jet angularity +HNAME(mc_eff_det_hfl_pt); // Detector level HF candidate pT +HNAME(mc_eff_det_hfl_mass); // Detector level HF candidate mass +HNAME(mc_eff_det_hfl_eta); // Detector level HF candidate eta +HNAME(mc_eff_det_hfl_phi); // Detector level HF candidate phi +HNAME(mc_eff_par_jet_pt); // Particle level jet pT +HNAME(mc_eff_par_jet_eta); // Particle level jet eta +HNAME(mc_eff_par_jet_phi); // Particle level jet phi +HNAME(mc_eff_par_jet_ang); // Particle level jet angularity +HNAME(mc_eff_par_hfl_pt); // Particle level HF candidate pT +// HNAME(mc_eff_par_hfl_mass); // Particle level HF candidate mass -> PDG value +HNAME(mc_eff_par_hfl_eta); // Particle level HF candidate eta +HNAME(mc_eff_par_hfl_phi); // Particle level HF candidate phi +/* +// Monte Carlo Data Matching (analyseMonteCarlo) +*/ +HNAME(mc_col); // Collision Counter +HNAME(mc_jet); // Jet Counter +HNAME(mc_det_jet_pt); // Detector level jet pT +HNAME(mc_det_jet_eta); // Detector level jet eta +HNAME(mc_det_jet_phi); // Detector level jet phi +HNAME(mc_det_jet_ang); // Detector level jet angularity +HNAME(mc_det_hfl_pt); // Detector level HF candidate pT (matched) +HNAME(mc_det_hfl_mass); // Detector level HF candidate mass (matched) +HNAME(mc_det_hfl_eta); // Detector level HF candidate eta (matched) +HNAME(mc_det_hfl_phi); // Detector level HF candidate phi (matched) +// HNAME(mc_par_jet_pt); // Particle level jet pT +// HNAME(mc_par_jet_eta); // Particle level jet eta +// HNAME(mc_par_jet_phi); // Particle level jet phi +// HNAME(mc_par_jet_ang); // Particle level jet angularity +// HNAME(mc_par_hfl_pt); // Particle level HF candidate pT +// HNAME(mc_par_hfl_mass); // Particle level HF candidate mass -> PDG value +// HNAME(mc_par_hfl_eta); // Particle level HF candidate eta +// HNAME(mc_par_hfl_phi); // Particle level HF candidate phi + +#undef HNAME +} // namespace histnames + +consteval float get_val_from_bin(int bin) +{ + return static_cast(bin) - 0.5f; +} + +enum BIN_EX_COLCNTR { + AllCollisions = 1, + Sel8ZCut = 2 +}; + +enum BIN_EX_JETCNTR { + ChargedJets = 1 +}; +enum BIN_MC_COLCNTR { + All = 1, + ZCut = 2, + Matched = 3, + MatchedSel8ZCut = 4 +}; + +enum BIN_MC_JETCNTR { + DetectorLevelJetInMCCollision = 1, + ParticleLevelJetInMCCollision = 2, + DetectorLevelJetWithMatchedCandidate = 3, + ParticleLevelJetWithMatchedCandidate = 4 +}; + struct JetD0AngSubstructure { - /** - * Histogram registry - * - * Contains: - * - Event and track histograms - * - Jet kinematic distributions - * - D0–jet substructure observables - */ - HistogramRegistry registry{"registry", - {{"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{2, 0., 2.}}}}, - {"h_jet_counter", ";# of D^{0} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}}, - {"h_d0_jet_projection", ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_jet_distance_vs_projection", ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}}, - {"h_d0_jet_distance", ";#DeltaR_{D^{0},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_jet_pt", ";p_{T,D^{0} jet};dN/dp_{T,D^{0} jet}", {HistType::kTH1F, {{200, 0., 10.}}}}, - {"h_d0_jet_eta", ";#eta_{T,D^{0} jet};dN/d#eta_{D^{0} jet}", {HistType::kTH1F, {{250, -5., 5.}}}}, - {"h_d0_jet_phi", ";#phi_{T,D^{0} jet};dN/d#phi_{D^{0} jet}", {HistType::kTH1F, {{250, -10., 10.}}}}, - {"h_d0_mass", ";m_{D^{0}});dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}}, - {"h_d0_eta", ";#eta_{D^{0}});dN/d#eta_{D_{}}", {HistType::kTH1F, {{250, -5., 5.}}}}, - {"h_d0_phi", ";#phi_{D^{0}});dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}}, - {"h_d0_ang", ";#lambda_{#kappa}^{#alpha};counts", {HistType::kTH1F, {{100, 0., 1.}}}}}}; + + // Output table producer + Produces ObjJetTable; + Produces mcdJetTable; + Produces mcpJetTable; + Produces matchJetTable; + // + using JetChargedTableD0 = soa::Join; + // MC Matching Tables + using JetD0MCDTable = soa::Join; + using JetD0MCPTable = soa::Join; + // Type Aliases for the constituent jets in the table (to be used for angularity calculation) + using JetD0MCDTableConstituent = JetD0MCDTable::iterator; + using JetD0MCPTableConstituent = JetD0MCPTable::iterator; + // Slices for access to proper HF MCD jet collision that is associated to MCCollision + PresliceUnsorted collisionsPerMCCollisionPreslice = aod::jmccollisionlb::mcCollisionId; + Preslice d0MCDJetsPerEXPCollisionPreslice = aod::jet::collisionId; + Preslice d0MCPJetsPerMCCollisionPreslice = aod::jet::mcCollisionId; + + HistogramRegistry registryD0{"registryD0", {}}; // Configurables Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; @@ -119,48 +339,155 @@ struct JetD0AngSubstructure { std::vector eventSelectionBits; int trackSelection = -1; - // Output table producer - Produces ObjJetTable; - - float angularity; - void init(o2::framework::InitContext&) { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + addHistograms(); } + void addHistograms() + { + /* + // Experimental Data Histograms + */ + // Definitions + registryD0.add(histnames::ex_col, "N_{coll};", {HistType::kTH1F, {{2, 0., 2.}}}); + registryD0.add(histnames::ex_jet, "N_{jet};", {HistType::kTH1F, {{1, 0., 1.0}}}); + registryD0.add(histnames::ex_jet_pt, ";p_{T,jet};dN/dp_{T,jet}", {HistType::kTH1F, {{1000, 0., 50.}}}); + registryD0.add(histnames::ex_jet_eta, ";#eta_{jet};dN/d#eta_{jet}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::ex_jet_phi, ";#phi_{jet};dN/d#phi_{jet}", {HistType::kTH1F, {{250, -10., 10.}}}); + registryD0.add(histnames::ex_jet_ang, ";angularity_{jet};dN/d(angularity)", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::ex_jet_dist, ";#DeltaR_{D^{0},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::ex_jet_dist_proj, ";#DeltaR_{D^{0},jet};z^{D^{0},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}); + registryD0.add(histnames::ex_jet_proj, ";z^{D^{0},jet}_{||};dN/dz^{D^{0},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::ex_hfl_pt, ";p_{T,D^{0}};dN/dp_{T,D^{0}}", {HistType::kTH1F, {{200, 0., 10.}}}); + registryD0.add(histnames::ex_hfl_mass, ";m_{D^{0}};dN/dm_{D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::ex_hfl_eta, ";#eta_{D^{0}};dN/d#eta_{D^{0}}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::ex_hfl_phi, ";#phi_{D^{0}};dN/d#phi_{D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}); + // Labels + auto expCollisionCounter = registryD0.get(HIST(histnames::ex_col)); + expCollisionCounter->GetXaxis()->SetBinLabel(BIN_EX_COLCNTR::AllCollisions, "all"); + expCollisionCounter->GetXaxis()->SetBinLabel(BIN_EX_COLCNTR::Sel8ZCut, "sel8 + zcut"); + + auto expJetCounter = registryD0.get(HIST(histnames::ex_jet)); + expJetCounter->GetXaxis()->SetBinLabel(BIN_EX_JETCNTR::ChargedJets, "Charged jets with D0"); + /* + // Monte Carlo Data Efficiency Histograms + */ + // Definitions + registryD0.add(histnames::mc_eff_col, "N_{coll};", {HistType::kTH1F, {{4, 0., 4.0}}}); + registryD0.add(histnames::mc_eff_jet, "N_{jet};", {HistType::kTH1F, {{4, 0., 4.0}}}); + registryD0.add(histnames::mc_eff_det_jet_pt, ";p_{T,det jet};dN/dp_{T,det jet}", {HistType::kTH1F, {{1000, 0., 50.}}}); + registryD0.add(histnames::mc_eff_det_jet_eta, ";#eta_{det jet};dN/d#eta_{det jet}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::mc_eff_det_jet_phi, ";#phi_{det jet};dN/d#phi_{det jet}", {HistType::kTH1F, {{250, -10., 10.}}}); + registryD0.add(histnames::mc_eff_det_jet_ang, ";angularity_{det jet};dN/d(angularity)", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::mc_eff_det_hfl_pt, ";p_{T,det D^{0}};dN/dp_{T,det D^{0}}", {HistType::kTH1F, {{200, 0., 10.}}}); + registryD0.add(histnames::mc_eff_det_hfl_mass, ";m_{det D^{0}};dN/dm_{det D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::mc_eff_det_hfl_eta, ";#eta_{det D^{0}};dN/d#eta_{det D^{0}}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::mc_eff_det_hfl_phi, ";#phi_{det D^{0}};dN/d#phi_{det D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}); + // Labels + auto mcCollisionCounter = registryD0.get(HIST(histnames::mc_eff_col)); + mcCollisionCounter->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::All, "mccollisions"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::ZCut, "z_cut"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::Matched, "collisions"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::MatchedSel8ZCut, "sel8"); + + auto jetCounter = registryD0.get(HIST(histnames::mc_eff_jet)); + jetCounter->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::ParticleLevelJetInMCCollision, "particle level"); + jetCounter->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::DetectorLevelJetInMCCollision, "detector level"); + jetCounter->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::DetectorLevelJetWithMatchedCandidate, "particle matched jets"); + jetCounter->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::ParticleLevelJetWithMatchedCandidate, "detector matched jets"); + /* + // Monte Carlo Data Histograms + */ + registryD0.add(histnames::mc_col, "N_{coll};", {HistType::kTH1F, {{4, 0., 4.0}}}); + registryD0.add(histnames::mc_jet, "N_{jet};", {HistType::kTH1F, {{4, 0., 4.0}}}); + registryD0.add(histnames::mc_det_jet_pt, ";p_{T,det jet};dN/dp_{T,det jet}", {HistType::kTH1F, {{1000, 0., 50.}}}); + registryD0.add(histnames::mc_det_jet_eta, ";#eta_{det jet};dN/d#eta_{det jet}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::mc_det_jet_phi, ";#phi_{det jet};dN/d#phi_{det jet}", {HistType::kTH1F, {{250, -10., 10.}}}); + registryD0.add(histnames::mc_det_jet_ang, ";angularity_{det jet};dN/d(angularity)", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::mc_det_hfl_pt, ";p_{T,det D^{0}};dN/dp_{T,det D^{0}}", {HistType::kTH1F, {{200, 0., 10.}}}); + registryD0.add(histnames::mc_det_hfl_mass, ";m_{det D^{0}};dN/dm_{det D^{0}}", {HistType::kTH1F, {{1000, 0., 10.}}}); + registryD0.add(histnames::mc_det_hfl_eta, ";#eta_{det D^{0}};dN/d#eta_{det D^{0}}", {HistType::kTH1F, {{250, -5., 5.}}}); + registryD0.add(histnames::mc_det_hfl_phi, ";#phi_{det D^{0}};dN/d#phi_{det D^{0}}", {HistType::kTH1F, {{250, -10., 10.}}}); + // Labels + auto mcCollisionCounter2 = registryD0.get(HIST(histnames::mc_col)); + mcCollisionCounter2->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::All, "mccollisions"); + mcCollisionCounter2->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::ZCut, "z_cut"); + mcCollisionCounter2->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::Matched, "collisions"); + mcCollisionCounter2->GetXaxis()->SetBinLabel(BIN_MC_COLCNTR::MatchedSel8ZCut, "sel8"); + + auto jetCounter2 = registryD0.get(HIST(histnames::mc_jet)); + jetCounter2->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::DetectorLevelJetInMCCollision, "detector level"); + jetCounter2->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::ParticleLevelJetInMCCollision, "particle level"); + jetCounter2->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::DetectorLevelJetWithMatchedCandidate, "particle matched jets"); + jetCounter2->GetXaxis()->SetBinLabel(BIN_MC_JETCNTR::ParticleLevelJetWithMatchedCandidate, "detector matched jets"); + }; + template - void jetCalculateAngularity(T const& jet, U const& /*tracks*/) + float jetCalculateAngularityEXP(T const& jet, U const& /*tracks*/) { - angularity = 0.0; + float tAngularity = 0.0; for (auto& constituent : jet.template tracks_as()) { - angularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent) / (jet.r() / 100.f), alpha); + tAngularity += std::pow(constituent.pt(), kappa) * std::pow(jetutilities::deltaR(jet, constituent) / (jet.r() / 100.f), alpha); + } + tAngularity /= std::pow(jet.pt(), kappa); + return tAngularity; + } + + template + float jetCalculateAngularityMCD( + JetTableMCDConstituent const& jet, + aod::JetTracks const& tracks) + { + float a = 0.f; + for (auto id : jet.tracksIds()) { + auto trk = tracks.iteratorAt(id); + a += std::pow(trk.pt(), kappa) * + std::pow(jetutilities::deltaR(jet, trk) / (jet.r() / 100.f), alpha); + } + return a / std::pow(jet.pt(), kappa); + } + + template + float jetCalculateAngularityMCP( + JetTableMCPConstituent const& jet, + aod::JetParticles const& particles) + { + float a = 0.f; + for (auto id : jet.tracksIds()) { + auto p = particles.iteratorAt(id); + a += std::pow(p.pt(), kappa) * + std::pow(jetutilities::deltaR(jet, p) / (jet.r() / 100.f), alpha); } - angularity /= std::pow(jet.pt(), kappa); + return a / std::pow(jet.pt(), kappa); } - void processDataChargedSubstructure(aod::JetCollision const& collision, - soa::Join const& jets, - aod::CandidatesD0Data const&, aod::JetTracks const& tracks) + template + void analyseDataChargedSubstructure(aod::JetCollision const& collision, + JetChargedTable const& jets, + CandidatesTable const& /*candidates*/, + aod::JetTracks const& tracks, + HistogramRegistry& registry) { // apply event selection and fill histograms for sanity check - registry.fill(HIST("h_collision_counter"), 0.5); + registry.fill(HIST(histnames::ex_col), get_val_from_bin(BIN_EX_COLCNTR::AllCollisions)); if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { return; } - registry.fill(HIST("h_collision_counter"), 1.5); + registry.fill(HIST(histnames::ex_col), get_val_from_bin(BIN_EX_COLCNTR::Sel8ZCut)); // Loop over jets containing D0 candidates for (const auto& jet : jets) { // number of charged jets with D0 - registry.fill(HIST("h_jet_counter"), 0.5); + registry.fill(HIST(histnames::ex_jet), get_val_from_bin(BIN_EX_JETCNTR::ChargedJets)); // obtaining jet 3-vector TVector3 jetVector(jet.px(), jet.py(), jet.pz()); // Loop over D0 candidates associated to the jet - for (const auto& d0Candidate : jet.candidates_as()) { + for (const auto& d0Candidate : jet.template candidates_as()) { // obtaining jet 3-vector TVector3 d0Vector(d0Candidate.px(), d0Candidate.py(), d0Candidate.pz()); @@ -170,24 +497,35 @@ struct JetD0AngSubstructure { // calculating angular distance in eta-phi plane double axisDistance = jetutilities::deltaR(jet, d0Candidate); - jetCalculateAngularity(jet, tracks); + float angularity = jetCalculateAngularityEXP(jet, tracks); // filling histograms - registry.fill(HIST("h_d0_jet_projection"), zParallel); - registry.fill(HIST("h_d0_jet_distance_vs_projection"), axisDistance, zParallel); - registry.fill(HIST("h_d0_jet_distance"), axisDistance); - registry.fill(HIST("h_d0_jet_pt"), jet.pt()); - registry.fill(HIST("h_d0_jet_eta"), jet.eta()); - registry.fill(HIST("h_d0_jet_phi"), jet.phi()); - registry.fill(HIST("h_d0_mass"), d0Candidate.m()); - registry.fill(HIST("h_d0_eta"), d0Candidate.eta()); - registry.fill(HIST("h_d0_phi"), d0Candidate.phi()); - registry.fill(HIST("h_d0_ang"), angularity); // add more axis + registry.fill(HIST(histnames::ex_jet_proj), zParallel); + registry.fill(HIST(histnames::ex_jet_dist_proj), axisDistance, zParallel); + registry.fill(HIST(histnames::ex_jet_dist), axisDistance); + registry.fill(HIST(histnames::ex_jet_pt), jet.pt()); + registry.fill(HIST(histnames::ex_jet_eta), jet.eta()); + registry.fill(HIST(histnames::ex_jet_phi), jet.phi()); + registry.fill(HIST(histnames::ex_jet_ang), angularity); + registry.fill(HIST(histnames::ex_hfl_pt), d0Candidate.pt()); + registry.fill(HIST(histnames::ex_hfl_mass), d0Candidate.m()); + registry.fill(HIST(histnames::ex_hfl_eta), d0Candidate.eta()); + registry.fill(HIST(histnames::ex_hfl_phi), d0Candidate.phi()); // add more axis // filling table ObjJetTable(axisDistance, - jet.pt(), jet.eta(), jet.phi(), jet.tracks_as().size(), angularity, - d0Candidate.pt(), d0Candidate.eta(), d0Candidate.phi(), d0Candidate.m(), d0Candidate.y(), d0Candidate.mlScores()[0], d0Candidate.mlScores()[1], d0Candidate.mlScores()[2]); + jet.pt(), + jet.eta(), + jet.phi(), + jet.template tracks_as().size(), angularity, + d0Candidate.pt(), + d0Candidate.eta(), + d0Candidate.phi(), + d0Candidate.m(), + d0Candidate.y(), + d0Candidate.mlScores()[0], + d0Candidate.mlScores()[1], + d0Candidate.mlScores()[2]); break; // get out of candidates' loop after first HF particle is found in jet } // end of D0 candidates loop @@ -195,7 +533,353 @@ struct JetD0AngSubstructure { } // end of jets loop } // end of process function - PROCESS_SWITCH(JetD0AngSubstructure, processDataChargedSubstructure, "charged HF jet substructure", false); + + template + void analyseMonteCarloEfficiency(MCDJetsPerMCCollissionPreslice const& jetmcdpreslice, + MCPJetsPerMCCollissionPreslice const& jetmcppreslice, + aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + JetTableMCD const& mcdjets, + JetTableMCP const& mcpjets, + CandidatesMCD const& /*mcdCandidates*/, + CandidatesMCP const& /*mcpCandidates*/, + aod::JetTracks const& tracks, + aod::JetParticles const& particles, + HistogramRegistry& registry) + { + for (const auto& mccollision : mccollisions) { + + registry.fill(HIST(histnames::mc_eff_col), get_val_from_bin(BIN_MC_COLCNTR::All)); + // skip collisions outside of |z| < vertexZCut + if (std::abs(mccollision.posZ()) > vertexZCut) { + continue; + } + registry.fill(HIST(histnames::mc_eff_col), get_val_from_bin(BIN_MC_COLCNTR::ZCut)); + + // reconstructed collisions associated to same mccollision + const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mccollision.globalIndex()); + for (const auto& collision : collisionsPerMCCollision) { + + registry.fill(HIST(histnames::mc_eff_col), get_val_from_bin(BIN_MC_COLCNTR::Matched)); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { + continue; + } + registry.fill(HIST(histnames::mc_eff_col), get_val_from_bin(BIN_MC_COLCNTR::MatchedSel8ZCut)); + + // d0 detector level jets associated to the current same collision + const auto d0mcdJetsPerCollision = mcdjets.sliceBy(jetmcdpreslice, collision.globalIndex()); + for (const auto& mcdjet : d0mcdJetsPerCollision) { + + registry.fill(HIST(histnames::mc_eff_jet), get_val_from_bin(BIN_MC_JETCNTR::DetectorLevelJetInMCCollision)); + + // obtain leading HF candidate in jet + auto mcdd0cand = mcdjet.template candidates_first_as(); + + if (mcdjet.has_matchedJetCand()) { + registry.fill(HIST(histnames::mc_eff_jet), get_val_from_bin(BIN_MC_JETCNTR::DetectorLevelJetWithMatchedCandidate)); + } + + // reflection information for storage: D0 = +1, D0bar = -1, neither = 0 + int matchedFrom = 0; + int decayChannel = o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK; + int selectedAs = 0; + + if (mcdd0cand.flagMcMatchRec() == decayChannel) { // matched to D0 on truth level + matchedFrom = 1; + } else if (mcdd0cand.flagMcMatchRec() == -decayChannel) { // matched to D0bar on truth level + matchedFrom = -1; + } + // bitwise AND operation: Checks whether BIT(i) is set, regardless of other bits + if (mcdd0cand.candidateSelFlag() & BIT(0)) { // CandidateSelFlag == BIT(0) -> selected as D0 + selectedAs = 1; + } else if (mcdd0cand.candidateSelFlag() & BIT(1)) { // CandidateSelFlag == BIT(1) -> selected as D0bar + selectedAs = -1; + } + + float angularity = jetCalculateAngularityMCD(mcdjet, tracks); + registry.fill(HIST(histnames::mc_eff_det_jet_pt), mcdjet.pt()); + registry.fill(HIST(histnames::mc_eff_det_jet_eta), mcdjet.eta()); + registry.fill(HIST(histnames::mc_eff_det_jet_phi), mcdjet.phi()); + registry.fill(HIST(histnames::mc_eff_det_jet_ang), angularity); + // Particle Histgrams + registry.fill(HIST(histnames::mc_eff_det_hfl_pt), mcdd0cand.pt()); + registry.fill(HIST(histnames::mc_eff_det_hfl_mass), mcdd0cand.m()); + registry.fill(HIST(histnames::mc_eff_det_hfl_eta), mcdd0cand.eta()); + registry.fill(HIST(histnames::mc_eff_det_hfl_phi), mcdd0cand.phi()); + + mcdJetTable(jetutilities::deltaR(mcdjet, mcdd0cand), + mcdjet.pt(), + mcdjet.eta(), + mcdjet.phi(), + mcdjet.template tracks_as().size(), // detector level jet + angularity, + mcdd0cand.pt(), + mcdd0cand.eta(), + mcdd0cand.phi(), + mcdd0cand.m(), + mcdd0cand.y(), + (mcdd0cand.originMcRec() == RecoDecay::OriginType::Prompt), // detector level D0 candidate + mcdjet.has_matchedJetCand(), + mcdd0cand.mlScores()[0], + mcdd0cand.mlScores()[1], + mcdd0cand.mlScores()[2], // // Machine Learning PID scores: background, prompt, non-prompt + matchedFrom, selectedAs); // D0 = +1, D0bar = -1, neither = 0 + } + } // end of reconstructed collisions loop (detector level mc collisions matching with real collisions) + + // d0 particle level jets associated to same mccollision + const auto d0mcpJetsPerMCCollision = mcpjets.sliceBy(jetmcppreslice, mccollision.globalIndex()); + for (const auto& mcpjet : d0mcpJetsPerMCCollision) { + + registry.fill(HIST(histnames::mc_eff_jet), get_val_from_bin(BIN_MC_JETCNTR::ParticleLevelJetInMCCollision)); + + // obtain leading HF particle in jet + auto mcpd0cand = mcpjet.template candidates_first_as(); + + if (mcpjet.has_matchedJetCand()) { + registry.fill(HIST(histnames::mc_eff_jet), get_val_from_bin(BIN_MC_JETCNTR::ParticleLevelJetWithMatchedCandidate)); + } + float angularity = jetCalculateAngularityMCP(mcpjet, particles); + + registry.fill(HIST(histnames::mc_eff_par_jet_pt), mcpjet.pt()); + registry.fill(HIST(histnames::mc_eff_par_jet_eta), mcpjet.eta()); + registry.fill(HIST(histnames::mc_eff_par_jet_phi), mcpjet.phi()); + registry.fill(HIST(histnames::mc_eff_par_jet_ang), angularity); + // Particle Histgrams + registry.fill(HIST(histnames::mc_eff_par_hfl_pt), mcpd0cand.pt()); + + registry.fill(HIST(histnames::mc_eff_par_hfl_eta), mcpd0cand.eta()); + registry.fill(HIST(histnames::mc_eff_par_hfl_phi), mcpd0cand.phi()); + // store data in MC detector level table (calculate angular distance in eta-phi plane on the fly) + mcpJetTable(jetutilities::deltaR(mcpjet, mcpd0cand), + mcpjet.pt(), + mcpjet.eta(), + mcpjet.phi(), + mcpjet.template tracks_as().size(), // particle level jet + angularity, + mcpd0cand.pt(), + mcpd0cand.eta(), + mcpd0cand.phi(), + // mcpd0cand.m(), + mcpd0cand.y(), + (mcpd0cand.originMcGen() == RecoDecay::OriginType::Prompt), // particle level D0 + mcpjet.has_matchedJetCand()); + } // End of particle level jets loop, related to detector level collisions and jets. + } + } + + template + void analyseMonteCarlo(MCPJetsPerMCCollissionPreslice jetmcpreslice, + aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + JetTableMCD const& /*mcdjets*/, + JetTableMCP const& mcpjets, + CandidatesMCD const& /*mcdCandidates*/, + CandidatesMCP const& /*mcpCandidates*/, + aod::JetTracks const& jettracks, + aod::JetParticles const& jetparticles, + HistogramRegistry& registry) + { + for (const auto& mccollision : mccollisions) { + registry.fill(HIST(histnames::mc_col), get_val_from_bin(BIN_MC_COLCNTR::All)); + // skip collisions outside of |z| < vertexZCut + if (std::abs(mccollision.posZ()) > vertexZCut) { + continue; + } + registry.fill(HIST(histnames::mc_col), get_val_from_bin(BIN_MC_COLCNTR::ZCut)); + + // hf particle level jets associated to same mccollision + const JetD0MCPTable mcpJetsPerMCCollision = mcpjets.sliceBy(jetmcpreslice, mccollision.globalIndex()); + for (const auto& mcpjet : mcpJetsPerMCCollision) { + + registry.fill(HIST(histnames::mc_jet), get_val_from_bin(BIN_MC_JETCNTR::ParticleLevelJetInMCCollision)); + + // obtain leading HF particle in jet + auto mcpcand = mcpjet.template candidates_first_as(); + + if (mcpjet.has_matchedJetCand()) { + registry.fill(HIST(histnames::mc_jet), get_val_from_bin(BIN_MC_JETCNTR::ParticleLevelJetWithMatchedCandidate)); + + // loop over detector level matched to current particle level + for (const auto& mcdjet : mcpjet.template matchedJetCand_as()) { + registry.fill(HIST(histnames::mc_jet), get_val_from_bin(BIN_MC_JETCNTR::DetectorLevelJetWithMatchedCandidate)); + + // apply collision sel8 selection on detector level jet's collision + const auto& collision = collisions.iteratorAt(mcdjet.collisionId()); + registry.fill(HIST(histnames::mc_col), get_val_from_bin(BIN_MC_COLCNTR::Matched)); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) { + continue; + } + registry.fill(HIST(histnames::mc_col), get_val_from_bin(BIN_MC_COLCNTR::MatchedSel8ZCut)); + + // obtain leading HF candidate in jet + auto mcdcand = mcdjet.template candidates_first_as(); + + // reflection information for storage: HF = +1, HFbar = -1, neither = 0 + int matchedFrom = 0; + int decayChannel = 0; + if (jethfutilities::isD0Table()) { + decayChannel = o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK; + } else if (jethfutilities::isLcTable()) { + decayChannel = o2::hf_decay::hf_cand_3prong::DecayChannelMain::LcToPKPi; + } + int selectedAs = 0; + + if (mcdcand.flagMcMatchRec() == decayChannel) { // matched to HF on truth level + matchedFrom = 1; + } else if (mcdcand.flagMcMatchRec() == -decayChannel) { // matched to HFbar on truth level + matchedFrom = -1; + } + // bitwise AND operation: Checks whether BIT(i) is set, regardless of other bits + if (mcdcand.candidateSelFlag() & BIT(0)) { // CandidateSelFlag == BIT(0) -> selected as HF + selectedAs = 1; + } else if (mcdcand.candidateSelFlag() & BIT(1)) { // CandidateSelFlag == BIT(1) -> selected as HFbar + selectedAs = -1; + } + + float mcpAngularity = jetCalculateAngularityMCP(mcpjet, jetparticles); + float mcdAngularity = jetCalculateAngularityMCD(mcdjet, jettracks); + + // filling histograms + // registry.fill(HIST(histnames::h_exp_d0_jet_projection), zParallel); + // registry.fill(HIST(histnames::h_exp_d0_jet_distance_vs_projection), axisDistance, zParallel); + // registry.fill(HIST(histnames::h_exp_d0_jet_distance), axisDistance); + // registry.fill(HIST(histnames::h_exp_d0_jet_pt), jet.pt()); + // registry.fill(HIST(histnames::h_exp_d0_jet_eta), jet.eta()); + // registry.fill(HIST(histnames::h_exp_d0_jet_phi), jet.phi()); + // registry.fill(HIST(histnames::h_exp_d0_jet_ang), angularity); + // Jet Histograms + + registry.fill(HIST(histnames::mc_det_jet_pt), mcdjet.pt()); + registry.fill(HIST(histnames::mc_det_jet_eta), mcdjet.eta()); + registry.fill(HIST(histnames::mc_det_jet_phi), mcdjet.phi()); + registry.fill(HIST(histnames::mc_det_jet_ang), mcdAngularity); + // Particle Histgrams + registry.fill(HIST(histnames::mc_det_hfl_pt), mcdcand.pt()); + registry.fill(HIST(histnames::mc_det_hfl_mass), mcdcand.m()); + registry.fill(HIST(histnames::mc_det_hfl_eta), mcdcand.eta()); + registry.fill(HIST(histnames::mc_det_hfl_phi), mcdcand.phi()); // add more axis + + // store matched particle and detector level data in one single table (calculate angular distance in eta-phi plane on the fly) + matchJetTable(jetutilities::deltaR(mcpjet, mcpcand), + mcpjet.pt(), + mcpjet.eta(), + mcpjet.phi(), + mcpjet.template tracks_as().size(), // particle level jet + mcpAngularity, + mcpcand.pt(), + mcpcand.eta(), + mcpcand.phi(), + mcpcand.y(), + (mcpcand.originMcGen() == RecoDecay::OriginType::Prompt), // particle level HF + jetutilities::deltaR(mcdjet, mcdcand), + mcdjet.pt(), + mcdjet.eta(), + mcdjet.phi(), + mcdjet.template tracks_as().size(), // detector level jet + mcdAngularity, + mcdcand.pt(), + mcdcand.eta(), + mcdcand.phi(), + mcdcand.m(), + mcdcand.y(), + (mcdcand.originMcRec() == RecoDecay::OriginType::Prompt), // detector level HF + mcdcand.mlScores()[0], + mcdcand.mlScores()[1], + mcdcand.mlScores()[2], // Machine Learning PID scores: background, prompt, non-prompt + matchedFrom, + selectedAs); // HF = +1, HFbar = -1, neither = 0 + } + } else { + // store matched particle and detector level data in one single table (calculate angular distance in eta-phi plane on the fly) + float mcpAngularity = jetCalculateAngularityMCP(mcpjet, jetparticles); + // float mcpAngularity = 0.; + matchJetTable(jetutilities::deltaR(mcpjet, mcpcand), + mcpjet.pt(), + mcpjet.eta(), + mcpjet.phi(), + mcpjet.template tracks_as().size(), // particle level jet + mcpAngularity, + + mcpcand.pt(), + mcpcand.eta(), + mcpcand.phi(), + mcpcand.y(), + (mcpcand.originMcGen() == RecoDecay::OriginType::Prompt), // particle level HF + -2, -2, -2, -2, -2, // detector level jet + -2, -2, -2, -2, -2, -2, -2, // detector level HF + -2, -2, -2, // Machine Learning PID scores: background, prompt, non-prompt + -2, -2); // HF = +1, HFbar = -1, neither = 0 + } + } // end of mcpjets loop + } // end of mccollisions loop + }; + + void processataChargedSubstructureD0(aod::JetCollision const& collision, + JetChargedTableD0 const& jets, + aod::CandidatesD0Data const& candidates, + aod::JetTracks const& tracks) + { + analyseDataChargedSubstructure(collision, jets, candidates, tracks, registryD0); + } + + PROCESS_SWITCH(JetD0AngSubstructure, processataChargedSubstructureD0, "charged HF jet substructure", false); + + void processMonteCarloEfficiencyD0(aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + JetD0MCDTable const& mcdjets, + JetD0MCPTable const& mcpjets, + aod::CandidatesD0MCD const& mcdCandidates, + aod::CandidatesD0MCP const& mcpCandidates, + aod::JetTracks const& jettracks, + aod::JetParticles const& jetparticles) + { + analyseMonteCarloEfficiency, + Preslice, + JetD0MCDTable, + JetD0MCPTable, + aod::CandidatesD0MCD, + aod::CandidatesD0MCP>( + d0MCDJetsPerEXPCollisionPreslice, d0MCPJetsPerMCCollisionPreslice, mccollisions, collisions, mcdjets, + mcpjets, mcdCandidates, mcpCandidates, jettracks, jetparticles, registryD0); + } + + PROCESS_SWITCH(JetD0AngSubstructure, processMonteCarloEfficiencyD0, "non-matched and matched MC D0 and jets", false); + + void processMonteCarloD0(aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + JetD0MCDTable const& mcdjets, + JetD0MCPTable const& mcpjets, + aod::CandidatesD0MCD const& mcdCandidates, + aod::CandidatesD0MCP const& mcpCandidates, + aod::JetTracks const& jettracks, + aod::JetParticles const& jetparticles) + { + analyseMonteCarlo, + JetD0MCDTable, + JetD0MCPTable, + aod::CandidatesD0MCD, + aod::CandidatesD0MCP>( + d0MCPJetsPerMCCollisionPreslice, mccollisions, collisions, mcdjets, + mcpjets, mcdCandidates, mcpCandidates, jettracks, jetparticles, registryD0); + } + + PROCESS_SWITCH(JetD0AngSubstructure, processMonteCarloD0, "Store all simulated D0 jets information with matched candidate (if any found)", false); }; // Workflow definition -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-d0-ang-substructure"})}; } +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask( + cfgc, + TaskName{"jet-d0-ang-substructure"})}; +}