diff --git a/PWGJE/Tasks/jetDsSpecSubs.cxx b/PWGJE/Tasks/jetDsSpecSubs.cxx index b4f33806754..cff0187f5d4 100644 --- a/PWGJE/Tasks/jetDsSpecSubs.cxx +++ b/PWGJE/Tasks/jetDsSpecSubs.cxx @@ -18,6 +18,8 @@ #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" +#include "Common/Core/RecoDecay.h" + #include #include #include @@ -72,12 +74,12 @@ struct JetDsSpecSubs { using DsDataJets = soa::Join; using DsMCDJets = soa::Join; - // using DsMCPJets = 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; + Preslice dsMCPJetsPerMCCollisionPreslice = aod::jet::mcCollisionId; // Configurables Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; @@ -91,6 +93,7 @@ struct JetDsSpecSubs { std::vector eventSelectionBits; int trackSelection = -1; + static constexpr float epsilon = 1.e-6f; // Filters // Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); // Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; @@ -104,8 +107,6 @@ struct JetDsSpecSubs { { {"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}}, {"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.}}}}, @@ -132,11 +133,13 @@ struct JetDsSpecSubs { {"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}}}}, + {"hSparse_ds_data", ";m_{D_{S}};#it{p}_{T,D_{S}};#it{p}_{T,jet};z^{D_{S},jet}_{||};#DeltaR_{D_{S},jet};Origin(D_{S})", {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}}}}, + // MC general histograms {"McEffJet", "N_{jet};", {HistType::kTH1F, {{4, 0., 4.0}}}}, + {"McEffCol", "N_{collisions};", {HistType::kTH1F, {{4, 0., 4.0}}}}, + + // MC detector-level histograms {"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.}}}}, @@ -154,24 +157,23 @@ struct JetDsSpecSubs { {"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}}}}, + {"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}, {2, -0.5, 1.5}}}}, // 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}}}}, + {"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 @@ -199,6 +201,11 @@ struct JetDsSpecSubs { jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetInMCCollision, "detector level"); jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetWithMatchedCandidate, "particle matched jets"); jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::ParticleLevelJetWithMatchedCandidate, "detector matched jets"); + + auto hSparse_ds_mcd = registry.get(HIST("hSparse_ds_mcd")); + auto* axisOrigin = hSparse_ds_mcd->GetAxis(5); + axisOrigin->SetBinLabel(1, "Prompt"); + axisOrigin->SetBinLabel(2, "Non-prompt"); } //=============== // Lambda compute @@ -216,9 +223,9 @@ struct JetDsSpecSubs { const float dr = jetutilities::deltaR(jet, trk); sum += std::pow(trk.pt(), kappa) * std::pow(dr, alpha); } - const float jetR = jet.r() / 100.f; + const float jetRadius = jet.r() / 100.f; - const float denom = std::pow(jet.pt(), kappa) * std::pow(jetR, alpha); + const float denom = std::pow(jet.pt(), kappa) * std::pow(jetRadius, alpha); if (denom <= 0.f) { return -1.f; } @@ -298,6 +305,10 @@ struct JetDsSpecSubs { registry.fill(HIST("h_collision_counter_data"), 3.0); for (const auto& jet : jets) { + const float jetRadius = jet.r() / 100.f; + if (std::abs(jetRadius - jetR.value) > epsilon) { + continue; + } registry.fill(HIST("h_dsjet_counter_data"), 0.5); // DsChargedJets entries registry.fill(HIST("h_jet_pt_data"), jet.pt()); @@ -374,13 +385,19 @@ struct JetDsSpecSubs { // MC function //============== template + typename DsMCPJets, + typename DsCandidatesMCD, + typename DsCandidatesMCP> void analyseMonteCarloEfficiency(MCDJetsPerMCCollissionPreslice const& jetmcdpreslice, + MCPJetsPerMCCollissionPreslice const& jetmcppreslice, aod::JetMcCollisions const& mccollisions, aod::JetCollisionsMCD const& collisions, DsMCDJets const& mcdjets, - CandidatesMCD const& /*mcdCandidates*/, + DsMCPJets const& mcpjets, + DsCandidatesMCD const& /*mcdCandidates*/, + DsCandidatesMCP const& /*mcpCandidates*/, aod::JetTracks const& tracks) { for (const auto& mccollision : mccollisions) { @@ -412,12 +429,22 @@ struct JetDsSpecSubs { // Detector-level Ds-tagged jets associated with the current reconstructed collision const auto dsmcdJetsPerCollision = mcdjets.sliceBy(jetmcdpreslice, collision.globalIndex()); for (const auto& mcdjet : dsmcdJetsPerCollision) { - + const float mcdjetRadius = mcdjet.r() / 100.f; + if (std::abs(mcdjetRadius - jetR.value) > epsilon) { + continue; + } // 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(); + auto mcdDscand = mcdjet.template candidates_first_as(); + + // Check if it's prompt + int origin = 0; + + if (mcdDscand.originMcRec() != RecoDecay::OriginType::Prompt) { + origin = 1; + } // Check whether a matched particle-level jet exists if (mcdjet.has_matchedJetCand()) { @@ -453,9 +480,55 @@ struct JetDsSpecSubs { mcdDscand.pt(), mcdjet.pt(), mcd_zParallel, - mcd_deltaR); + mcd_deltaR, + origin); } } + // Particle level + const auto dsmcpJetsPerMCCollision = mcpjets.sliceBy(jetmcppreslice, mccollision.globalIndex()); + for (const auto& mcpjet : dsmcpJetsPerMCCollision) { + const float mcpjetRadius = mcpjet.r() / 100.f; + if (std::abs(mcpjetRadius - jetR.value) > epsilon) { + continue; + } + + registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::ParticleLevelJetInMCCollision)); + + // obtain leading HF particle in jet + auto mcpDscand = mcpjet.template candidates_first_as(); + + if (mcpjet.has_matchedJetCand()) { + registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::ParticleLevelJetWithMatchedCandidate)); + } + + TVector3 mcp_jetvector(mcpjet.px(), mcpjet.py(), mcpjet.pz()); + TVector3 mcp_candvector(mcpDscand.px(), mcpDscand.py(), mcpDscand.pz()); + + float mcp_zParallel = (mcp_jetvector * mcp_candvector) / (mcp_jetvector * mcp_jetvector); + // Axis distance Delta_R + float mcp_deltaR = jetutilities::deltaR(mcpjet, mcpDscand); + float mcp_lambda11 = computeLambda(mcpjet, tracks, 1.f, 1.f); + float mcp_lambda12 = computeLambda(mcpjet, tracks, 2.f, 1.f); + + // Particle-level Jet Histograms + registry.fill(HIST("h_jet_pt_mcp"), mcpjet.pt()); + registry.fill(HIST("h_jet_eta_mcp"), mcpjet.eta()); + registry.fill(HIST("h_jet_phi_mcp"), mcpjet.phi()); + registry.fill(HIST("h_ds_jet_projection_mcp"), mcp_zParallel); + registry.fill(HIST("h_ds_jet_lambda11_mcp"), mcp_lambda11); + registry.fill(HIST("h_ds_jet_lambda12_mcp"), mcp_lambda12); + // Particle-level Ds Histgrams + registry.fill(HIST("h_ds_jet_pt_mcp"), mcpDscand.pt()); + registry.fill(HIST("h_ds_jet_eta_mcp"), mcpDscand.eta()); + registry.fill(HIST("h_ds_jet_phi_mcp"), mcpDscand.phi()); + + // Main THnSparse: invariant mass, pT, z, and DeltaR + registry.fill(HIST("hSparse_ds_mcp"), + mcpDscand.pt(), + mcpjet.pt(), + mcp_zParallel, + mcp_deltaR); + } } } //============== @@ -465,18 +538,25 @@ struct JetDsSpecSubs { void processMonteCarloEfficiencyDs(aod::JetMcCollisions const& mccollisions, aod::JetCollisionsMCD const& collisions, DsMCDJets const& mcdjets, - aod::CandidatesDsMCD const& mcdCandidates, + DsMCPJets const& mcpjets, + DsCandidatesMCD const& mcdDscand, + DsCandidatesMCP const& mcpDscand, aod::JetTracks const& jettracks) { - analyseMonteCarloEfficiency, + Preslice, DsMCDJets, - aod::CandidatesDsMCD>(dsMCDJetsPerEXPCollisionPreslice, - mccollisions, - collisions, - mcdjets, - mcdCandidates, - jettracks); + DsMCPJets, + DsCandidatesMCD, + DsCandidatesMCP>(dsMCDJetsPerEXPCollisionPreslice, + dsMCPJetsPerMCCollisionPreslice, + mccollisions, + collisions, + mcdjets, + mcpjets, + mcdDscand, + mcpDscand, + jettracks); } PROCESS_SWITCH(JetDsSpecSubs, processMonteCarloEfficiencyDs, "Non-matched and matched MC Ds and jets", false); };