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 7 commits into
base: master
Choose a base branch
from
+161 −46
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
7 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
309a35b
fix error
zhangbiao-phy Feb 1, 2026
457ae60
fix comments from vit
zhangbiao-phy Feb 1, 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
18 changes: 14 additions & 4 deletions PWGHF/Core/SelectorCuts.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,13 +84,23 @@ static const std::vector<std::string> labelsRowsPid = {"ProtonInLcToPKPi", "Prot
namespace hf_presel_lightnuclei
{

static constexpr int NParticleRows = 3; // number of particles / rows
static constexpr int NVarCuts = 10; // number of cuts for each particles
static constexpr int NBetheBlochParams = 6; // number of parameters for Bethe-Bloch

// 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[NParticleRows][NVarCuts] = {{-4.f, 3.f, 5.f, 0.f, 100.f, 100.f, 0.83, 160.f, 1.f, 5.f},
{-4.f, 3.f, 5.f, 0.f, 100.f, 100.f, 0.83, 160.f, 1.f, 5.f},
{-4.f, 3.f, 5.f, 0.f, 100.f, 100.f, 0.83, 160.f, 1.f, 5.f}};
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[NParticleRows][NBetheBlochParams] = {{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
189 changes: 147 additions & 42 deletions PWGHF/TableProducer/trackIndexSkimCreator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@
#include <Framework/InitContext.h>
#include <Framework/Logger.h>
#include <Framework/runDataProcessing.h>
#include <MathUtils/BetheBlochAleph.h>
#include <ReconstructionDataFormats/Track.h>
#include <ReconstructionDataFormats/Vertex.h> // for PV refit

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], hf_presel_lightnuclei::NParticleRows, hf_presel_lightnuclei::NVarCuts, 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], hf_presel_lightnuclei::NParticleRows, hf_presel_lightnuclei::NBetheBlochParams, 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,59 +1662,149 @@ 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 itsPidNsigma = -999.f;
Copy link
Collaborator

@vkucera vkucera Feb 2, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

its is not a quantity, Pid is redundant, Nsigma is not a word.

Suggested change
float itsPidNsigma = -999.f;
float nSigmaIts = -999.f;


switch (lightnuclei) {
case ChannelsNucleiQA::Deuteron:
itsPidNsigma = track.itsNSigmaDe();
break;
case ChannelsNucleiQA::Triton:
itsPidNsigma = track.itsNSigmaTr();
break;
case ChannelsNucleiQA::Helium3:
itsPidNsigma = track.itsNSigmaHe();
break;
default:
LOG(fatal) << "Unhandled ChannelsNucleiQA " << static_cast<int>(lightnuclei);
}

// Load cuts for the selected species.
const float minItsNSigmaPid = config.selectionsLightNuclei->get(row, 0u);
const int minItsClusterSizes = config.selectionsLightNuclei->get(row, 1u);
const int minItsCluster = config.selectionsLightNuclei->get(row, 2u);
const int minItsIbCluster = config.selectionsLightNuclei->get(row, 3u);
const int minTpcCluster = config.selectionsLightNuclei->get(row, 4u);
const int minTpcRow = config.selectionsLightNuclei->get(row, 5u);
const float minTpcCrossedOverFound = config.selectionsLightNuclei->get(row, 6u);
const int maxTpcShared = config.selectionsLightNuclei->get(row, 7u);
const float maxTpcFracShared = config.selectionsLightNuclei->get(row, 8u);

if (nSigmaItsNuclei < minItsNSigmaPid) {
const float itsPidNsigmaMin = config.selectionsLightNuclei->get(row, 0u);
const float itsClusterSizeMin = config.selectionsLightNuclei->get(row, 1u);
const float itsClusterMin = config.selectionsLightNuclei->get(row, 2u);
const float itsIbClusterMin = config.selectionsLightNuclei->get(row, 3u);
const float tpcClusterMin = config.selectionsLightNuclei->get(row, 4u);
const float tpcCrossedRowsMin = config.selectionsLightNuclei->get(row, 5u);
const float tpcCrossedRowsOverFindMin = config.selectionsLightNuclei->get(row, 6u);
const float tpcSharedMax = config.selectionsLightNuclei->get(row, 7u);
const float tpcFracSharedMax = config.selectionsLightNuclei->get(row, 8u);

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

if (itsPidNsigma < itsPidNsigmaMin) {
return false;
}
if (track.itsClusterSizes() < static_cast<unsigned int>(minItsClusterSizes)) {
if (track.itsClusterSizes() < static_cast<unsigned int>(itsClusterSizeMin)) {
return false;
}
if (track.itsNCls() < minItsCluster) {
if (track.itsNCls() < itsClusterMin) {
return false;
}
if (track.itsNClsInnerBarrel() < minItsIbCluster) {
if (track.itsNClsInnerBarrel() < itsIbClusterMin) {
return false;
}
if (track.tpcNClsFound() < minTpcCluster) {
if (track.tpcNClsFound() < tpcClusterMin) {
return false;
}
if (track.tpcNClsCrossedRows() < minTpcRow) {
if (track.tpcNClsCrossedRows() < tpcCrossedRowsMin) {
return false;
}
if (track.tpcCrossedRowsOverFindableCls() < minTpcCrossedOverFound) {
if (track.tpcCrossedRowsOverFindableCls() < tpcCrossedRowsOverFindMin) {
return false;
}
if (maxTpcShared >= 0 && track.tpcNClsShared() > maxTpcShared) {
if (track.tpcNClsShared() > tpcSharedMax) {
return false;
}
if (track.tpcFractionSharedCls() > maxTpcFracShared) {
if (track.tpcFractionSharedCls() > tpcFracSharedMax) {
return false;
}

if (config.applyLightNucleiTpcPidBasedOnBB) {
const float tpcBbPidNsigma = getTPCnSigmaBB(track, lightnuclei);
if (std::abs(tpcBbPidNsigma) > tpcBbPidNsigmaMax) {
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)
Copy link
Collaborator

@vkucera vkucera Feb 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Follow naming conventions.

{
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.
double mass = 0.;
switch (lightnuclei) {
case ChannelsNucleiQA::Deuteron:
mass = MassDeuteron;
break;
case ChannelsNucleiQA::Triton:
mass = MassTriton;
break;
case ChannelsNucleiQA::Helium3:
mass = MassHelium3;
break;
default:
LOG(fatal) << "Unhandled ChannelsNucleiQA " << static_cast<int>(lightnuclei);
}

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 +1860,34 @@ 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) {
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();
ChannelsNucleiQA nucleiType;

switch (iDecay3P) {
case hf_cand_3prong::DecayType::CdToDeKPi:
nucleiType = ChannelsNucleiQA::Deuteron;
break;
case hf_cand_3prong::DecayType::CtToTrKPi:
nucleiType = ChannelsNucleiQA::Triton;
break;
case hf_cand_3prong::DecayType::ChToHeKPi:
nucleiType = ChannelsNucleiQA::Helium3;
break;
default:
LOG(fatal) << "Unhandled DecayType " << static_cast<int>(iDecay3P);
continue;
}

// 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