Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
12 changes: 6 additions & 6 deletions PWGCF/Femto/Core/femtoUtils.h
Original file line number Diff line number Diff line change
Expand Up @@ -167,14 +167,14 @@ concept HasEventShapeRow = requires(T row) {
template <typename T>
concept HasEventShape = HasEventShapeRow<std::decay_t<T>> || (requires { typename std::decay_t<T>::iterator; } && HasEventShapeRow<typename std::decay_t<T>::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);
Expand Down
96 changes: 60 additions & 36 deletions PWGCF/Femto/Core/kinkBuilder.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "Common/Core/RecoDecay.h"

#include <CommonConstants/MathConstants.h>
#include <CommonConstants/PhysicsConstants.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/Configurable.h>
#include <Framework/HistogramRegistry.h>
Expand Down Expand Up @@ -60,34 +61,39 @@ struct ConfKinkFilters : o2::framework::ConfigurableGroup {

// selections bits for all kinks
// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
#define KINK_DEFAULT_BITS \
o2::framework::Configurable<bool> passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."}; \
o2::framework::Configurable<std::vector<float>> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"}; \
o2::framework::Configurable<std::vector<float>> transRadMin{"transRadMin", {20.f}, "Minimum transverse radius (cm)"}; \
o2::framework::Configurable<std::vector<float>> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"}; \
o2::framework::Configurable<std::vector<float>> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"}; \
o2::framework::Configurable<std::vector<float>> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"}; \
o2::framework::Configurable<std::vector<float>> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"}; \
o2::framework::Configurable<std::vector<float>> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"};
#define KINK_DEFAULT_BITS \
o2::framework::Configurable<bool> passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."}; \
o2::framework::Configurable<std::vector<float>> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"}; \
o2::framework::Configurable<std::vector<float>> transRadMin{"transRadMin", {19.6f}, "Minimum transverse radius (cm)"}; \
o2::framework::Configurable<std::vector<float>> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"}; \
o2::framework::Configurable<std::vector<float>> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"}; \
o2::framework::Configurable<std::vector<float>> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"}; \
o2::framework::Configurable<std::vector<float>> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"}; \
o2::framework::Configurable<std::vector<float>> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"}; \
o2::framework::Configurable<std::vector<float>> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"}; \
o2::framework::Configurable<std::vector<float>> ptOriginalMin{"ptOriginalMin", {1.2f}, "Minimum original (not recalculated) pT of the mother (GeV/c)"}; \
o2::framework::Configurable<std::vector<float>> 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<std::vector<float>> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"};
o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"};
o2::framework::Configurable<std::vector<float>> chaDauTpcPion{"chaDauTpcPion", {5.f}, "Maximum |nsigma_Pion| TPC for charged daughter tracks"};
};

// derived selection bits for sigma plus
struct ConfSigmaPlusBits : o2::framework::ConfigurableGroup {
std::string prefix = std::string("SigmaPlusBits");
KINK_DEFAULT_BITS
o2::framework::Configurable<std::vector<float>> alphaAPMin{"alphaAPMin", {0.0f}, "Minimum Alpha_AP for SigmaPlus candidates"};
o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {1.0f}, "Maximum Alpha_AP for SigmaPlus candidates"};
o2::framework::Configurable<std::vector<float>> chaDauTpcProton{"chaDauTpcProton", {5.f}, "Maximum |nsigma_Proton| TPC for charged daughter tracks"};
o2::framework::Configurable<std::vector<float>> chaDauTofProton{"chaDauTofProton", {5.f}, "Maximum combined |nsigma_Proton| (TPC+TOF) for charged daughter tracks"};
o2::framework::Configurable<float> pidThres{"pidThres", 0.75f, "Momentum threshold for using TOF/combined pid for daughter tracks (GeV/c)"};
o2::framework::Configurable<std::vector<float>> chaDauTofProton{"chaDauTofProton", {}, "Maximum |nsigma_Proton| TOF for charged daughter tracks"};
o2::framework::Configurable<bool> requireTof{"requireTof", false, "If true, TOF PID is a minimal selection. If false, TOF PID is optional"};
o2::framework::Configurable<bool> 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
Expand Down Expand Up @@ -153,6 +159,9 @@ enum KinkSeles {
kQtAPMax,
kCosPointingAngleMin,

kPtOriginalMin,
kPtOriginalMax,

kKinkSelsMax
};

Expand All @@ -173,7 +182,9 @@ const std::unordered_map<KinkSeles, std::string> 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 {
Expand Down Expand Up @@ -231,9 +242,9 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) {
mMassSigmaPlusLowerLimit = filter.massMinSigmaPlus.value;
mMassSigmaPlusUpperLimit = filter.massMaxSigmaPlus.value;
mPidThreshold = config.pidThres.value;
this->addSelection(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);
Expand All @@ -247,6 +258,8 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
this->addSelection(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<SelectionHistName>(registry);

Expand Down Expand Up @@ -281,7 +294,7 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
};

template <typename T1, typename T2, typename T3>
void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& col)
void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& /*col*/)
{
std::array<float, 3> momMother = {kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()};
float kinkMomP = RecoDecay::p(momMother);
Expand All @@ -300,12 +313,11 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
float p2A = kinkDauP * kinkDauP;
mQtAp = (p2V0 > 0.f) ? std::sqrt(std::max(0.f, p2A - dp * dp / p2V0)) : 0.f;

std::array<float, 3> vMother = {kinkCand.xDecVtx() - col.posX(), kinkCand.yDecVtx() - col.posY(), kinkCand.zDecVtx() - col.posZ()};
std::array<float, 3> 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;
Expand All @@ -327,6 +339,9 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
this->evaluateObservable(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
Expand All @@ -345,14 +360,11 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
this->evaluateObservable(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);
}
}

Expand All @@ -368,9 +380,20 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
mKinkMotherEta = RecoDecay::eta(momMother);
mKinkMotherPhi = RecoDecay::phi(momMother);

// Recalculate pT using kinematic constraints
float ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2]);
mKinkMotherPt = (ptRecalc > 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 <typename T>
Expand Down Expand Up @@ -443,7 +466,7 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
float mMassSigmaUpperLimit = 1.25f;
float mMassSigmaPlusLowerLimit = 1.15f;
float mMassSigmaPlusUpperLimit = 1.25f;
float mPidThreshold = 0.75f;
bool mKeepTracksWithoutTof = true;

// kinematic filters
float mPtMin = 0.f;
Expand All @@ -454,7 +477,8 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
float mPhiMax = o2::constants::math::TwoPI;

// mother kinematic
float mKinkMotherPt = 0.f;
float mKinkMotherPt = 0.f; // recalculated pT (falls back to the original pT if the recalculation fails)
float mKinkMotherPtOriginal = 0.f; // pT of the mother track as reconstructed
float mKinkMotherEta = 0.f;
float mKinkMotherPhi = 0.f;

Expand All @@ -464,7 +488,6 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
float mCosPointingAngle = 0.f;
float mTransRadius = 0.f;
float mKinkDauEta = 0.f;
float mKinkDauP = 0.f;
float mKinkAngle = 0.f;
};

Expand Down Expand Up @@ -549,7 +572,8 @@ class KinkBuilder
}
}

if (mProduceSigmas || mProduceSigmaMasks || mProduceSigmaExtras || mProduceSigmaPlus || mProduceSigmaPlusMasks || mProduceSigmaPlusExtras) {
if (mProduceSigmas || mProduceLiteSigmas || mProduceSigmaMasks || mProduceSigmaExtras ||
mProduceSigmaPlus || mProduceLiteSigmaPlus || mProduceSigmaPlusMasks || mProduceSigmaPlusExtras) {
mFillAnyTable = true;
} else {
LOG(info) << "No tables configured, Selection object will not be configured...";
Expand Down
16 changes: 10 additions & 6 deletions PWGCF/Femto/Core/kinkHistManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -284,25 +284,27 @@ class KinkHistManager
std::map<trackhistmanager::TrackHist, std::vector<o2::framework::AxisSpec>> 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;
}
Expand Down Expand Up @@ -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;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "PWGCF/Femto/Core/femtoUtils.h"
#include "PWGCF/Femto/DataModel/FemtoTables.h"

#include <CommonConstants/PhysicsConstants.h>
#include <Framework/AnalysisDataModel.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/AnalysisTask.h>
Expand Down Expand Up @@ -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();
Expand Down
Loading