diff --git a/PWGCF/Femto/Core/dataTypes.h b/PWGCF/Femto/Core/dataTypes.h index 2fa1273dbcc..f8e349a77d3 100644 --- a/PWGCF/Femto/Core/dataTypes.h +++ b/PWGCF/Femto/Core/dataTypes.h @@ -35,6 +35,10 @@ using TrackType = uint16_t; using V0MaskType = uint16_t; using V0Type = uint8_t; +// datatypes for kinks +using KinkMaskType = uint32_t; +using KinkType = uint8_t; + // datatypes for two track resonances using TwoTrackResonanceMaskType = uint32_t; // two track resonance types diff --git a/PWGCF/Femto/Core/kinkBuilder.h b/PWGCF/Femto/Core/kinkBuilder.h new file mode 100644 index 00000000000..33a843dee0a --- /dev/null +++ b/PWGCF/Femto/Core/kinkBuilder.h @@ -0,0 +1,421 @@ +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file kinkBuilder.h +/// \brief kink builder +/// \author Anton Riedel, TU München, anton.riedel@cern.ch +/// \author Henrik Fribert, TU München, henrik.fribert@cern.ch + +#ifndef PWGCF_FEMTO_CORE_KINKBUILDER_H_ +#define PWGCF_FEMTO_CORE_KINKBUILDER_H_ + +#include "PWGCF/Femto/Core/baseSelection.h" +#include "PWGCF/Femto/Core/dataTypes.h" +#include "PWGCF/Femto/Core/femtoUtils.h" +#include "PWGCF/Femto/Core/modes.h" +#include "PWGCF/Femto/Core/selectionContainer.h" +#include "PWGCF/Femto/DataModel/FemtoTables.h" + +#include "Common/Core/RecoDecay.h" + +#include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisHelpers.h" +#include "Framework/Configurable.h" + +#include "fairlogger/Logger.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::femto +{ +namespace kinkbuilder +{ + +// filters applied in the producer task +struct ConfKinkFilters : o2::framework::ConfigurableGroup { + std::string prefix = std::string("KinkFilters"); + o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; + o2::framework::Configurable ptMax{"ptMax", 99.f, "Maximum pT"}; + o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; + o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; + o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum phi"}; + o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; + o2::framework::Configurable massMinSigma{"massMinSigma", 1.1f, "Minimum mass for Sigma hypothesis"}; + o2::framework::Configurable massMaxSigma{"massMaxSigma", 1.3f, "Maximum mass for Sigma hypothesis"}; +}; + +// selections bits for all kinks +#define KINK_DEFAULT_BITS \ + o2::framework::Configurable> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"}; \ + o2::framework::Configurable> transRadMin{"transRadMin", {0.2f}, "Minimum transverse radius (cm)"}; \ + o2::framework::Configurable> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"}; \ + o2::framework::Configurable> dauAbsEtaMax{"dauAbsEtaMax", {0.8f}, "Maximum absolute pseudorapidity for daughter track"}; \ + o2::framework::Configurable> dauDcaPvMin{"dauDcaPvMin", {0.0f}, "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> 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"}; + +// derived selection bits for sigma +struct ConfSigmaBits : o2::framework::ConfigurableGroup { + std::string prefix = std::string("SigmaBits"); + KINK_DEFAULT_BITS + o2::framework::Configurable> chaDauTpcPion{"chaDauTpcPion", {5.f}, "Maximum |nsigma_Pion| TPC for charged daughter tracks"}; +}; + +#undef KINK_DEFAULT_BITS + +// base selection for analysis task for kinks +#define KINK_DEFAULT_SELECTIONS(defaultMassMin, defaultMassMax, defaultPdgCode) \ + o2::framework::Configurable pdgCode{"pdgCode", defaultPdgCode, "Kink PDG code"}; \ + o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ + o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ + o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ + o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ + o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum phi"}; \ + o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ + o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Sigma"}; \ + o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Sigma"}; \ + o2::framework::Configurable mask{"mask", 0, "Bitmask for kink selection"}; + +// base selection for analysis task for sigmas +template +struct ConfSigmaSelection : o2::framework::ConfigurableGroup { + std::string prefix = Prefix; + KINK_DEFAULT_SELECTIONS(1.1, 1.3, 3112) + o2::framework::Configurable sign{"sign", 1, "Sign of the Sigma mother track (e.g. -1 for Sigma- or +1 for AntiSigma-)"}; +}; + +#undef KINK_DEFAULT_SELECTIONS + +constexpr const char PrefixSigmaSelection1[] = "SigmaSelection1"; +using ConfSigmaSelection1 = ConfSigmaSelection; + +/// The different selections for kinks +enum KinkSeles { + kKinkTopoDcaMax, + kTransRadMin, + kTransRadMax, + + kDauAbsEtaMax, + kDauDcaPvMin, + kMothDcaPvMax, + + kChaDaughTpcPion, + + kAlphaAPMax, + kQtAPMin, + kQtAPMax, + kCosPointingAngleMin, + + kKinkSelsMax +}; + +const char kinkSelsName[] = "Kink selection object"; +const std::unordered_map kinkSelsToStrings = { + {kKinkTopoDcaMax, "kinkTopoDcaMax"}, + {kTransRadMin, "transRadMin"}, + {kTransRadMax, "transRadMax"}, + {kDauAbsEtaMax, "dauAbsEtaMax"}, + {kDauDcaPvMin, "dauDcaPvMin"}, + {kMothDcaPvMax, "mothDcaPvMax"}, + {kChaDaughTpcPion, "chaDauTpcPion"}, + {kAlphaAPMax, "alphaAPMax"}, + {kQtAPMin, "qtAPMin"}, + {kQtAPMax, "qtAPMax"}, + {kCosPointingAngleMin, "cosPointingAngleMin"}}; + +/// \class KinkCuts +/// \brief Cut class to contain and execute all cuts applied to kinks +template +class KinkSelection : public BaseSelection +{ + public: + KinkSelection() {} + virtual ~KinkSelection() = default; + + template + void configure(T1& config, T2& filter) + { + mPtMin = filter.ptMin.value; + mPtMax = filter.ptMax.value; + mEtaMin = filter.etaMin.value; + mEtaMax = filter.etaMax.value; + mPhiMin = filter.phiMin.value; + mPhiMax = filter.phiMax.value; + + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { + mMassSigmaLowerLimit = filter.massMinSigma.value; + mMassSigmaUpperLimit = filter.massMaxSigma.value; + // Only add PID selection if we need it - will be checked at runtime + this->addSelection(config.chaDauTpcPion.value, kChaDaughTpcPion, limits::kAbsUpperLimit, true, true); + } + + this->addSelection(config.kinkTopoDcaMax.value, kKinkTopoDcaMax, limits::kUpperLimit, true, true); + this->addSelection(config.transRadMin.value, kTransRadMin, limits::kLowerLimit, true, true); + this->addSelection(config.transRadMax.value, kTransRadMax, limits::kUpperLimit, true, true); + this->addSelection(config.dauAbsEtaMax.value, kDauAbsEtaMax, limits::kAbsUpperLimit, true, true); + this->addSelection(config.dauDcaPvMin.value, kDauDcaPvMin, limits::kLowerLimit, true, true); + this->addSelection(config.mothDcaPvMax.value, kMothDcaPvMax, limits::kUpperLimit, true, true); + this->addSelection(config.alphaAPMax.value, kAlphaAPMax, limits::kUpperLimit, true, true); + this->addSelection(config.qtAPMin.value, kQtAPMin, limits::kLowerLimit, true, true); + this->addSelection(config.qtAPMax.value, kQtAPMax, limits::kUpperLimit, true, true); + this->addSelection(config.cosPointingAngleMin.value, kCosPointingAngleMin, limits::kLowerLimit, true, true); + }; + + template + void applySelections(T1 const& kinkCand, T2 const& /*tracks*/) + { + this->reset(); + // kink selections + std::array momMother = {kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}; + std::array momDaughter = {kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}; + + // Alpha_AP + std::array momMissing = {momMother[0] - momDaughter[0], momMother[1] - momDaughter[1], momMother[2] - momDaughter[2]}; + float lQlP = std::inner_product(momMother.begin(), momMother.end(), momDaughter.begin(), 0.f); + float lQlN = std::inner_product(momMother.begin(), momMother.end(), momMissing.begin(), 0.f); + float alphaAP = (lQlP + lQlN != 0.f) ? (lQlP - lQlN) / (lQlP + lQlN) : 0.f; + this->evaluateObservable(kAlphaAPMax, alphaAP); + + // qT_AP + float dp = lQlP; + float p2V0 = std::inner_product(momMother.begin(), momMother.end(), momMother.begin(), 0.f); + float p2A = std::inner_product(momDaughter.begin(), momDaughter.end(), momDaughter.begin(), 0.f); + float qtAP = std::sqrt(std::max(0.f, p2A - dp * dp / p2V0)); + this->evaluateObservable(kQtAPMin, qtAP); + this->evaluateObservable(kQtAPMax, qtAP); + + std::array vMother = {kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}; + float pMother = std::sqrt(std::inner_product(momMother.begin(), momMother.end(), momMother.begin(), 0.f)); + float vMotherNorm = std::sqrt(std::inner_product(vMother.begin(), vMother.end(), vMother.begin(), 0.f)); + float cosPointingAngle = (vMotherNorm > 0.f && pMother > 0.f) ? (std::inner_product(momMother.begin(), momMother.end(), vMother.begin(), 0.f)) / (pMother * vMotherNorm) : 0.f; + this->evaluateObservable(kCosPointingAngleMin, cosPointingAngle); + + this->evaluateObservable(kKinkTopoDcaMax, kinkCand.dcaKinkTopo()); + + // Compute transRadius + float transRadius = std::hypot(kinkCand.xDecVtx(), kinkCand.yDecVtx()); + this->evaluateObservable(kTransRadMin, transRadius); + this->evaluateObservable(kTransRadMax, transRadius); + + // Compute daughter eta + float pxDaug = kinkCand.pxDaug(); + float pyDaug = kinkCand.pyDaug(); + float pzDaug = kinkCand.pzDaug(); + float pDaug = std::sqrt(pxDaug * pxDaug + pyDaug * pyDaug + pzDaug * pzDaug); + float etaDaug = (pDaug > 0.f) ? 0.5f * std::log((pDaug + pzDaug) / (pDaug - pzDaug)) : 0.f; + this->evaluateObservable(kDauAbsEtaMax, std::fabs(etaDaug)); + + this->evaluateObservable(kDauDcaPvMin, std::abs(kinkCand.dcaDaugPv())); + this->evaluateObservable(kMothDcaPvMax, std::abs(kinkCand.dcaMothPv())); + + auto chaDaughter = kinkCand.template trackDaug_as(); + this->evaluateObservable(kChaDaughTpcPion, chaDaughter.tpcNSigmaPi()); + + this->assembleBitmask(); + }; + + template + bool checkFilters(const T& kink) const + { + // Compute mother eta and phi + float px = kink.pxMoth(); + float py = kink.pyMoth(); + float pz = kink.pzMoth(); + float pt = std::hypot(px, py); + float p = std::sqrt(px * px + py * py + pz * pz); + float eta = (p > 0.f) ? 0.5f * std::log((p + pz) / (p - pz)) : 0.f; + float phi = std::atan2(py, px); + + return ((pt > mPtMin && pt < mPtMax) && + (eta > mEtaMin && eta < mEtaMax) && + (phi > mPhiMin && phi < mPhiMax)); + } + + template + bool checkHypothesis(T const& kinkCand) const + { + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { + // Compute mass + float pxmoth = kinkCand.pxMoth(); + float pymoth = kinkCand.pyMoth(); + float pzmoth = kinkCand.pzMoth(); + float pxch = kinkCand.pxDaug(); + float pych = kinkCand.pyDaug(); + float pzch = kinkCand.pzDaug(); + float pxneut = pxmoth - pxch; + float pyneut = pymoth - pych; + float pzneut = pzmoth - pzch; + + float sigmaMass = RecoDecay::m( + std::array{std::array{pxch, pych, pzch}, std::array{pxneut, pyneut, pzneut}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron}); + + return (sigmaMass > mMassSigmaLowerLimit && sigmaMass < mMassSigmaUpperLimit); + } + return false; + } + + protected: + float mMassSigmaLowerLimit = 1.15f; + float mMassSigmaUpperLimit = 1.25f; + + // kinematic filters + float mPtMin = 0.f; + float mPtMax = 6.f; + float mEtaMin = -1.f; + float mEtaMax = 1.f; + float mPhiMin = 0.f; + float mPhiMax = o2::constants::math::TwoPI; +}; + +struct KinkBuilderProducts : o2::framework::ProducesGroup { + o2::framework::Produces producedSigmas; + o2::framework::Produces producedSigmaMasks; + o2::framework::Produces producedSigmaExtras; +}; + +struct ConfKinkTables : o2::framework::ConfigurableGroup { + std::string prefix = std::string("KinkTables"); + o2::framework::Configurable produceSigmas{"produceSigmas", -1, "Produce Sigmas (-1: auto; 0 off; 1 on)"}; + o2::framework::Configurable produceSigmaMasks{"produceSigmaMasks", -1, "Produce SigmaMasks (-1: auto; 0 off; 1 on)"}; + o2::framework::Configurable produceSigmaExtras{"produceSigmaExtras", -1, "Produce SigmaExtras (-1: auto; 0 off; 1 on)"}; +}; + +template +class KinkBuilder +{ + public: + KinkBuilder() {} + virtual ~KinkBuilder() = default; + + template + void init(T1& config, T2& filter, T3& table, T4& initContext) + { + mKinkSelection.configure(config, filter); + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { + LOG(info) << "Initialize femto Sigma builder..."; + mProduceSigmas = utils::enableTable("FSigmas_001", table.produceSigmas.value, initContext); + mProduceSigmaMasks = utils::enableTable("FSigmaMasks_001", table.produceSigmaMasks.value, initContext); + mProduceSigmaExtras = utils::enableTable("FSigmaExtras_001", table.produceSigmaExtras.value, initContext); + } + + if (mProduceSigmas || mProduceSigmaMasks || mProduceSigmaExtras) { + mFillAnyTable = true; + mKinkSelection.printSelections(kinkSelsName, kinkSelsToStrings); + } + } + + template + void fillKinks(T1& collisionProducts, T2& trackProducts, T3& kinkProducts, T4 const& kinks, T5 const& tracks, T6& trackBuilder, T7& indexMap) + { + if (!mFillAnyTable) { + return; + } + int64_t daughterIndex = 0; + for (const auto& kink : kinks) { + if (!mKinkSelection.checkFilters(kink)) { + continue; + } + mKinkSelection.applySelections(kink, tracks); + if (mKinkSelection.passesAllRequiredSelections() && mKinkSelection.checkHypothesis(kink)) { + auto daughter = kink.template trackDaug_as(); + daughterIndex = trackBuilder.template getDaughterIndex(daughter, trackProducts, collisionProducts, indexMap); + if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { + fillSigma(collisionProducts, kinkProducts, kink, daughterIndex); + } + } + } + } + + template + void fillSigma(T1& collisionProducts, T2& kinkProducts, T3 const& kink, int daughterIndex) + { + // Compute mass + float pxmoth = kink.pxMoth(); + float pymoth = kink.pyMoth(); + float pzmoth = kink.pzMoth(); + float pxch = kink.pxDaug(); + float pych = kink.pyDaug(); + float pzch = kink.pzDaug(); + float pxneut = pxmoth - pxch; + float pyneut = pymoth - pych; + float pzneut = pzmoth - pzch; + + float mass = RecoDecay::m( + std::array{std::array{pxch, pych, pzch}, std::array{pxneut, pyneut, pzneut}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron}); + + if (mProduceSigmas) { + // Compute mother eta and phi + float px = kink.pxMoth(); + float py = kink.pyMoth(); + float pz = kink.pzMoth(); + float pt = std::hypot(px, py); + float p = std::sqrt(px * px + py * py + pz * pz); + float eta = (p > 0.f) ? 0.5f * std::log((p + pz) / (p - pz)) : 0.f; + float phi = std::atan2(py, px); + + kinkProducts.producedSigmas(collisionProducts.producedCollision.lastIndex(), + kink.mothSign() * pt, + eta, + phi, + mass, + daughterIndex); + } + if (mProduceSigmaMasks) { + kinkProducts.producedSigmaMasks(mKinkSelection.getBitmask()); + } + if (mProduceSigmaExtras) { + // Compute kink angle and transRadius + float pMoth = std::sqrt(pxmoth * pxmoth + pymoth * pymoth + pzmoth * pzmoth); + float pDaug = std::sqrt(pxch * pxch + pych * pych + pzch * pzch); + float kinkAngle = 0.f; + if (pMoth > 0.f && pDaug > 0.f) { + float dotProduct = pxmoth * pxch + pymoth * pych + pzmoth * pzch; + float cosAngle = dotProduct / (pMoth * pDaug); + cosAngle = std::max(-1.0f, std::min(1.0f, cosAngle)); // Clamp + kinkAngle = std::acos(cosAngle); + } + + float transRadius = std::hypot(kink.xDecVtx(), kink.yDecVtx()); + + kinkProducts.producedSigmaExtras( + kinkAngle, + kink.dcaDaugPv(), + kink.dcaMothPv(), + kink.xDecVtx(), + kink.yDecVtx(), + kink.zDecVtx(), + transRadius); + } + } + bool fillAnyTable() { return mFillAnyTable; } + + private: + KinkSelection mKinkSelection; + bool mFillAnyTable = false; + bool mProduceSigmas = false; + bool mProduceSigmaMasks = false; + bool mProduceSigmaExtras = false; +}; +} // namespace kinkbuilder +} // namespace o2::analysis::femto +#endif // PWGCF_FEMTO_CORE_KINKBUILDER_H_ diff --git a/PWGCF/Femto/Core/kinkHistManager.h b/PWGCF/Femto/Core/kinkHistManager.h new file mode 100644 index 00000000000..8371a78d5e3 --- /dev/null +++ b/PWGCF/Femto/Core/kinkHistManager.h @@ -0,0 +1,289 @@ +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file kinkHistManager.h +/// \brief histogram manager for kink histograms +/// \author Anton Riedel, TU München, anton.riedel@cern.ch +/// \author Henrik Fribert, TU München, henrik.fribert@cern.ch + +#ifndef PWGCF_FEMTO_CORE_KINKHISTMANAGER_H_ +#define PWGCF_FEMTO_CORE_KINKHISTMANAGER_H_ + +#include "PWGCF/Femto/Core/histManager.h" +#include "PWGCF/Femto/Core/modes.h" +#include "PWGCF/Femto/Core/trackHistManager.h" + +#include "CommonConstants/MathConstants.h" +#include "Framework/Configurable.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" + +#include +#include +#include +#include +#include + +namespace o2::analysis::femto +{ +namespace kinkhistmanager +{ +// enum for kink histograms +enum KinkHist { + // analysis + kPt, + kEta, + kPhi, + kMass, + kSign, + // qa variables + kKinkAngle, + kDcaMothToPV, + kDcaDaugToPV, + kDecayVtxX, + kDecayVtxY, + kDecayVtxZ, + kDecayVtx, + kTransRadius, + // 2d qa + kPtVsEta, + kPtVsPhi, + kPhiVsEta, + kPtVsKinkAngle, + kPtVsDecayRadius, + kKinkHistLast +}; + +#define KINK_DEFAULT_BINNING(defaultMassMin, defaultMassMax) \ + o2::framework::ConfigurableAxis pt{"pt", {{600, 0, 6}}, "Pt"}; \ + o2::framework::ConfigurableAxis eta{"eta", {{300, -1.5, 1.5}}, "Eta"}; \ + o2::framework::ConfigurableAxis phi{"phi", {{720, 0, 1.f * o2::constants::math::TwoPI}}, "Phi"}; \ + o2::framework::ConfigurableAxis mass{"mass", {{200, defaultMassMin, defaultMassMax}}, "Mass"}; \ + o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; + +template +struct ConfSigmaBinning : o2::framework::ConfigurableGroup { + std::string prefix = Prefix; + KINK_DEFAULT_BINNING(1.1, 1.3) +}; +#undef KINK_DEFAULT_BINNING + +constexpr const char PrefixSigmaBinning1[] = "SigmaBinning1"; +using ConfSigmaBinning1 = ConfSigmaBinning; + +template +struct ConfKinkQaBinning : o2::framework::ConfigurableGroup { + std::string prefix = Prefix; + o2::framework::ConfigurableAxis kinkAngle{"kinkAngle", {{100, 0, 3.15}}, "Kink Angle (rad)"}; + o2::framework::ConfigurableAxis dcaMothToPV{"dcaMothToPV", {{150, 0, 1.5}}, "Mother DCA to PV (cm)"}; + o2::framework::ConfigurableAxis dcaDaugToPV{"dcaDaugToPV", {{1000, 0, 100}}, "Daughter DCA to PV (cm)"}; + o2::framework::ConfigurableAxis decayVertex{"decayVertex", {{100, 0, 100}}, "Decay vertex position (cm)"}; + o2::framework::ConfigurableAxis transRadius{"transRadius", {{100, 0, 100}}, "Transverse radius (cm)"}; +}; + +constexpr const char PrefixSigmaQaBinning1[] = "SigmaQaBinning1"; +using ConfSigmaQaBinning1 = ConfKinkQaBinning; + +// must be in sync with enum KinkHist +// the enum gives the correct index in the array +constexpr std::array, kKinkHistLast> HistTable = { + {{kPt, o2::framework::kTH1F, "hPt", "Transverse Momentum; p_{T} (GeV/#it{c}); Entries"}, + {kEta, o2::framework::kTH1F, "hEta", "Pseudorapidity; #eta; Entries"}, + {kPhi, o2::framework::kTH1F, "hPhi", "Azimuthal angle; #varphi; Entries"}, + {kMass, o2::framework::kTH1F, "hMass", "Invariant Mass; m_{Inv} (GeV/#it{c}^{2}); Entries"}, + {kSign, o2::framework::kTH1F, "hSign", "Sign; sign; Entries"}, + {kKinkAngle, o2::framework::kTH1F, "hKinkAngle", "Kink Angle; Angle (rad); Entries"}, + {kDcaMothToPV, o2::framework::kTH1F, "hDcaMothToPV", "Mother DCA to PV; DCA (cm); Entries"}, + {kDcaDaugToPV, o2::framework::kTH1F, "hDcaDaugToPV", "Daughter DCA to PV; DCA (cm); Entries"}, + {kDecayVtxX, o2::framework::kTH1F, "hDecayVtxX", "Decay Vertex X; x (cm); Entries"}, + {kDecayVtxY, o2::framework::kTH1F, "hDecayVtxY", "Decay Vertex Y; y (cm); Entries"}, + {kDecayVtxZ, o2::framework::kTH1F, "hDecayVtxZ", "Decay Vertex Z; z (cm); Entries"}, + {kDecayVtx, o2::framework::kTH1F, "hDecayVtx", "Decay Distance from PV; r (cm); Entries"}, + {kTransRadius, o2::framework::kTH1F, "hTransRadius", "Transverse Decay Radius; r_{xy} (cm); Entries"}, + {kPtVsEta, o2::framework::kTH2F, "hPtVsEta", "p_{T} vs #eta; p_{T} (GeV/#it{c}); #eta"}, + {kPtVsPhi, o2::framework::kTH2F, "hPtVsPhi", "p_{T} vs #varphi; p_{T} (GeV/#it{c}); #varphi"}, + {kPhiVsEta, o2::framework::kTH2F, "hPhiVsEta", "#varphi vs #eta; #varphi; #eta"}, + {kPtVsKinkAngle, o2::framework::kTH2F, "hPtVsKinkAngle", "p_{T} vs kink angle; p_{T} (GeV/#it{c}); kink angle (rad)"}, + {kPtVsDecayRadius, o2::framework::kTH2F, "hPtVsDecayRadius", "p_{T} vs transverse decay radius; p_{T} (GeV/#it{c}); r_{xy} (cm)"}}}; + +template +auto makeKinkHistSpecMap(const T& confBinningAnalysis) +{ + return std::map>{ + {kPt, {confBinningAnalysis.pt}}, + {kEta, {confBinningAnalysis.eta}}, + {kPhi, {confBinningAnalysis.phi}}, + {kMass, {confBinningAnalysis.mass}}, + {kSign, {confBinningAnalysis.sign}}}; +} + +template +std::map> makeKinkQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) +{ + return std::map>{ + {kPt, {confBinningAnalysis.pt}}, + {kEta, {confBinningAnalysis.eta}}, + {kPhi, {confBinningAnalysis.phi}}, + {kMass, {confBinningAnalysis.mass}}, + {kSign, {confBinningAnalysis.sign}}, + {kKinkAngle, {confBinningQa.kinkAngle}}, + {kDcaMothToPV, {confBinningQa.dcaMothToPV}}, + {kDcaDaugToPV, {confBinningQa.dcaDaugToPV}}, + {kDecayVtxX, {confBinningQa.decayVertex}}, + {kDecayVtxY, {confBinningQa.decayVertex}}, + {kDecayVtxZ, {confBinningQa.decayVertex}}, + {kDecayVtx, {confBinningQa.decayVertex}}, + {kTransRadius, {confBinningQa.transRadius}}, + {kPtVsEta, {confBinningAnalysis.pt, confBinningAnalysis.eta}}, + {kPtVsPhi, {confBinningAnalysis.pt, confBinningAnalysis.phi}}, + {kPhiVsEta, {confBinningAnalysis.phi, confBinningAnalysis.eta}}, + {kPtVsKinkAngle, {confBinningAnalysis.pt, confBinningQa.kinkAngle}}, + {kPtVsDecayRadius, {confBinningAnalysis.pt, confBinningQa.transRadius}}}; +} + +constexpr char PrefixSigmaQa[] = "SigmaQA/"; +constexpr char PrefixSigma[] = "Sigma/"; +constexpr char PrefixKinkChaDaughter[] = "KinkChaDau/"; +constexpr char PrefixKinkChaDaughterQa[] = "KinkChaDauQa/"; + +constexpr std::string_view AnalysisDir = "Kinematics/"; +constexpr std::string_view QaDir = "QA/"; + +/// \class KinkHistManager +/// \brief Class for histogramming event properties +// template +template +class KinkHistManager +{ + public: + /// Destructor + virtual ~KinkHistManager() = default; + + /// Initializes histograms for the task + /// \param registry Histogram registry to be passed + /// + void init(o2::framework::HistogramRegistry* registry, + std::map> KinkSpecs, + std::map> ChaDauSpecs) + { + mHistogramRegistry = registry; + mChaDauManager.init(registry, ChaDauSpecs); + + if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { + std::string analysisDir = std::string(kinkPrefix) + std::string(AnalysisDir); + mHistogramRegistry->add(analysisDir + GetHistNamev2(kPt, HistTable), GetHistDesc(kPt, HistTable), GetHistType(kPt, HistTable), {KinkSpecs[kPt]}); + mHistogramRegistry->add(analysisDir + GetHistNamev2(kEta, HistTable), GetHistDesc(kEta, HistTable), GetHistType(kEta, HistTable), {KinkSpecs[kEta]}); + mHistogramRegistry->add(analysisDir + GetHistNamev2(kPhi, HistTable), GetHistDesc(kPhi, HistTable), GetHistType(kPhi, HistTable), {KinkSpecs[kPhi]}); + mHistogramRegistry->add(analysisDir + GetHistNamev2(kMass, HistTable), GetHistDesc(kMass, HistTable), GetHistType(kMass, HistTable), {KinkSpecs[kMass]}); + mHistogramRegistry->add(analysisDir + GetHistNamev2(kSign, HistTable), GetHistDesc(kSign, HistTable), GetHistType(kSign, HistTable), {KinkSpecs[kSign]}); + } + + if constexpr (isFlagSet(mode, modes::Mode::kQa)) { + std::string qaDir = std::string(kinkPrefix) + std::string(QaDir); + + // Basic kinematic histograms + mHistogramRegistry->add(qaDir + GetHistNamev2(kPt, HistTable), GetHistDesc(kPt, HistTable), GetHistType(kPt, HistTable), {KinkSpecs[kPt]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kEta, HistTable), GetHistDesc(kEta, HistTable), GetHistType(kEta, HistTable), {KinkSpecs[kEta]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kPhi, HistTable), GetHistDesc(kPhi, HistTable), GetHistType(kPhi, HistTable), {KinkSpecs[kPhi]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kMass, HistTable), GetHistDesc(kMass, HistTable), GetHistType(kMass, HistTable), {KinkSpecs[kMass]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kSign, HistTable), GetHistDesc(kSign, HistTable), GetHistType(kSign, HistTable), {KinkSpecs[kSign]}); + + // Kink-specific QA histograms + mHistogramRegistry->add(qaDir + GetHistNamev2(kKinkAngle, HistTable), GetHistDesc(kKinkAngle, HistTable), GetHistType(kKinkAngle, HistTable), {KinkSpecs[kKinkAngle]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDcaMothToPV, HistTable), GetHistDesc(kDcaMothToPV, HistTable), GetHistType(kDcaMothToPV, HistTable), {KinkSpecs[kDcaMothToPV]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDcaDaugToPV, HistTable), GetHistDesc(kDcaDaugToPV, HistTable), GetHistType(kDcaDaugToPV, HistTable), {KinkSpecs[kDcaDaugToPV]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDecayVtxX, HistTable), GetHistDesc(kDecayVtxX, HistTable), GetHistType(kDecayVtxX, HistTable), {KinkSpecs[kDecayVtxX]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDecayVtxY, HistTable), GetHistDesc(kDecayVtxY, HistTable), GetHistType(kDecayVtxY, HistTable), {KinkSpecs[kDecayVtxY]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDecayVtxZ, HistTable), GetHistDesc(kDecayVtxZ, HistTable), GetHistType(kDecayVtxZ, HistTable), {KinkSpecs[kDecayVtxZ]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kDecayVtx, HistTable), GetHistDesc(kDecayVtx, HistTable), GetHistType(kDecayVtx, HistTable), {KinkSpecs[kDecayVtx]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kTransRadius, HistTable), GetHistDesc(kTransRadius, HistTable), GetHistType(kTransRadius, HistTable), {KinkSpecs[kTransRadius]}); + + // 2D QA histograms + mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsEta, HistTable), GetHistDesc(kPtVsEta, HistTable), GetHistType(kPtVsEta, HistTable), {KinkSpecs[kPtVsEta]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsPhi, HistTable), GetHistDesc(kPtVsPhi, HistTable), GetHistType(kPtVsPhi, HistTable), {KinkSpecs[kPtVsPhi]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kPhiVsEta, HistTable), GetHistDesc(kPhiVsEta, HistTable), GetHistType(kPhiVsEta, HistTable), {KinkSpecs[kPhiVsEta]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsKinkAngle, HistTable), GetHistDesc(kPtVsKinkAngle, HistTable), GetHistType(kPtVsKinkAngle, HistTable), {KinkSpecs[kPtVsKinkAngle]}); + mHistogramRegistry->add(qaDir + GetHistNamev2(kPtVsDecayRadius, HistTable), GetHistDesc(kPtVsDecayRadius, HistTable), GetHistType(kPtVsDecayRadius, HistTable), {KinkSpecs[kPtVsDecayRadius]}); + } + } + + /// Fill histograms for kink candidates + /// \param kinkcandidate Kink candidate to fill histograms for + template + void fill(T const& kinkcandidate) + { + if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(AnalysisDir) + HIST(GetHistName(kPt, HistTable)), kinkcandidate.pt()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(AnalysisDir) + HIST(GetHistName(kEta, HistTable)), kinkcandidate.eta()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(AnalysisDir) + HIST(GetHistName(kPhi, HistTable)), kinkcandidate.phi()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(AnalysisDir) + HIST(GetHistName(kMass, HistTable)), kinkcandidate.mass()); + + if constexpr (isEqual(kink, modes::Kink::kSigma)) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(AnalysisDir) + HIST(GetHistName(kSign, HistTable)), kinkcandidate.sign()); + } + } + + if constexpr (isFlagSet(mode, modes::Mode::kQa)) { + // Basic kinematic histograms + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPt, HistTable)), kinkcandidate.pt()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kEta, HistTable)), kinkcandidate.eta()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPhi, HistTable)), kinkcandidate.phi()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kMass, HistTable)), kinkcandidate.mass()); + + if constexpr (isEqual(kink, modes::Kink::kSigma)) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kSign, HistTable)), kinkcandidate.sign()); + } + + // Kink-specific QA histograms + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kKinkAngle, HistTable)), kinkcandidate.kinkAngle()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDcaMothToPV, HistTable)), kinkcandidate.dcaMothToPV()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDcaDaugToPV, HistTable)), kinkcandidate.dcaDaugToPV()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDecayVtxX, HistTable)), kinkcandidate.decayVtxX()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDecayVtxY, HistTable)), kinkcandidate.decayVtxY()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDecayVtxZ, HistTable)), kinkcandidate.decayVtxZ()); + + // Calculate decay distance from PV + float decayDistance = std::sqrt(kinkcandidate.decayVtxX() * kinkcandidate.decayVtxX() + + kinkcandidate.decayVtxY() * kinkcandidate.decayVtxY() + + kinkcandidate.decayVtxZ() * kinkcandidate.decayVtxZ()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kDecayVtx, HistTable)), decayDistance); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kTransRadius, HistTable)), kinkcandidate.transRadius()); + + // 2D QA histograms + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPtVsEta, HistTable)), kinkcandidate.pt(), kinkcandidate.eta()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPtVsPhi, HistTable)), kinkcandidate.pt(), kinkcandidate.phi()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPhiVsEta, HistTable)), kinkcandidate.phi(), kinkcandidate.eta()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPtVsKinkAngle, HistTable)), kinkcandidate.pt(), kinkcandidate.kinkAngle()); + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(QaDir) + HIST(GetHistName(kPtVsDecayRadius, HistTable)), kinkcandidate.pt(), kinkcandidate.transRadius()); + } + } + + /// Fill histograms for kink candidates - overload with track table argument + /// \param kinkcandidate Kink candidate to fill histograms for + /// \param tracks Track table for daughter access + template + void fill(T1 const& kinkcandidate, T2 const& /*tracks*/) + { + auto chaDaughter = kinkcandidate.template chaDau_as(); + mChaDauManager.fill(chaDaughter); + fill(kinkcandidate); + } + + private: + o2::framework::HistogramRegistry* mHistogramRegistry; + trackhistmanager::TrackHistManager mChaDauManager; +}; +}; // namespace kinkhistmanager +}; // namespace o2::analysis::femto +#endif // PWGCF_FEMTO_CORE_KINKHISTMANAGER_H_ diff --git a/PWGCF/Femto/Core/modes.h b/PWGCF/Femto/Core/modes.h index dccd9ac8559..012442bb3f1 100644 --- a/PWGCF/Femto/Core/modes.h +++ b/PWGCF/Femto/Core/modes.h @@ -69,7 +69,8 @@ enum class Track : o2::aod::femtodatatypes::TrackType { kPrimaryTrack, kV0Daughter, kCascadeBachelor, - kResonanceDaughter + kResonanceDaughter, + kKinkDaughter }; enum class V0 : o2::aod::femtodatatypes::V0Type { @@ -78,6 +79,10 @@ enum class V0 : o2::aod::femtodatatypes::V0Type { kK0short }; +enum class Kink : o2::aod::femtodatatypes::KinkType { + kSigma +}; + enum class Cascade : o2::aod::femtodatatypes::CascadeType { kXi, kOmega @@ -95,14 +100,16 @@ enum class Pairs : o2::aod::femtodatatypes::PairType { kTrackTrack, kTrackV0, kTrackResonance, - kTrackCascade + kTrackCascade, + kTrackKink }; enum class TrackPairs : o2::aod::femtodatatypes::PairType { kTrackTrack, kTrackPosDaughter, kTrackNegDaughter, - kTrackBachelor + kTrackBachelor, + kTrackChaDaughter }; }; // namespace modes diff --git a/PWGCF/Femto/Core/partitions.h b/PWGCF/Femto/Core/partitions.h index adbac4d2213..ffbdf6fd214 100644 --- a/PWGCF/Femto/Core/partitions.h +++ b/PWGCF/Femto/Core/partitions.h @@ -110,4 +110,16 @@ (femtobase::stored::mass < selection.massMax) && \ ncheckbit(femtocascades::mask, selection.mask) +#define MAKE_SIGMA_PARTITION(selection) \ + ifnode(selection.sign.node() > 0, femtobase::stored::signedPt > 0.f, femtobase::stored::signedPt < 0.f) && \ + (nabs(femtobase::stored::signedPt) > selection.ptMin) && \ + (nabs(femtobase::stored::signedPt) < selection.ptMax) && \ + (femtobase::stored::eta > selection.etaMin) && \ + (femtobase::stored::eta < selection.etaMax) && \ + (femtobase::stored::phi > selection.phiMin) && \ + (femtobase::stored::phi < selection.phiMax) && \ + (femtobase::stored::mass > selection.massMin) && \ + (femtobase::stored::mass < selection.massMax) && \ + ncheckbit(femtokinks::mask, selection.mask) + #endif // PWGCF_FEMTO_CORE_PARTITIONS_H_ diff --git a/PWGCF/Femto/Core/trackHistManager.h b/PWGCF/Femto/Core/trackHistManager.h index 3f22459ac5b..677e88e5487 100644 --- a/PWGCF/Femto/Core/trackHistManager.h +++ b/PWGCF/Femto/Core/trackHistManager.h @@ -129,6 +129,7 @@ constexpr const char PrefixV0NegDauBinning[] = "V0NegDauBinning"; constexpr const char PrefixCascadePosDauBinning[] = "CascadePosDauBinning"; constexpr const char PrefixCascadeNegDauBinning[] = "CascadeNegDauBinning"; constexpr const char PrefixCascadeBachelorBinning[] = "CascadeBachelorBinning"; +constexpr const char PrefixKinkChaDauBinning[] = "KinkChaDauBinning"; using ConfTrackBinning1 = ConfTrackBinning; using ConfTrackBinning2 = ConfTrackBinning; @@ -136,9 +137,11 @@ using ConfResonancePosDauBinning = ConfTrackBinning; using ConfV0PosDauBinning = ConfTrackBinning; using ConfV0NegDauBinning = ConfTrackBinning; +using ConfKinkChaDauBinning = ConfTrackBinning; using ConfCascadePosDauBinning = ConfTrackBinning; using ConfCascadeNegDauBinning = ConfTrackBinning; using ConfCascadeBachelorBinning = ConfTrackBinning; +using ConfKinkChaDauBinning = ConfTrackBinning; template struct ConfTrackQaBinning : o2::framework::ConfigurableGroup { @@ -203,6 +206,7 @@ constexpr const char PrefixV0NegDauQaBinning[] = "V0NegDauQaBinning"; constexpr const char PrefixCascadePosDauQaBinning[] = "CascadePosDauQaBinning"; constexpr const char PrefixCascadeNegDauQaBinning[] = "CascadeNegDauQaBinning"; constexpr const char PrefixCascadeBachelorQaBinning[] = "CascadeBachelorQaBinning"; +constexpr const char PrefixKinkChaDauQaBinning[] = "KinkChaDauQaBinning"; using ConfTrackQaBinning1 = ConfTrackQaBinning; using ConfTrackQaBinning2 = ConfTrackQaBinning; @@ -213,6 +217,7 @@ using ConfV0NegDauQaBinning = ConfTrackQaBinning; using ConfCascadePosDauQaBinning = ConfTrackQaBinning; using ConfCascadeNegDauQaBinning = ConfTrackQaBinning; using ConfCascadeBachelorQaBinning = ConfTrackQaBinning; +using ConfKinkChaDauQaBinning = ConfTrackQaBinning; // must be in sync with enum TrackVariables // the enum gives the correct index in the array @@ -374,6 +379,9 @@ constexpr char PrefixCascadePosDaughterQa[] = "CascadePosDauQa/"; constexpr char PrefixCascadeNegDaughterQa[] = "CascadeNegDauQa/"; constexpr char PrefixCascadeBachelorQa[] = "CascadeBachelorQa/"; +constexpr char PrefixKinkChaDaughter[] = "KinkChaDau/"; +constexpr char PrefixKinkChaDaughterQa[] = "KinkChaDauQa/"; + constexpr std::string_view AnalysisDir = "Kinematics/"; constexpr std::string_view QaDir = "QA/"; constexpr std::string_view PidDir = "PID/"; diff --git a/PWGCF/Femto/DataModel/FemtoTables.h b/PWGCF/Femto/DataModel/FemtoTables.h index 1708ba66ba7..c7bc8334c52 100644 --- a/PWGCF/Femto/DataModel/FemtoTables.h +++ b/PWGCF/Femto/DataModel/FemtoTables.h @@ -488,6 +488,57 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FK0shortExtras_001, "FK0SHORTEXTRA", 1, //! k using FK0shortExtras = FK0shortExtras_001; +namespace femtokinks +{ +// columns for bit masks +DECLARE_SOA_COLUMN(Mask, mask, femtodatatypes::KinkMaskType); //! Bitmask for kink selections + +// columns for debug information +DECLARE_SOA_COLUMN(KinkAngle, kinkAngle, float); //! Kink angle between mother and charged daughter at decay vertex +DECLARE_SOA_COLUMN(DcaMothToPV, dcaMothToPV, float); //! DCA of the mother track to the primary vertex +DECLARE_SOA_COLUMN(DcaDaugToPV, dcaDaugToPV, float); //! DCA of the charged daughter track to the primary vertex +DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float); //! x coordinate of decay vertex (relative to PV) +DECLARE_SOA_COLUMN(DecayVtxY, decayVtxY, float); //! y coordinate of decay vertex (relative to PV) +DECLARE_SOA_COLUMN(DecayVtxZ, decayVtxZ, float); //! z coordinate of decay vertex (relative to PV) +DECLARE_SOA_COLUMN(TransRadius, transRadius, float); //! Transverse decay radius from PV + +// id column for charged daughter track +DECLARE_SOA_INDEX_COLUMN_FULL(ChaDau, chaDau, int32_t, FTracks, "_ChaDau"); //! +} // namespace femtokinks + +// table for basic sigma minus information +DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmas_001, "FSIGMA", 1, + o2::soa::Index<>, + femtobase::stored::CollisionId, // use sign to differentiate between sigma minus (-1) and anti sigma minus (+1) + femtobase::stored::SignedPt, + femtobase::stored::Eta, + femtobase::stored::Phi, + femtobase::stored::Mass, + femtokinks::ChaDauId, + femtobase::dynamic::Sign, + femtobase::dynamic::Pt, + femtobase::dynamic::P, + femtobase::dynamic::Px, + femtobase::dynamic::Py, + femtobase::dynamic::Pz, + femtobase::dynamic::Theta); +using FSigmas = FSigmas_001; + +DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmaMasks_001, "FSIGMAMASKS", 1, + femtokinks::Mask); +using FSigmaMasks = FSigmaMasks_001; + +DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmaExtras_001, "FSIGMAEXTRAS", 1, + femtokinks::KinkAngle, + femtokinks::DcaDaugToPV, + femtokinks::DcaMothToPV, + femtokinks::DecayVtxX, + femtokinks::DecayVtxY, + femtokinks::DecayVtxZ, + femtokinks::TransRadius); + +using FSigmaExtras = FSigmaExtras_001; + namespace femtocascades { // columns for cascade bit masks diff --git a/PWGCF/Femto/TableProducer/femtoProducer.cxx b/PWGCF/Femto/TableProducer/femtoProducer.cxx index 205ca34a986..c5140fb5855 100644 --- a/PWGCF/Femto/TableProducer/femtoProducer.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducer.cxx @@ -15,10 +15,12 @@ #include "PWGCF/Femto/Core/cascadeBuilder.h" #include "PWGCF/Femto/Core/collisionBuilder.h" +#include "PWGCF/Femto/Core/kinkBuilder.h" #include "PWGCF/Femto/Core/modes.h" #include "PWGCF/Femto/Core/trackBuilder.h" #include "PWGCF/Femto/Core/twoTrackResonanceBuilder.h" #include "PWGCF/Femto/Core/v0Builder.h" +#include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "Common/DataModel/Centrality.h" @@ -70,6 +72,8 @@ using Run3PpVzeros = V0Datas; using Run3PpCascades = CascDatas; +using Run3PpKinks = KinkCands; + } // namespace consumeddata } // namespace o2::analysis::femto @@ -116,6 +120,13 @@ struct FemtoProducer { cascadebuilder::ConfOmegaBits confOmegaBits; cascadebuilder::CascadeBuilder omegaBuilder; + // kink builder + kinkbuilder::KinkBuilderProducts kinkBuilderProducts; + kinkbuilder::ConfKinkTables confKinkTables; + kinkbuilder::ConfKinkFilters confKinkFilters; + kinkbuilder::ConfSigmaBits confSigmaBits; + kinkbuilder::KinkBuilder sigmaBuilder; + // resonance daughter filters and partitions twotrackresonancebuilder::ConfTwoTrackResonanceDaughterFilters confResonanceDaughterFilters; // caching and preslicing @@ -191,6 +202,9 @@ struct FemtoProducer { lambdaBuilder.init(confLambdaBits, confV0Filters, confV0Tables, context); antilambdaBuilder.init(confLambdaBits, confV0Filters, confV0Tables, context); + // configure kink builder + sigmaBuilder.init(confSigmaBits, confKinkFilters, confKinkTables, context); + // cascade selections xiBuilder.init(confXiBits, confCascadeFilters, confCascadeTables, context); omegaBuilder.init(confOmegaBits, confCascadeFilters, confCascadeTables, context); @@ -201,11 +215,14 @@ struct FemtoProducer { kstar0Builder.init(confKstar0Bits, confKstarFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); kstar0barBuilder.init(confKstar0Bits, confKstarFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); - if ((xiBuilder.fillAnyTable() || omegaBuilder.fillAnyTable()) && !doprocessTracksV0sCascadesRun3pp) { - LOG(fatal) << "At least one cascade tabel is enabled, but wrong process function is enabled. Breaking..."; + if ((xiBuilder.fillAnyTable() || omegaBuilder.fillAnyTable()) && (!doprocessTracksV0sCascadesRun3pp && !doprocessTracksV0sCascadesKinksRun3pp)) { + LOG(fatal) << "At least one cascade table is enabled, but wrong process function is enabled. Breaking..."; + } + if ((lambdaBuilder.fillAnyTable() || antilambdaBuilder.fillAnyTable() || k0shortBuilder.fillAnyTable()) && (!doprocessTracksV0sCascadesRun3pp && !doprocessTracksV0sRun3pp && !doprocessTracksV0sCascadesKinksRun3pp)) { + LOG(info) << "At least one v0 table is enabled, but wrong process function is enabled. Breaking..."; } - if ((lambdaBuilder.fillAnyTable() || antilambdaBuilder.fillAnyTable() || k0shortBuilder.fillAnyTable()) && (!doprocessTracksV0sCascadesRun3pp && !doprocessTracksV0sRun3pp)) { - LOG(info) << "At least one v0 tabel is enbaled, but wrong process function is enabled. Breaking..."; + if (sigmaBuilder.fillAnyTable() && !doprocessTracksKinksRun3pp && !doprocessTracksV0sCascadesKinksRun3pp) { + LOG(fatal) << "At least one kink table is enabled, but wrong process function is enabled. Breaking..."; } } @@ -244,6 +261,14 @@ struct FemtoProducer { k0shortBuilder.fillV0s(collisionBuilderProducts, trackBuilderProducts, v0builderProducts, v0s, tracks, trackBuilder, indexMapTracks); } + // add kinks + template + void processTracksKinks(T1 const& col, T2 const& bcs, T3 const& tracks, T4 const& tracksWithItsPid, T5 const& kinks) + { + processTracks(col, bcs, tracks, tracksWithItsPid); + sigmaBuilder.fillKinks(collisionBuilderProducts, trackBuilderProducts, kinkBuilderProducts, kinks, tracks, trackBuilder, indexMapTracks); + } + // add cascades template void processTracksV0sCascades(T1 const& col, T2 const& bcs, T3 const& tracks, T4 const& tracksWithItsPid, T5 const& v0s, T6 const& cascades) @@ -255,6 +280,18 @@ struct FemtoProducer { cascades, tracks, col, trackBuilder, indexMapTracks); } + // add kinks + template + void processTracksV0sCascadesKinks(T1 const& col, T2 const& bcs, T3 const& tracks, T4 const& tracksWithItsPid, T5 const& v0s, T6 const& cascades, T7 const& kinks) + { + processTracksV0s(col, bcs, tracks, tracksWithItsPid, v0s); + sigmaBuilder.fillKinks(collisionBuilderProducts, trackBuilderProducts, kinkBuilderProducts, kinks, tracks, trackBuilder, indexMapTracks); + xiBuilder.fillCascades(collisionBuilderProducts, trackBuilderProducts, cascadeBuilderProducts, + cascades, tracks, col, trackBuilder, indexMapTracks); + omegaBuilder.fillCascades(collisionBuilderProducts, trackBuilderProducts, cascadeBuilderProducts, + cascades, tracks, col, trackBuilder, indexMapTracks); + } + // proccess functions void processTracksRun3pp(consumeddata::Run3PpCollisions::iterator const& col, BCsWithTimestamps const& bcs, @@ -280,7 +317,20 @@ struct FemtoProducer { }; PROCESS_SWITCH(FemtoProducer, processTracksV0sRun3pp, "Process tracks and v0s", false); - // process tracks, v0s and casacades + // process tracks and kinks + void processTracksKinksRun3pp(consumeddata::Run3PpCollisions::iterator const& col, + BCsWithTimestamps const& bcs, + consumeddata::Run3FullPidTracks const& tracks, + consumeddata::Run3PpKinks const& kinks) + { + // its pid information is generated dynamically, so we need to add it here + auto tracksWithItsPid = o2::soa::Attach(tracks); + processTracksKinks(col, bcs, tracks, tracksWithItsPid, kinks); + } + PROCESS_SWITCH(FemtoProducer, processTracksKinksRun3pp, "Process tracks and kinks", false); + + // process tracks, v0s and cascades void processTracksV0sCascadesRun3pp(consumeddata::Run3PpCollisions::iterator const& col, BCsWithTimestamps const& bcs, consumeddata::Run3FullPidTracks const& tracks, @@ -293,6 +343,21 @@ struct FemtoProducer { processTracksV0sCascades(col, bcs, tracks, tracksWithItsPid, v0s, cascades); } PROCESS_SWITCH(FemtoProducer, processTracksV0sCascadesRun3pp, "Provide Tracks, V0s and Cascades for Run3", false); + + // process tracks, v0s, cascades and kinks + void processTracksV0sCascadesKinksRun3pp(consumeddata::Run3PpCollisions::iterator const& col, + BCsWithTimestamps const& bcs, + consumeddata::Run3FullPidTracks const& tracks, + consumeddata::Run3PpVzeros const& v0s, + consumeddata::Run3PpCascades const& cascades, + consumeddata::Run3PpKinks const& kinks) + { + // its pid information is generated dynamically, so we need to add it here + auto tracksWithItsPid = o2::soa::Attach(tracks); + processTracksV0sCascadesKinks(col, bcs, tracks, tracksWithItsPid, v0s, cascades, kinks); + } + PROCESS_SWITCH(FemtoProducer, processTracksV0sCascadesKinksRun3pp, "Provide Tracks, V0s and Cascades for Run3", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto/Tasks/CMakeLists.txt b/PWGCF/Femto/Tasks/CMakeLists.txt index eb16171c62b..7ed493cc309 100644 --- a/PWGCF/Femto/Tasks/CMakeLists.txt +++ b/PWGCF/Femto/Tasks/CMakeLists.txt @@ -24,6 +24,11 @@ o2physics_add_dpl_workflow(femto-v0-qa PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(femto-kink-qa + SOURCES femtoKinkQa.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(femto-cascade-qa SOURCES femtoCascadeQa.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGCF/Femto/Tasks/femtoKinkQa.cxx b/PWGCF/Femto/Tasks/femtoKinkQa.cxx new file mode 100644 index 00000000000..1fdb4e58693 --- /dev/null +++ b/PWGCF/Femto/Tasks/femtoKinkQa.cxx @@ -0,0 +1,120 @@ +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file femtoKinkQa.cxx +/// \brief QA task for kinks +/// \author Anton Riedel, TU München, anton.riedel@cern.ch +/// \author Henrik Fribert, TU München, henrik.fribert@cern.ch + +#include "PWGCF/Femto/Core/collisionBuilder.h" +#include "PWGCF/Femto/Core/collisionHistManager.h" +#include "PWGCF/Femto/Core/kinkBuilder.h" +#include "PWGCF/Femto/Core/kinkHistManager.h" +#include "PWGCF/Femto/Core/modes.h" +#include "PWGCF/Femto/Core/partitions.h" +#include "PWGCF/Femto/Core/trackHistManager.h" +#include "PWGCF/Femto/DataModel/FemtoTables.h" +#include "PWGLF/DataModel/LFKinkDecayTables.h" + +#include "Framework/ASoA.h" +#include "Framework/AnalysisHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/Configurable.h" +#include "Framework/Expressions.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/InitContext.h" +#include "Framework/OutputObjHeader.h" +#include "Framework/runDataProcessing.h" + +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::analysis::femto; + +struct FemtoKinkQa { + + // setup for collisions + collisionbuilder::ConfCollisionFilters collisionSelection; + Filter collisionFilter = MAKE_COLLISION_FILTER(collisionSelection); + + colhistmanager::CollisionHistManager colHistManager; + colhistmanager::ConfCollisionBinning confCollisionBinning; + + // using Collisions = o2::soa::Join; + using Collisions = FCols; + using Collision = Collisions::iterator; + + using FilteredCollisions = o2::soa::Filtered; + using FilteredCollision = FilteredCollisions::iterator; + + // Define kink/sigma tables (joining tables for comprehensive information) + using Sigmas = o2::soa::Join; + using Tracks = o2::soa::Join; + + SliceCache cache; + + // setup for sigmas + kinkbuilder::ConfSigmaSelection1 confSigmaSelection; + + Partition sigmaPartition = MAKE_SIGMA_PARTITION(confSigmaSelection); + Preslice perColSigmas = aod::femtobase::stored::collisionId; + + kinkhistmanager::ConfSigmaBinning1 confSigmaBinning; + kinkhistmanager::ConfSigmaQaBinning1 confSigmaQaBinning; + kinkhistmanager::KinkHistManager< + kinkhistmanager::PrefixSigmaQa, + trackhistmanager::PrefixKinkChaDaughterQa, + modes::Mode::kAnalysis_Qa, + modes::Kink::kSigma> + sigmaHistManager; + + // setup for daughters + trackhistmanager::ConfKinkChaDauBinning confKinkChaDaughterBinning; + trackhistmanager::ConfKinkChaDauQaBinning confKinkChaDaughterQaBinning; + + HistogramRegistry hRegistry{"FemtoKinkQa", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext&) + { + auto sigmaHistSpec = kinkhistmanager::makeKinkQaHistSpecMap(confSigmaBinning, confSigmaQaBinning); + auto chaDauHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confKinkChaDaughterBinning, confKinkChaDaughterQaBinning); + + sigmaHistManager.init(&hRegistry, sigmaHistSpec, chaDauHistSpec); + + auto collisionHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); + colHistManager.init(&hRegistry, collisionHistSpec); + }; + + // Process function for sigma particles from femto tables + void processSigma(FilteredCollision const& col, Sigmas const& /*sigmas*/, Tracks const& tracks) + { + colHistManager.fill(col); + auto sigmaSlice = sigmaPartition->sliceByCached(femtobase::stored::collisionId, col.globalIndex(), cache); + for (auto const& sigma : sigmaSlice) { + sigmaHistManager.fill(sigma, tracks); + } + } + PROCESS_SWITCH(FemtoKinkQa, processSigma, "Process sigmas", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc), + }; + return workflow; +}