diff --git a/PWGHF/Core/SelectorCuts.h b/PWGHF/Core/SelectorCuts.h index 0351681ce23..c051aa238b5 100644 --- a/PWGHF/Core/SelectorCuts.h +++ b/PWGHF/Core/SelectorCuts.h @@ -84,13 +84,23 @@ static const std::vector 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 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 labelsCutsTrack = {"nSigmaMinIts", "minItsClusterSizes", "minItsCluster", "minItsIbCluster", "minTpcCluster", "minTpcRow", "minTpcCrossedOverFound", "maxTpcShared", "maxTpcFracShared", "maxTPCnSigmaBB"}; static const std::vector 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 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..c2ab79f9aeb 100644 --- a/PWGHF/TableProducer/trackIndexSkimCreator.cxx +++ b/PWGHF/TableProducer/trackIndexSkimCreator.cxx @@ -65,6 +65,7 @@ #include #include #include +#include #include #include // for PV refit @@ -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], 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> 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 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,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 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 itsPidNsigma = -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(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(minItsClusterSizes)) { + if (track.itsClusterSizes() < static_cast(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 + 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. + 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(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(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 +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(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) {