Skip to content

[PWGDQ] Adding MC mother information #14485

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ppoava wants to merge 1 commit into
base: master
Choose a base branch
from
+104 −1
Open

[PWGDQ] Adding MC mother information #14485

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
105 changes: 104 additions & 1 deletion PWGDQ/Tasks/qaMatching.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
#include <Math/ProbFunc.h>

#include <algorithm>
#include <iostream>

Check failure on line 33 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[include-iostream] 

Do not include iostream. Use O2 logging instead.
#include <limits>
#include <map>
#include <memory>
Expand Down Expand Up @@ -151,6 +151,10 @@
Configurable<std::string> fConfigGrpMagPath{"grpmagPath-", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
} fConfigCCDB;

struct : ConfigurableGroup {
Configurable<bool> fCreatePdgMomHistograms{"cfgCreatePdgMomHistograms", false, "create matching characteristics plots with particle mom PDG codes"};
} fConfigQAs;

/// Variables for histograms configuration
Configurable<int> fNCandidatesMax{"nCandidatesMax", 5, ""};

Expand Down Expand Up @@ -297,12 +301,20 @@
struct EfficiencyPlotter {
o2::framework::HistPtr p_num;
o2::framework::HistPtr p_den;
o2::framework::HistPtr p_pdg_num;
o2::framework::HistPtr p_pdg_den;
o2::framework::HistPtr pt_num;
o2::framework::HistPtr pt_den;
o2::framework::HistPtr pt_pdg_num;
o2::framework::HistPtr pt_pdg_den;
o2::framework::HistPtr phi_num;
o2::framework::HistPtr phi_den;
o2::framework::HistPtr phi_pdg_num;
o2::framework::HistPtr phi_pdg_den;
o2::framework::HistPtr eta_num;
o2::framework::HistPtr eta_den;
o2::framework::HistPtr eta_pdg_num;
o2::framework::HistPtr eta_pdg_den;

EfficiencyPlotter(std::string path, std::string title,
HistogramRegistry& registry)
Expand All @@ -311,6 +323,7 @@
AxisSpec pTAxis = {100, 0, 10, "p_{T} (GeV/c)"};
AxisSpec etaAxis = {100, -4, -2, "#eta"};
AxisSpec phiAxis = {90, -180, 180, "#phi (degrees)"};
AxisSpec motherPDGAxis{1201, -600.5, 600.5, "Direct mother PDG"};

std::string histName;
std::string histTitle;
Expand All @@ -324,6 +337,14 @@
histTitle = title + " vs. p - den";
p_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pAxis}});

histName = path + "p_pdg_num";
histTitle = title + " vs. p vs pdg ID - num";
p_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}});

histName = path + "p_pdg_den";
histTitle = title + " vs. p vs pdg ID - den";
p_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}});

// pT dependence
histName = path + "pt_num";
histTitle = title + " vs. p_{T} - num";
Expand All @@ -333,6 +354,14 @@
histTitle = title + " vs. p_{T} - den";
pt_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pTAxis}});

histName = path + "pt_pdg_num";
histTitle = title + " vs. p_{T} vs pdg ID - num";
pt_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}});

histName = path + "pt_pdg_den";
histTitle = title + " vs. p_{T} vs pdg ID - den";
pt_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}});

// eta dependence
histName = path + "eta_num";
histTitle = title + " vs. #eta - num";
Expand All @@ -342,6 +371,14 @@
histTitle = title + " vs. #eta - den";
eta_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {etaAxis}});

histName = path + "eta_pdg_num";
histTitle = title + " vs. #eta vs pdg ID - num";
eta_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}});

histName = path + "eta_pdg_den";
histTitle = title + " vs. #eta vs pdg ID - den";
eta_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}});

// phi dependence
histName = path + "phi_num";
histTitle = title + " vs. #phi - num";
Expand All @@ -350,6 +387,14 @@
histName = path + "phi_den";
histTitle = title + " vs. #phi - den";
phi_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {phiAxis}});

histName = path + "phi_pdg_num";
histTitle = title + " vs. #phi vs pdg ID - num";
phi_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}});

histName = path + "phi_pdg_den";
histTitle = title + " vs. #phi vs pdg ID - den";
phi_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}});
}

template <class T>
Expand All @@ -368,6 +413,24 @@
std::get<std::shared_ptr<TH1>>(phi_num)->Fill(phi);
}
}

// Study the PDG origin of particles and their effect on the purity score
template <class T>
void Fill(const T& track, int pdgCode, bool passed)
{
double phi = track.phi() * 180 / TMath::Pi();
std::get<std::shared_ptr<TH2>>(p_pdg_den)->Fill(track.p(), pdgCode);
std::get<std::shared_ptr<TH2>>(pt_pdg_den)->Fill(track.pt(), pdgCode);
std::get<std::shared_ptr<TH2>>(eta_pdg_den)->Fill(track.eta(), pdgCode);
std::get<std::shared_ptr<TH2>>(phi_pdg_den)->Fill(phi, pdgCode);

if (passed) {
std::get<std::shared_ptr<TH2>>(p_pdg_num)->Fill(track.p(), pdgCode);
std::get<std::shared_ptr<TH2>>(pt_pdg_num)->Fill(track.pt(), pdgCode);
std::get<std::shared_ptr<TH2>>(eta_pdg_num)->Fill(track.eta(), pdgCode);
std::get<std::shared_ptr<TH2>>(phi_pdg_num)->Fill(phi, pdgCode);
}
}
};

struct MatchRankingHistos {
Expand Down Expand Up @@ -1221,8 +1284,8 @@
// extrapolation with MCH tools
auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex);
double pMCH = mchTrackAtVertex.getP();
double px = pMCH * sin(M_PI / 2 - atan(mftTrack.tgl())) * cos(mftTrack.phi());

Check failure on line 1287 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double py = pMCH * sin(M_PI / 2 - atan(mftTrack.tgl())) * sin(mftTrack.phi());

Check failure on line 1288 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
double sign = mchTrack.sign();

Expand Down Expand Up @@ -1413,7 +1476,7 @@
continue;
// get the index associated to the MC particle
auto muonMcParticle = muonTrack.mcParticle();
if (std::abs(muonMcParticle.pdgCode()) != 13)

Check failure on line 1479 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
continue;

int64_t muonMcTrackIndex = muonMcParticle.globalIndex();
Expand Down Expand Up @@ -1470,7 +1533,7 @@
if (mchMcParticle.globalIndex() != mftMcParticle.globalIndex())
return false;

if (std::abs(mchMcParticle.pdgCode()) != 13)

Check failure on line 1536 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
return false;

return true;
Expand Down Expand Up @@ -1520,7 +1583,7 @@
}

if (result == kMatchTypeUndefined) {
std::cout << std::format("[GetMatchType] isPaired={} isMuon={} decayRanking={} result={}",

Check failure on line 1586 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
isPaired, isMuon, decayRanking, static_cast<int>(result))
<< std::endl;
}
Expand Down Expand Up @@ -2040,7 +2103,7 @@
// ====================================
// Matching purity
if (verbose)
std::cout << std::format(" Filling matching purity plots with score cut {}", matchingScoreCut) << std::endl;

Check failure on line 2106 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
for (auto [mchIndex, globalTracksVector] : matchingCandidates) {
if (globalTracksVector.size() < 1)
continue;
Expand All @@ -2066,10 +2129,26 @@
// check if the matching candidate is a true one
bool isTrueMatch = IsTrueGlobalMatching(muonTrack, matchablePairs);

if (verbose)
// ---- MC ancestry ----
auto motherParticles = GetMotherParticles(muonTrack);
int motherPDG = 0;
if (motherParticles.size() > 1) {
motherPDG = motherParticles[1].first;
}

if (verbose) {
std::cout << std::format(" MCH track #{} -> Muon track #{}, isTrueMatch={}", mchIndex, globalTracksVector[0].globalTrackId, isTrueMatch) << std::endl;

Check failure on line 2140 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
std::cout << " MC ancestry (pdg): ";

Check failure on line 2141 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
for (auto const& [pdg, idx] : motherParticles) {
std::cout << "(" << pdg << ") ";

Check failure on line 2143 in PWGDQ/Tasks/qaMatching.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
}
std::cout << std::endl;
}
// fill matching purity plots
plotter->fMatchingPurityPlotter.Fill(mchTrack, isTrueMatch);
if (fConfigQAs.fCreatePdgMomHistograms) {
plotter->fMatchingPurityPlotter.Fill(mchTrack, motherPDG, isTrueMatch);
}
}

// ====================================
Expand Down Expand Up @@ -2105,10 +2184,34 @@
}
}

// ---- MC ancestry ----
auto motherParticles = GetMotherParticles(mchTrack);
int motherPDG = 0;
if (motherParticles.size() > 1) {
motherPDG = motherParticles[1].first;
}

if (verbose) {
std::cout << " MC ancestry (pdg): ";
for (auto const& [pdg, idx] : motherParticles) {
std::cout << "(" << pdg << ") ";
}
std::cout << std::endl;
}

// fill matching efficiency plots
plotter->fPairingEfficiencyPlotter.Fill(mchTrack, goodMatchFound);
if (fConfigQAs.fCreatePdgMomHistograms) {
plotter->fPairingEfficiencyPlotter.Fill(mchTrack, motherPDG, goodMatchFound);
}
plotter->fMatchingEfficiencyPlotter.Fill(mchTrack, (goodMatchFound && isTrueMatch));
if (fConfigQAs.fCreatePdgMomHistograms) {
plotter->fMatchingEfficiencyPlotter.Fill(mchTrack, motherPDG, (goodMatchFound && isTrueMatch));
}
plotter->fFakeMatchingEfficiencyPlotter.Fill(mchTrack, (goodMatchFound && !isTrueMatch));
if (fConfigQAs.fCreatePdgMomHistograms) {
plotter->fFakeMatchingEfficiencyPlotter.Fill(mchTrack, motherPDG, (goodMatchFound && !isTrueMatch));
}
}
}

Expand Down
Loading