diff --git a/PWGJE/Tasks/jetDsSpecSubs.cxx b/PWGJE/Tasks/jetDsSpecSubs.cxx index 59fbf1ff4fb..c5e45f0a899 100644 --- a/PWGJE/Tasks/jetDsSpecSubs.cxx +++ b/PWGJE/Tasks/jetDsSpecSubs.cxx @@ -13,6 +13,7 @@ /// \brief Ds-tagged jet analysis with substructure histogram outputs /// \author Monalisa Melo , Universidade de São Paulo +#include "PWGHF/Core/DecayChannels.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGJE/Core/JetDerivedDataUtilities.h" @@ -22,18 +23,19 @@ #include "PWGJE/DataModel/JetSubstructure.h" #include "Common/Core/RecoDecay.h" +#include "Common/DataModel/TrackSelectionTables.h" #include #include #include #include #include +#include #include #include #include #include -#include #include #include @@ -45,6 +47,27 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +consteval float getValFromBin(int bin) +{ + return static_cast(bin) - 0.5f; +} + +enum BinExpColCntr { AllCollisions = 1, + Sel8ZCut = 2 }; + +enum BinExpJetCntr { ChargedJets = 1 }; +enum BinMCColCntr { All = 1, + ZCut = 2, + Matched = 3, + MatchedSel8ZCut = 4 +}; + +enum BinMCJetCntr { DetectorLevelJetInMCCollision = 1, + ParticleLevelJetInMCCollision = 2, + DetectorLevelJetWithMatchedCandidate = 3, + ParticleLevelJetWithMatchedCandidate = 4 +}; + struct JetDsSpecSubs { //================== @@ -56,8 +79,13 @@ struct JetDsSpecSubs { using DsCandidatesMCP = aod::CandidatesDsMCP; using DsDataJets = soa::Join; - using DsMCDJets = soa::Join; - using DsMCPJets = soa::Join; + using DsMCDJets = soa::Join; + // using DsMCPJets = soa::Join; + + // Slices for access to proper HF MCD jet collision that is associated to MCCollision + PresliceUnsorted collisionsPerMCCollisionPreslice = aod::jmccollisionlb::mcCollisionId; + Preslice dsMCDJetsPerEXPCollisionPreslice = aod::jet::collisionId; + // Preslice dsMCPJetsPerMCCollisionPreslice = aod::jet::mcCollisionId; // Configurables Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; @@ -72,8 +100,8 @@ struct JetDsSpecSubs { 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; + // Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); + // Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; //============= // Histograms @@ -82,37 +110,77 @@ struct JetDsSpecSubs { 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_collision_counter_data", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}}, + {"h_collision_counter_mcd", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}}, + {"h_collision_counter_mcp", ";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.}}}}}}; - + // Data histograms + {"h_dsjet_counter_data", ";type;counts", {HistType::kTH1F, {{3, 0., 3.}}}}, + + {"h_jet_pt_data", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_jet_eta_data", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_jet_phi_data", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + + {"h_ds_mass_data", ";m_{D_{S}} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{1000, 0., 6.}}}}, + {"h_ds_pt_data", ";#it{p}_{T,D_{S}} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + {"h_ds_eta_data", ";#eta_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_phi_data", ";#phi_{D_{S}};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + + {"h_ds_jet_pt_data", ";#it{p}_{T,D_{S} jet} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + {"h_ds_jet_eta_data", ";#eta_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_jet_phi_data", ";#phi_{D_{S} jet};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + + {"h_ds_jet_projection_data", ";z^{D_{S},jet}_{||};entries", {HistType::kTH1F, {{1000, 0., 2.}}}}, + {"h_ds_jet_distance_data", ";#DeltaR_{D_{S},jet};entries", {HistType::kTH1F, {{1000, 0., 1.}}}}, + {"h_ds_jet_mass_data", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}}, + {"h_ds_jet_lambda11_data", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"h_ds_jet_lambda12_data", ";#lambda_{2}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"hSparse_ds_data", ";m_{D_{S}};#it{p}_{T,D_{S}};#it{p}_{T,jet};z^{D_{S},jet}_{||};#DeltaR_{D_{S},jet}", {HistType::kTHnSparseF, {{60, 1.7, 2.1}, {60, 0., 100.}, {60, 0., 100.}, {60, 0., 2.}, {60, 0., 1.0}}}}, + + // MC detector-level histograms + {"McEffCol", "N_{collisions};", {HistType::kTH1F, {{4, 0., 4.0}}}}, + {"McEffJet", "N_{jet};", {HistType::kTH1F, {{4, 0., 4.0}}}}, + {"h_jet_pt_mcd", "detector-level jet pT;#it{p}_{T,jet}^{det} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_jet_eta_mcd", "detector-level jet #eta;#eta_{jet}^{det};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_jet_phi_mcd", "detector-level jet #phi;#phi_{jet}^{det};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + + {"h_ds_mass_mcd", ";m_{D_{S}}^{rec} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{1000, 0., 6.}}}}, + {"h_ds_pt_mcd", ";#it{p}_{T,D_{S}}^{det} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + {"h_ds_eta_mcd", ";#eta_{D_{S}}^{det};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_phi_mcd", ";#phi_{D_{S}}^{det};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + + {"h_ds_jet_pt_mcd", ";#it{p}_{T,D_{S} jet}^{det} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + {"h_ds_jet_eta_mcd", ";#eta_{D_{S} jet}^{det};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + {"h_ds_jet_phi_mcd", ";#phi_{D_{S} jet}^{det};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + {"h_ds_jet_projection_mcd", ";z^{D_{S},jet}_{||, det};entries", {HistType::kTH1F, {{1000, 0., 2.}}}}, + {"h_ds_jet_distance_mcd", ";#DeltaR_{D_{S},jet}^{det};entries", {HistType::kTH1F, {{1000, 0., 1.}}}}, + {"h_ds_jet_mass_mcd", ";m_{jet}^{ch, det} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}}, + {"h_ds_jet_lambda11_mcd", ";#lambda_{1}^{1, det};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"h_ds_jet_lambda12_mcd", ";#lambda_{2}^{1, det};entries", {HistType::kTH1F, {{200, 0., 1.0}}}}, + {"hSparse_ds_mcd", ";m_{D_{S}}^{rec};#it{p}_{T,D_{S}}^{det};#it{p}_{T,jet}^{det};z^{D_{S},jet}_{||, det};#DeltaR_{D_{S},jet}^{det}", {HistType::kTHnSparseF, {{60, 1.7, 2.1}, {60, 0., 100.}, {60, 0., 100.}, {60, 0., 2.}, {60, 0., 1.0}}}}, + + // MC particle-level histograms + // {"h_dsjet_counter_mcp", ";type;counts", {HistType::kTH1F, {{3, 0., 3.}}}}, + // {"h_jet_pt_mcp", "particle-level jet pT;#it{p}_{T,jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + // {"h_jet_eta_mcp", "particle-level jet #eta;#eta_{jet}^{part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + // {"h_jet_phi_mcp", "particle-level jet #phi;#phi_{jet}^{part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + + // {"h_ds_pt_mcp", ";#it{p}_{T,D_{S}}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + // {"h_ds_eta_mcp", ";#eta_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + // {"h_ds_phi_mcp", ";#phi_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + + // {"h_ds_jet_pt_mcp", ";#it{p}_{T,D_{S} jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}}, + // {"h_ds_jet_eta_mcp", ";#eta_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}}, + // {"h_ds_jet_phi_mcp", ";#phi_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}}, + // {"h_ds_jet_projection_mcp", ";z^{D_{S},jet}_{||, part};entries", {HistType::kTH1F, {{1000, 0., 2.}}}}, + // {"h_ds_jet_distance_mcp", ";#DeltaR_{D_{S},jet}^{part};entries", {HistType::kTH1F, {{1000, 0., 1.}}}}, + // {"hSparse_ds_mcp", ";#it{p}_{T,D_{S}}^{part};#it{p}_{T,jet}^{part};z^{D_{S},jet}_{||, part};#DeltaR_{D_{S},jet}^{part}", {HistType::kTHnSparseF, {{60, 0., 100.}, {60, 0., 100.}, {60, 0., 2.}, {60, 0., 1.0}}}}, + }}; //======== // INIT //======== @@ -122,21 +190,30 @@ struct JetDsSpecSubs { 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"); + auto hData = registry.get(HIST("h_dsjet_counter_data")); + hData->GetXaxis()->SetBinLabel(1, "Ds-jet entries"); + hData->GetXaxis()->SetBinLabel(2, "Ds candidates"); + hData->GetXaxis()->SetBinLabel(3, "Ds jets with >=1 cand."); + + // Labels + auto mcCollisionCounter = registry.get(HIST("McEffCol")); + mcCollisionCounter->GetXaxis()->SetBinLabel(BinMCColCntr::All, "mccollisions"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BinMCColCntr::ZCut, "z_cut"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BinMCColCntr::Matched, "collisions"); + mcCollisionCounter->GetXaxis()->SetBinLabel(BinMCColCntr::MatchedSel8ZCut, "sel8"); + + auto jetCounter = registry.get(HIST("McEffJet")); + jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::ParticleLevelJetInMCCollision, "particle level"); + jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetInMCCollision, "detector level"); + jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetWithMatchedCandidate, "particle matched jets"); + jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::ParticleLevelJetWithMatchedCandidate, "detector matched jets"); } - //=============== // Lambda compute //=============== template - float computeLambda(JET const& jet, - TRACKS const& tracks, - float alpha, - float kappa) + float computeLambda(JET const& jet, TRACKS const& tracks, float alpha, float kappa) { if (jet.pt() <= 0.f) { return -1.f; @@ -189,8 +266,7 @@ struct JetDsSpecSubs { // 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); @@ -203,7 +279,7 @@ struct JetDsSpecSubs { for (auto const& track : tracks) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; + return; } registry.fill(HIST("h_track_pt"), track.pt()); @@ -211,213 +287,207 @@ struct JetDsSpecSubs { registry.fill(HIST("h_track_phi"), track.phi()); } } - PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "collision QA", false); - //=============== - // DATA analysis - //=============== + //============== + // DATA process + //============== - template - void analyseData(JETS const& jets, - CANDS const&) + void processDataChargedSubstructure(aod::JetCollision const& collision, soa::Join const& jets, + aod::CandidatesDsData const&, aod::JetTracks const&) { - 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()); + registry.fill(HIST("h_collision_counter_data"), 2.0); - 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); - } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || + !(std::abs(collision.posZ()) < vertexZCut)) { + return; } - } - //============== - // MCD analysis - //============== + registry.fill(HIST("h_collision_counter_data"), 3.0); - template - void analyseMCD(JETS const& jets) - { for (const auto& jet : jets) { + registry.fill(HIST("h_dsjet_counter_data"), 0.5); // DsChargedJets entries - auto jetTracks = - jet.template tracks_as(); + registry.fill(HIST("h_jet_pt_data"), jet.pt()); + registry.fill(HIST("h_jet_eta_data"), jet.eta()); + registry.fill(HIST("h_jet_phi_data"), jet.phi()); - const float lambda11 = computeLambda(jet, jetTracks, 1.f, 1.f); + auto jetTracks = jet.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); TVector3 jetVector(jet.px(), jet.py(), jet.pz()); - for (const auto& ds : jet.template candidates_as()) { + int nDsInJet = 0; - TVector3 dsVector(ds.px(), ds.py(), ds.pz()); + // Loop over Ds candidates (particle level) + for (const auto& dsCandidate : jet.candidates_as()) { + ++nDsInJet; + registry.fill(HIST("h_dsjet_counter_data"), 1.5); // Ds candidates associated with the jet - const float deltaR = jetutilities::deltaR(jet, ds); + TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz()); + // Axis distance Delta_R + const float deltaR = jetutilities::deltaR(jet, dsCandidate); + // zParallel defined as longitudinal momentum fraction along the jet axis 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 = computeJetMass(jetTracks); - - TVector3 jetVector(jet.px(), jet.py(), jet.pz()); - - for (const auto& ds : jet.template candidates_as()) { + // --- Ds-level observables --- + registry.fill(HIST("h_ds_mass_data"), dsCandidate.m()); + registry.fill(HIST("h_ds_pt_data"), dsCandidate.pt()); + registry.fill(HIST("h_ds_eta_data"), dsCandidate.eta()); + registry.fill(HIST("h_ds_phi_data"), dsCandidate.phi()); - TVector3 dsVector(ds.px(), ds.py(), ds.pz()); + registry.fill(HIST("h_ds_jet_projection_data"), zParallel); - const float deltaR = jetutilities::deltaR(jet, ds); + registry.fill(HIST("h_ds_jet_distance_data"), deltaR); - const float zParallel = (jetVector * dsVector) / (jetVector * jetVector); - - registry.fill(HIST("hSparse_ds"), - ds.pt(), + // Main THnSparse: invariant mass, pT, z, and DeltaR + registry.fill(HIST("hSparse_ds_data"), + dsCandidate.m(), + dsCandidate.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); + // Jet-level quantities (filled once per jet containing at least one Ds) + if (nDsInJet > 0) { + + registry.fill(HIST("h_dsjet_counter_data"), 2.5); // Ds jets with at least one associated candidate + + // Jet properties + registry.fill(HIST("h_ds_jet_pt_data"), jet.pt()); + registry.fill(HIST("h_ds_jet_eta_data"), jet.eta()); + registry.fill(HIST("h_ds_jet_phi_data"), jet.phi()); + // Jet Mass + registry.fill(HIST("h_ds_jet_mass_data"), mjet); + + // Jet substructure observables + if (lambda11 >= 0.f) { + registry.fill(HIST("h_ds_jet_lambda11_data"), lambda11); + } + if (lambda12 >= 0.f) { + registry.fill(HIST("h_ds_jet_lambda12_data"), lambda12); + } } } } + PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "Data charged jets", false); //============== - // DATA process + // MC function //============== - - void processDataChargedSubstructure( - aod::JetCollision const& collision, - soa::Filtered const& jets, - DsCandidatesData const& candidates, - aod::JetTracks const&) + template + void analyseMonteCarloEfficiency(MCDJetsPerMCCollissionPreslice const& jetmcdpreslice, + aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + DsMCDJets const& mcdjets, + CandidatesMCD const& /*mcdCandidates*/, + aod::JetTracks const& tracks) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; - } - - analyseData(jets, candidates); - } - - PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "data charged jets", false); - - //============== - // MCD process - //============== + for (const auto& mccollision : mccollisions) { + // Count all generated MC collisions + registry.fill(HIST("McEffCol"), getValFromBin(BinMCColCntr::All)); - void processMCDChargedSubstructure( - aod::JetCollision const& collision, - soa::Filtered const& jets, - aod::JetTracks const&) - { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; + // Apply MC vertex selection + if (std::abs(mccollision.posZ()) > vertexZCut) { + continue; + } + // MC collisions passing z_cut selection + registry.fill(HIST("McEffCol"), getValFromBin(BinMCColCntr::ZCut)); + + // Reconstructed collisions associated to this mccollision + const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mccollision.globalIndex()); + for (const auto& collision : collisionsPerMCCollision) { + + // Successfully matched reconstructed collision + registry.fill(HIST("McEffCol"), getValFromBin(BinMCColCntr::Matched)); + + // Apply standard event selection and vertex cut + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || + !(std::abs(collision.posZ()) < vertexZCut)) { + continue; + } + // Matched collision passing analysis selections + registry.fill(HIST("McEffCol"), getValFromBin(BinMCColCntr::MatchedSel8ZCut)); + + // Detector-level Ds-tagged jets associated with the current reconstructed collision + const auto dsmcdJetsPerCollision = mcdjets.sliceBy(jetmcdpreslice, collision.globalIndex()); + for (const auto& mcdjet : dsmcdJetsPerCollision) { + + // Detector-level jet found in a matched collision + registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::DetectorLevelJetInMCCollision)); + + // Leading Ds candidate associated to the jet + auto mcdDscand = mcdjet.template candidates_first_as(); + + // Check whether a matched particle-level jet exists + if (mcdjet.has_matchedJetCand()) { + registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::DetectorLevelJetWithMatchedCandidate)); + } + + // Compute jet-substructure observables + TVector3 mcd_jetvector(mcdjet.px(), mcdjet.py(), mcdjet.pz()); + TVector3 mcd_candvector(mcdDscand.px(), mcdDscand.py(), mcdDscand.pz()); + + float mcd_zParallel = (mcd_jetvector * mcd_candvector) / (mcd_jetvector * mcd_jetvector); + // Axis distance Delta_R + float mcd_deltaR = jetutilities::deltaR(mcdjet, mcdDscand); + float mcd_lambda11 = computeLambda(mcdjet, tracks, 1.f, 1.f); + float mcd_lambda12 = computeLambda(mcdjet, tracks, 2.f, 1.f); + + // Detector-level Jet Histograms + registry.fill(HIST("h_jet_pt_mcd"), mcdjet.pt()); + registry.fill(HIST("h_jet_eta_mcd"), mcdjet.eta()); + registry.fill(HIST("h_jet_phi_mcd"), mcdjet.phi()); + registry.fill(HIST("h_ds_jet_projection_mcd"), mcd_zParallel); + registry.fill(HIST("h_ds_jet_lambda11_mcd"), mcd_lambda11); + registry.fill(HIST("h_ds_jet_lambda12_mcd"), mcd_lambda12); + // Detector-level Ds Histgrams + registry.fill(HIST("h_ds_jet_pt_mcd"), mcdDscand.pt()); + registry.fill(HIST("h_ds_jet_mass_mcd"), mcdDscand.m()); + registry.fill(HIST("h_ds_jet_eta_mcd"), mcdDscand.eta()); + registry.fill(HIST("h_ds_jet_phi_mcd"), mcdDscand.phi()); + + // Main THnSparse: invariant mass, pT, z, and DeltaR + registry.fill(HIST("hSparse_ds_mcd"), + mcdDscand.m(), + mcdDscand.pt(), + mcdjet.pt(), + mcd_zParallel, + mcd_deltaR); + } + } } - - analyseMCD(jets); } - - PROCESS_SWITCH(JetDsSpecSubs, processMCDChargedSubstructure, "MC detector level", false); - //============== - // MCP process + // MC process //============== - void processMCPChargedSubstructure( - aod::JetCollision const& collision, - soa::Filtered const& jets, - aod::JetTracks const&) + void processMonteCarloEfficiencyDs(aod::JetMcCollisions const& mccollisions, + aod::JetCollisionsMCD const& collisions, + DsMCDJets const& mcdjets, + aod::CandidatesDsMCD const& mcdCandidates, + aod::JetTracks const& jettracks) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; - } - analyseMCP(jets); + analyseMonteCarloEfficiency, + DsMCDJets, + aod::CandidatesDsMCD>(dsMCDJetsPerEXPCollisionPreslice, + mccollisions, + collisions, + mcdjets, + mcdCandidates, + jettracks); } - - PROCESS_SWITCH(JetDsSpecSubs, processMCPChargedSubstructure, "MC particle level", false); + PROCESS_SWITCH(JetDsSpecSubs, processMonteCarloEfficiencyDs, "Non-matched and matched MC Ds and jets", false); }; - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)};