diff --git a/Common/Core/TrackSelectorPID.h b/Common/Core/TrackSelectorPID.h index 6c0267c3bab..ccb6ae05c4b 100644 --- a/Common/Core/TrackSelectorPID.h +++ b/Common/Core/TrackSelectorPID.h @@ -31,8 +31,8 @@ class TrackSelectorPID /// Standard constructor with PDG code initialisation explicit TrackSelectorPID(int pdg) - : mPdg(std::abs(pdg)) { + setPDG(pdg); } /// Default destructor @@ -46,7 +46,36 @@ class TrackSelectorPID PIDAccepted }; - void setPDG(int pdg) { mPdg = std::abs(pdg); } + void setPDG(int pdg) + { + 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 @@ -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 + bool isValidTrackBayesPID(const T& track) + { + auto pt = track.pt(); + return (mPtBayesMin <= pt && pt <= mPtBayesMax); + } + + /// Bayesian maximum probability algorithm. + template + 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 + 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] @@ -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_ diff --git a/PWGHF/TableProducer/HFLcCandidateSelector.cxx b/PWGHF/TableProducer/HFLcCandidateSelector.cxx index 808c877cdd3..0a9ebac436c 100644 --- a/PWGHF/TableProducer/HFLcCandidateSelector.cxx +++ b/PWGHF/TableProducer/HFLcCandidateSelector.cxx @@ -33,7 +33,7 @@ struct HFLcCandidateSelector { Configurable d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"}; Configurable d_pTCandMax{"d_pTCandMax", 36., "Upper bound of candidate pT"}; - Configurable d_FilterPID{"d_FilterPID", true, "Bool to use or not the PID at filtering level"}; + Configurable d_FilterPID{"d_FilterPID", true, "Bool to use or not the PID based on nSigma cut at filtering level"}; // TPC Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 0.1, "Lower bound of track pT for TPC PID"}; Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 1., "Upper bound of track pT for TPC PID"}; @@ -45,6 +45,10 @@ struct HFLcCandidateSelector { Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 2.5, "Upper bound of track pT for TOF PID"}; Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; Configurable d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 5., "Nsigma cut on TOF combined with TPC"}; + // Bayesian + Configurable d_BayesPID{"d_BayesPID", true, "Bool to use or not the PID based on Bayesian probability cut at filtering level"}; + Configurable d_pidBayesMinpT{"d_pidBayesMinpT", 0., "Lower bound of track pT for Bayesian PID"}; + Configurable d_pidBayesMaxpT{"d_pidBayesMaxpT", 100, "Upper bound of track pT for Bayesian PID"}; // topological cuts Configurable> pTBins{"pTBins", std::vector{hf_cuts_lc_topkpi::pTBins_v}, "pT bin limits"}; Configurable> cuts{"Lc_to_p_K_pi_cuts", {hf_cuts_lc_topkpi::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Lc candidate selection per pT bin"}; @@ -132,7 +136,9 @@ struct HFLcCandidateSelector { return true; } - void process(aod::HfCandProng3 const& candidates, aod::BigTracksPID const&) + using TrksPID = soa::Join; + + void process(aod::HfCandProng3 const& candidates, TrksPID const&) { TrackSelectorPID selectorPion(kPiPlus); selectorPion.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); @@ -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); @@ -160,9 +167,9 @@ struct HFLcCandidateSelector { continue; } - auto trackPos1 = candidate.index0_as(); // positive daughter (negative for the antiparticles) - auto trackNeg = candidate.index1_as(); // negative daughter (positive for the antiparticles) - auto trackPos2 = candidate.index2_as(); // positive daughter (negative for the antiparticles) + auto trackPos1 = candidate.index0_as(); // positive daughter (negative for the antiparticles) + auto trackNeg = candidate.index1_as(); // negative daughter (positive for the antiparticles) + auto trackPos2 = candidate.index2_as(); // positive daughter (negative for the antiparticles) /* // daughter track validity selection @@ -192,6 +199,8 @@ struct HFLcCandidateSelector { auto pidLcpKpi = -1; auto pidLcpiKp = -1; + auto pidBayesLcpKpi = -1; + auto pidBayesLcpiKp = -1; if (!d_FilterPID) { // PID non applied @@ -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 } diff --git a/PWGHF/Tasks/taskLc.cxx b/PWGHF/Tasks/taskLc.cxx index a79c30973bb..c263e026f48 100644 --- a/PWGHF/Tasks/taskLc.cxx +++ b/PWGHF/Tasks/taskLc.cxx @@ -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.}}}}}}; @@ -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.}}}); @@ -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()); } } @@ -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()); } }