From 4e00e4b0ecb4984bcbcd852a0c6c91d6ed909bb5 Mon Sep 17 00:00:00 2001 From: Biao Zhang Date: Wed, 15 Feb 2023 01:16:09 +0800 Subject: [PATCH 1/2] Add the option for using TPC splines PID calibration --- EventFiltering/PWGHF/HFFilter.cxx | 49 ++++---- EventFiltering/PWGHF/HFFilterHelpers.h | 148 +++++++++++++++++-------- 2 files changed, 128 insertions(+), 69 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilter.cxx b/EventFiltering/PWGHF/HFFilter.cxx index 48ad70f210b..b4e1d93d1f2 100644 --- a/EventFiltering/PWGHF/HFFilter.cxx +++ b/EventFiltering/PWGHF/HFFilter.cxx @@ -113,8 +113,12 @@ struct HfFilter { // Main struct for HF triggers Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable ccdbPathTPC{"ccdbPathTPC", "Users/i/iarsene/Calib/TPCpostCalib", "base path to the ccdb object"}; - // TPC PID post calibrations - Configurable computeTPCPostCalib{"computeTPCPostCalib", false, "If true, compute TPC post-calibrated n-sigmas"}; + // TPC PID calibrations + Configurable setTPCCalib{"setTPCCalib", 0, "0 is not use re-calibrations, 1 is compute TPC post-calibrated n-sigmas, 2 is using TPC Spline"}; + Configurable ccdbBBProton{"ccdbBBProton", "Users/l/lserksny/PIDProton", "Path to the CCDB ocject for proton BB param"}; + Configurable ccdbBBAntiProton{"ccdbBBAntiProton", "Users/l/lserksny/PIDAntiProton", "Path to the CCDB ocject for antiproton BB param"}; + Configurable ccdbBBPion{"ccdbBBPion", "Users/l/lserksny/PIDPion", "Path to the CCDB ocject for Pion BB param"}; + Configurable ccdbBBAntiPion{"ccdbBBAntiPion", "Users/l/lserksny/PIDAntiPion", "Path to the CCDB ocject for antiPion BB param"}; // parameter for Optimisation Tree Configurable applyOptimisation{"applyOptimisation", false, "Flag to enable or disable optimisation"}; @@ -140,12 +144,11 @@ struct HfFilter { // Main struct for HF triggers std::shared_ptr hGammaSelected, hGammaEtaBefore, hGammaEtaAfter; std::shared_ptr hGammaArmPodBefore, hGammaArmPodAfter; - // Histogram of TPC postcalibration map for pion and proton - TH3F* hMapPionMean = nullptr; - TH3F* hMapPionSigma = nullptr; - TH3F* hMapProtonMean = nullptr; - TH3F* hMapProtonSigma = nullptr; - + // Histograms of TPC calibration for pion and proton + std::array hMapPion = {nullptr, nullptr}; + std::array hMapProton = {nullptr, nullptr}; + std::array, 2> hBBProton{}; + std::array, 2> hBBPion{}; // ONNX std::array, kNCharmParticles> sessionML = {nullptr, nullptr, nullptr, nullptr, nullptr}; std::array>, kNCharmParticles> inputShapesML{}; @@ -275,20 +278,26 @@ struct HfFilter { // Main struct for HF triggers o2::base::Propagator::initFieldFromGRP(grpo); // needed for TPC PID postcalibrations - if (computeTPCPostCalib) { + if (setTPCCalib == 1) { auto calibList = ccdb->getForTimeStamp(ccdbPathTPC.value, bc.timestamp()); if (!calibList) { LOG(fatal) << "Can not find the TPC Post Calibration object!"; } - hMapPionMean = (TH3F*)calibList->FindObject("mean_map_pion"); - hMapPionSigma = (TH3F*)calibList->FindObject("sigma_map_pion"); - hMapProtonMean = (TH3F*)calibList->FindObject("mean_map_proton"); - hMapProtonSigma = (TH3F*)calibList->FindObject("sigma_map_proton"); + hMapPion[0] = (TH3F*)calibList->FindObject("mean_map_pion"); + hMapPion[1] = (TH3F*)calibList->FindObject("sigma_map_pion"); + hMapProton[0] = (TH3F*)calibList->FindObject("mean_map_proton"); + hMapProton[1] = (TH3F*)calibList->FindObject("sigma_map_proton"); - if (!hMapPionMean || !hMapPionSigma || !hMapProtonMean || !hMapProtonSigma) { + if (!hMapPion[0] || !hMapPion[1] || !hMapProton[0] || !hMapProton[1]) { LOG(fatal) << "Can not find histograms!"; } + } else if (setTPCCalib > 1) { + + hBBProton[0] = setValuesBB(ccdbApi, bc, ccdbBBProton); + hBBProton[1] = setValuesBB(ccdbApi, bc, ccdbBBAntiProton); + hBBPion[0] = setValuesBB(ccdbApi, bc, ccdbBBPion); + hBBPion[1] = setValuesBB(ccdbApi, bc, ccdbBBAntiPion); } currentRun = bc.runNumber(); @@ -311,7 +320,7 @@ struct HfFilter { // Main struct for HF triggers auto trackPos = cand2Prong.prong0_as(); // positive daughter auto trackNeg = cand2Prong.prong1_as(); // negative daughter - auto preselD0 = isDzeroPreselected(trackPos, trackNeg, nsigmaTPCPionKaonDzero, nsigmaTOFPionKaonDzero, computeTPCPostCalib, hMapPionMean, hMapPionSigma); + auto preselD0 = isDzeroPreselected(trackPos, trackNeg, nsigmaTPCPionKaonDzero, nsigmaTOFPionKaonDzero, setTPCCalib, hMapPion, hBBPion); if (!preselD0) { continue; } @@ -465,7 +474,7 @@ struct HfFilter { // Main struct for HF triggers // 2-prong femto if (!keepEvent[kFemto2P] && isCharmTagged && track.collisionId() == thisCollId) { - bool isProton = isSelectedProton4Femto(track, trackParThird, femtoMinProtonPt, femtoMaxNsigmaProton, femtoProtonOnlyTOF, computeTPCPostCalib, hMapProtonMean, hMapProtonSigma, activateQA, hProtonTPCPID, hProtonTOFPID); + bool isProton = isSelectedProton4Femto(track, trackParThird, femtoMinProtonPt, femtoMaxNsigmaProton, femtoProtonOnlyTOF, setTPCCalib, hMapProton, hBBProton, activateQA, hProtonTPCPID, hProtonTOFPID); if (isProton) { float relativeMomentum = computeRelativeMomentum(pVecThird, pVec2Prong, massD0); if (applyOptimisation) { @@ -543,13 +552,13 @@ struct HfFilter { // Main struct for HF triggers } if (is3Prong[0]) { // D+ preselections - is3Prong[0] = isDplusPreselected(trackSecond, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapPionMean, hMapPionSigma); + is3Prong[0] = isDplusPreselected(trackSecond, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapPion, hBBPion); } if (is3Prong[1]) { // Ds preselections - is3Prong[1] = isDsPreselected(pVecFirst, pVecThird, pVecSecond, trackSecond, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapPionMean, hMapPionSigma); + is3Prong[1] = isDsPreselected(pVecFirst, pVecThird, pVecSecond, trackSecond, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapPion, hBBPion); } if (is3Prong[2] || is3Prong[3]) { // charm baryon preselections - auto presel = isCharmBaryonPreselected(trackFirst, trackThird, trackSecond, nsigmaTPCProtonLc, nsigmaTOFProtonLc, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapProtonMean, hMapProtonSigma, hMapPionMean, hMapPionSigma); + auto presel = isCharmBaryonPreselected(trackFirst, trackThird, trackSecond, nsigmaTPCProtonLc, nsigmaTOFProtonLc, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapProton, hBBProton, hMapPion, hBBPion); if (is3Prong[2]) { is3Prong[2] = presel; } @@ -696,7 +705,7 @@ struct HfFilter { // Main struct for HF triggers } // end beauty selection // 3-prong femto - bool isProton = isSelectedProton4Femto(track, trackParFourth, femtoMinProtonPt, femtoMaxNsigmaProton, femtoProtonOnlyTOF, computeTPCPostCalib, hMapProtonMean, hMapProtonSigma, activateQA, hProtonTPCPID, hProtonTOFPID); + bool isProton = isSelectedProton4Femto(track, trackParFourth, femtoMinProtonPt, femtoMaxNsigmaProton, femtoProtonOnlyTOF, setTPCCalib, hMapProton, hBBProton, activateQA, hProtonTPCPID, hProtonTOFPID); if (isProton && track.collisionId() == thisCollId) { for (int iHypo{0}; iHypo < kNCharmParticles - 1 && !keepEvent[kFemto3P]; ++iHypo) { if (isCharmTagged[iHypo]) { diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index 871c26f6c14..5ef496f01fd 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -29,6 +29,7 @@ #include "DataFormatsParameters/GRPObject.h" #include "DetectorsBase/Propagator.h" #include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsTPC/BetheBlochAleph.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" @@ -141,6 +142,44 @@ static const AxisSpec bdtAxis{100, 0.f, 1.f}; static const std::array massAxisC = {AxisSpec{100, 1.65f, 2.05f}, AxisSpec{100, 1.65f, 2.05f}, AxisSpec{100, 1.75f, 2.15f}, AxisSpec{100, 2.05f, 2.45f}, AxisSpec{100, 2.25f, 2.65f}, AxisSpec{100, 1.98f, 2.08f}, AxisSpec{100, 1.98f, 2.08f}, AxisSpec{100, 2.08f, 2.18f}}; static const std::array massAxisB = {AxisSpec{100, 5.0f, 5.6f}, AxisSpec{100, 5.0f, 5.6f}, AxisSpec{100, 5.0f, 5.6f}, AxisSpec{100, 5.0f, 5.6f}, AxisSpec{100, 5.3f, 5.9f}, AxisSpec{100, 5.3f, 5.9f}}; +/// load the TPC spline from the CCDB +/// \param ccdbApi is Api for CCDB +/// \param bunchCrossing is the timestamp of bunchcrossing for the run number +/// \param ccdbPath is the path on CCDB +/// \return a vector include parameters for BetheBloch formula +std::vector setValuesBB(o2::ccdb::CcdbApi& ccdbApi, aod::BCsWithTimestamps::iterator const& bunchCrossing, const std::string ccdbPath) +{ + map metadata; + auto hSpline = ccdbApi.retrieveFromTFileAny(ccdbPath, metadata, bunchCrossing.timestamp()); + + if (!hSpline) { + LOG(fatal) << "File from CCDB in path " << ccdbPath << " was not found for run " << bunchCrossing.runNumber(); + } + + TAxis* axis = hSpline->GetXaxis(); + std::vector parsBB{static_cast(hSpline->GetBinContent(axis->FindBin("bb1"))), + static_cast(hSpline->GetBinContent(axis->FindBin("bb2"))), + static_cast(hSpline->GetBinContent(axis->FindBin("bb3"))), + static_cast(hSpline->GetBinContent(axis->FindBin("bb4"))), + static_cast(hSpline->GetBinContent(axis->FindBin("bb5"))), + static_cast(hSpline->GetBinContent(axis->FindBin("Resolution")))}; + return parsBB; +} + +/// Update the TPC PID baesd on the spline of particles +/// \param track is a track parameter +/// \param mMassPar is the mass of particles +/// \param paraBetheBloch vector for the parameters of BetheBloch formula +/// \return updated nsigma value for TPC PID +template +double updateTPCPID(const T& track, const float mMassPar, const std::vector paraBetheBloch) +{ + auto bgScaling = 1 / mMassPar; + double expBethe = tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScaling), paraBetheBloch[0], paraBetheBloch[1], paraBetheBloch[2], paraBetheBloch[3], paraBetheBloch[4]); + double expSigma = expBethe * paraBetheBloch[5]; + return static_cast((track.tpcSignal() - expBethe) / expSigma); +} + /// Single-track cuts for bachelor track of beauty candidates /// \param trackPar is a track parameter /// \param dca is the 2d array with dcaXY and dcaZ of the track @@ -191,14 +230,14 @@ int isSelectedTrackForBeauty(const T1& trackPar, const T2& dca, const float& pTM /// \param femtoMaxNsigmaProton max Nsigma for proton candidates /// \param femtoProtonOnlyTOF flag to activate PID selection with TOF only /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapProtonMean map of nSigma mean calibrations -/// \param hMapProtonSigma map of nSigma sigma calibrations +/// \param hMapProton map of nSigma mean and sigma calibrations for proton +/// \param hSplineProton spline of proton and anti-proton calibrations /// \param activateQA flag to activate the filling of QA histos /// \param hProtonTPCPID histo with NsigmaTPC vs. p /// \param hProtonTOFPID histo with NsigmaTOF vs. p /// \return true if track passes all cuts template -bool isSelectedProton4Femto(const T1& track, const T2& trackPar, const float& femtoMinProtonPt, const float& femtoMaxNsigmaProton, const bool femtoProtonOnlyTOF, const bool computeTPCPostCalib, H3 hMapProtonMean, H3 hMapProtonSigma, const int& activateQA, H2 hProtonTPCPID, H2 hProtonTOFPID) +bool isSelectedProton4Femto(const T1& track, const T2& trackPar, const float& femtoMinProtonPt, const float& femtoMaxNsigmaProton, const bool femtoProtonOnlyTOF, const int setTPCCalib, H3 hMapProton, const std::array, 2>& hSplineProton, const int& activateQA, H2 hProtonTPCPID, H2 hProtonTOFPID) { if (trackPar.getPt() < femtoMinProtonPt) { return false; @@ -217,8 +256,10 @@ bool isSelectedProton4Femto(const T1& track, const T2& trackPar, const float& fe float NSigmaTOF = track.tofNSigmaPr(); float NSigma; - if (computeTPCPostCalib) { - NSigmaTPC = getTPCPostCalib(hMapProtonMean, hMapProtonSigma, track, kPr); + if (setTPCCalib == 1) { + NSigmaTPC = getTPCPostCalib(hMapProton, track, kPr); + } else if (setTPCCalib == 2) { + NSigmaTPC = updateTPCPID(track, massProton, hSplineProton[0]); } if (femtoProtonOnlyTOF) { @@ -244,17 +285,19 @@ bool isSelectedProton4Femto(const T1& track, const T2& trackPar, const float& fe /// \param nsigmaTPCProtonLc max NsigmaTPC for proton candidates /// \param nsigmaTOFProtonLc max NsigmaTOF for proton candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapProtonMean map of nSigma mean calibrations -/// \param hMapProtonSigma map of nSigma sigma calibrations +/// \param hMapProton map of nSigma mean and sigma calibrations for proton +/// \param hSplineProton spline of proton and anti-proton calibrations /// \return true if track passes all cuts template -bool isSelectedProton4CharmBaryons(const T& track, const float& nsigmaTPCProtonLc, const float& nsigmaTOFProtonLc, const bool computeTPCPostCalib, H3 hMapProtonMean, H3 hMapProtonSigma) +bool isSelectedProton4CharmBaryons(const T& track, const float& nsigmaTPCProtonLc, const float& nsigmaTOFProtonLc, const int setTPCCalib, H3 hMapProton, const std::array, 2>& hSplineProton) { float NSigmaTPC = track.tpcNSigmaPr(); float NSigmaTOF = track.tofNSigmaPr(); - if (computeTPCPostCalib) { - NSigmaTPC = getTPCPostCalib(hMapProtonMean, hMapProtonSigma, track, kPr); + if (setTPCCalib == 1) { + NSigmaTPC = getTPCPostCalib(hMapProton, track, kPr); + } else if (setTPCCalib == 2) { + NSigmaTPC = updateTPCPID(track, massProton, hSplineProton[0]); } if (std::abs(NSigmaTPC) > nsigmaTPCProtonLc) { @@ -272,17 +315,19 @@ bool isSelectedProton4CharmBaryons(const T& track, const float& nsigmaTPCProtonL /// \param nsigmaTPCKaon3Prong max NsigmaTPC for kaon candidates /// \param nsigmaTOFKaon3Prong max NsigmaTOF for kaon candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapPionMean map of nSigma mean calibrations -/// \param hMapPionSigma map of nSigma sigma calibrations +/// \param hMapPion map of nSigma mean and sigma calibrations for pion +/// \param hSplinePion spline of pion and anti-pion calibrations /// \return true if track passes all cuts template -bool isSelectedKaon4Charm3Prong(const T& track, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const bool computeTPCPostCalib, H3 hMapPionMean, H3 hMapPionSigma) +bool isSelectedKaon4Charm3Prong(const T& track, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const int setTPCCalib, H3 hMapPion, const std::array, 2>& hSplinePion) { float NSigmaTPC = track.tpcNSigmaKa(); float NSigmaTOF = track.tofNSigmaKa(); - if (computeTPCPostCalib) { - NSigmaTPC = getTPCPostCalib(hMapPionMean, hMapPionSigma, track, kKa); // use pion correction map for kaon for the moment + if (setTPCCalib == 1) { + NSigmaTPC = getTPCPostCalib(hMapPion, track, kKa); // use pion correction map for kaon for the moment + } else if (setTPCCalib == 2) { + NSigmaTPC = updateTPCPID(track, massK, hSplinePion[0]); } if (std::abs(NSigmaTPC) > nsigmaTPCKaon3Prong) { @@ -300,16 +345,16 @@ bool isSelectedKaon4Charm3Prong(const T& track, const float& nsigmaTPCKaon3Prong /// \param nsigmaTPCKaon3Prong max NsigmaTPC for kaon candidates /// \param nsigmaTOFKaon3Prong max NsigmaTOF for kaon candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapPionMean map of nSigma mean calibrations -/// \param hMapPionSigma map of nSigma sigma calibrations +/// \param hMapPion map of nSigma mean and sigma calibrations for pion +/// \param hSplinePion spline of pion and anti-pion calibrations /// \return BIT(0) for Kpipi template -int8_t isDplusPreselected(const T& trackOppositeCharge, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const bool computeTPCPostCalib, H3 hMapPionMean, H3 hMapPionSigma) +int8_t isDplusPreselected(const T& trackOppositeCharge, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const int setTPCCalib, H3 hMapPion, std::array, 2>& hSplinePion) { int8_t retValue = 0; // check PID of opposite charge track - if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapPionMean, hMapPionSigma)) { + if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapPion, hSplinePion)) { return retValue; } @@ -325,16 +370,16 @@ int8_t isDplusPreselected(const T& trackOppositeCharge, const float& nsigmaTPCKa /// \param nsigmaTPCKaon3Prong max NsigmaTPC for kaon candidates /// \param nsigmaTOFKaon3Prong max NsigmaTOF for kaon candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapPionMean map of nSigma mean calibrations -/// \param hMapPionSigma map of nSigma sigma calibrations +/// \param hMapPion map of nSigma mean and sigma calibrations for pion +/// \param hSplinePion spline of pion and anti-pion calibrations /// \return BIT(0) for KKpi, BIT(1) for piKK template -int8_t isDsPreselected(const P& pTrackSameChargeFirst, const P& pTrackSameChargeSecond, const P& pTrackOppositeCharge, const T& trackOppositeCharge, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const bool computeTPCPostCalib, H3 hMapPionMean, H3 hMapPionSigma) +int8_t isDsPreselected(const P& pTrackSameChargeFirst, const P& pTrackSameChargeSecond, const P& pTrackOppositeCharge, const T& trackOppositeCharge, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const int setTPCCalib, H3 hMapPion, const std::array, 2>& hSplinePion) { int8_t retValue = 0; // check PID of opposite charge track - if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapPionMean, hMapPionSigma)) { + if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapPion, hSplinePion)) { return retValue; } @@ -361,24 +406,24 @@ int8_t isDsPreselected(const P& pTrackSameChargeFirst, const P& pTrackSameCharge /// \param nsigmaTPCKaon3Prong max NsigmaTPC for kaon candidates /// \param nsigmaTOFKaon3Prong max NsigmaTOF for kaon candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapProtonMean map of nSigma mean calibrations for protons -/// \param hMapProtonSigma map of nSigma sigma calibrations for protons -/// \param hMapPionMean map of nSigma mean calibrations for pions -/// \param hMapPionSigma map of nSigma sigma calibrations for pions +/// \param hMapProton map of nSigma mean and sigma calibrations for proton +/// \param hSplineProton spline of proton and anti-proton calibrations +/// \param hMapPion map of nSigma mean and sigma calibrations for pion +/// \param hSplinePion spline of pion and anti-pion calibrations /// \return BIT(0) for pKpi, BIT(1) for piKp template -int8_t isCharmBaryonPreselected(const T& trackSameChargeFirst, const T& trackSameChargeSecond, const T& trackOppositeCharge, const float& nsigmaTPCProtonLc, const float& nsigmaTOFProtonLc, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const bool computeTPCPostCalib, H3 hMapProtonMean, H3 hMapProtonSigma, H3 hMapPionMean, H3 hMapPionSigma) +int8_t isCharmBaryonPreselected(const T& trackSameChargeFirst, const T& trackSameChargeSecond, const T& trackOppositeCharge, const float& nsigmaTPCProtonLc, const float& nsigmaTOFProtonLc, const float& nsigmaTPCKaon3Prong, const float& nsigmaTOFKaon3Prong, const int setTPCCalib, H3 hMapProton, const std::array, 2>& hSplineProton, H3 hMapPion, const std::array, 2>& hSplinePion) { int8_t retValue = 0; // check PID of opposite charge track - if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, computeTPCPostCalib, hMapPionMean, hMapPionSigma)) { + if (!isSelectedKaon4Charm3Prong(trackOppositeCharge, nsigmaTPCKaon3Prong, nsigmaTOFKaon3Prong, setTPCCalib, hMapPion, hSplinePion)) { return retValue; } - if (isSelectedProton4CharmBaryons(trackSameChargeFirst, nsigmaTPCProtonLc, nsigmaTOFProtonLc, computeTPCPostCalib, hMapProtonMean, hMapProtonSigma)) { + if (isSelectedProton4CharmBaryons(trackSameChargeFirst, nsigmaTPCProtonLc, nsigmaTOFProtonLc, setTPCCalib, hMapProton, hSplineProton)) { retValue |= BIT(0); } - if (isSelectedProton4CharmBaryons(trackSameChargeSecond, nsigmaTPCProtonLc, nsigmaTOFProtonLc, computeTPCPostCalib, hMapProtonMean, hMapProtonSigma)) { + if (isSelectedProton4CharmBaryons(trackSameChargeSecond, nsigmaTPCProtonLc, nsigmaTOFProtonLc, setTPCCalib, hMapProton, hSplineProton)) { retValue |= BIT(1); } @@ -391,11 +436,11 @@ int8_t isCharmBaryonPreselected(const T& trackSameChargeFirst, const T& trackSam /// \param nsigmaTPCPionKaonDzero max NsigmaTPC for pion/kaon candidates /// \param nsigmaTOFPionKaonDzero max NsigmaTOF for pion/kaon candidates /// \param computeTPCPostCalib flag to activate TPC PID postcalibrations -/// \param hMapPionMean map of nSigma mean calibrations -/// \param hMapPionSigma map of nSigma sigma calibrations +/// \param hMapPion map of nSigma mean and sigma calibrations for pion +/// \param hSplinePion spline of pion and anti-pion calibrations /// \return BIT(0) for D0, BIT(1) for D0bar template -int8_t isDzeroPreselected(const T& trackPos, const T& trackNeg, const float& nsigmaTPCPionKaonDzero, const float& nsigmaTOFPionKaonDzero, const bool computeTPCPostCalib, H3 hMapPionMean, H3 hMapPionSigma) +int8_t isDzeroPreselected(const T& trackPos, const T& trackNeg, const float& nsigmaTPCPionKaonDzero, const float& nsigmaTOFPionKaonDzero, const int setTPCCalib, H3 hMapPion, const std::array, 2>& hSplinePion) { int8_t retValue = 0; @@ -409,11 +454,16 @@ int8_t isDzeroPreselected(const T& trackPos, const T& trackNeg, const float& nsi float NSigmaKaTPCNeg = trackNeg.tpcNSigmaKa(); float NSigmaKaTOFNeg = trackNeg.tofNSigmaKa(); - if (computeTPCPostCalib) { - NSigmaPiTPCPos = getTPCPostCalib(hMapPionMean, hMapPionSigma, trackPos, kPi); - NSigmaPiTPCNeg = getTPCPostCalib(hMapPionMean, hMapPionSigma, trackNeg, kPi); - NSigmaKaTPCPos = getTPCPostCalib(hMapPionMean, hMapPionSigma, trackPos, kKa); // use pion correction map for kaon for the moment - NSigmaKaTPCNeg = getTPCPostCalib(hMapPionMean, hMapPionSigma, trackNeg, kKa); // use pion correction map for kaon for the moment + if (setTPCCalib == 1) { + NSigmaPiTPCPos = getTPCPostCalib(hMapPion, trackPos, kPi); + NSigmaPiTPCNeg = getTPCPostCalib(hMapPion, trackNeg, kPi); + NSigmaKaTPCPos = getTPCPostCalib(hMapPion, trackPos, kKa); // use pion correction map for kaon for the moment + NSigmaKaTPCNeg = getTPCPostCalib(hMapPion, trackNeg, kKa); // use pion correction map for kaon for the moment + } else if (setTPCCalib == 2) { + NSigmaPiTPCPos = updateTPCPID(trackPos, massPi, hSplinePion[0]); + NSigmaPiTPCNeg = updateTPCPID(trackNeg, massPi, hSplinePion[1]); + NSigmaKaTPCPos = updateTPCPID(trackPos, massK, hSplinePion[0]); // use pion correction for kaon for the moment + NSigmaKaTPCNeg = updateTPCPID(trackNeg, massK, hSplinePion[1]); // use pion correction for kaon for the moment } if ((std::abs(NSigmaPiTPCPos) <= nsigmaTPCPionKaonDzero && (!trackPos.hasTOF() || std::abs(NSigmaPiTOFPos) <= nsigmaTOFPionKaonDzero)) && (std::abs(NSigmaKaTPCNeg) <= nsigmaTPCPionKaonDzero && (!trackNeg.hasTOF() || std::abs(NSigmaKaTOFNeg) <= nsigmaTOFPionKaonDzero))) { @@ -819,8 +869,8 @@ std::array PredictONNX(std::vector& inputFeatures, std::shared_ptr -float getTPCPostCalib(const TH3F* hCalibMean, const TH3F* hCalibSigma, const T& track, const int pidSpecies) +template +float getTPCPostCalib(const array& hCalibMap, const T& track, const int pidSpecies) { auto tpcNCls = track.tpcNClsFound(); auto tpcPin = track.tpcInnerParam(); @@ -836,18 +886,18 @@ float getTPCPostCalib(const TH3F* hCalibMean, const TH3F* hCalibSigma, const T& } else { LOG(fatal) << "Wrong PID Species be selected, please check!"; } - auto binTPCNCls = hCalibMean->GetXaxis()->FindBin(tpcNCls); + auto binTPCNCls = hCalibMap[0]->GetXaxis()->FindBin(tpcNCls); binTPCNCls = (binTPCNCls == 0 ? 1 : binTPCNCls); - binTPCNCls = std::min(hCalibMean->GetXaxis()->GetNbins(), binTPCNCls); - auto binPin = hCalibMean->GetYaxis()->FindBin(tpcPin); + binTPCNCls = std::min(hCalibMap[0]->GetXaxis()->GetNbins(), binTPCNCls); + auto binPin = hCalibMap[0]->GetYaxis()->FindBin(tpcPin); binPin = (binPin == 0 ? 1 : binPin); - binPin = std::min(hCalibMean->GetYaxis()->GetNbins(), binPin); - auto binEta = hCalibMean->GetZaxis()->FindBin(eta); + binPin = std::min(hCalibMap[0]->GetYaxis()->GetNbins(), binPin); + auto binEta = hCalibMap[0]->GetZaxis()->FindBin(eta); binEta = (binEta == 0 ? 1 : binEta); - binEta = std::min(hCalibMean->GetZaxis()->GetNbins(), binEta); + binEta = std::min(hCalibMap[0]->GetZaxis()->GetNbins(), binEta); - auto mean = hCalibMean->GetBinContent(binTPCNCls, binPin, binEta); - auto width = hCalibSigma->GetBinContent(binTPCNCls, binPin, binEta); + auto mean = hCalibMap[0]->GetBinContent(binTPCNCls, binPin, binEta); + auto width = hCalibMap[1]->GetBinContent(binTPCNCls, binPin, binEta); return (tpcNSigma - mean) / width; } From 873ec8bf2697aa12d2baae9525dff016ca4596b0 Mon Sep 17 00:00:00 2001 From: Biao Zhang Date: Wed, 15 Feb 2023 01:33:11 +0800 Subject: [PATCH 2/2] Change to proper name --- EventFiltering/PWGHF/HFFilterHelpers.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index 5ef496f01fd..637941b77b1 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -172,7 +172,7 @@ std::vector setValuesBB(o2::ccdb::CcdbApi& ccdbApi, aod::BCsWithTimestam /// \param paraBetheBloch vector for the parameters of BetheBloch formula /// \return updated nsigma value for TPC PID template -double updateTPCPID(const T& track, const float mMassPar, const std::vector paraBetheBloch) +double getTPCSplineCalib(const T& track, const float mMassPar, const std::vector paraBetheBloch) { auto bgScaling = 1 / mMassPar; double expBethe = tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScaling), paraBetheBloch[0], paraBetheBloch[1], paraBetheBloch[2], paraBetheBloch[3], paraBetheBloch[4]); @@ -259,7 +259,7 @@ bool isSelectedProton4Femto(const T1& track, const T2& trackPar, const float& fe if (setTPCCalib == 1) { NSigmaTPC = getTPCPostCalib(hMapProton, track, kPr); } else if (setTPCCalib == 2) { - NSigmaTPC = updateTPCPID(track, massProton, hSplineProton[0]); + NSigmaTPC = getTPCSplineCalib(track, massProton, hSplineProton[0]); } if (femtoProtonOnlyTOF) { @@ -297,7 +297,7 @@ bool isSelectedProton4CharmBaryons(const T& track, const float& nsigmaTPCProtonL if (setTPCCalib == 1) { NSigmaTPC = getTPCPostCalib(hMapProton, track, kPr); } else if (setTPCCalib == 2) { - NSigmaTPC = updateTPCPID(track, massProton, hSplineProton[0]); + NSigmaTPC = getTPCSplineCalib(track, massProton, hSplineProton[0]); } if (std::abs(NSigmaTPC) > nsigmaTPCProtonLc) { @@ -327,7 +327,7 @@ bool isSelectedKaon4Charm3Prong(const T& track, const float& nsigmaTPCKaon3Prong if (setTPCCalib == 1) { NSigmaTPC = getTPCPostCalib(hMapPion, track, kKa); // use pion correction map for kaon for the moment } else if (setTPCCalib == 2) { - NSigmaTPC = updateTPCPID(track, massK, hSplinePion[0]); + NSigmaTPC = getTPCSplineCalib(track, massK, hSplinePion[0]); } if (std::abs(NSigmaTPC) > nsigmaTPCKaon3Prong) { @@ -460,10 +460,10 @@ int8_t isDzeroPreselected(const T& trackPos, const T& trackNeg, const float& nsi NSigmaKaTPCPos = getTPCPostCalib(hMapPion, trackPos, kKa); // use pion correction map for kaon for the moment NSigmaKaTPCNeg = getTPCPostCalib(hMapPion, trackNeg, kKa); // use pion correction map for kaon for the moment } else if (setTPCCalib == 2) { - NSigmaPiTPCPos = updateTPCPID(trackPos, massPi, hSplinePion[0]); - NSigmaPiTPCNeg = updateTPCPID(trackNeg, massPi, hSplinePion[1]); - NSigmaKaTPCPos = updateTPCPID(trackPos, massK, hSplinePion[0]); // use pion correction for kaon for the moment - NSigmaKaTPCNeg = updateTPCPID(trackNeg, massK, hSplinePion[1]); // use pion correction for kaon for the moment + NSigmaPiTPCPos = getTPCSplineCalib(trackPos, massPi, hSplinePion[0]); + NSigmaPiTPCNeg = getTPCSplineCalib(trackNeg, massPi, hSplinePion[1]); + NSigmaKaTPCPos = getTPCSplineCalib(trackPos, massK, hSplinePion[0]); // use pion correction for kaon for the moment + NSigmaKaTPCNeg = getTPCSplineCalib(trackNeg, massK, hSplinePion[1]); // use pion correction for kaon for the moment } if ((std::abs(NSigmaPiTPCPos) <= nsigmaTPCPionKaonDzero && (!trackPos.hasTOF() || std::abs(NSigmaPiTOFPos) <= nsigmaTOFPionKaonDzero)) && (std::abs(NSigmaKaTPCNeg) <= nsigmaTPCPionKaonDzero && (!trackNeg.hasTOF() || std::abs(NSigmaKaTOFNeg) <= nsigmaTOFPionKaonDzero))) {