diff --git a/PWGHF/Core/SelectorCuts.h b/PWGHF/Core/SelectorCuts.h index 0351681ce23..6f1e3044ad9 100644 --- a/PWGHF/Core/SelectorCuts.h +++ b/PWGHF/Core/SelectorCuts.h @@ -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 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 labelsCutsTrack = {"nSigmaMinIts", "minItsClusterSizes", "minItsCluster", "minItsIbCluster", "minTpcCluster", "minTpcRow", "minTpcCrossedOverFound", "maxTpcShared", "maxTpcFracShared", "maxTPCnSigmaBB"}; static const std::vector 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 labelsBetheBlochParams = {"p0", "p1", "p2", "p3", "p4", "resolution"}; + } // namespace hf_presel_lightnuclei namespace hf_cuts_bdt_multiclass diff --git a/PWGHF/TableProducer/trackIndexSkimCreator.cxx b/PWGHF/TableProducer/trackIndexSkimCreator.cxx index 396e8dea819..61be9a3a589 100644 --- a/PWGHF/TableProducer/trackIndexSkimCreator.cxx +++ b/PWGHF/TableProducer/trackIndexSkimCreator.cxx @@ -45,6 +45,7 @@ #include "Common/DataModel/TrackSelectionTables.h" #include "Tools/ML/MlResponse.h" +#include "MathUtils/BetheBlochAleph.h" #include // for PV refit #include #include @@ -1293,7 +1294,7 @@ struct HfTrackIndexSkimCreator { Configurable> 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> binsPtCtToTrKPi{"binsPtCtToTrKPi", std::vector{hf_cuts_presel_3prong::vecBinsPt}, "pT bin limits for Ct->tKpi pT-dependent cuts"}; - Configurable> 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> 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> binsPtChToHeKPi{"binsPtChToHeKPi", std::vector{hf_cuts_presel_3prong::vecBinsPt}, "pT bin limits for Ch->heKpi pT-dependent cuts"}; Configurable> 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"}; @@ -1302,7 +1303,10 @@ struct HfTrackIndexSkimCreator { Configurable> 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> 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> 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> 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 applyProtonPidForLcToPKPi{"applyProtonPidForLcToPKPi", false, "Apply proton PID for Lc->pKpi"}; @@ -1311,6 +1315,8 @@ struct HfTrackIndexSkimCreator { Configurable applyDeuteronPidForCdToDeKPi{"applyDeuteronPidForCdToDeKPi", false, "Require deuteron PID for Cd->deKpi"}; Configurable applyTritonPidForCtToTrKPi{"applyTritonPidForCtToTrKPi", false, "Require triton PID for Ct->tKpi"}; Configurable applyHeliumPidForChToHeKPi{"applyHeliumPidForChToHeKPi", false, "Require helium3 PID for Ch->heKpi"}; + Configurable applyLightNucleiTpcPidBasedOnBB{"applyLightNucleiTpcPidBasedOnBB", false, "Apply TPC PID for light nuclei using Bethe–Bloch parameterization"}; + // lightnuclei track selection for charmnuclei Configurable applyNucleiSelcForCharmNuclei{"applyNucleiSelcForCharmNuclei", false, "Require track selection for charm nuclei"}; // ML models for triggers @@ -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)"}; @@ -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 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(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); @@ -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; } @@ -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 + float getTPCnSigmaBB(const TrackType& track, ChannelsNucleiQA lightnuclei) + { + if (!track.hasTPC()) { + return -999.f; + } + + const int row = static_cast(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(charge) * static_cast(rigidity) / mass; + const double expBethe = common::BetheBlochAleph(x, bb0, bb1, bb2, bb3, bb4); + const double expSigma = expBethe * static_cast(relRes); + + if (expSigma <= 0.) { + return -999.f; + } + + return static_cast((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 @@ -1762,27 +1853,20 @@ 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 = + (iDecay3P == hf_cand_3prong::DecayType::CdToDeKPi) ? ChannelsNucleiQA::Deuteron : (iDecay3P == hf_cand_3prong::DecayType::CtToTrKPi) ? ChannelsNucleiQA::Triton + : ChannelsNucleiQA::Helium3; // 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) {