diff --git a/PWGCF/Femto/Core/femtoUtils.h b/PWGCF/Femto/Core/femtoUtils.h index 643ddffe654..d7a4a888e74 100644 --- a/PWGCF/Femto/Core/femtoUtils.h +++ b/PWGCF/Femto/Core/femtoUtils.h @@ -167,14 +167,14 @@ concept HasEventShapeRow = requires(T row) { template concept HasEventShape = HasEventShapeRow> || (requires { typename std::decay_t::iterator; } && HasEventShapeRow::iterator>); -/// Recalculate pT for Kinks (Sigmas) using kinematic constraints -inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter) +/// Recalculate the pT of a kink mother from its direction and the charged daughter momentum +inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter, + float massMother, float massChargedDaughter, float massNeutralDaughter) { float almost0 = 1e-6f; - // Particle masses in GeV/c^2 - auto massPion = o2::constants::physics::MassPionCharged; - auto massNeutron = o2::constants::physics::MassNeutron; - auto massSigmaMinus = o2::constants::physics::MassSigmaMinus; + const float massPion = massChargedDaughter; + const float massNeutron = massNeutralDaughter; + const float massSigmaMinus = massMother; // Calculate mother momentum and direction versor float pMother = std::sqrt(pxMother * pxMother + pyMother * pyMother + pzMother * pzMother); diff --git a/PWGCF/Femto/Core/kinkBuilder.h b/PWGCF/Femto/Core/kinkBuilder.h index 009be2f6dd4..7dc4fb1b560 100644 --- a/PWGCF/Femto/Core/kinkBuilder.h +++ b/PWGCF/Femto/Core/kinkBuilder.h @@ -27,6 +27,7 @@ #include "Common/Core/RecoDecay.h" #include +#include #include #include #include @@ -60,24 +61,26 @@ struct ConfKinkFilters : o2::framework::ConfigurableGroup { // selections bits for all kinks // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define KINK_DEFAULT_BITS \ - o2::framework::Configurable passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."}; \ - o2::framework::Configurable> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"}; \ - o2::framework::Configurable> transRadMin{"transRadMin", {20.f}, "Minimum transverse radius (cm)"}; \ - o2::framework::Configurable> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"}; \ - o2::framework::Configurable> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"}; \ - o2::framework::Configurable> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"}; \ - o2::framework::Configurable> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"}; \ - o2::framework::Configurable> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"}; \ - o2::framework::Configurable> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"}; \ - o2::framework::Configurable> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"}; \ - o2::framework::Configurable> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"}; \ - o2::framework::Configurable> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"}; +#define KINK_DEFAULT_BITS \ + o2::framework::Configurable passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."}; \ + o2::framework::Configurable> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"}; \ + o2::framework::Configurable> transRadMin{"transRadMin", {19.6f}, "Minimum transverse radius (cm)"}; \ + o2::framework::Configurable> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"}; \ + o2::framework::Configurable> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"}; \ + o2::framework::Configurable> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"}; \ + o2::framework::Configurable> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"}; \ + o2::framework::Configurable> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"}; \ + o2::framework::Configurable> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"}; \ + o2::framework::Configurable> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"}; \ + o2::framework::Configurable> ptOriginalMin{"ptOriginalMin", {1.2f}, "Minimum original (not recalculated) pT of the mother (GeV/c)"}; \ + o2::framework::Configurable> ptOriginalMax{"ptOriginalMax", {10.f}, "Maximum original (not recalculated) pT of the mother (GeV/c)"}; // derived selection bits for sigma struct ConfSigmaBits : o2::framework::ConfigurableGroup { std::string prefix = std::string("SigmaBits"); KINK_DEFAULT_BITS + o2::framework::Configurable> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"}; + o2::framework::Configurable> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"}; o2::framework::Configurable> chaDauTpcPion{"chaDauTpcPion", {5.f}, "Maximum |nsigma_Pion| TPC for charged daughter tracks"}; }; @@ -85,9 +88,12 @@ struct ConfSigmaBits : o2::framework::ConfigurableGroup { struct ConfSigmaPlusBits : o2::framework::ConfigurableGroup { std::string prefix = std::string("SigmaPlusBits"); KINK_DEFAULT_BITS + o2::framework::Configurable> alphaAPMin{"alphaAPMin", {0.0f}, "Minimum Alpha_AP for SigmaPlus candidates"}; + o2::framework::Configurable> alphaAPMax{"alphaAPMax", {1.0f}, "Maximum Alpha_AP for SigmaPlus candidates"}; o2::framework::Configurable> chaDauTpcProton{"chaDauTpcProton", {5.f}, "Maximum |nsigma_Proton| TPC for charged daughter tracks"}; - o2::framework::Configurable> chaDauTofProton{"chaDauTofProton", {5.f}, "Maximum combined |nsigma_Proton| (TPC+TOF) for charged daughter tracks"}; - o2::framework::Configurable pidThres{"pidThres", 0.75f, "Momentum threshold for using TOF/combined pid for daughter tracks (GeV/c)"}; + o2::framework::Configurable> chaDauTofProton{"chaDauTofProton", {}, "Maximum |nsigma_Proton| TOF for charged daughter tracks"}; + o2::framework::Configurable requireTof{"requireTof", false, "If true, TOF PID is a minimal selection. If false, TOF PID is optional"}; + o2::framework::Configurable keepTracksWithoutTof{"keepTracksWithoutTof", true, "If true, the bit mask for the TOF selection will be true for all limits if the daughter track has no TOF"}; }; #undef KINK_DEFAULT_BITS @@ -153,6 +159,9 @@ enum KinkSeles { kQtAPMax, kCosPointingAngleMin, + kPtOriginalMin, + kPtOriginalMax, + kKinkSelsMax }; @@ -173,7 +182,9 @@ const std::unordered_map kinkSelectionNames = { {kAlphaAPMax, "alphaAPMax"}, {kQtAPMin, "qtAPMin"}, {kQtAPMax, "qtAPMax"}, - {kCosPointingAngleMin, "cosPointingAngleMin"}}; + {kCosPointingAngleMin, "cosPointingAngleMin"}, + {kPtOriginalMin, "ptOriginalMin"}, + {kPtOriginalMax, "ptOriginalMax"}}; /// enum for all kink pre-filters (evaluated in checkFilters, before the selection bitmask) enum KinkFilters { @@ -231,9 +242,9 @@ class KinkSelection : public baseselection::BaseSelectionaddSelection(kChaDaughTpcProton, kinkSelectionNames.at(kChaDaughTpcProton), config.chaDauTpcProton.value, limits::kAbsUpperLimit, false, false, true); - this->addSelection(kChaDaughTofProton, kinkSelectionNames.at(kChaDaughTofProton), config.chaDauTofProton.value, limits::kUpperLimit, false, false, true); + mKeepTracksWithoutTof = config.keepTracksWithoutTof.value; + this->addSelection(kChaDaughTpcProton, kinkSelectionNames.at(kChaDaughTpcProton), config.chaDauTpcProton.value, limits::kAbsUpperLimit, true, true, false); + this->addSelection(kChaDaughTofProton, kinkSelectionNames.at(kChaDaughTofProton), config.chaDauTofProton.value, limits::kAbsUpperLimit, true, config.requireTof.value, false); } this->addSelection(kKinkTopoDcaMax, kinkSelectionNames.at(kKinkTopoDcaMax), config.kinkTopoDcaMax.value, limits::kUpperLimit, true, true, false); @@ -247,6 +258,8 @@ class KinkSelection : public baseselection::BaseSelectionaddSelection(kQtAPMin, kinkSelectionNames.at(kQtAPMin), config.qtAPMin.value, limits::kLowerLimit, true, true, false); this->addSelection(kQtAPMax, kinkSelectionNames.at(kQtAPMax), config.qtAPMax.value, limits::kUpperLimit, true, true, false); this->addSelection(kCosPointingAngleMin, kinkSelectionNames.at(kCosPointingAngleMin), config.cosPointingAngleMin.value, limits::kLowerLimit, true, true, false); + this->addSelection(kPtOriginalMin, kinkSelectionNames.at(kPtOriginalMin), config.ptOriginalMin.value, limits::kLowerLimit, true, true, false); + this->addSelection(kPtOriginalMax, kinkSelectionNames.at(kPtOriginalMax), config.ptOriginalMax.value, limits::kUpperLimit, true, true, false); this->setupSelectionHistogram(registry); @@ -281,7 +294,7 @@ class KinkSelection : public baseselection::BaseSelection - void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& col) + void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& /*col*/) { std::array momMother = {kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}; float kinkMomP = RecoDecay::p(momMother); @@ -300,12 +313,11 @@ class KinkSelection : public baseselection::BaseSelection 0.f) ? std::sqrt(std::max(0.f, p2A - dp * dp / p2V0)) : 0.f; - std::array vMother = {kinkCand.xDecVtx() - col.posX(), kinkCand.yDecVtx() - col.posY(), kinkCand.zDecVtx() - col.posZ()}; + std::array vMother = {kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}; float vMotherNorm = std::sqrt(std::inner_product(vMother.begin(), vMother.end(), vMother.begin(), 0.f)); mCosPointingAngle = (vMotherNorm > 0.f && kinkMomP > 0.f) ? (std::inner_product(momMother.begin(), momMother.end(), vMother.begin(), 0.f)) / (kinkMomP * vMotherNorm) : 0.f; mTransRadius = std::hypot(kinkCand.xDecVtx(), kinkCand.yDecVtx()); - mKinkDauP = kinkDauP; mKinkDauEta = RecoDecay::eta(momDaughter); mKinkAngle = 0.f; @@ -327,6 +339,9 @@ class KinkSelection : public baseselection::BaseSelectionevaluateObservable(kQtAPMin, mQtAp); this->evaluateObservable(kQtAPMax, mQtAp); this->evaluateObservable(kCosPointingAngleMin, mCosPointingAngle); + // the stored pT is the recalculated one, the original pT is only available as selection bits + this->evaluateObservable(kPtOriginalMin, mKinkMotherPtOriginal); + this->evaluateObservable(kPtOriginalMax, mKinkMotherPtOriginal); this->evaluateObservable(kKinkTopoDcaMax, kinkCand.dcaKinkTopo()); // Compute transRadius @@ -345,14 +360,11 @@ class KinkSelection : public baseselection::BaseSelectionevaluateObservable(kChaDaughTpcPion, chaDaughter.tpcNSigmaPi()); } if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) { - if (mKinkDauP < mPidThreshold) { - this->evaluateObservable(kChaDaughTpcProton, chaDaughter.tpcNSigmaPr()); + this->evaluateObservable(kChaDaughTpcProton, chaDaughter.tpcNSigmaPr()); + if (chaDaughter.hasTOF()) { + this->evaluateObservable(kChaDaughTofProton, chaDaughter.tofNSigmaPr()); } else { - if (chaDaughter.hasTOF()) { - this->evaluateObservable(kChaDaughTofProton, std::abs(chaDaughter.tofNSigmaPr())); - } else { - this->evaluateObservable(kChaDaughTofProton, 999.f); - } + this->evaluateObservable(kChaDaughTofProton, mKeepTracksWithoutTof ? 0.f : 999.f); } } @@ -368,9 +380,20 @@ class KinkSelection : public baseselection::BaseSelection 0.f) ? ptRecalc : std::hypot(momMother[0], momMother[1]); + // Recalculate pT using kinematic constraints of the decay channel + float ptRecalc = -999.f; + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { + // Sigma- -> pi- n + ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2], + o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron); + } + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) { + // Sigma+ -> p pi0 + ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2], + o2::constants::physics::MassSigmaPlus, o2::constants::physics::MassProton, o2::constants::physics::MassPionNeutral); + } + mKinkMotherPtOriginal = std::hypot(momMother[0], momMother[1]); + mKinkMotherPt = (ptRecalc > 0.f) ? ptRecalc : mKinkMotherPtOriginal; } template @@ -443,7 +466,7 @@ class KinkSelection : public baseselection::BaseSelection> const& ChaDauSpecs) { mHistogramRegistry = registry; - mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value; + mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value); int chaDauPdgCodeAbs = 0; int chaDauCharge = 0; const int absCharge = 1; - if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) { + if (mPdgCode == PDG_t::kSigmaMinus) { if (ConfKinkSelection.sign.value < 0) { chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); chaDauCharge = -1; } else { + mPdgCode = -1 * mPdgCode; // anti-Sigma- is positively charged and has negative pdg code chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); chaDauCharge = 1; } - } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) { + } else if (mPdgCode == PDG_t::kSigmaPlus) { if (ConfKinkSelection.sign.value > 0) { chaDauPdgCodeAbs = std::abs(PDG_t::kProton); chaDauCharge = 1; } else { + mPdgCode = -1 * mPdgCode; // anti-Sigma+ is negatively charged and has negative pdg code chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); chaDauCharge = -1; } @@ -332,26 +334,28 @@ class KinkHistManager T3 const& ConfChaDauBinningQa) { mHistogramRegistry = registry; - mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value; + mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value); this->enableOptionalHistograms(ConfKinkBinningQa); int chaDauPdgCodeAbs = 0; int chaDauCharge = 0; const int absCharge = 1; - if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) { + if (mPdgCode == PDG_t::kSigmaMinus) { if (ConfKinkSelection.sign.value < 0) { chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); chaDauCharge = -1; } else { + mPdgCode = -1 * mPdgCode; // anti-Sigma- is positively charged and has negative pdg code chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); chaDauCharge = 1; } - } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) { + } else if (mPdgCode == PDG_t::kSigmaPlus) { if (ConfKinkSelection.sign.value > 0) { chaDauPdgCodeAbs = std::abs(PDG_t::kProton); chaDauCharge = 1; } else { + mPdgCode = -1 * mPdgCode; // anti-Sigma+ is negatively charged and has negative pdg code chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); chaDauCharge = -1; } diff --git a/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx b/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx index 8f61ef26d16..b97853e6829 100644 --- a/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx @@ -16,6 +16,7 @@ #include "PWGCF/Femto/Core/femtoUtils.h" #include "PWGCF/Femto/DataModel/FemtoTables.h" +#include #include #include #include @@ -73,7 +74,9 @@ struct FemtoProducerKinkPtConverter { float pyMoth = sigma.pt() * std::sin(sigma.phi()); float pzMoth = sigma.pt() * std::sinh(sigma.eta()); - float ptRecalc = utils::calcPtnew(pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug); + // Sigma- -> pi- n + float ptRecalc = utils::calcPtnew(pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug, + o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron); ROOT::Math::PtEtaPhiMVector recalcVec(ptRecalc, sigma.eta(), sigma.phi(), sigma.mass()); float ptFrom4Vec = recalcVec.Pt();