Skip to content

[PWGHF] Implement TPC pid for light nuclei based on Bethe-Bloch parametrization in HF track skimming task #14729

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

Draft
zhangbiao-phy wants to merge 5 commits into
base: master
Choose a base branch
from
+122 −24
Draft

[PWGHF] Implement TPC pid for light nuclei based on Bethe-Bloch parametrization in HF track skimming task #14729

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
027e597
enabled TPC pid for light nuclei based on Bethe-Bloch parametrization
zhangbiao-phy Jan 31, 2026
590e563
Add PID QA plots
zhangbiao-phy Jan 31, 2026
5200e72
fix typo
zhangbiao-phy Jan 31, 2026
129256a
fix warning
zhangbiao-phy Jan 31, 2026
01c6f00
fix warning
zhangbiao-phy Jan 31, 2026
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
14 changes: 10 additions & 4 deletions PWGHF/Core/SelectorCuts.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,12 +85,18 @@ namespace hf_presel_lightnuclei
{

// default values for the track cuts for lightnuclei in the track-index-skim-creator
constexpr float CutsTrackQuality[3][9] = {{-4, 3, 5., 0., 100, 100, 0.83, 160., 1.},
{-4, 3, 5., 0., 100, 100, 0.83, 160., 1.},
{-4, 3, 5., 0., 100, 100, 0.83, 160., 1.}};
static const std::vector<std::string> labelsCutsTrack = {"nSigmaMinIts", "minItsClusterSizes", "minItsCluster", "minItsIbCluster", "minTpcCluster", "minTpcRow", "minTpcCrossedOverFound", "maxTpcShared", "maxTpcFracShared"};
constexpr float CutsTrackQuality[3][10] = {{-4, 3, 5., 0., 100, 100, 0.83, 160., 1., 5},
{-4, 3, 5., 0., 100, 100, 0.83, 160., 1., 5},
{-4, 3, 5., 0., 100, 100, 0.83, 160., 1., 5}};
static const std::vector<std::string> labelsCutsTrack = {"nSigmaMinIts", "minItsClusterSizes", "minItsCluster", "minItsIbCluster", "minTpcCluster", "minTpcRow", "minTpcCrossedOverFound", "maxTpcShared", "maxTpcFracShared", "maxTPCnSigmaBB"};
static const std::vector<std::string> labelsRowsNucleiType = {"Deutron", "Triton", "Helium3"};

constexpr float BetheBlochParams[3][6] = {{5.39302, 7.859534, 0.004048, 2.323197, 1.609307, 0.09},
{5.39302, 7.859534, 0.004048, 2.323197, 1.609307, 0.09},
{-126.55736, -0.858569, 1.11164, 1.21032, 2.656374, 0.09}};

static const std::vector<std::string> labelsBetheBlochParams = {"p0", "p1", "p2", "p3", "p4", "resolution"};

} // namespace hf_presel_lightnuclei

namespace hf_cuts_bdt_multiclass
Expand Down
132 changes: 112 additions & 20 deletions PWGHF/TableProducer/trackIndexSkimCreator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
#include "Common/DataModel/TrackSelectionTables.h"
#include "Tools/ML/MlResponse.h"

#include "MathUtils/BetheBlochAleph.h"
#include <CCDB/BasicCCDBManager.h> // for PV refit
#include <CCDB/CcdbApi.h>
#include <CommonConstants/PhysicsConstants.h>
Expand Down Expand Up @@ -1293,7 +1294,7 @@ struct HfTrackIndexSkimCreator {
Configurable<LabeledArray<double>> cutsCdToDeKPi{"cutsCdToDeKPi", {hf_cuts_presel_3prong::Cuts[0], hf_cuts_presel_3prong::NBinsPt, hf_cuts_presel_3prong::NCutVars, hf_cuts_presel_3prong::labelsPt, hf_cuts_presel_3prong::labelsCutVar}, "Cd->deKpi selections per pT bin"};
// Ct cuts
Configurable<std::vector<double>> binsPtCtToTrKPi{"binsPtCtToTrKPi", std::vector<double>{hf_cuts_presel_3prong::vecBinsPt}, "pT bin limits for Ct->tKpi pT-dependent cuts"};
Configurable<LabeledArray<double>> cutsCtToTrKPi{"cutsCdToTrKPi", {hf_cuts_presel_3prong::Cuts[0], hf_cuts_presel_3prong::NBinsPt, hf_cuts_presel_3prong::NCutVars, hf_cuts_presel_3prong::labelsPt, hf_cuts_presel_3prong::labelsCutVar}, "Ct->tKpi selections per pT bin"};
Configurable<LabeledArray<double>> cutsCtToTrKPi{"cutsCtToTrKPi", {hf_cuts_presel_3prong::Cuts[0], hf_cuts_presel_3prong::NBinsPt, hf_cuts_presel_3prong::NCutVars, hf_cuts_presel_3prong::labelsPt, hf_cuts_presel_3prong::labelsCutVar}, "Ct->tKpi selections per pT bin"};
// Ch cuts
Configurable<std::vector<double>> binsPtChToHeKPi{"binsPtChToHeKPi", std::vector<double>{hf_cuts_presel_3prong::vecBinsPt}, "pT bin limits for Ch->heKpi pT-dependent cuts"};
Configurable<LabeledArray<double>> cutsChToHeKPi{"cutsChToHeKPi", {hf_cuts_presel_3prong::Cuts[0], hf_cuts_presel_3prong::NBinsPt, hf_cuts_presel_3prong::NCutVars, hf_cuts_presel_3prong::labelsPt, hf_cuts_presel_3prong::labelsCutVar}, "Ch->heKpi selections per pT bin"};
Expand All @@ -1302,7 +1303,10 @@ struct HfTrackIndexSkimCreator {
Configurable<LabeledArray<double>> cutsDstarToD0Pi{"cutsDstarToD0Pi", {hf_cuts_presel_dstar::Cuts[0], hf_cuts_presel_dstar::NBinsPt, hf_cuts_presel_dstar::NCutVars, hf_cuts_presel_dstar::labelsPt, hf_cuts_presel_dstar::labelsCutVar}, "D*+->D0pi selections per pT bin"};

// CharmNuclei track selection
Configurable<LabeledArray<float>> selectionsLightNuclei{"selectionsLightNuclei", {hf_presel_lightnuclei::CutsTrackQuality[0], 3, 9, hf_presel_lightnuclei::labelsRowsNucleiType, hf_presel_lightnuclei::labelsCutsTrack}, "nuclei track selections for deuteron / triton / helium applied if proper process function enabled"};
Configurable<LabeledArray<float>> selectionsLightNuclei{"selectionsLightNuclei", {hf_presel_lightnuclei::CutsTrackQuality[0], 3, 10, hf_presel_lightnuclei::labelsRowsNucleiType, hf_presel_lightnuclei::labelsCutsTrack}, "nuclei track selections for deuteron / triton / helium applied if proper process function enabled"};
Configurable<LabeledArray<float>> tpcPidBBParamsLightNuclei{"tpcPidBBParamsLightNuclei", {hf_presel_lightnuclei::BetheBlochParams[0], 3, 6, hf_presel_lightnuclei::labelsRowsNucleiType, hf_presel_lightnuclei::labelsBetheBlochParams},
"TPC PID Bethe–Bloch parameter configurations for light nuclei "
"(deuteron, triton, helium-3), used in BB-based PID when enabled"};

// proton PID selections for Lc and Xic
Configurable<bool> applyProtonPidForLcToPKPi{"applyProtonPidForLcToPKPi", false, "Apply proton PID for Lc->pKpi"};
Expand All @@ -1311,6 +1315,8 @@ struct HfTrackIndexSkimCreator {
Configurable<bool> applyDeuteronPidForCdToDeKPi{"applyDeuteronPidForCdToDeKPi", false, "Require deuteron PID for Cd->deKpi"};
Configurable<bool> applyTritonPidForCtToTrKPi{"applyTritonPidForCtToTrKPi", false, "Require triton PID for Ct->tKpi"};
Configurable<bool> applyHeliumPidForChToHeKPi{"applyHeliumPidForChToHeKPi", false, "Require helium3 PID for Ch->heKpi"};
Configurable<bool> applyLightNucleiTpcPidBasedOnBB{"applyLightNucleiTpcPidBasedOnBB", false, "Apply TPC PID for light nuclei using Bethe–Bloch parameterization"};

// lightnuclei track selection for charmnuclei
Configurable<bool> applyNucleiSelcForCharmNuclei{"applyNucleiSelcForCharmNuclei", false, "Require track selection for charm nuclei"};
// ML models for triggers
Expand Down Expand Up @@ -1477,7 +1483,9 @@ struct HfTrackIndexSkimCreator {
registry.add("hMassCdToDeKPi", "C Deuteron candidates;inv. mass (De K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{500, 0., 5.}}});
registry.add("hMassCtToTrKPi", "C Triton candidates;inv. mass (Tr K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{500, 0., 5.}}});
registry.add("hMassChToHeKPi", "C Helium3 candidates;inv. mass (He3 K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{500, 0., 5.}}});

if (config.applyNucleiSelcForCharmNuclei && config.applyLightNucleiTpcPidBasedOnBB) {
registry.add("hTPCSignalsLightNuclei", "Light Nuclei TPC signal (a.u.)", {HistType::kTH2D, {{2000, -10., 10.}, {1000, 0., 2000.}}});
}
// needed for PV refitting
if (doprocess2And3ProngsWithPvRefit || doprocess2And3ProngsWithPvRefitWithPidForHfFiltersBdt) {
const AxisSpec axisCollisionX{100, -20.f, 20.f, "X (cm)"};
Expand Down Expand Up @@ -1654,18 +1662,36 @@ struct HfTrackIndexSkimCreator {
}

/// Apply track-quality (ITS/TPC) + optional ITS-PID preselection for light-nucleus daughters used in charm-nuclei 3-prong channels (Cd/Ct/Ch).
/// \tparam TrackType Track/ASoA row type providing ITS/TPC quality accessors.
/// \tparam TrackType Track providing ITS/TPC quality accessors.
/// \param track Daughter track to be tested (either prong0 or prong2).
/// \param lightnuclei Species selector 0: Deuteron, 1: Triton, 2: Helium3.
/// \param nSigmaItsPid ITS nσ value for the selected light nucleus hypothesis
/// \param lightnuclei Species selector: 0=Deuteron, 1=Triton, 2=Helium3.
/// \return true if the track passes all enabled selections.
template <typename TrackType>
bool applyTrackSelectionForCharmNuclei(const TrackType& track,
ChannelsNucleiQA lightnuclei,
float nSigmaItsPid)
ChannelsNucleiQA lightnuclei)
{
// Row index in the selection table: 0 (De), 1 (Tr), 2 (He3)
const int row = static_cast<int>(lightnuclei);
const float nSigmaItsNuclei = nSigmaItsPid;
if (row < 0 || row >= NChannelsLightNucleiPid) {
return false;
}

float nSigmaItsNuclei = -999.f;

switch (lightnuclei) {
case ChannelsNucleiQA::Deuteron:
nSigmaItsNuclei = track.itsNSigmaDe();
break;
case ChannelsNucleiQA::Triton:
nSigmaItsNuclei = track.itsNSigmaTr();
break;
case ChannelsNucleiQA::Helium3:
nSigmaItsNuclei = track.itsNSigmaHe();
break;
default:
return false;
}

// Load cuts for the selected species.
const float minItsNSigmaPid = config.selectionsLightNuclei->get(row, 0u);
const int minItsClusterSizes = config.selectionsLightNuclei->get(row, 1u);
Expand All @@ -1677,6 +1703,9 @@ struct HfTrackIndexSkimCreator {
const int maxTpcShared = config.selectionsLightNuclei->get(row, 7u);
const float maxTpcFracShared = config.selectionsLightNuclei->get(row, 8u);

// Optional: BB-based TPC nσ selection (only if enabled)
const float maxTPCnSigmaBB = config.selectionsLightNuclei->get(row, 9u);

if (nSigmaItsNuclei < minItsNSigmaPid) {
return false;
}
Expand Down Expand Up @@ -1704,9 +1733,71 @@ struct HfTrackIndexSkimCreator {
if (track.tpcFractionSharedCls() > maxTpcFracShared) {
return false;
}

if (config.applyLightNucleiTpcPidBasedOnBB) {
const float nSigmaTpcNuclei = getTPCnSigmaBB(track, lightnuclei);
if (nSigmaTpcNuclei < -999.f) { // invalid marker
return false;
}
// Correct inequality: reject if |nσ| exceeds allowed window
if (std::abs(nSigmaTpcNuclei) > maxTPCnSigmaBB) {
return false;
}
}

return true;
}

/// Compute TPC nσ for light nuclei (De/Tr/He3) using a Bethe–Bloch parameter configuration (BB-based PID).
///
/// \tparam TrackType Track/ASoA row type providing TPC accessors.
/// \param track Track to be tested.
/// \param lightnuclei Species selector: 0=Deuteron, 1=Triton, 2=Helium3.
/// \return TPC nσ for the chosen nucleus hypothesis (or -999 if not applicable).
template <typename TrackType>
float getTPCnSigmaBB(const TrackType& track, ChannelsNucleiQA lightnuclei)
{
if (!track.hasTPC()) {
return -999.f;
}

const int row = static_cast<int>(lightnuclei);
if (row < 0 || row >= NChannelsLightNucleiPid) {
return -999.f;
}

// Columns: [0..4] BB params, [5] relative resolution (sigma/mean)
const double bb0 = config.tpcPidBBParamsLightNuclei->get(row, 0u);
const double bb1 = config.tpcPidBBParamsLightNuclei->get(row, 1u);
const double bb2 = config.tpcPidBBParamsLightNuclei->get(row, 2u);
const double bb3 = config.tpcPidBBParamsLightNuclei->get(row, 3u);
const double bb4 = config.tpcPidBBParamsLightNuclei->get(row, 4u);
const double relRes = config.tpcPidBBParamsLightNuclei->get(row, 5u);

if (relRes <= 0.f) {
return -999.f;
}

// Mass/charge hypothesis for the selected nucleus.
const double mass =
(lightnuclei == ChannelsNucleiQA::Deuteron) ? MassDeuteron : (lightnuclei == ChannelsNucleiQA::Triton) ? MassTriton
: MassHelium3;

const int charge = (lightnuclei == ChannelsNucleiQA::Helium3) ? 2 : 1;

const float rigidity = track.tpcInnerParam(); // p / |q| note: here we didn't apply rigidity correction

const double x = static_cast<double>(charge) * static_cast<double>(rigidity) / mass;
const double expBethe = common::BetheBlochAleph(x, bb0, bb1, bb2, bb3, bb4);
const double expSigma = expBethe * static_cast<double>(relRes);

if (expSigma <= 0.) {
return -999.f;
}

return static_cast<float>((track.tpcSignal() - expBethe) / expSigma);
}

/// Method to perform selections on difference from nominal mass for phi decay
/// \param binPt pt bin for the cuts
/// \param pVecTrack0 is the momentum array of the first daughter track
Expand Down Expand Up @@ -1762,27 +1853,28 @@ struct HfTrackIndexSkimCreator {
iDecay3P == hf_cand_3prong::DecayType::CtToTrKPi ||
iDecay3P == hf_cand_3prong::DecayType::ChToHeKPi)) {

ChannelsNucleiQA nucleiType = ChannelsNucleiQA::Deuteron;
float nSigmaItsNuclei0 = track0.itsNSigmaDe();
float nSigmaItsNuclei2 = track2.itsNSigmaDe();

if (iDecay3P == hf_cand_3prong::DecayType::CtToTrKPi) {
ChannelsNucleiQA nucleiType;
if (iDecay3P == hf_cand_3prong::DecayType::CdToDeKPi) {
nucleiType = ChannelsNucleiQA::Deuteron;
} else if (iDecay3P == hf_cand_3prong::DecayType::CtToTrKPi) {
nucleiType = ChannelsNucleiQA::Triton;
nSigmaItsNuclei0 = track0.itsNSigmaTr();
nSigmaItsNuclei2 = track2.itsNSigmaTr();
} else if (iDecay3P == hf_cand_3prong::DecayType::ChToHeKPi) {
nucleiType = ChannelsNucleiQA::Helium3;
nSigmaItsNuclei0 = track0.itsNSigmaHe();
nSigmaItsNuclei2 = track2.itsNSigmaHe();
} else {
LOGF(fatal, "Unhandled DecayType = %d", static_cast<int>(iDecay3P));
}

// hypo0: nucleus on track0
if (!applyTrackSelectionForCharmNuclei(track0, nucleiType, nSigmaItsNuclei0)) {
if (!applyTrackSelectionForCharmNuclei(track0, nucleiType)) {
CLRBIT(whichHypo[iDecay3P], 0);
} else {
registry.fill(HIST("hTPCSignalsLightNuclei"), track0.tpcInnerParam() * track0.sign(), track0.tpcSignal());
}
// hypo1: nucleus on track2
if (!applyTrackSelectionForCharmNuclei(track2, nucleiType, nSigmaItsNuclei2)) {
if (!applyTrackSelectionForCharmNuclei(track2, nucleiType)) {
CLRBIT(whichHypo[iDecay3P], 1);
} else {
registry.fill(HIST("hTPCSignalsLightNuclei"), track2.tpcInnerParam() * track2.sign(), track2.tpcSignal());
}

if (whichHypo[iDecay3P] == 0) {
Expand Down
Loading