diff --git a/PWGJE/Tasks/jetDsSpecSubs.cxx b/PWGJE/Tasks/jetDsSpecSubs.cxx index 12253f021c6..59fbf1ff4fb 100644 --- a/PWGJE/Tasks/jetDsSpecSubs.cxx +++ b/PWGJE/Tasks/jetDsSpecSubs.cxx @@ -13,10 +13,15 @@ /// \brief Ds-tagged jet analysis with substructure histogram outputs /// \author Monalisa Melo , Universidade de São Paulo +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGJE/Core/JetDerivedDataUtilities.h" #include "PWGJE/Core/JetUtilities.h" #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" +#include "PWGJE/DataModel/JetSubstructure.h" + +#include "Common/Core/RecoDecay.h" #include #include @@ -28,6 +33,7 @@ #include #include +#include #include #include @@ -40,38 +46,20 @@ using namespace o2::framework; using namespace o2::framework::expressions; struct JetDsSpecSubs { - 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, {{200, 0., 200.}}}}, - {"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, {{200, 1.7, 2.1}, {200, 0., 100.}, {200, 0., 2.}, {200, 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.}}}}, - {"h_ds_jet_mass", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}}, - {"h_ds_jet_lambda11", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, - {"h_ds_jet_lambda12", ";#lambda_{2}^{1} (thrust);entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, - {"h_ds_jet_girth", ";g (#equiv #lambda_{1}^{1}R);entries", {HistType::kTH1F, {{200, 0., 0.5}}}}, - {"h2_dsjet_pt_lambda11", ";#it{p}_{T,jet} (GeV/#it{c});#lambda_{1}^{1}", {HistType::kTH2F, {{100, 0., 100.}, {200, 0., 1.0}}}}, - {"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}}}}, - }}; + //================== + // Type definitions + //================== + + using DsCandidatesData = aod::CandidatesDsData; + using DsCandidatesMCD = aod::CandidatesDsMCD; + using DsCandidatesMCP = aod::CandidatesDsMCP; + + using DsDataJets = soa::Join; + using DsMCDJets = soa::Join; + using DsMCPJets = soa::Join; + + // 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"}; @@ -79,30 +67,101 @@ struct JetDsSpecSubs { Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; + // internals std::vector eventSelectionBits; int trackSelection = -1; + // Filters + Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); + Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; + + //============= + // Histograms + //============= + + HistogramRegistry registry{ + "registry", + { + + {"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}}, + {"h_collision_counter", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}}, + + {"h_track_pt", ";#it{p}_{T,track};entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_track_eta", ";#eta_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"h_track_phi", ";#varphi_{track};entries", {HistType::kTH1F, {{80, -1., 7.}}}}, + + {"h_jet_counter", ";type;counts", {HistType::kTH1F, {{2, 0., 2.}}}}, + + {"h_jet_pt", ";#it{p}_{T,jet};entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_jet_eta", ";#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"h_jet_phi", ";#phi_{jet};entries", {HistType::kTH1F, {{80, -1., 7.}}}}, + + {"h_ds_mass", ";m_{D_{S}};entries", {HistType::kTH1F, {{400, 1.7, 2.2}}}}, + {"h_ds_pt", ";p_{T,D_{S}};entries", {HistType::kTH1F, {{200, 0., 100.}}}}, + {"h_ds_eta", ";#eta_{D_{S}};entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"h_ds_phi", ";#phi_{D_{S}};entries", {HistType::kTH1F, {{100, -1., 7.}}}}, + + {"h_ds_jet_projection", ";z_{||}^{D_{S},jet};entries", {HistType::kTH1F, {{200, 0., 2.}}}}, + {"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};entries", {HistType::kTH1F, {{200, 0., 1.}}}}, + {"h_ds_jet_lambda11", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.}}}}, + {"h_ds_jet_lambda12", ";#lambda_{2}^{1};entries", {HistType::kTH1F, {{200, 0., 1.}}}}, + {"h_ds_jet_mass", ";m_{jet};entries", {HistType::kTH1F, {{200, 0., 50.}}}}, + + {"hSparse_ds", ";m_{D_{S}};p_{T,D_{S}};p_{T,jet};z_{||};#DeltaR", {HistType::kTHnSparseF, {{200, 1.7, 2.2}, {200, 0., 100.}, {200, 0., 100.}, {200, 0., 2.}, {200, 0., 1.}}}}, + + {"h2_response_jet_pt", ";p_{T}^{det};p_{T}^{part}", {HistType::kTH2F, {{200, 0., 100.}, {200, 0., 100.}}}}, + + {"h2_response_lambda11", ";#lambda_{1}^{1,det};#lambda_{1}^{1,part}", {HistType::kTH2F, {{200, 0., 1.}, {200, 0., 1.}}}}}}; + + //======== + // INIT + //======== + + void init(InitContext const&) + { + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + + auto h = registry.get(HIST("h_jet_counter")); + + h->GetXaxis()->SetBinLabel(1, "All jets"); + h->GetXaxis()->SetBinLabel(2, "Ds-tagged jets"); + } + + //=============== + // Lambda compute + //=============== + template - float computeLambda(JET const& jet, TRACKS const& tracks, float a, float k) + float computeLambda(JET const& jet, + TRACKS const& tracks, + float alpha, + float kappa) { if (jet.pt() <= 0.f) { return -1.f; } + float sum = 0.f; for (auto const& trk : tracks) { const float dr = jetutilities::deltaR(jet, trk); - sum += std::pow(trk.pt(), k) * std::pow(dr, a); + sum += std::pow(trk.pt(), kappa) * std::pow(dr, alpha); } const float jetR = jet.r() / 100.f; - const float denom = std::pow(jet.pt(), k) * std::pow(jetR, a); + + const float denom = std::pow(jet.pt(), kappa) * std::pow(jetR, alpha); if (denom <= 0.f) { return -1.f; } return sum / denom; } + //================= + // Jet Mass compute + //================= + template - float computeJetMassFromTracksMass(TRACKS const& tracks) + float computeJetMass(TRACKS const& tracks) { double sumPx = 0.0, sumPy = 0.0, sumPz = 0.0, sumE = 0.0; @@ -126,22 +185,15 @@ struct JetDsSpecSubs { return (m2 > 0.0) ? static_cast(std::sqrt(m2)) : 0.f; } - void init(o2::framework::InitContext&) - { - eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); - trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); - - auto h = registry.get(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); - Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; + //============== + // Collision QA + //============== - void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks) + void processCollisions(aod::JetCollision const& collision, + aod::JetTracks const& tracks) { registry.fill(HIST("h_collisions"), 0.5); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } @@ -149,121 +201,223 @@ struct JetDsSpecSubs { registry.fill(HIST("h_collisions"), 1.5); 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(JetDsSpecSubs, processCollisions, "process JE collisions", false); - void processDataCharged(soa::Filtered::iterator const& collision, - soa::Filtered const& jets) + PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "collision QA", false); + + //=============== + // DATA analysis + //=============== + + template + void analyseData(JETS const& jets, + CANDS const&) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; - } for (const auto& jet : jets) { + + registry.fill(HIST("h_jet_counter"), 0.5); + registry.fill(HIST("h_jet_pt"), jet.pt()); registry.fill(HIST("h_jet_eta"), jet.eta()); registry.fill(HIST("h_jet_phi"), jet.phi()); + + auto jetTracks = jet.template tracks_as(); + + const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f); + + const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f); + + const float mjet = computeJetMass(jetTracks); + + bool hasDs = false; + + TVector3 jetVector(jet.px(), jet.py(), jet.pz()); + + for (const auto& ds : jet.template candidates_as()) { + + hasDs = true; + + TVector3 dsVector(ds.px(), ds.py(), ds.pz()); + + const float deltaR = jetutilities::deltaR(jet, ds); + + const float zParallel = (jetVector * dsVector) / (jetVector * jetVector); + + registry.fill(HIST("h_ds_mass"), ds.m()); + registry.fill(HIST("h_ds_pt"), ds.pt()); + registry.fill(HIST("h_ds_eta"), ds.eta()); + registry.fill(HIST("h_ds_phi"), ds.phi()); + + registry.fill(HIST("h_ds_jet_projection"), zParallel); + + registry.fill(HIST("h_ds_jet_distance"), deltaR); + + registry.fill(HIST("hSparse_ds"), + ds.m(), + ds.pt(), + jet.pt(), + zParallel, + deltaR); + } + + if (hasDs) { + + registry.fill(HIST("h_jet_counter"), 1.5); + + registry.fill(HIST("h_ds_jet_lambda11"), lambda11); + registry.fill(HIST("h_ds_jet_lambda12"), lambda12); + registry.fill(HIST("h_ds_jet_mass"), mjet); + } } } - PROCESS_SWITCH(JetDsSpecSubs, processDataCharged, "charged jets in data", false); - void processDataChargedSubstructure(aod::JetCollision const& collision, - soa::Join const& jets, - aod::CandidatesDsData const&, - aod::JetTracks const&) - { - registry.fill(HIST("h_collision_counter"), 2.0); - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || - !(std::abs(collision.posZ()) < vertexZCut)) { - return; - } - registry.fill(HIST("h_collision_counter"), 3.0); + //============== + // MCD analysis + //============== + template + void analyseMCD(JETS const& jets) + { for (const auto& jet : jets) { - registry.fill(HIST("h_jet_counter"), 0.5); + auto jetTracks = + jet.template tracks_as(); - bool hasDsCandidate = false; + const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f); + + const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f); + + const float mjet = computeJetMass(jetTracks); TVector3 jetVector(jet.px(), jet.py(), jet.pz()); - // Compute jet-level quantities once (independent of Ds candidates) - auto jetTracks = jet.tracks_as(); + for (const auto& ds : jet.template candidates_as()) { + + TVector3 dsVector(ds.px(), ds.py(), ds.pz()); + + const float deltaR = jetutilities::deltaR(jet, ds); + + const float zParallel = (jetVector * dsVector) / (jetVector * jetVector); + + registry.fill(HIST("hSparse_ds"), + ds.m(), + ds.pt(), + jet.pt(), + zParallel, + deltaR); + + registry.fill(HIST("h_ds_jet_lambda11"), lambda11); + registry.fill(HIST("h_ds_jet_lambda12"), lambda12); + registry.fill(HIST("h_ds_jet_mass"), mjet); + } + } + } + + //================== + // MC particle level + //================== + + template + void analyseMCP(JETS const& jets) + { + for (const auto& jet : jets) { + + auto jetTracks = jet.template tracks_as(); 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 jetR = jet.r() / 100.f; - const float girth = (lambda11 >= 0.f) ? (lambda11 * jetR) : -1.f; + const float lambda12 = computeLambda(jet, jetTracks, 2.f, 1.f); - // Loop over Ds candidates (particle level) - for (const auto& dsCandidate : jet.candidates_as()) { + const float mjet = computeJetMass(jetTracks); - hasDsCandidate = true; + TVector3 jetVector(jet.px(), jet.py(), jet.pz()); - TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz()); + for (const auto& ds : jet.template candidates_as()) { - // zParallel defined as longitudinal momentum fraction along the jet axis - const double zParallel = (jetVector * dsVector) / (jetVector * jetVector); - const float axisDistance = jetutilities::deltaR(jet, dsCandidate); + TVector3 dsVector(ds.px(), ds.py(), ds.pz()); - // --- 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); + const float deltaR = jetutilities::deltaR(jet, ds); - registry.fill(HIST("h_ds_mass"), dsCandidate.m()); - registry.fill(HIST("h_ds_eta"), dsCandidate.eta()); - registry.fill(HIST("h_ds_phi"), dsCandidate.phi()); + const float zParallel = (jetVector * dsVector) / (jetVector * jetVector); - // Main THnSparse: invariant mass, pT, z, and DeltaR registry.fill(HIST("hSparse_ds"), - dsCandidate.m(), - dsCandidate.pt(), + ds.pt(), + jet.pt(), zParallel, - axisDistance); + deltaR); + + registry.fill(HIST("h_ds_jet_lambda11"), lambda11); + registry.fill(HIST("h_ds_jet_lambda12"), lambda12); + registry.fill(HIST("h_ds_jet_mass"), mjet); } + } + } - // Jet-level quantities (filled once per jet containing at least one Ds) - if (hasDsCandidate) { + //============== + // DATA process + //============== - registry.fill(HIST("h_jet_counter"), 1.5); + void processDataChargedSubstructure( + aod::JetCollision const& collision, + soa::Filtered const& jets, + DsCandidatesData const& candidates, + aod::JetTracks const&) + { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } - // 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()); + analyseData(jets, candidates); + } - // 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); - } + PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "data charged jets", false); - if (lambda12 >= 0.f) { - registry.fill(HIST("h_ds_jet_lambda12"), lambda12); - registry.fill(HIST("h2_dsjet_pt_lambda12"), jet.pt(), lambda12); - } + //============== + // MCD process + //============== - registry.fill(HIST("h_ds_jet_mass"), mjet); - registry.fill(HIST("h2_dsjet_pt_mass"), jet.pt(), mjet); + void processMCDChargedSubstructure( + aod::JetCollision const& collision, + soa::Filtered const& jets, + aod::JetTracks const&) + { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } - if (girth >= 0.f) { - registry.fill(HIST("h_ds_jet_girth"), girth); - registry.fill(HIST("h2_dsjet_pt_girth"), jet.pt(), girth); - } - } + analyseMCD(jets); + } + + PROCESS_SWITCH(JetDsSpecSubs, processMCDChargedSubstructure, "MC detector level", false); + + //============== + // MCP process + //============== + + void processMCPChargedSubstructure( + aod::JetCollision const& collision, + soa::Filtered const& jets, + aod::JetTracks const&) + { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; } + + analyseMCP(jets); } - PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false); + + PROCESS_SWITCH(JetDsSpecSubs, processMCPChargedSubstructure, "MC particle level", false); }; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)};