Skip to content

PWGHF: Adding Bayesian PID selections in the Lc->pKpi selector task #977

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

Merged
vkucera merged 2 commits into from
Jul 19, 2022
Merged

PWGHF: Adding Bayesian PID selections in the Lc->pKpi selector task #977

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
05b3d6d
PWGHF: Adding Bayesian PID selections in the Lc->pKpi selector task
Jul 12, 2022
39edf9c
Merge branch 'master' into master
vkucera Jul 15, 2022
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
80 changes: 78 additions & 2 deletions Common/Core/TrackSelectorPID.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@ class TrackSelectorPID

/// Standard constructor with PDG code initialisation
explicit TrackSelectorPID(int pdg)
: mPdg(std::abs(pdg))
{
setPDG(pdg);
}

/// Default destructor
Expand All @@ -46,7 +46,36 @@ class TrackSelectorPID
PIDAccepted
};

void setPDG(int pdg) { mPdg = std::abs(pdg); }
void setPDG(int pdg)
Comment thread
vkucera marked this conversation as resolved.
{
mPdg = std::abs(pdg);
switch (mPdg) {
case kElectron: {
mSpecies = track::PID::Electron;
break;
}
case kMuonMinus: {
mSpecies = track::PID::Muon;
break;
}
case kPiPlus: {
mSpecies = track::PID::Pion;
break;
}
case kKPlus: {
mSpecies = track::PID::Kaon;
break;
}
case kProton: {
mSpecies = track::PID::Proton;
break;
}
default: {
LOGF(error, "ERROR: Species not implemented for PDG %d", mPdg);
assert(false);
}
}
}

// TPC

Expand Down Expand Up @@ -477,8 +506,51 @@ class TrackSelectorPID
return isSelRICH || isSelTOF;
}

/// Set pT range where Bayes PID is applicable.
void setRangePtBayes(float ptMin, float ptMax)
{
mPtBayesMin = ptMin;
mPtBayesMax = ptMax;
}

/// Checks if track is OK for Bayesian PID.
/// \param track track
/// \return true if track is OK for Bayesian PID
template <typename T>
bool isValidTrackBayesPID(const T& track)
{
auto pt = track.pt();
return (mPtBayesMin <= pt && pt <= mPtBayesMax);
}

/// Bayesian maximum probability algorithm.
template <typename T>
bool isSelectedTrackBayesPID(const T& track)
{
// Get index of the most probable species for a given track.
return track.bayesID() == mSpecies;
}

/// Returns status of Bayesian PID selection for a given track.
/// \param track track
/// \return Bayesian selection status (see TrackSelectorPID::Status)
template <typename T>
int getStatusTrackBayesPID(const T& track)
{
if (isValidTrackBayesPID(track)) {
if (isSelectedTrackBayesPID(track)) {
return Status::PIDAccepted; // accepted
} else {
return Status::PIDRejected; // rejected
}
} else {
return Status::PIDNotApplicable; // PID not applicable
}
}

private:
uint mPdg = kPiPlus; ///< PDG code of the expected particle
uint mSpecies = track::PID::Pion; ///< Expected species of the track

// TPC
float mPtTPCMin = 0.; ///< minimum pT for TPC PID [GeV/c]
Expand All @@ -503,6 +575,10 @@ class TrackSelectorPID
float mNSigmaRICHMax = 3.; ///< maximum number of RICH σ
float mNSigmaRICHMinCondTOF = -1000.; ///< minimum number of RICH σ if combined with TOF
float mNSigmaRICHMaxCondTOF = 1000.; ///< maximum number of RICH σ if combined with TOF

// Bayesian
float mPtBayesMin = 0.; ///< minimum pT for Bayesian PID [GeV/c]
float mPtBayesMax = 100.; ///< maximum pT for Bayesian PID [GeV/c]
};

#endif // O2_ANALYSIS_TRACKSELECTORPID_H_
59 changes: 52 additions & 7 deletions PWGHF/TableProducer/HFLcCandidateSelector.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ struct HFLcCandidateSelector {

Configurable<double> d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"};
Configurable<double> d_pTCandMax{"d_pTCandMax", 36., "Upper bound of candidate pT"};
Configurable<bool> d_FilterPID{"d_FilterPID", true, "Bool to use or not the PID at filtering level"};
Configurable<bool> d_FilterPID{"d_FilterPID", true, "Bool to use or not the PID based on nSigma cut at filtering level"};
// TPC
Configurable<double> d_pidTPCMinpT{"d_pidTPCMinpT", 0.1, "Lower bound of track pT for TPC PID"};
Configurable<double> d_pidTPCMaxpT{"d_pidTPCMaxpT", 1., "Upper bound of track pT for TPC PID"};
Expand All @@ -45,6 +45,10 @@ struct HFLcCandidateSelector {
Configurable<double> d_pidTOFMaxpT{"d_pidTOFMaxpT", 2.5, "Upper bound of track pT for TOF PID"};
Configurable<double> d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"};
Configurable<double> d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 5., "Nsigma cut on TOF combined with TPC"};
// Bayesian
Configurable<bool> d_BayesPID{"d_BayesPID", true, "Bool to use or not the PID based on Bayesian probability cut at filtering level"};
Configurable<double> d_pidBayesMinpT{"d_pidBayesMinpT", 0., "Lower bound of track pT for Bayesian PID"};
Configurable<double> d_pidBayesMaxpT{"d_pidBayesMaxpT", 100, "Upper bound of track pT for Bayesian PID"};
// topological cuts
Configurable<std::vector<double>> pTBins{"pTBins", std::vector<double>{hf_cuts_lc_topkpi::pTBins_v}, "pT bin limits"};
Configurable<LabeledArray<double>> cuts{"Lc_to_p_K_pi_cuts", {hf_cuts_lc_topkpi::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Lc candidate selection per pT bin"};
Expand Down Expand Up @@ -132,7 +136,9 @@ struct HFLcCandidateSelector {
return true;
}

void process(aod::HfCandProng3 const& candidates, aod::BigTracksPID const&)
using TrksPID = soa::Join<aod::BigTracksPID, aod::pidBayesPi, aod::pidBayesKa, aod::pidBayesPr, aod::pidBayes>;

void process(aod::HfCandProng3 const& candidates, TrksPID const&)
{
TrackSelectorPID selectorPion(kPiPlus);
selectorPion.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT);
Expand All @@ -141,6 +147,7 @@ struct HFLcCandidateSelector {
selectorPion.setRangePtTOF(d_pidTOFMinpT, d_pidTOFMaxpT);
selectorPion.setRangeNSigmaTOF(-d_nSigmaTOF, d_nSigmaTOF);
selectorPion.setRangeNSigmaTOFCondTPC(-d_nSigmaTOFCombined, d_nSigmaTOFCombined);
selectorPion.setRangePtBayes(d_pidBayesMinpT, d_pidBayesMaxpT);

TrackSelectorPID selectorKaon(selectorPion);
selectorKaon.setPDG(kKPlus);
Expand All @@ -160,9 +167,9 @@ struct HFLcCandidateSelector {
continue;
}

auto trackPos1 = candidate.index0_as<aod::BigTracksPID>(); // positive daughter (negative for the antiparticles)
auto trackNeg = candidate.index1_as<aod::BigTracksPID>(); // negative daughter (positive for the antiparticles)
auto trackPos2 = candidate.index2_as<aod::BigTracksPID>(); // positive daughter (negative for the antiparticles)
auto trackPos1 = candidate.index0_as<TrksPID>(); // positive daughter (negative for the antiparticles)
auto trackNeg = candidate.index1_as<TrksPID>(); // negative daughter (positive for the antiparticles)
auto trackPos2 = candidate.index2_as<TrksPID>(); // positive daughter (negative for the antiparticles)

/*
// daughter track validity selection
Expand Down Expand Up @@ -192,6 +199,8 @@ struct HFLcCandidateSelector {

auto pidLcpKpi = -1;
auto pidLcpiKp = -1;
auto pidBayesLcpKpi = -1;
auto pidBayesLcpiKp = -1;

if (!d_FilterPID) {
// PID non applied
Expand Down Expand Up @@ -225,15 +234,51 @@ struct HFLcCandidateSelector {
}
}

if (!d_BayesPID) {
// PID non applied
pidBayesLcpKpi = 1;
pidBayesLcpiKp = 1;
} else {
int pidBayesTrackPos1Proton = selectorProton.getStatusTrackBayesPID(trackPos1);
int pidBayesTrackPos2Proton = selectorProton.getStatusTrackBayesPID(trackPos2);
int pidBayesTrackPos1Pion = selectorPion.getStatusTrackBayesPID(trackPos1);
int pidBayesTrackPos2Pion = selectorPion.getStatusTrackBayesPID(trackPos2);
int pidBayesTrackNegKaon = selectorKaon.getStatusTrackBayesPID(trackNeg);

if (pidBayesTrackPos1Proton == TrackSelectorPID::Status::PIDAccepted &&
pidBayesTrackNegKaon == TrackSelectorPID::Status::PIDAccepted &&
pidBayesTrackPos2Pion == TrackSelectorPID::Status::PIDAccepted) {
pidBayesLcpKpi = 1; // accept LcpKpi
} else if (pidBayesTrackPos1Proton == TrackSelectorPID::Status::PIDRejected ||
pidBayesTrackNegKaon == TrackSelectorPID::Status::PIDRejected ||
pidBayesTrackPos2Pion == TrackSelectorPID::Status::PIDRejected) {
pidBayesLcpKpi = 0; // exclude LcpKpi
}
if (pidBayesTrackPos2Proton == TrackSelectorPID::Status::PIDAccepted &&
pidBayesTrackNegKaon == TrackSelectorPID::Status::PIDAccepted &&
pidBayesTrackPos1Pion == TrackSelectorPID::Status::PIDAccepted) {
pidBayesLcpiKp = 1; // accept LcpiKp
} else if (pidBayesTrackPos1Pion == TrackSelectorPID::Status::PIDRejected ||
pidBayesTrackNegKaon == TrackSelectorPID::Status::PIDRejected ||
pidBayesTrackPos2Proton == TrackSelectorPID::Status::PIDRejected) {
pidBayesLcpiKp = 0; // exclude LcpiKp
}
}

if (pidLcpKpi == 0 && pidLcpiKp == 0) {
hfSelLcCandidate(statusLcpKpi, statusLcpiKp);
continue;
}

if ((pidLcpKpi == -1 || pidLcpKpi == 1) && topolLcpKpi) {
if (pidBayesLcpKpi == 0 && pidBayesLcpiKp == 0) {
hfSelLcCandidate(statusLcpKpi, statusLcpiKp);
continue;
}

if ((pidLcpKpi == -1 || pidLcpKpi == 1) && (pidBayesLcpKpi == -1 || pidBayesLcpKpi == 1) && topolLcpKpi) {
statusLcpKpi = 1; // identified as LcpKpi
}
if ((pidLcpiKp == -1 || pidLcpiKp == 1) && topolLcpiKp) {
if ((pidLcpiKp == -1 || pidLcpiKp == 1) && (pidBayesLcpiKp == -1 || pidBayesLcpiKp == 1) && topolLcpiKp) {
statusLcpiKp = 1; // identified as LcpiKp
}

Expand Down
12 changes: 10 additions & 2 deletions PWGHF/Tasks/taskLc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,10 @@ struct TaskLc {
{"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}},
{"hCPARecSig", "3-prong candidates (matched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
{"hCPARecBg", "3-prong candidates (unmatched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
{"hDecLengthRecSig", "3-prong candidates (matched);decay length (cm);entries", {HistType::kTH1F, {{200, 0, 1.0}}}},
{"hDecLengthRecBg", "3-prong candidates (unmatched);decay length (cm);entries", {HistType::kTH1F, {{200, 0, 1.0}}}},
{"hDca2RecSig", "3-prong candidates (matched);prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH1F, {{100, 0, 0.5}}}},
{"hDca2RecBg", "3-prong candidates (unmatched);prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH1F, {{100, 0, 0.5}}}},
{"hEtaRecSig", "3-prong candidates (matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}},
{"hEtaRecBg", "3-prong candidates (unmatched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}},
{"hEtaGen", "MC particles (matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}}};
Expand All @@ -66,7 +70,7 @@ struct TaskLc {
registry.add("hSelectionStatus", "3-prong candidates;selection status;entries", {HistType::kTH2F, {{5, -0.5, 4.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
registry.add("hImpParErr", "3-prong candidates;impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
registry.add("hDecLenErr", "3-prong candidates;decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
registry.add("hDCA2", "3-prong candidates;prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
registry.add("hDca2", "3-prong candidates;prong Chi2PCA to sec. vertex (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});

registry.add("hPtRecSig", "3-prong candidates (matched);#it{p}_{T}^{rec.} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}});
registry.add("hPtRecSigPrompt", "3-prong candidates (matched, prompt);#it{p}_{T}^{rec.} (GeV/#it{c});entries", {HistType::kTH1F, {{360, 0., 36.}}});
Expand Down Expand Up @@ -126,7 +130,7 @@ struct TaskLc {
registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1(), candidate.pt());
registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter2(), candidate.pt());
registry.fill(HIST("hDecLenErr"), candidate.errorDecayLength(), candidate.pt());
registry.fill(HIST("hDecLenErr"), candidate.chi2PCA(), candidate.pt());
registry.fill(HIST("hDca2"), candidate.chi2PCA(), candidate.pt());
}
}

Expand Down Expand Up @@ -156,10 +160,14 @@ struct TaskLc {
registry.fill(HIST("hPtRecSigNonPrompt"), ptRec); // rec. level pT, non-prompt
}
registry.fill(HIST("hCPARecSig"), candidate.cpa());
registry.fill(HIST("hDecLengthRecSig"), candidate.decayLength());
registry.fill(HIST("hDca2RecSig"), candidate.chi2PCA());
registry.fill(HIST("hEtaRecSig"), candidate.eta());
} else {
registry.fill(HIST("hPtRecBg"), candidate.pt());
registry.fill(HIST("hCPARecBg"), candidate.cpa());
registry.fill(HIST("hDecLengthRecBg"), candidate.decayLength());
registry.fill(HIST("hDca2RecBg"), candidate.chi2PCA());
registry.fill(HIST("hEtaRecBg"), candidate.eta());
}
}
Expand Down