Skip to content

Commit dc440ec

Browse files
monamelopmonamelop
authored
[PWGJE] Add MC efficiency histograms for Ds-jet analysis (#16957)
1 parent 577ba1e commit dc440ec

1 file changed

Lines changed: 119 additions & 38 deletions

File tree

‎PWGJE/Tasks/jetDsSpecSubs.cxx‎

Lines changed: 119 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@
1818
#include "PWGJE/DataModel/Jet.h"
1919
#include "PWGJE/DataModel/JetReducedData.h"
2020

21+
#include "Common/Core/RecoDecay.h"
22+
2123
#include <Framework/ASoA.h>
2224
#include <Framework/AnalysisDataModel.h>
2325
#include <Framework/AnalysisTask.h>
@@ -29,6 +31,7 @@
2931
#include <Framework/runDataProcessing.h>
3032

3133
#include <TH1.h>
34+
#include <THnSparse.h>
3235
#include <TVector3.h>
3336

3437
#include <cmath>
@@ -72,12 +75,12 @@ struct JetDsSpecSubs {
7275

7376
using DsDataJets = soa::Join<aod::DsChargedJets, aod::DsChargedJetConstituents>;
7477
using DsMCDJets = soa::Join<aod::DsChargedMCDetectorLevelJets, aod::DsChargedMCDetectorLevelJetConstituents, aod::DsChargedMCDetectorLevelJetsMatchedToDsChargedMCParticleLevelJets>;
75-
// using DsMCPJets = soa::Join<aod::DsChargedMCParticleLevelJets, aod::DsChargedMCParticleLevelJetConstituents, aod::DsChargedMCDetectorLevelJetsMatchedToDsChargedMCParticleLevelJets>;
78+
using DsMCPJets = soa::Join<aod::DsChargedMCParticleLevelJets, aod::DsChargedMCParticleLevelJetConstituents, aod::DsChargedMCDetectorLevelJetsMatchedToDsChargedMCParticleLevelJets>;
7679

7780
// Slices for access to proper HF MCD jet collision that is associated to MCCollision
7881
PresliceUnsorted<aod::JetCollisionsMCD> collisionsPerMCCollisionPreslice = aod::jmccollisionlb::mcCollisionId;
7982
Preslice<DsMCDJets> dsMCDJetsPerEXPCollisionPreslice = aod::jet::collisionId;
80-
// Preslice<DsMCPJets> dsMCPJetsPerMCCollisionPreslice = aod::jet::mcCollisionId;
83+
Preslice<DsMCPJets> dsMCPJetsPerMCCollisionPreslice = aod::jet::mcCollisionId;
8184

8285
// Configurables
8386
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
@@ -91,6 +94,7 @@ struct JetDsSpecSubs {
9194
std::vector<int> eventSelectionBits;
9295
int trackSelection = -1;
9396

97+
static constexpr float epsilon = 1.e-6f;
9498
// Filters
9599
// Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
96100
// Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut;
@@ -104,8 +108,6 @@ struct JetDsSpecSubs {
104108
{
105109
{"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}},
106110
{"h_collision_counter_data", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}},
107-
{"h_collision_counter_mcd", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}},
108-
{"h_collision_counter_mcp", ";event counter;entries", {HistType::kTH1F, {{10, 0., 10.}}}},
109111

110112
{"h_track_pt", ";#it{p}_{T,track};entries", {HistType::kTH1F, {{200, 0., 200.}}}},
111113
{"h_track_eta", ";#eta_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}},
@@ -132,11 +134,13 @@ struct JetDsSpecSubs {
132134
{"h_ds_jet_mass_data", ";m_{jet}^{ch} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}},
133135
{"h_ds_jet_lambda11_data", ";#lambda_{1}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
134136
{"h_ds_jet_lambda12_data", ";#lambda_{2}^{1};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
135-
{"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}}}},
137+
{"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}}}},
136138

137-
// MC detector-level histograms
138-
{"McEffCol", "N_{collisions};", {HistType::kTH1F, {{4, 0., 4.0}}}},
139+
// MC general histograms
139140
{"McEffJet", "N_{jet};", {HistType::kTH1F, {{4, 0., 4.0}}}},
141+
{"McEffCol", "N_{collisions};", {HistType::kTH1F, {{4, 0., 4.0}}}},
142+
143+
// MC detector-level histograms
140144
{"h_jet_pt_mcd", "detector-level jet pT;#it{p}_{T,jet}^{det} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
141145
{"h_jet_eta_mcd", "detector-level jet #eta;#eta_{jet}^{det};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
142146
{"h_jet_phi_mcd", "detector-level jet #phi;#phi_{jet}^{det};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
@@ -154,24 +158,23 @@ struct JetDsSpecSubs {
154158
{"h_ds_jet_mass_mcd", ";m_{jet}^{ch, det} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{200, 0., 50.}}}},
155159
{"h_ds_jet_lambda11_mcd", ";#lambda_{1}^{1, det};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
156160
{"h_ds_jet_lambda12_mcd", ";#lambda_{2}^{1, det};entries", {HistType::kTH1F, {{200, 0., 1.0}}}},
157-
{"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}}}},
161+
{"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}}}},
158162

159163
// MC particle-level histograms
160-
// {"h_dsjet_counter_mcp", ";type;counts", {HistType::kTH1F, {{3, 0., 3.}}}},
161-
// {"h_jet_pt_mcp", "particle-level jet pT;#it{p}_{T,jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
162-
// {"h_jet_eta_mcp", "particle-level jet #eta;#eta_{jet}^{part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
163-
// {"h_jet_phi_mcp", "particle-level jet #phi;#phi_{jet}^{part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
164-
165-
// {"h_ds_pt_mcp", ";#it{p}_{T,D_{S}}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}},
166-
// {"h_ds_eta_mcp", ";#eta_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
167-
// {"h_ds_phi_mcp", ";#phi_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
168-
169-
// {"h_ds_jet_pt_mcp", ";#it{p}_{T,D_{S} jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}},
170-
// {"h_ds_jet_eta_mcp", ";#eta_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
171-
// {"h_ds_jet_phi_mcp", ";#phi_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
172-
// {"h_ds_jet_projection_mcp", ";z^{D_{S},jet}_{||, part};entries", {HistType::kTH1F, {{1000, 0., 2.}}}},
173-
// {"h_ds_jet_distance_mcp", ";#DeltaR_{D_{S},jet}^{part};entries", {HistType::kTH1F, {{1000, 0., 1.}}}},
174-
// {"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}}}},
164+
{"h_jet_pt_mcp", "particle-level jet pT;#it{p}_{T,jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
165+
{"h_jet_eta_mcp", "particle-level jet #eta;#eta_{jet}^{part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
166+
{"h_jet_phi_mcp", "particle-level jet #phi;#phi_{jet}^{part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
167+
168+
{"h_ds_pt_mcp", ";#it{p}_{T,D_{S}}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}},
169+
{"h_ds_eta_mcp", ";#eta_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
170+
{"h_ds_phi_mcp", ";#phi_{D_{S}}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
171+
172+
{"h_ds_jet_pt_mcp", ";#it{p}_{T,D_{S} jet}^{part} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 100.}}}},
173+
{"h_ds_jet_eta_mcp", ";#eta_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 1.}}}},
174+
{"h_ds_jet_phi_mcp", ";#phi_{D_{S} jet}^{part};entries", {HistType::kTH1F, {{250, -1., 7.}}}},
175+
{"h_ds_jet_projection_mcp", ";z^{D_{S},jet}_{||, part};entries", {HistType::kTH1F, {{1000, 0., 2.}}}},
176+
{"h_ds_jet_distance_mcp", ";#DeltaR_{D_{S},jet}^{part};entries", {HistType::kTH1F, {{1000, 0., 1.}}}},
177+
{"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}}}},
175178
}};
176179
//========
177180
// INIT
@@ -199,6 +202,11 @@ struct JetDsSpecSubs {
199202
jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetInMCCollision, "detector level");
200203
jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::DetectorLevelJetWithMatchedCandidate, "particle matched jets");
201204
jetCounter->GetXaxis()->SetBinLabel(BinMCJetCntr::ParticleLevelJetWithMatchedCandidate, "detector matched jets");
205+
206+
auto hSparse_ds_mcd = registry.get<THnSparse>(HIST("hSparse_ds_mcd"));
207+
auto* axisOrigin = hSparse_ds_mcd->GetAxis(5);
208+
axisOrigin->SetBinLabel(1, "Prompt");
209+
axisOrigin->SetBinLabel(2, "Non-prompt");
202210
}
203211
//===============
204212
// Lambda compute
@@ -216,9 +224,9 @@ struct JetDsSpecSubs {
216224
const float dr = jetutilities::deltaR(jet, trk);
217225
sum += std::pow(trk.pt(), kappa) * std::pow(dr, alpha);
218226
}
219-
const float jetR = jet.r() / 100.f;
227+
const float jetRadius = jet.r() / 100.f;
220228

221-
const float denom = std::pow(jet.pt(), kappa) * std::pow(jetR, alpha);
229+
const float denom = std::pow(jet.pt(), kappa) * std::pow(jetRadius, alpha);
222230
if (denom <= 0.f) {
223231
return -1.f;
224232
}
@@ -298,6 +306,10 @@ struct JetDsSpecSubs {
298306
registry.fill(HIST("h_collision_counter_data"), 3.0);
299307

300308
for (const auto& jet : jets) {
309+
const float jetRadius = jet.r() / 100.f;
310+
if (std::abs(jetRadius - jetR.value) > epsilon) {
311+
continue;
312+
}
301313
registry.fill(HIST("h_dsjet_counter_data"), 0.5); // DsChargedJets entries
302314

303315
registry.fill(HIST("h_jet_pt_data"), jet.pt());
@@ -374,13 +386,19 @@ struct JetDsSpecSubs {
374386
// MC function
375387
//==============
376388
template <typename MCDJetsPerMCCollissionPreslice,
389+
typename MCPJetsPerMCCollissionPreslice,
377390
typename DsMCDJets,
378-
typename CandidatesMCD>
391+
typename DsMCPJets,
392+
typename DsCandidatesMCD,
393+
typename DsCandidatesMCP>
379394
void analyseMonteCarloEfficiency(MCDJetsPerMCCollissionPreslice const& jetmcdpreslice,
395+
MCPJetsPerMCCollissionPreslice const& jetmcppreslice,
380396
aod::JetMcCollisions const& mccollisions,
381397
aod::JetCollisionsMCD const& collisions,
382398
DsMCDJets const& mcdjets,
383-
CandidatesMCD const& /*mcdCandidates*/,
399+
DsMCPJets const& mcpjets,
400+
DsCandidatesMCD const& /*mcdCandidates*/,
401+
DsCandidatesMCP const& /*mcpCandidates*/,
384402
aod::JetTracks const& tracks)
385403
{
386404
for (const auto& mccollision : mccollisions) {
@@ -412,12 +430,22 @@ struct JetDsSpecSubs {
412430
// Detector-level Ds-tagged jets associated with the current reconstructed collision
413431
const auto dsmcdJetsPerCollision = mcdjets.sliceBy(jetmcdpreslice, collision.globalIndex());
414432
for (const auto& mcdjet : dsmcdJetsPerCollision) {
415-
433+
const float mcdjetRadius = mcdjet.r() / 100.f;
434+
if (std::abs(mcdjetRadius - jetR.value) > epsilon) {
435+
continue;
436+
}
416437
// Detector-level jet found in a matched collision
417438
registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::DetectorLevelJetInMCCollision));
418439

419440
// Leading Ds candidate associated to the jet
420-
auto mcdDscand = mcdjet.template candidates_first_as<CandidatesMCD>();
441+
auto mcdDscand = mcdjet.template candidates_first_as<DsCandidatesMCD>();
442+
443+
// Check if it's prompt
444+
int origin = 0;
445+
446+
if (mcdDscand.originMcRec() != RecoDecay::OriginType::Prompt) {
447+
origin = 1;
448+
}
421449

422450
// Check whether a matched particle-level jet exists
423451
if (mcdjet.has_matchedJetCand()) {
@@ -453,9 +481,55 @@ struct JetDsSpecSubs {
453481
mcdDscand.pt(),
454482
mcdjet.pt(),
455483
mcd_zParallel,
456-
mcd_deltaR);
484+
mcd_deltaR,
485+
origin);
457486
}
458487
}
488+
// Particle level
489+
const auto dsmcpJetsPerMCCollision = mcpjets.sliceBy(jetmcppreslice, mccollision.globalIndex());
490+
for (const auto& mcpjet : dsmcpJetsPerMCCollision) {
491+
const float mcpjetRadius = mcpjet.r() / 100.f;
492+
if (std::abs(mcpjetRadius - jetR.value) > epsilon) {
493+
continue;
494+
}
495+
496+
registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::ParticleLevelJetInMCCollision));
497+
498+
// obtain leading HF particle in jet
499+
auto mcpDscand = mcpjet.template candidates_first_as<DsCandidatesMCP>();
500+
501+
if (mcpjet.has_matchedJetCand()) {
502+
registry.fill(HIST("McEffJet"), getValFromBin(BinMCJetCntr::ParticleLevelJetWithMatchedCandidate));
503+
}
504+
505+
TVector3 mcp_jetvector(mcpjet.px(), mcpjet.py(), mcpjet.pz());
506+
TVector3 mcp_candvector(mcpDscand.px(), mcpDscand.py(), mcpDscand.pz());
507+
508+
float mcp_zParallel = (mcp_jetvector * mcp_candvector) / (mcp_jetvector * mcp_jetvector);
509+
// Axis distance Delta_R
510+
float mcp_deltaR = jetutilities::deltaR(mcpjet, mcpDscand);
511+
float mcp_lambda11 = computeLambda(mcpjet, tracks, 1.f, 1.f);
512+
float mcp_lambda12 = computeLambda(mcpjet, tracks, 2.f, 1.f);
513+
514+
// Particle-level Jet Histograms
515+
registry.fill(HIST("h_jet_pt_mcp"), mcpjet.pt());
516+
registry.fill(HIST("h_jet_eta_mcp"), mcpjet.eta());
517+
registry.fill(HIST("h_jet_phi_mcp"), mcpjet.phi());
518+
registry.fill(HIST("h_ds_jet_projection_mcp"), mcp_zParallel);
519+
registry.fill(HIST("h_ds_jet_lambda11_mcp"), mcp_lambda11);
520+
registry.fill(HIST("h_ds_jet_lambda12_mcp"), mcp_lambda12);
521+
// Particle-level Ds Histgrams
522+
registry.fill(HIST("h_ds_jet_pt_mcp"), mcpDscand.pt());
523+
registry.fill(HIST("h_ds_jet_eta_mcp"), mcpDscand.eta());
524+
registry.fill(HIST("h_ds_jet_phi_mcp"), mcpDscand.phi());
525+
526+
// Main THnSparse: invariant mass, pT, z, and DeltaR
527+
registry.fill(HIST("hSparse_ds_mcp"),
528+
mcpDscand.pt(),
529+
mcpjet.pt(),
530+
mcp_zParallel,
531+
mcp_deltaR);
532+
}
459533
}
460534
}
461535
//==============
@@ -465,18 +539,25 @@ struct JetDsSpecSubs {
465539
void processMonteCarloEfficiencyDs(aod::JetMcCollisions const& mccollisions,
466540
aod::JetCollisionsMCD const& collisions,
467541
DsMCDJets const& mcdjets,
468-
aod::CandidatesDsMCD const& mcdCandidates,
542+
DsMCPJets const& mcpjets,
543+
DsCandidatesMCD const& mcdDscand,
544+
DsCandidatesMCP const& mcpDscand,
469545
aod::JetTracks const& jettracks)
470546
{
471-
472547
analyseMonteCarloEfficiency<Preslice<DsMCDJets>,
548+
Preslice<DsMCPJets>,
473549
DsMCDJets,
474-
aod::CandidatesDsMCD>(dsMCDJetsPerEXPCollisionPreslice,
475-
mccollisions,
476-
collisions,
477-
mcdjets,
478-
mcdCandidates,
479-
jettracks);
550+
DsMCPJets,
551+
DsCandidatesMCD,
552+
DsCandidatesMCP>(dsMCDJetsPerEXPCollisionPreslice,
553+
dsMCPJetsPerMCCollisionPreslice,
554+
mccollisions,
555+
collisions,
556+
mcdjets,
557+
mcpjets,
558+
mcdDscand,
559+
mcpDscand,
560+
jettracks);
480561
}
481562
PROCESS_SWITCH(JetDsSpecSubs, processMonteCarloEfficiencyDs, "Non-matched and matched MC Ds and jets", false);
482563
};

0 commit comments

Comments
 (0)