diff --git a/PWGLF/DataModel/LFResonanceTables.h b/PWGLF/DataModel/LFResonanceTables.h index 79ea6bc9821..1de673d77c4 100644 --- a/PWGLF/DataModel/LFResonanceTables.h +++ b/PWGLF/DataModel/LFResonanceTables.h @@ -34,6 +34,7 @@ #include +#include #include #include #include @@ -45,7 +46,7 @@ namespace resocollision { enum { kECbegin = 0, - kINEL = 1, + kINEL, kINEL10, kINELg0, kINELg010, @@ -79,7 +80,7 @@ DECLARE_SOA_COLUMN(IsTriggerTVX, isTriggerTVX, bool); //! TriggerTVX DECLARE_SOA_COLUMN(IsInSel8, isInSel8, bool); //! InSel8 DECLARE_SOA_COLUMN(IsInAfterAllCuts, isInAfterAllCuts, bool); //! InAfterAllCuts DECLARE_SOA_COLUMN(ImpactParameter, impactParameter, float); //! ImpactParameter -DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity +DECLARE_SOA_COLUMN(MCMultiplicity, mcMultiplicity, float); //! MC Multiplicity, o2-linter: disable=name/o2-column (pre-existing public column name kept for schema and API compatibility) } // namespace resocollision DECLARE_SOA_TABLE(ResoCollisions, "AOD", "RESOCOLLISION", @@ -163,32 +164,47 @@ struct ResoTrackFlags { return (flags & mask) == mask; } }; +// These macros build framework expression nodes and cannot be replaced by +// ordinary constexpr functions without changing their public DSL API. +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requireTrackFlag(mask) ((o2::aod::resodaughter::trackFlags & o2::aod::resodaughter::mask) == o2::aod::resodaughter::mask) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requirePassedITSRefit() requireTrackFlag(ResoTrackFlags::kPassedITSRefit) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requirePassedTPCRefit() requireTrackFlag(ResoTrackFlags::kPassedTPCRefit) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requireGlobalTrack() requireTrackFlag(ResoTrackFlags::kIsGlobalTrack) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requireGlobalTrackWoDCA() requireTrackFlag(ResoTrackFlags::kIsGlobalTrackWoDCA) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requirePrimaryTrack() requireTrackFlag(ResoTrackFlags::kIsPrimaryTrack) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requirePVContributor() requireTrackFlag(ResoTrackFlags::kIsPVContributor) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requireHasTOF() requireTrackFlag(ResoTrackFlags::kHasTOF) +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define requireSign() requireTrackFlag(ResoTrackFlags::kSign) -#define DECLARE_DYN_TRKSEL_COLUMN(name, getter, mask) \ - DECLARE_SOA_DYNAMIC_COLUMN(name, getter, [](ResoTrackFlags::flagtype flags) -> bool { return ResoTrackFlags::checkFlag(flags, mask); }); +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define DECLARE_DYN_TRKSEL_COLUMN(_Name_, _Getter_, _Mask_) \ + DECLARE_SOA_DYNAMIC_COLUMN(_Name_, _Getter_, [](ResoTrackFlags::flagtype flags) -> bool { return ResoTrackFlags::checkFlag(flags, _Mask_); }); DECLARE_SOA_INDEX_COLUMN(ResoCollision, resoCollision); DECLARE_SOA_INDEX_COLUMN(ResoCollisionDF, resoCollisionDF); -DECLARE_SOA_INDEX_COLUMN_FULL(Track, track, int, Tracks, "_Trk"); //! -DECLARE_SOA_INDEX_COLUMN_FULL(V0, v0, int, V0s, "_V0"); //! -DECLARE_SOA_INDEX_COLUMN_FULL(Cascade, cascade, int, Cascades, "_Cas"); //! -DECLARE_SOA_COLUMN(Pt, pt, float); //! p_t (GeV/c) -DECLARE_SOA_COLUMN(Px, px, float); //! p_x (GeV/c) -DECLARE_SOA_COLUMN(Py, py, float); //! p_y (GeV/c) -DECLARE_SOA_COLUMN(Pz, pz, float); //! p_z (GeV/c) -DECLARE_SOA_COLUMN(PartType, partType, uint8_t); //! Type of the particle, according to resodaughter::ParticleType -DECLARE_SOA_COLUMN(TempFitVar, tempFitVar, float); //! Observable for the template fitting (Track: DCA_xy, V0: CPA) -DECLARE_SOA_COLUMN(Indices, indices, int[2]); //! Field for the track indices to remove auto-correlations +DECLARE_SOA_INDEX_COLUMN_FULL(Track, track, int, Tracks, "_Trk"); //! Soft link to the original track +DECLARE_SOA_INDEX_COLUMN_FULL(V0, v0, int, V0s, "_V0"); //! Soft link to the original V0 +DECLARE_SOA_INDEX_COLUMN_FULL(Cascade, cascade, int, Cascades, "_Cas"); //! Soft link to the original cascade +DECLARE_SOA_COLUMN(Pt, pt, float); //! p_t (GeV/c) +DECLARE_SOA_COLUMN(Px, px, float); //! p_x (GeV/c) +DECLARE_SOA_COLUMN(Py, py, float); //! p_y (GeV/c) +DECLARE_SOA_COLUMN(Pz, pz, float); //! p_z (GeV/c) +DECLARE_SOA_COLUMN(PartType, partType, uint8_t); //! Type of the particle, according to resodaughter::ParticleType +DECLARE_SOA_COLUMN(TempFitVar, tempFitVar, float); //! Observable for the template fitting (Track: DCA_xy, V0: CPA) +// Fixed-size C arrays are part of the existing persistent AOD schema. +// NOLINTNEXTLINE(modernize-avoid-c-arrays) +DECLARE_SOA_COLUMN(Indices, indices, int[2]); //! Field for the track indices to remove auto-correlations +// NOLINTNEXTLINE(modernize-avoid-c-arrays) DECLARE_SOA_COLUMN(CascadeIndices, cascadeIndices, int[3]); //! Field for the track indices to remove auto-correlations (ordered: positive, negative, bachelor) DECLARE_SOA_COLUMN(TpcNClsCrossedRows, tpcNClsCrossedRows, uint8_t); //! Number of TPC crossed rows DECLARE_SOA_COLUMN(TpcNClsFound, tpcNClsFound, uint8_t); //! Number of TPC clusters found @@ -215,7 +231,7 @@ DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float); DECLARE_SOA_COLUMN(DecayVtxY, decayVtxY, float); //! Y position of the decay vertex DECLARE_SOA_COLUMN(DecayVtxZ, decayVtxZ, float); //! Z position of the decay vertex DECLARE_SOA_COLUMN(Alpha, alpha, float); //! Alpha of the decay vertex -DECLARE_SOA_COLUMN(QtArm, qtarm, float); //! Armenteros Qt of the decay vertex +DECLARE_SOA_COLUMN(QtArm, qtarm, float); //! Armenteros Qt of the decay vertex, o2-linter: disable=name/o2-column (pre-existing public column name kept for schema and API compatibility) DECLARE_SOA_COLUMN(TpcSignal10, tpcSignal10, int16_t); //! TPC signal of the track x10 DECLARE_SOA_COLUMN(DaughterTPCNSigmaPosPi10, daughterTPCNSigmaPosPi10, int8_t); //! TPC PID x10 of the positive daughter as Pion DECLARE_SOA_COLUMN(DaughterTPCNSigmaPosKa10, daughterTPCNSigmaPosKa10, int8_t); //! TPC PID x10 of the positive daughter as Kaon @@ -250,9 +266,10 @@ DECLARE_SOA_COLUMN(DaughterPDG1, daughterPDG1, int); //! PDG code of the first D DECLARE_SOA_COLUMN(DaughterPDG2, daughterPDG2, int); //! PDG code of the second Daughter particle DECLARE_SOA_COLUMN(DaughterID1, daughterID1, int); //! Id of the first Daughter particle DECLARE_SOA_COLUMN(DaughterID2, daughterID2, int); //! Id of the second Daughter particle -DECLARE_SOA_COLUMN(SiblingIds, siblingIds, int[2]); //! Index of the particles with the same mother -DECLARE_SOA_COLUMN(BachTrkID, bachTrkID, int); //! Id of the bach track from cascade -DECLARE_SOA_COLUMN(V0ID, v0ID, int); //! Id of the V0 from cascade +// NOLINTNEXTLINE(modernize-avoid-c-arrays) -- persistent fixed-size AOD column +DECLARE_SOA_COLUMN(SiblingIds, siblingIds, int[2]); //! Index of the particles with the same mother +DECLARE_SOA_COLUMN(BachTrkID, bachTrkID, int); //! Id of the bach track from cascade +DECLARE_SOA_COLUMN(V0ID, v0ID, int); //! Id of the V0 from cascade // Dynamic columns // DCA_xy x10,000 DECLARE_SOA_DYNAMIC_COLUMN(DcaXY, dcaXY, @@ -358,6 +375,10 @@ namespace resomicrodaughter /// @brief Save TPC & TOF nSigma info with 8-bit variable struct PidNSigma { + static constexpr double MinNSigma = 1.5; + static constexpr double NSigmaStep = 0.2; + static constexpr uint8_t MaxNSigmaCode = 10; + uint8_t flag; /// @brief Constructor: Convert TPC & TOF values and save @@ -372,29 +393,34 @@ struct PidNSigma { static uint8_t encodeNSigma(float nSigma) { const float x = std::abs(nSigma); - if (x <= 1.5) + if (x <= MinNSigma) { return 0; // Return 0 when absolute nSigma is smaller than 1.5 - float t = (x - 1.5) / 0.2; + } + float t = (x - MinNSigma) / NSigmaStep; int encoded = static_cast(std::ceil(t)); // (1.5,1.7]->1, ..., (3.3,3.5]->10 - if (encoded < 1) + if (encoded < 1) { encoded = 1; - if (encoded > 10) - encoded = 10; + } + if (encoded > MaxNSigmaCode) { + encoded = MaxNSigmaCode; + } return static_cast(encoded); } /// @brief Decode 0~10 value to original 1.5~3.5 sigma range static float decodeNSigma(uint8_t encoded) { - if (encoded == 0) - return 1.5; - if (encoded > 10) - encoded = 10; - return 1.5 + static_cast(encoded) * 0.2; + if (encoded == 0) { + return MinNSigma; + } + if (encoded > MaxNSigmaCode) { + encoded = MaxNSigmaCode; + } + return MinNSigma + static_cast(encoded) * NSigmaStep; } /// @brief Check if TOF info is available - bool hasTOF() const + [[nodiscard]] bool hasTOF() const { return (flag & 0x0F) != 0x0F; // Check if lower 4 bits are not all 1 } @@ -430,12 +456,16 @@ DECLARE_SOA_DYNAMIC_COLUMN(HasTOF, hasTOF, /// @brief DCAxy & DCAz selection flag struct ResoMicroTrackSelFlag { + static constexpr double DCAEncodingStep = 0.1; + static constexpr uint8_t MaxRegularDCAFlag = 14; + static constexpr uint8_t OverflowDCAFlag = 15; + uint8_t flag; // Flag for DCAxy & DCAz selection (8-bit variable) /// @brief Default constructor ResoMicroTrackSelFlag() + : flag(0x00) { - flag = 0x00; } /// @brief Constructor: Convert DCAxy/DCAz and save (default 1~15 values) @@ -450,13 +480,16 @@ struct ResoMicroTrackSelFlag { static uint8_t encodeDCA(float DCA) { float x = std::fabs(DCA); - if (x < 0.1) + if (x < DCAEncodingStep) { return 0; - int encoded = static_cast(std::ceil((x - 0.1) / 0.1)); // (0.1, 0.2] -> 1, ..., (1.4, 1.5] -> 14 - if (encoded < 1) + } + int encoded = static_cast(std::ceil((x - DCAEncodingStep) / DCAEncodingStep)); // (0.1, 0.2] -> 1, ..., (1.4, 1.5] -> 14 + if (encoded < 1) { encoded = 1; - if (encoded > 14) - encoded = 15; + } + if (encoded > MaxRegularDCAFlag) { + encoded = OverflowDCAFlag; + } return static_cast(encoded); } @@ -467,13 +500,13 @@ struct ResoMicroTrackSelFlag { } /// @brief Get DCAxy value - uint8_t getDCAxyFlag() const + [[nodiscard]] uint8_t getDCAxyFlag() const { return (flag >> 4) & 0x0F; // Extract upper 4 bits } /// @brief Get DCAz value - uint8_t getDCAzFlag() const + [[nodiscard]] uint8_t getDCAzFlag() const { return flag & 0x0F; // Extract lower 4 bits } @@ -507,6 +540,167 @@ struct ResoMicroTrackSelFlag { DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float px, float py) -> float { return RecoDecay::sqrtSumOfSquares(px, py); }); } // namespace resomicrodaughter +// Ultra-micro track representation. The momentum components are quantised +// to 1 MeV/c and only one (pion/kaon/proton) PID flag is retained. +namespace resoultramicrodaughter +{ +/// @brief Compact absolute TPC/TOF n-sigma values for one configured species. +/// Code 0 represents |n-sigma| < 1.5. Codes 1..14 represent lower-inclusive +/// 0.25-sigma-wide bins [1.5, 1.75), ..., [4.75, 5.0) and decode to their +/// lower edges. Code 15 represents |n-sigma| >= 5. Missing TOF information +/// is carried independently by resodaughter::TrackFlags::kHasTOF. +struct PidNSigma { + static constexpr float MinFineNSigma = 1.5f; + static constexpr float Step = 0.25f; + static constexpr float MaxNSigma = 5.f; + static constexpr uint8_t MaxRegularCode = 14; + static constexpr uint8_t AboveRangeCode = 15; + + uint8_t flag; + + PidNSigma(float tpcNSigma, float tofNSigma, bool hasTOF) + { + const uint8_t tpcEncoded = encodeNSigma(tpcNSigma); + const uint8_t tofEncoded = hasTOF && std::isfinite(tofNSigma) ? encodeNSigma(tofNSigma) : AboveRangeCode; + flag = (tpcEncoded << 4) | tofEncoded; + } + + static uint8_t encodeNSigma(float nSigma) + { + const float value = std::abs(nSigma); + if (!std::isfinite(value)) { + return AboveRangeCode; + } + if (value < MinFineNSigma) { + return 0; + } + if (value >= MaxNSigma) { + return AboveRangeCode; + } + const int encoded = 1 + static_cast(std::floor((value - MinFineNSigma) / Step)); + return static_cast(std::clamp(encoded, 1, static_cast(MaxRegularCode))); + } + + static float decodeNSigma(uint8_t encoded) + { + const uint8_t code = encoded & 0x0F; + if (code == 0) { + return 0.f; + } + return code == AboveRangeCode ? MaxNSigma : MinFineNSigma + static_cast(code - 1) * Step; + } + + static float getTPCNSigma(uint8_t encoded) + { + return decodeNSigma((encoded >> 4) & 0x0F); + } + + static float getTOFNSigma(uint8_t encoded, bool hasTOF) + { + return hasTOF ? decodeNSigma(encoded & 0x0F) : NAN; + } + + operator uint8_t() const { return flag; } +}; + +/// @brief Compact absolute DCAxy/DCAz values into two four-bit fields. +/// Codes 0..14 represent lower-inclusive 0.01 cm-wide bins [0.00, 0.01), +/// ..., [0.14, 0.15) cm and decode to their lower edges. Code 15 directly +/// represents |DCA| >= 0.15 cm and decodes to the 0.15 cm marker; no additional +/// tight-DCA or overflow flag is stored. +struct DCAEncoding { + static constexpr float MaxDCA = 0.15f; + static constexpr float Step = 0.01f; + static constexpr uint8_t MaxRegularCode = 14; + static constexpr uint8_t AboveRangeCode = 15; + + uint8_t flag = 0; + + DCAEncoding() = default; + DCAEncoding(float dcaXY, float dcaZ) + : flag(static_cast((encodeDCA(dcaXY) << 4) | encodeDCA(dcaZ))) + { + } + + static bool isValid(float dca) + { + return std::isfinite(dca); + } + + static uint8_t encodeDCA(float dca) + { + const float value = std::abs(dca); + if (!std::isfinite(value)) { + return AboveRangeCode; + } + if (value >= MaxDCA) { + return AboveRangeCode; + } + const int encoded = static_cast(std::floor(value / Step)); + return static_cast(std::clamp(encoded, 0, static_cast(MaxRegularCode))); + } + + static float decodeDCA(uint8_t encoded) + { + const uint8_t code = encoded & 0x0F; + return code == AboveRangeCode ? MaxDCA : static_cast(code) * Step; + } + + static float decodeDCAxy(uint8_t encoded) + { + return decodeDCA((encoded >> 4) & 0x0F); + } + + static float decodeDCAz(uint8_t encoded) + { + return decodeDCA(encoded & 0x0F); + } + + operator uint8_t() const { return flag; } +}; + +DECLARE_SOA_COLUMN(PidNSigmaFlag, pidNSigmaFlag, uint8_t); //! TPC/TOF PID flag for the configured species +DECLARE_SOA_COLUMN(TrackSelectionFlags, trackSelectionFlags, uint8_t); //! Packed absolute DCAxy/DCAz values +DECLARE_SOA_COLUMN(Px1000, px1000, int16_t); //! p_x x 1000 (GeV/c) +DECLARE_SOA_COLUMN(Py1000, py1000, int16_t); //! p_y x 1000 (GeV/c) +DECLARE_SOA_COLUMN(Pz1000, pz1000, int16_t); //! p_z x 1000 (GeV/c) + +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, + [](int16_t px1000) { return static_cast(px1000) / 1000.f; }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, + [](int16_t py1000) { return static_cast(py1000) / 1000.f; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, + [](int16_t pz1000) { return static_cast(pz1000) / 1000.f; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, + [](int16_t px1000, int16_t py1000) { + const float px = static_cast(px1000) / 1000.f; + const float py = static_cast(py1000) / 1000.f; + return RecoDecay::sqrtSumOfSquares(px, py); + }); +DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, + [](int16_t px1000, int16_t py1000, int16_t pz1000) { + return RecoDecay::eta(std::array{static_cast(px1000) / 1000.f, + static_cast(py1000) / 1000.f, + static_cast(pz1000) / 1000.f}); + }); +DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, + [](int16_t px1000, int16_t py1000) { + return RecoDecay::phi(static_cast(px1000) / 1000.f, + static_cast(py1000) / 1000.f); + }); +DECLARE_SOA_DYNAMIC_COLUMN(TpcNSigma, tpcNSigma, + [](uint8_t pidNSigmaFlag) { return PidNSigma::getTPCNSigma(pidNSigmaFlag); }); +DECLARE_SOA_DYNAMIC_COLUMN(TofNSigma, tofNSigma, + [](uint8_t pidNSigmaFlag, uint8_t trackFlags) -> float { + const bool hasTOF = resodaughter::ResoTrackFlags::checkFlag(trackFlags, resodaughter::ResoTrackFlags::kHasTOF); + return PidNSigma::getTOFNSigma(pidNSigmaFlag, hasTOF); + }); +DECLARE_SOA_DYNAMIC_COLUMN(DcaXY, dcaXY, + [](uint8_t trackSelectionFlags) { return DCAEncoding::decodeDCAxy(trackSelectionFlags); }); +DECLARE_SOA_DYNAMIC_COLUMN(DcaZ, dcaZ, + [](uint8_t trackSelectionFlags) { return DCAEncoding::decodeDCAz(trackSelectionFlags); }); +} // namespace resoultramicrodaughter + DECLARE_SOA_TABLE(ResoTracks, "AOD", "RESOTRACK", o2::soa::Index<>, resodaughter::ResoCollisionId, @@ -582,6 +776,40 @@ DECLARE_SOA_TABLE(ResoMicroTrackTracks, "AOD", "RESOMICROTRACKTRACK", resodaughter::TrackId); using ResoMicroTrackTrack = ResoMicroTrackTracks::iterator; +DECLARE_SOA_TABLE(ResoUltraMicroTracks, "AOD", "RESOULTRAMTRK", + o2::soa::Index<>, + resodaughter::ResoCollisionId, + resoultramicrodaughter::Px1000, + resoultramicrodaughter::Py1000, + resoultramicrodaughter::Pz1000, + resoultramicrodaughter::PidNSigmaFlag, + resoultramicrodaughter::TrackSelectionFlags, + resodaughter::TrackFlags, + // Dynamic columns + resoultramicrodaughter::Px, + resoultramicrodaughter::Py, + resoultramicrodaughter::Pz, + resoultramicrodaughter::Pt, + resoultramicrodaughter::Eta, + resoultramicrodaughter::Phi, + resoultramicrodaughter::TpcNSigma, + resoultramicrodaughter::TofNSigma, + resoultramicrodaughter::DcaXY, + resoultramicrodaughter::DcaZ, + resodaughter::PassedITSRefit, + resodaughter::PassedTPCRefit, + resodaughter::IsGlobalTrackWoDCA, + resodaughter::IsGlobalTrack, + resodaughter::IsPrimaryTrack, + resodaughter::IsPVContributor, + resodaughter::HasTOF, + resodaughter::Sign); +using ResoUltraMicroTrack = ResoUltraMicroTracks::iterator; + +DECLARE_SOA_TABLE(ResoUltraMicroTrackTracks, "AOD", "RESOULTRAMTRKID", + resodaughter::TrackId); +using ResoUltraMicroTrackTrack = ResoUltraMicroTrackTracks::iterator; + // For DF mixing study DECLARE_SOA_TABLE(ResoTrackDFs, "AOD", "RESOTRACKDF", o2::soa::Index<>, diff --git a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx index 17e1048c06e..500e6219c01 100644 --- a/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx +++ b/PWGLF/TableProducer/Resonances/resonanceModuleInitializer.cxx @@ -12,21 +12,24 @@ /// \file resonanceModuleInitializer.cxx /// \brief Initializes variables for the resonance candidate producers /// -/// \author Bong-Hwi Lim +/// \author Bong-Hwi Lim , Minjae Kim +/// \since Aug.18 2026 #include "PWGLF/DataModel/LFResonanceTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/mcCentrality.h" #include "PWGLF/Utils/collisionCuts.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" -#include "Common/Core/EventPlaneHelper.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/Qvectors.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" #include #include +#include #include #include #include @@ -41,9 +44,12 @@ #include #include +#include #include #include +#include #include +#include #include #include @@ -62,77 +68,76 @@ using namespace o2::aod::rctsel; * for resonance studies. It handles event selection, centrality estimation, * and QA histogram filling. */ +// Framework Service members are populated by the analysis framework. +// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) struct ResonanceModuleInitializer { - int mRunNumber; ///< Run number for the current data - int multEstimator; ///< Multiplicity estimator type - float dBz; ///< Magnetic field value - float centrality; ///< Centrality value for the event + static constexpr double BzOverrideThreshold = -990.; + static constexpr double MinimumNonzeroBz = 1.e-5; + static constexpr double MinimumChargedParticleCharge = 3.; // ROOT particle charge is stored in units of e/3 + static constexpr int MCCentralityRecoEstimator = 0; + static constexpr int MCCentralityGeneratorEstimator = 1; + static constexpr int MCCentralityImpactParameterEstimator = 2; + static constexpr float MCVertexZMax = 10.f; + + int mRunNumber = 0; ///< Run number for the current data + int multEstimator = 0; ///< Multiplicity estimator type + float dBz = 0.f; ///< Magnetic field value + float centrality = 0.f; ///< Centrality value for the event Service ccdb; ///< CCDB manager service Service pdg; ///< PDG database service - EventPlaneHelper helperEP; ///< Helper for event plane calculations - Produces resoCollisions; ///< Output table for resonance collisions - Produces resoCollisionColls; ///< Output table for collision references - Produces resoMCCollisions; ///< Output table for MC resonance collisions - Produces resoSpheroCollisions; ///< Output table for spherocity - Produces resoEvtPlCollisions; ///< Output table for event plane + Produces resoCollisions; ///< Output table for resonance collisions + Produces resoCollisionColls; ///< Output table for collision references + Produces resoMCCollisions; ///< Output table for MC resonance collisions // CCDB options - Configurable ccdbURL{"ccdbURL", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - - Configurable cfgFatalWhenNull{"cfgFatalWhenNull", true, "Fatal when null on ccdb access"}; - Configurable cfgBypassCollIndexFill{"cfgBypassCollIndexFill", false, "Bypass collision index fill"}; - - // Configurables - Configurable dBzInput{"dBzInput", -999, "bz field, -999 is automatic"}; - Configurable cfgFillQA{"cfgFillQA", false, "Fill QA histograms"}; - Configurable cfgBypassCCDB{"cfgBypassCCDB", true, "Bypass loading CCDB part to save CPU time and memory"}; // will be affected to b_z value. - Configurable cfgMultName{"cfgMultName", "FT0M", "The name of multiplicity estimator"}; - Configurable cfgCentralityMC{"cfgCentralityMC", 0, "Centrality estimator for MC (0: Reco, 1: MC, 2: impact parameter)"}; - - // EventCorrection for MC - ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 0.01, 0.1, 1.0, 5.0, 10., 15., 20., 30., 40., 50., 70., 100.0, 105.}, "Binning of the centrality axis"}; - ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -20, -15, -10, -7, -5, -3, -2, -1, 0, 1, 2, 3, 5, 7, 10, 15, 20}, "Mixing bins - z-vertex"}; + struct : ConfigurableGroup { + Configurable ccdbURL{"ccdbURL", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable cfgFatalWhenNull{"cfgFatalWhenNull", true, "Fatal when null on ccdb access"}; + Configurable cfgBypassCollIndexFill{"cfgBypassCollIndexFill", false, "Unsupported in the modular workflow; must remain false"}; + } CCDB; + + // General event options + struct : ConfigurableGroup { + Configurable dBzInput{"dBzInput", -999, "bz field, -999 is automatic"}; + Configurable cfgFillQA{"cfgFillQA", true, "Fill QA histograms"}; + Configurable cfgBypassCCDB{"cfgBypassCCDB", true, "Bypass loading CCDB part to save CPU time and memory"}; // will be affected to b_z value. + Configurable cfgMultName{"cfgMultName", "FT0M", "The name of multiplicity estimator"}; + Configurable cfgCentralityMC{"cfgCentralityMC", 0, "Centrality estimator for MC (0: Reco, 1: MC, 2: impact parameter)"}; + ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 0.01, 0.1, 1.0, 5.0, 10., 15., 20., 30., 40., 50., 70., 100.0, 105.}, "Binning of the centrality axis"}; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -20, -15, -10, -7, -5, -3, -2, -1, 0, 1, 2, 3, 5, 7, 10, 15, 20}, "Mixing bins - z-vertex"}; + } EventConfig; /// Event cuts o2::analysis::CollisonCuts colCuts; - Configurable cfgEvtZvtx{"cfgEvtZvtx", 10.f, "Evt sel: Max. z-Vertex (cm)"}; - Configurable cfgEvtOccupancyInTimeRange{"cfgEvtOccupancyInTimeRange", -1, "Evt sel: maximum track occupancy"}; - Configurable cfgEvtTriggerCheck{"cfgEvtTriggerCheck", false, "Evt sel: check for trigger"}; - Configurable cfgEvtOfflineCheck{"cfgEvtOfflineCheck", true, "Evt sel: check for offline selection"}; - Configurable cfgEvtTriggerTVXSel{"cfgEvtTriggerTVXSel", false, "Evt sel: triggerTVX selection (MB)"}; - Configurable cfgEvtTFBorderCut{"cfgEvtTFBorderCut", false, "Evt sel: apply TF border cut"}; - Configurable cfgEvtUseITSTPCvertex{"cfgEvtUseITSTPCvertex", false, "Evt sel: use at lease on ITS-TPC track for vertexing"}; - Configurable cfgEvtCollInTimeRangeNarrow{"cfgEvtCollInTimeRangeNarrow", false, "Evt sel: apply NoCollInTimeRangeNarrow"}; - Configurable cfgEvtZvertexTimedifference{"cfgEvtZvertexTimedifference", false, "Evt sel: apply Z-vertex time difference"}; - Configurable cfgEvtPileupRejection{"cfgEvtPileupRejection", false, "Evt sel: apply pileup rejection"}; - Configurable cfgEvtNoITSROBorderCut{"cfgEvtNoITSROBorderCut", false, "Evt sel: apply NoITSRO border cut"}; - Configurable cfgEvtRun2AliEventCuts{"cfgEvtRun2AliEventCuts", true, "Evt sel: apply Run2 AliEventCuts"}; - Configurable cfgEvtRun2INELgtZERO{"cfgEvtRun2INELgtZERO", false, "Evt sel: apply Run2 INELgtZERO"}; - Configurable cfgEvtUseRCTFlagChecker{"cfgEvtUseRCTFlagChecker", false, "Evt sel: use RCT flag checker"}; - Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; - Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; - Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", false, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; + struct : ConfigurableGroup { + Configurable cfgEvtZvtx{"cfgEvtZvtx", 10.f, "Evt sel: Max. z-Vertex (cm)"}; + Configurable cfgEvtOccupancyInTimeRange{"cfgEvtOccupancyInTimeRange", -1, "Evt sel: maximum track occupancy"}; + Configurable cfgEvtTriggerCheck{"cfgEvtTriggerCheck", false, "Evt sel: check for trigger"}; + Configurable cfgEvtOfflineCheck{"cfgEvtOfflineCheck", true, "Evt sel: check for offline selection"}; + Configurable cfgEvtTriggerTVXSel{"cfgEvtTriggerTVXSel", false, "Evt sel: triggerTVX selection (MB)"}; + Configurable cfgEvtTFBorderCut{"cfgEvtTFBorderCut", false, "Evt sel: apply TF border cut"}; + Configurable cfgEvtUseITSTPCvertex{"cfgEvtUseITSTPCvertex", false, "Evt sel: use at lease on ITS-TPC track for vertexing"}; + Configurable cfgEvtCollInTimeRangeNarrow{"cfgEvtCollInTimeRangeNarrow", false, "Evt sel: apply NoCollInTimeRangeNarrow"}; + Configurable cfgEvtZvertexTimedifference{"cfgEvtZvertexTimedifference", false, "Evt sel: apply Z-vertex time difference"}; + Configurable cfgEvtPileupRejection{"cfgEvtPileupRejection", false, "Evt sel: apply pileup rejection"}; + Configurable cfgEvtNoITSROBorderCut{"cfgEvtNoITSROBorderCut", false, "Evt sel: apply NoITSRO border cut"}; + Configurable cfgEvtRun2AliEventCuts{"cfgEvtRun2AliEventCuts", true, "Evt sel: apply Run2 AliEventCuts"}; + Configurable cfgEvtRun2INELgtZERO{"cfgEvtRun2INELgtZERO", false, "Evt sel: apply Run2 INELgtZERO"}; + Configurable cfgEvtUseRCTFlagChecker{"cfgEvtUseRCTFlagChecker", false, "Evt sel: use RCT flag checker"}; + Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; + Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; + Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", false, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; + } EventCuts; RCTFlagsChecker rctChecker; - // Spherocity configuration - Configurable cfgTrackSphMin{"cfgTrackSphMin", 10, "Number of tracks for Spherocity Calculation"}; - Configurable cfgTrackSphDef{"cfgTrackSphDef", 0, "Spherocity Definition: |pT| = 1 -> 0, otherwise -> 1"}; - - // Qvector configuration - Configurable cfgEvtPl{"cfgEvtPl", 40500, "Configuration of three subsystems for the event plane and its resolution, 10000*RefA + 100*RefB + S, where FT0C:0, FT0A:1, FT0M:2, FV0A:3, BPos:5, BNeg:6"}; - - int evtPlRefAId = static_cast(cfgEvtPl / 10000); - int evtPlRefBId = static_cast((cfgEvtPl - evtPlRefAId * 10000) / 100); - int evtPlDetId = cfgEvtPl - evtPlRefAId * 10000 - evtPlRefBId * 100; - HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Filter collisionFilter = nabs(aod::collision::posZ) < cfgEvtZvtx; + Filter collisionFilter = nabs(aod::collision::posZ) < EventCuts.cfgEvtZvtx; /** * @brief Initializes the task @@ -146,14 +151,14 @@ struct ResonanceModuleInitializer { centrality = 0; // Determine the multiplicity estimator based on the configuration multEstimator = 0; - if (cfgMultName.value == "FT0M") { + if (EventConfig.cfgMultName.value == "FT0M") { multEstimator = 0; - } else if (cfgMultName.value == "FT0C") { + } else if (EventConfig.cfgMultName.value == "FT0C") { multEstimator = 1; - } else if (cfgMultName.value == "FT0A") { + } else if (EventConfig.cfgMultName.value == "FT0A") { multEstimator = 2; } - LOGF(info, "Mult estimator: %d, %s", multEstimator, cfgMultName.value.c_str()); + LOGF(info, "Mult estimator: %d, %s", multEstimator, EventConfig.cfgMultName.value.c_str()); // Ensure that only one process type is active at a time if (doprocessRun3 && doprocessRun2) { @@ -162,39 +167,42 @@ struct ResonanceModuleInitializer { if (doprocessRun2MC && doprocessRun3MC) { LOG(fatal) << "You cannot run both Run2 and Run3 MC processes at the same time"; } + if (CCDB.cfgBypassCollIndexFill) { + LOG(fatal) << "cfgBypassCollIndexFill is incompatible with ResonanceDaughterInitializer"; + } // Initialize event selection cuts based on the process type if (doprocessRun2) { - colCuts.setCuts(cfgEvtZvtx, cfgEvtTriggerCheck, cfgEvtOfflineCheck, false); + colCuts.setCuts(EventCuts.cfgEvtZvtx, EventCuts.cfgEvtTriggerCheck, EventCuts.cfgEvtOfflineCheck, false); } else if (doprocessRun3) { - colCuts.setCuts(cfgEvtZvtx, cfgEvtTriggerCheck, cfgEvtOfflineCheck, true, false, cfgEvtOccupancyInTimeRange); + colCuts.setCuts(EventCuts.cfgEvtZvtx, EventCuts.cfgEvtTriggerCheck, EventCuts.cfgEvtOfflineCheck, true, false, EventCuts.cfgEvtOccupancyInTimeRange); } colCuts.init(&qaRegistry); - colCuts.setTriggerTVX(cfgEvtTriggerTVXSel); - colCuts.setApplyTFBorderCut(cfgEvtTFBorderCut); - colCuts.setApplyITSTPCvertex(cfgEvtUseITSTPCvertex); - colCuts.setApplyCollInTimeRangeNarrow(cfgEvtCollInTimeRangeNarrow); - colCuts.setApplyZvertexTimedifference(cfgEvtZvertexTimedifference); - colCuts.setApplyPileupRejection(cfgEvtPileupRejection); - colCuts.setApplyNoITSROBorderCut(cfgEvtNoITSROBorderCut); - colCuts.setApplyRun2AliEventCuts(cfgEvtRun2AliEventCuts); - colCuts.setApplyRun2INELgtZERO(cfgEvtRun2INELgtZERO); - - rctChecker.init(cfgEvtRCTFlagCheckerLabel, cfgEvtRCTFlagCheckerZDCCheck, cfgEvtRCTFlagCheckerLimitAcceptAsBad); + colCuts.setTriggerTVX(EventCuts.cfgEvtTriggerTVXSel); + colCuts.setApplyTFBorderCut(EventCuts.cfgEvtTFBorderCut); + colCuts.setApplyITSTPCvertex(EventCuts.cfgEvtUseITSTPCvertex); + colCuts.setApplyCollInTimeRangeNarrow(EventCuts.cfgEvtCollInTimeRangeNarrow); + colCuts.setApplyZvertexTimedifference(EventCuts.cfgEvtZvertexTimedifference); + colCuts.setApplyPileupRejection(EventCuts.cfgEvtPileupRejection); + colCuts.setApplyNoITSROBorderCut(EventCuts.cfgEvtNoITSROBorderCut); + colCuts.setApplyRun2AliEventCuts(EventCuts.cfgEvtRun2AliEventCuts); + colCuts.setApplyRun2INELgtZERO(EventCuts.cfgEvtRun2INELgtZERO); + + rctChecker.init(EventCuts.cfgEvtRCTFlagCheckerLabel, EventCuts.cfgEvtRCTFlagCheckerZDCCheck, EventCuts.cfgEvtRCTFlagCheckerLimitAcceptAsBad); // Configure CCDB access if not bypassed - if (!cfgBypassCCDB) { - ccdb->setURL(ccdbURL.value); + if (!EventConfig.cfgBypassCCDB) { + ccdb->setURL(CCDB.ccdbURL.value); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(cfgFatalWhenNull); + ccdb->setFatalWhenNull(CCDB.cfgFatalWhenNull); uint64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); ccdb->setCreatedNotAfter(now); // TODO must become global parameter from the train creation time } // Initialize QA histograms if required - if (doprocessRun3MC || doprocessRun2MC) { - AxisSpec centAxis = {binsCent, "Centrality (%)"}; + if (EventConfig.cfgFillQA && (doprocessRun3MC || doprocessRun2MC)) { + AxisSpec centAxis = {EventConfig.binsCent, "Centrality (%)"}; AxisSpec idxMCAxis = {26, -0.5, 25.5, "Index"}; qaRegistry.add("Event/hMCEventIndices", "hMCEventIndices", kTH2D, {centAxis, idxMCAxis}); } @@ -207,18 +215,18 @@ struct ResonanceModuleInitializer { */ void initCCDB(aod::BCsWithTimestamps::iterator const& bc) // Simple copy from LambdaKzeroFinder.cxx { - if (cfgBypassCCDB) + if (EventConfig.cfgBypassCCDB) { return; + } if (mRunNumber == bc.runNumber()) { return; } // In case override, don't proceed, please - no CCDB access required - if (dBzInput > -990) { - dBz = dBzInput; - ; + if (EventConfig.dBzInput > BzOverrideThreshold) { + dBz = EventConfig.dBzInput; o2::parameters::GRPMagField grpmag; - if (std::fabs(dBz) > 1e-5) { + if (std::fabs(dBz) > MinimumNonzeroBz) { grpmag.setL3Current(30000.f / (dBz / 5.0f)); } o2::base::Propagator::initFieldFromGRP(&grpmag); @@ -227,17 +235,17 @@ struct ResonanceModuleInitializer { } auto run3grpTimestamp = bc.timestamp(); - o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(grpPath, run3grpTimestamp); - o2::parameters::GRPMagField* grpmag = 0x0; + auto* grpo = ccdb->getForTimeStamp(CCDB.grpPath, run3grpTimestamp); + o2::parameters::GRPMagField* grpmag = nullptr; if (grpo) { o2::base::Propagator::initFieldFromGRP(grpo); // Fetch magnetic field from ccdb for current collision dBz = grpo->getNominalL3Field(); LOG(info) << "Retrieved GRP for timestamp " << run3grpTimestamp << " with magnetic field of " << dBz << " kZG"; } else { - grpmag = ccdb->getForTimeStamp(grpmagPath, run3grpTimestamp); + grpmag = ccdb->getForTimeStamp(CCDB.grpmagPath, run3grpTimestamp); if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grpTimestamp; + LOG(fatal) << "Got nullptr from CCDB for path " << CCDB.grpmagPath << " of object GRPMagField and " << CCDB.grpPath << " of object GRPObject for timestamp " << run3grpTimestamp; } o2::base::Propagator::initFieldFromGRP(grpmag); // Fetch magnetic field from ccdb for current collision @@ -260,13 +268,15 @@ struct ResonanceModuleInitializer { bool isTrueINEL0(MCPart const& mcparts) { for (auto const& mcparticle : mcparts) { - if (!mcparticle.isPhysicalPrimary()) + if (!mcparticle.isPhysicalPrimary()) { continue; + } auto p = pdg->GetParticle(mcparticle.pdgCode()); if (p != nullptr) { - if (std::abs(p->Charge()) >= 3) { - if (std::abs(mcparticle.eta()) < 1) + if (std::abs(p->Charge()) >= MinimumChargedParticleCharge) { + if (std::abs(mcparticle.eta()) < 1) { return true; + } } } } @@ -278,23 +288,23 @@ struct ResonanceModuleInitializer { * * @tparam ResoColl Type of resonance collision * @tparam isMC Boolean indicating if it's MC - * @param ResoEvents Resonance events + * @param resoEvents Resonance events * @return Centrality value */ template - float centEst(ResoColl ResoEvents) + float centEst(ResoColl const& resoEvents) { float returnValue = -999.0; switch (multEstimator) { case 0: - returnValue = ResoEvents.centFT0M(); + returnValue = resoEvents.centFT0M(); break; case 1: if constexpr (isMC) { LOG(fatal) << "CentFT0C is not available for MC"; return returnValue; } else { - returnValue = ResoEvents.centFT0C(); + returnValue = resoEvents.centFT0C(); break; } case 2: @@ -302,11 +312,11 @@ struct ResonanceModuleInitializer { LOG(fatal) << "CentFT0A is not available for MC"; return returnValue; } else { - returnValue = ResoEvents.centFT0A(); + returnValue = resoEvents.centFT0A(); break; } default: - returnValue = ResoEvents.centFT0M(); + returnValue = resoEvents.centFT0M(); break; } return returnValue; @@ -314,97 +324,6 @@ struct ResonanceModuleInitializer { using GenMCCollisions = soa::Join; float centEstMC(const GenMCCollisions::iterator& collision) { return centEst(collision); } - /** - * @brief Computes the spherocity of an event - * - * @tparam T Type of the tracks - * @param tracks All tracks - * @param nTracksMin Minimum number of tracks - * @param spdef Spherocity definition - * @return Spherocity value - */ - template - float computeSpherocity(T const& tracks, int nTracksMin, int spdef) - { - // if number of tracks is not enough for spherocity estimation. - int ntrks = tracks.size(); - if (ntrks < nTracksMin) - return -99.; - - // start computing spherocity - - float ptSum = 0.; - for (auto const& track : tracks) { - if (cfgFillQA) { - qaRegistry.fill(HIST("Phi"), track.phi()); - } - if (spdef == 0) { - ptSum += 1.; - } else { - ptSum += track.pt(); - } - } - - float tempSph = 1.; - for (int i = 0; i < 360 / 0.1; ++i) { - float sum = 0., pt = 0.; - float phiparm = (PI * i * 0.1) / 180.; - float nx = std::cos(phiparm); - float ny = std::sin(phiparm); - for (auto const& trk : tracks) { - pt = trk.pt(); - if (spdef == 0) { - pt = 1.; - } - float phi = trk.phi(); - float px = pt * std::cos(phi); - float py = pt * std::sin(phi); - // sum += pt * abs(sin(phiparm - phi)); - sum += std::abs(px * ny - py * nx); - } - float sph = std::pow((sum / ptSum), 2); - if (sph < tempSph) - tempSph = sph; - } - - return std::pow(PIHalf, 2) * tempSph; - } - - /** - * @brief Gets the event plane - * - * @tparam ResoColl Type of resonance collision - * @param ResoEvents Resonance events - * @return Event plane value - */ - template - float getEvtPl(ResoColl ResoEvents) - { - float returnValue = -999.0; - if (ResoEvents.qvecAmp()[evtPlDetId] > 1e-8) - returnValue = helperEP.GetEventPlane(ResoEvents.qvecRe()[evtPlDetId * 4 + 3], ResoEvents.qvecIm()[evtPlDetId * 4 + 3], 2); - return returnValue; - } - - /** - * @brief Gets the event plane resolution - * - * @tparam ResoColl Type of resonance collision - * @param ResoEvents Resonance events - * @param a First index - * @param b Second index - * @return Event plane resolution - */ - template - float getEvtPlRes(ResoColl ResoEvents, int a, int b) - { - float returnValue = -999.0; - if (ResoEvents.qvecAmp()[a] < 1e-8 || ResoEvents.qvecAmp()[b] < 1e-8) - return returnValue; - returnValue = helperEP.GetResolution(helperEP.GetEventPlane(ResoEvents.qvecRe()[a * 4 + 3], ResoEvents.qvecIm()[a * 4 + 3], 2), helperEP.GetEventPlane(ResoEvents.qvecRe()[b * 4 + 3], ResoEvents.qvecIm()[b * 4 + 3], 2), 2); - return returnValue; - } - /** * @brief Fills MC particles * @@ -459,65 +378,82 @@ struct ResonanceModuleInitializer { const auto& mcColg = mccol.template mcCollision_as(); float mcCent = 999.0; if constexpr (isRun2) { - if (cfgCentralityMC == 0) { + if (EventConfig.cfgCentralityMC == MCCentralityRecoEstimator) { mcCent = mccol.centRun2V0M(); } else { mcCent = mcColg.impactParameter(); } } else { - if (cfgCentralityMC == 0) { + if (EventConfig.cfgCentralityMC == MCCentralityRecoEstimator) { mcCent = centEst(mccol); - } else if (cfgCentralityMC == 1) { + } else if (EventConfig.cfgCentralityMC == MCCentralityGeneratorEstimator) { mcCent = centEstMC(mcColg); - } else if (cfgCentralityMC == 2) { + } else if (EventConfig.cfgCentralityMC == MCCentralityImpactParameterEstimator) { mcCent = mcColg.impactParameter(); } } - bool inVtx10 = (std::abs(mcColg.posZ()) > 10.) ? false : true; + const bool inVtx10 = !(std::abs(mcColg.posZ()) > MCVertexZMax); bool isTrueINELgt0 = isTrueINEL0(mcparts); bool isTriggerTVX = mccol.selection_bit(aod::evsel::kIsTriggerTVX); bool isSel8 = mccol.sel8(); - bool isSelected = colCuts.isSelected(mccol); + bool isSelected = colCuts.isSelected(mccol, EventConfig.cfgFillQA); resoMCCollisions(inVtx10, isTrueINELgt0, isTriggerTVX, isSel8, isSelected, mcCent, -1.0f); - // QA for Trigger efficiency - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL); - if (inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL10); - if (isTrueINELgt0) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINELg0); - if (inVtx10 && isTrueINELgt0) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINELg010); - - // TVX MB trigger - if (isTriggerTVX) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrig); - if (isTriggerTVX && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrig10); - if (isTriggerTVX && isTrueINELgt0) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrigINELg0); - if (isTriggerTVX && isTrueINELgt0 && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrigINELg010); - - // Sel8 event selection - if (isSel8) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8); - if (isSel8 && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel810); - if (isSel8 && isTrueINELgt0) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8INELg0); - if (isSel8 && isTrueINELgt0 && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8INELg010); - - // CollisionCuts selection - if (isSelected) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCuts); - if (isSelected && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCuts10); - if (isSelected && isTrueINELgt0) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCutsINELg0); - if (isSelected && isTrueINELgt0 && inVtx10) - qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCutsINELg010); + if (EventConfig.cfgFillQA) { + // QA for trigger efficiency + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL); + if (inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINEL10); + } + if (isTrueINELgt0) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINELg0); + } + if (inVtx10 && isTrueINELgt0) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kINELg010); + } + + // TVX MB trigger + if (isTriggerTVX) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrig); + } + if (isTriggerTVX && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrig10); + } + if (isTriggerTVX && isTrueINELgt0) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrigINELg0); + } + if (isTriggerTVX && isTrueINELgt0 && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kTrigINELg010); + } + + // Sel8 event selection + if (isSel8) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8); + } + if (isSel8 && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel810); + } + if (isSel8 && isTrueINELgt0) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8INELg0); + } + if (isSel8 && isTrueINELgt0 && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kSel8INELg010); + } + + // CollisionCuts selection + if (isSelected) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCuts); + } + if (isSelected && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCuts10); + } + if (isSelected && isTrueINELgt0) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCutsINELg0); + } + if (isSelected && isTrueINELgt0 && inVtx10) { + qaRegistry.fill(HIST("Event/hMCEventIndices"), mcCent, aod::resocollision::kAllCutsINELg010); + } + } } /** @@ -542,17 +478,20 @@ struct ResonanceModuleInitializer { auto bc = collision.bc_as(); initCCDB(bc); // Default event selection - if (!colCuts.isSelected(collision)) + if (!colCuts.isSelected(collision, EventConfig.cfgFillQA)) { return; - if (cfgEvtUseRCTFlagChecker && !rctChecker(collision)) + } + if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) { return; - colCuts.fillQA(collision); + } + if (EventConfig.cfgFillQA) { + colCuts.fillQA(collision); + } + const bool isRecINELgt0 = collision.isInelGt0(); centrality = centEst(collision); - resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); - if (!cfgBypassCollIndexFill) { - resoCollisionColls(collision.globalIndex()); - } + resoCollisions(collision.multNTracksPV(), collision.multNTracksPVeta1(), collision.multNTracksPVetaHalf(), collision.posX(), collision.posY(), collision.posZ(), centEst(collision), dBz, isRecINELgt0); + resoCollisionColls(collision.globalIndex()); } PROCESS_SWITCH(ResonanceModuleInitializer, processRun3, "Default process for RUN3", false); @@ -567,15 +506,16 @@ struct ResonanceModuleInitializer { { // auto bc = collision.bc_as(); // Default event selection - if (!colCuts.isSelected(collision)) + if (!colCuts.isSelected(collision, EventConfig.cfgFillQA)) { return; - colCuts.fillQARun2(collision); + } + if (EventConfig.cfgFillQA) { + colCuts.fillQARun2(collision); + } centrality = collision.centRun2V0M(); resoCollisions(0, 0, 0, collision.posX(), collision.posY(), collision.posZ(), centrality, dBz, 0); - if (!cfgBypassCollIndexFill) { - resoCollisionColls(collision.globalIndex()); - } + resoCollisionColls(collision.globalIndex()); } PROCESS_SWITCH(ResonanceModuleInitializer, processRun2, "process for RUN2", false); @@ -589,8 +529,9 @@ struct ResonanceModuleInitializer { void processRun3MC(soa::Filtered::iterator const& collision, aod::McParticles const& mcParticles, GenMCCollisions const&) { - if (cfgEvtUseRCTFlagChecker && !rctChecker(collision)) + if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) { return; + } fillMCCollision(collision, mcParticles); } PROCESS_SWITCH(ResonanceModuleInitializer, processRun3MC, "process MC for RUN3", false); @@ -607,31 +548,6 @@ struct ResonanceModuleInitializer { fillMCCollision(collision, mcParticles); } PROCESS_SWITCH(ResonanceModuleInitializer, processRun2MC, "process MC for RUN2", false); - - /** - * @brief Processes Spherocity - * - * @param collision Collision data - * @param tracks Track data - */ - void processSpherocity(soa::Filtered::iterator const& /*collision*/, aod::ResoTrackCandidates const& tracks) - { - float spherocity = computeSpherocity(tracks, cfgTrackSphMin, cfgTrackSphDef); - resoSpheroCollisions(spherocity); - } - PROCESS_SWITCH(ResonanceModuleInitializer, processSpherocity, "process Spherocity", false); - - /** - * @brief Processes Event Plane - * - * @param collision Collision data with Qvectors - * @param tracks Track data - */ - void processEventPlane(soa::Filtered>::iterator const& collision) - { - resoEvtPlCollisions(getEvtPl(collision), getEvtPlRes(collision, evtPlDetId, evtPlRefAId), getEvtPlRes(collision, evtPlDetId, evtPlRefBId), getEvtPlRes(collision, evtPlRefAId, evtPlRefBId)); - } - PROCESS_SWITCH(ResonanceModuleInitializer, processEventPlane, "process Event Plane", false); }; /** @@ -641,42 +557,89 @@ struct ResonanceModuleInitializer { * It applies daughter selection criteria and fills QA histograms for daughter properties. */ struct ResonanceDaughterInitializer { - SliceCache cache; - Produces reso2trks; ///< Output table for resonance tracks - Produces resoTrackTracks; ///< Output table for resonance track tracks - Produces reso2microtrks; ///< Output table for resonance microtracks - Produces resoMicroTrackTracks; ///< Output table for resonance microtrack tracks - Produces reso2mctracks; ///< Output table for MC resonance tracks - Produces reso2v0s; ///< Output table for resonance V0s - Produces resoV0V0s; ///< Output table for resonance V0-V0s - Produces reso2mcv0s; ///< Output table for MC resonance V0s - Produces reso2cascades; ///< Output table for resonance cascades - Produces resoCascadeCascades; ///< Output table for resonance cascade-cascades - Produces reso2mccascades; ///< Output table for MC resonance cascades - - // Configurables + enum class UltraMicroPidSpecies : uint8_t { + Pion, + Kaon, + Proton + }; + + static constexpr int TrackSelectionNone = 0; + static constexpr int TrackSelectionGlobal = 1; + static constexpr int TrackSelectionGlobalWoPtEta = 2; + static constexpr int TrackSelectionGlobalWoDCA = 3; + static constexpr int TrackSelectionQuality = 4; + static constexpr int TrackSelectionInAcceptance = 5; + static constexpr float MomentumQuantizationScale = 1000.f; + static constexpr std::size_t StoredMCRelationCount = 2; + + UltraMicroPidSpecies ultraMicroPidSpecies = UltraMicroPidSpecies::Pion; + bool warnedUltraMicroMomentumRange = false; + + Preslice> tracksPerCollision = aod::track::collisionId; + Preslice> tracksMCPerCollision = aod::track::collisionId; + Preslice v0sPerCollision = aod::v0data::collisionId; + Preslice v0sMCPerCollision = aod::v0data::collisionId; + Preslice cascadesPerCollision = aod::cascdata::collisionId; + Preslice cascadesMCPerCollision = aod::cascdata::collisionId; + Produces reso2trks; ///< Output table for resonance tracks + Produces resoTrackTracks; ///< Output table for original track row IDs + Produces reso2microtrks; ///< Output table for resonance microtracks + Produces resoMicroTrackTracks; ///< Output table for original microtrack row IDs + Produces reso2ultramicrotrks; ///< Output table for resonance ultra-microtracks + Produces resoUltraMicroTrackTracks; ///< Output table for original ultra-microtrack row IDs + Produces reso2mctracks; ///< Output table for MC resonance tracks + Produces reso2v0s; ///< Output table for resonance V0s + Produces resoV0V0s; ///< Output table for original V0 row IDs + Produces reso2mcv0s; ///< Output table for MC resonance V0s + Produces reso2cascades; ///< Output table for resonance cascades + Produces resoCascadeCascades; ///< Output table for original cascade row IDs + Produces reso2mccascades; ///< Output table for MC resonance cascades + + // General daughter output options Configurable cfgFillQA{"cfgFillQA", false, "Fill QA histograms"}; - Configurable cfgFillMicroTracks{"cfgFillMicroTracks", false, "Fill micro tracks"}; - Configurable cfgBypassTrackFill{"cfgBypassTrackFill", true, "Bypass track fill"}; - Configurable cfgBypassTrackIndexFill{"cfgBypassTrackIndexFill", false, "Bypass track index fill"}; - - // Configurables for tracks - Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 2.0, "Track DCAr cut to PV Maximum"}; - Configurable cMinDCArToPVcut{"cMinDCArToPVcut", 0.0, "Track DCAr cut to PV Minimum"}; - Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; - Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; - Configurable pidnSigmaPreSelectionCut{"pidnSigmaPreSelectionCut", 5.0f, "TPC and TOF PID cut (loose, improve performance)"}; - Configurable trackSelection{"trackSelection", 1, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; - - // Configurables for V0s - Configurable cMinV0Radius{"cMinV0Radius", 0.0, "Minimum V0 radius from PV"}; - Configurable cMaxV0Radius{"cMaxV0Radius", 200.0, "Maximum V0 radius from PV"}; - Configurable cMinV0CosPA{"cMinV0CosPA", 0.995, "Minimum V0 CosPA to PV"}; - - // Configurables for cascades - Configurable cMinCascRadius{"cMinCascRadius", 0.0, "Minimum Cascade radius from PV"}; - Configurable cMaxCascRadius{"cMaxCascRadius", 200.0, "Maximum Cascade radius from PV"}; - Configurable cMinCascCosPA{"cMinCascCosPA", 0.97, "Minimum Cascade CosPA to PV"}; + Configurable cfgDetailTrackQA{"cfgDetailTrackQA", false, "Fill detailed QA histograms for enabled track output tables"}; + + // Track pre-selection and DCA cuts + struct : ConfigurableGroup { + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable cfgCutMinPt{"cfgCutMinPt", 0.1f, "Minimum pT for tracks (GeV/c)"}; + Configurable cfgCutMaxPt{"cfgCutMaxPt", 32.767f, "Maximum pT for tracks (GeV/c)"}; + Configurable pidnSigmaPreSelectionCut{"pidnSigmaPreSelectionCut", 5.0f, "TPC PID cut (loose, improve performance)"}; + Configurable mincrossedrows{"mincrossedrows", 70, "Minimum crossed rows for V0 daughter tracks"}; + Configurable trackSelection{"trackSelection", 3, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; + Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 2.0, "Track DCAr cut to PV Maximum"}; + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; + Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; + Configurable cfgApplyTightDCAPtDepSelection{"cfgApplyTightDCAPtDepSelection", true, "Apply the pT-dependent tight DCA selection"}; + Configurable cfgTightDCAOffset{"cfgTightDCAOffset", 0.004f, "Constant term of the tight DCA threshold (cm)"}; + Configurable cfgTightDCAPtCoefficient{"cfgTightDCAPtCoefficient", 0.013f, "Coefficient of the pT-dependent tight DCA threshold"}; + Configurable cfgTightDCAPtPower{"cfgTightDCAPtPower", 1.f, "Power in tight DCA = offset + coefficient / pT^power"}; + } TrackCuts; + + // V0 and V0-daughter cuts + struct : ConfigurableGroup { + Configurable cMinV0PosDCArToPVcut{"cMinV0PosDCArToPVcut", 0.05f, "V0 Positive Track DCAr cut to PV Minimum"}; + Configurable cMinV0NegDCArToPVcut{"cMinV0NegDCArToPVcut", 0.05f, "V0 Negative Track DCAr cut to PV Minimum"}; + Configurable cMinV0Radius{"cMinV0Radius", 0.0, "Minimum V0 radius from PV"}; + Configurable cMaxV0Radius{"cMaxV0Radius", 200.0, "Maximum V0 radius from PV"}; + Configurable cMinV0CosPA{"cMinV0CosPA", 0.995, "Minimum V0 CosPA to PV"}; + } V0Cuts; + + // Cascade and cascade-daughter cuts + struct : ConfigurableGroup { + Configurable cfgMinCrossedRowsCascBach{"cfgMinCrossedRowsCascBach", 70, "min crossed rows for bachelor track from cascade"}; + Configurable cMinCascBachDCArToPVcut{"cMinCascBachDCArToPVcut", 0.05f, "Cascade Bachelor Track DCAr cut to PV Minimum"}; + Configurable cMaxCascBachDCArToPVcut{"cMaxCascBachDCArToPVcut", 999.0f, "Cascade Bachelor Track DCAr cut to PV Maximum"}; + Configurable cMaxCascDCAV0Daughters{"cMaxCascDCAV0Daughters", 1.6, "Cascade DCA between V0 daughters Maximum"}; + Configurable cMaxCascDCACascDaughters{"cMaxCascDCACascDaughters", 1.6, "Cascade DCA between Casc daughters Maximum"}; + Configurable cMinCascV0CosPA{"cMinCascV0CosPA", 0.97, "Minimum Cascade V0 CosPA to PV"}; + Configurable cMaxCascV0Radius{"cMaxCascV0Radius", 200.0, "Maximum Cascade V0 radius from PV"}; + Configurable cMinCascV0Radius{"cMinCascV0Radius", 0.0, "Minimum Cascade V0 radius from PV"}; + Configurable cMinCascRadius{"cMinCascRadius", 0.0, "Minimum Cascade radius from PV"}; + Configurable cMaxCascRadius{"cMaxCascRadius", 200.0, "Maximum Cascade radius from PV"}; + Configurable cMinCascCosPA{"cMinCascCosPA", 0.97, "Minimum Cascade CosPA to PV"}; + Configurable cCascMassResol{"cCascMassResol", 999, "Cascade mass resolution"}; + } CascadeCuts; // Derived dataset selections struct : ConfigurableGroup { @@ -686,22 +649,57 @@ struct ResonanceDaughterInitializer { Configurable cfgFillPionMicroTracks{"cfgFillPionMicroTracks", false, "Fill pion micro tracks"}; Configurable cfgFillKaonMicroTracks{"cfgFillKaonMicroTracks", false, "Fill kaon micro tracks"}; Configurable cfgFillProtonMicroTracks{"cfgFillProtonMicroTracks", false, "Fill proton micro tracks"}; + Configurable cfgFillPionUltraMicroTracks{"cfgFillPionUltraMicroTracks", true, "Fill pion ultra micro tracks"}; + Configurable cfgFillKaonUltraMicroTracks{"cfgFillKaonUltraMicroTracks", false, "Fill kaon ultra micro tracks"}; + Configurable cfgFillProtonUltraMicroTracks{"cfgFillProtonUltraMicroTracks", false, "Fill proton ultra micro tracks"}; + Configurable cfgFillK0s{"cfgFillK0s", false, "Fill K0s"}; + Configurable cfgFillLambda0{"cfgFillLambda0", false, "Fill Lambda0"}; + Configurable cfgBypassNoPairV0s{"cfgBypassNoPairV0s", false, "In a *WithPairGate process, bypass track fill if no V0 passes the configured selections"}; + Configurable cfgBypassNoPairCascades{"cfgBypassNoPairCascades", true, "In a *WithPairGate process, bypass track fill if no cascade passes the configured selections"}; + Configurable cfgFillMicroTracks{"cfgFillMicroTracks", false, "Fill micro tracks"}; + Configurable cfgFillUltraMicroTracks{"cfgFillUltraMicroTracks", false, "Fill ultra micro tracks"}; + Configurable cfgBypassTrackFill{"cfgBypassTrackFill", false, "Bypass the full ResoTracks table fill"}; + Configurable cfgBypassTrackIndexFill{"cfgBypassTrackIndexFill", false, "Bypass original daughter ID table fill"}; } FilterForDerivedTables; - // Filters - Filter dcaXYFilter = nabs(aod::track::dcaXY) < cMaxDCArToPVcut && nabs(aod::track::dcaXY) > cMinDCArToPVcut; - Filter dcaZFilter = nabs(aod::track::dcaZ) < cMaxDCAzToPVcut && nabs(aod::track::dcaZ) > cMinDCAzToPVcut; - Preslice perMcCollision = aod::mcparticle::mcCollisionId; + // Secondary selections for K0s and Lambda0 + struct : ConfigurableGroup { + Configurable cfgSecondaryRequire{"cfgSecondaryRequire", false, "Secondary cuts on/off"}; + Configurable cfgSecondaryArmenterosCut{"cfgSecondaryArmenterosCut", false, "cut on Armenteros-Podolanski graph"}; + Configurable cfgSecondaryCrossMassHypothesisCut{"cfgSecondaryCrossMassHypothesisCut", false, "Apply cut based on the lambda mass hypothesis"}; + Configurable cfgByPassDauPIDSelection{"cfgByPassDauPIDSelection", true, "Bypass TPC PID preselection for V0 daughters"}; + Configurable cfgSecondaryDauDCAMax{"cfgSecondaryDauDCAMax", 0.2, "Maximum DCA Secondary daughters to PV"}; + Configurable cfgSecondaryDauPosDCAtoPVMin{"cfgSecondaryDauPosDCAtoPVMin", 0.0, "Minimum DCA Secondary positive daughters to PV"}; + Configurable cfgSecondaryDauNegDCAtoPVMin{"cfgSecondaryDauNegDCAtoPVMin", 0.0, "Minimum DCA Secondary negative daughters to PV"}; + Configurable cfgSecondaryPtMin{"cfgSecondaryPtMin", 0.f, "Minimum transverse momentum of Secondary"}; + Configurable cfgSecondaryRapidityMax{"cfgSecondaryRapidityMax", 0.5, "Maximum rapidity of Secondary"}; + Configurable cfgSecondaryRadiusMin{"cfgSecondaryRadiusMin", 0.0, "Minimum transverse radius of Secondary"}; + Configurable cfgSecondaryRadiusMax{"cfgSecondaryRadiusMax", 999.9, "Maximum transverse radius of Secondary"}; + Configurable cfgSecondaryCosPAMin{"cfgSecondaryCosPAMin", 0.998, "Mininum cosine pointing angle of Secondary"}; + Configurable cfgSecondaryDCAtoPVMax{"cfgSecondaryDCAtoPVMax", 0.4, "Maximum DCA Secondary to PV"}; + Configurable cfgSecondaryProperLifetimeMax{"cfgSecondaryProperLifetimeMax", 20., "Maximum Secondary Lifetime"}; + Configurable cfgSecondaryparamArmenterosCut{"cfgSecondaryparamArmenterosCut", 0.2, "parameter for Armenteros Cut"}; + Configurable cfgSecondaryMassWindow{"cfgSecondaryMassWindow", 0.03, "Secondary inv mass selection window"}; + Configurable cfgSecondaryCrossMassCutWindow{"cfgSecondaryCrossMassCutWindow", 0.05, "Secondary inv mass selection window with (anti)lambda hypothesis"}; + } SecondaryCuts; // Track selection filter based on configuration - Filter trackFilter = (trackSelection.node() == 0) || - ((trackSelection.node() == 1) && requireGlobalTrackInFilter()) || // kGlobalTrack = kQualityTracks | kPrimaryTracks | kInAcceptanceTracks - ((trackSelection.node() == 2) && requireGlobalTrackWoPtEtaInFilter()) || // kGlobalTrackWoPtEta = kQualityTracks | kPrimaryTracks - ((trackSelection.node() == 3) && requireGlobalTrackWoDCAInFilter()) || // kGlobalTrackWoDCA = kQualityTracks | kInAcceptanceTracks - ((trackSelection.node() == 4) && requireQualityTracksInFilter()) || // kQualityTracks = kQualityTracksITS | kQualityTracksTPC - ((trackSelection.node() == 5) && requireTrackCutInFilter(TrackSelectionFlags::kInAcceptanceTracks)); // kInAcceptanceTracks = kPtRange | kEtaRange + Filter trackFilter = (TrackCuts.trackSelection.node() == TrackSelectionNone) || + ((TrackCuts.trackSelection.node() == TrackSelectionGlobal) && requireGlobalTrackInFilter()) || // kGlobalTrack = kQualityTracks | kPrimaryTracks | kInAcceptanceTracks + ((TrackCuts.trackSelection.node() == TrackSelectionGlobalWoPtEta) && requireGlobalTrackWoPtEtaInFilter()) || // kGlobalTrackWoPtEta = kQualityTracks | kPrimaryTracks + ((TrackCuts.trackSelection.node() == TrackSelectionGlobalWoDCA) && requireGlobalTrackWoDCAInFilter()) || // kGlobalTrackWoDCA = kQualityTracks | kInAcceptanceTracks + ((TrackCuts.trackSelection.node() == TrackSelectionQuality) && requireQualityTracksInFilter()) || // kQualityTracks = kQualityTracksITS | kQualityTracksTPC + ((TrackCuts.trackSelection.node() == TrackSelectionInAcceptance) && requireTrackCutInFilter(TrackSelectionFlags::kInAcceptanceTracks)); // kInAcceptanceTracks = kPtRange | kEtaRange + Filter trackKinematicsFilter = nabs(aod::track::eta) < TrackCuts.cfgCutEta && + aod::track::pt >= TrackCuts.cfgCutMinPt && + aod::track::pt <= TrackCuts.cfgCutMaxPt; HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject}; + // The daughter task needs this row-wise mapping back to the original collision. + // Keep ResonanceModuleInitializer::cfgBypassCollIndexFill disabled and enable + // the matching Run 2/Run 3 base event process for MC workflows. + using ResoCollisionWithIndex = soa::Join; + /** * @brief Initializes the task * @@ -709,23 +707,162 @@ struct ResonanceDaughterInitializer { */ void init(InitContext&) { + const bool processTrackDataEnabled = doprocessData || doprocessDataWithPairGate; + const bool processTrackMCEnabled = doprocessMC || doprocessMCWithPairGate; + const int enabledTrackProcesses = static_cast(doprocessData) + + static_cast(doprocessDataWithPairGate) + + static_cast(doprocessMC) + + static_cast(doprocessMCWithPairGate); + + if (enabledTrackProcesses > 1) { + LOGF(fatal, "Only one track process can be enabled in ResonanceDaughterInitializer"); + } + if ((doprocessData || doprocessMC) && + (FilterForDerivedTables.cfgBypassNoPairV0s || FilterForDerivedTables.cfgBypassNoPairCascades)) { + LOGF(warn, "Pair-gate options are ignored by processData/processMC; enable the matching *WithPairGate process to apply them"); + } + if (doprocessDataWithPairGate && FilterForDerivedTables.cfgBypassNoPairV0s && !doprocessV0Data) { + LOGF(fatal, "cfgBypassNoPairV0s requires processV0Data so an accepted V0 is written for every retained collision"); + } + if (doprocessDataWithPairGate && FilterForDerivedTables.cfgBypassNoPairCascades && !doprocessCascData) { + LOGF(fatal, "cfgBypassNoPairCascades requires processCascData so an accepted cascade is written for every retained collision"); + } + + if (!std::isfinite(TrackCuts.cfgCutMinPt.value) || + !std::isfinite(TrackCuts.cfgCutMaxPt.value) || + TrackCuts.cfgCutMinPt.value < 0.f || + TrackCuts.cfgCutMaxPt.value <= TrackCuts.cfgCutMinPt.value) { + LOGF(fatal, "Track pT limits must be finite and satisfy 0 <= cfgCutMinPt < cfgCutMaxPt"); + } + if (!std::isfinite(TrackCuts.cfgCutEta.value) || TrackCuts.cfgCutEta.value <= 0.f) { + LOGF(fatal, "cfgCutEta must be finite and positive"); + } + if (TrackCuts.trackSelection.value < TrackSelectionNone || TrackCuts.trackSelection.value > TrackSelectionInAcceptance) { + LOGF(fatal, "trackSelection must be in [0, 5]"); + } + if (!std::isfinite(TrackCuts.cMaxDCArToPVcut.value) || + !std::isfinite(TrackCuts.cMaxDCAzToPVcut.value) || + !std::isfinite(TrackCuts.cMinDCAzToPVcut.value) || + TrackCuts.cMaxDCArToPVcut.value < 0. || + TrackCuts.cMinDCAzToPVcut.value < 0. || + TrackCuts.cMaxDCAzToPVcut.value < TrackCuts.cMinDCAzToPVcut.value) { + LOGF(fatal, "Track DCA limits must be finite and satisfy 0 <= cMinDCAzToPVcut <= cMaxDCAzToPVcut and cMaxDCArToPVcut >= 0"); + } + if (!std::isfinite(TrackCuts.pidnSigmaPreSelectionCut.value) || + TrackCuts.pidnSigmaPreSelectionCut.value < 0.f) { + LOGF(fatal, "pidnSigmaPreSelectionCut must be finite and non-negative"); + } + + if (TrackCuts.cfgApplyTightDCAPtDepSelection.value && + (!std::isfinite(TrackCuts.cfgTightDCAOffset.value) || + !std::isfinite(TrackCuts.cfgTightDCAPtCoefficient.value) || + !std::isfinite(TrackCuts.cfgTightDCAPtPower.value) || + TrackCuts.cfgTightDCAOffset.value < 0.f || + TrackCuts.cfgTightDCAPtCoefficient.value < 0.f || + TrackCuts.cfgTightDCAPtPower.value < 0.f)) { + LOGF(fatal, "Tight-DCA offset, pT coefficient, and power must be finite and non-negative"); + } + + if (FilterForDerivedTables.cfgFillUltraMicroTracks) { + constexpr float MaxQuantizedMomentum = static_cast(std::numeric_limits::max()) / MomentumQuantizationScale; + const float maxLongitudinalMomentum = TrackCuts.cfgCutMaxPt.value * std::sinh(TrackCuts.cfgCutEta.value); + if (TrackCuts.cfgCutMaxPt.value > MaxQuantizedMomentum || + !std::isfinite(maxLongitudinalMomentum) || + maxLongitudinalMomentum > MaxQuantizedMomentum) { + LOGF(fatal, "Ultra-micro cfgCutMaxPt/cfgCutEta allow momentum components beyond the int16_t quantization range"); + } + int enabledUltraMicroSpecies = 0; + enabledUltraMicroSpecies += FilterForDerivedTables.cfgFillPionUltraMicroTracks ? 1 : 0; + enabledUltraMicroSpecies += FilterForDerivedTables.cfgFillKaonUltraMicroTracks ? 1 : 0; + enabledUltraMicroSpecies += FilterForDerivedTables.cfgFillProtonUltraMicroTracks ? 1 : 0; + if (enabledUltraMicroSpecies != 1) { + LOGF(fatal, "Exactly one pion/kaon/proton PID species must be enabled when filling ultra-micro tracks"); + } + if (TrackCuts.pidnSigmaPreSelectionCut.value > o2::aod::resoultramicrodaughter::PidNSigma::MaxNSigma) { + LOGF(fatal, "Ultra-micro PID encoding requires pidnSigmaPreSelectionCut <= 5"); + } + if (FilterForDerivedTables.cfgFillKaonUltraMicroTracks) { + ultraMicroPidSpecies = UltraMicroPidSpecies::Kaon; + } else if (FilterForDerivedTables.cfgFillProtonUltraMicroTracks) { + ultraMicroPidSpecies = UltraMicroPidSpecies::Proton; + } else { + ultraMicroPidSpecies = UltraMicroPidSpecies::Pion; + } + } + if (cfgFillQA) { - AxisSpec idxAxis = {8, -0.5, 7.5, "Index"}; + AxisSpec idxAxis = {8, 0.0, 8.0, "Index"}; AxisSpec ptAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"}; + // The DCA maps cover the full configured pT range in 0.1 GeV/c steps. + constexpr float DcaPtBinWidth = 0.1f; + const int maxDcaPtBin = static_cast(std::ceil(TrackCuts.cfgCutMaxPt.value / DcaPtBinWidth)); + const int nDcaPtBins = maxDcaPtBin + 1; + const float dcaPtAxisHalfBin = 0.5f * DcaPtBinWidth; + AxisSpec dcaPtAxis = {nDcaPtBins, + -dcaPtAxisHalfBin, + maxDcaPtBin * DcaPtBinWidth + dcaPtAxisHalfBin, + "#it{p}_{T} (GeV/#it{c})"}; AxisSpec etaAxis = {100, -1.0f, 1.0f, "#eta"}; AxisSpec phiAxis = {100, 0.0f, TwoPI, "#phi"}; - - qaRegistry.add("QA/hGoodTrackIndices", "hGoodTrackIndices", kTH1D, {idxAxis}); - if (doprocessMC) { - qaRegistry.add("QA/hGoodMCTrackIndices", "hGoodMCTrackIndices", kTH1D, {idxAxis}); + // Keep the configured signed DCA limits at bin centres so tracks exactly + // on an accepted selection boundary are not moved to overflow. + constexpr int NDcaBins = 201; + constexpr float DcaAxisPaddingFraction = 1.f / 200.f; + const auto configuredDcaXYMax = static_cast(TrackCuts.cMaxDCArToPVcut.value); + const auto configuredDcaZMax = static_cast(TrackCuts.cMaxDCAzToPVcut.value); + const float dcaXYAxisHalfRange = configuredDcaXYMax > 0.f ? configuredDcaXYMax * (1.f + DcaAxisPaddingFraction) : 1.e-4f; + const float dcaZAxisHalfRange = configuredDcaZMax > 0.f ? configuredDcaZMax * (1.f + DcaAxisPaddingFraction) : 1.e-4f; + AxisSpec dcaXYAxis = {NDcaBins, -dcaXYAxisHalfRange, dcaXYAxisHalfRange, "DCA_{xy} (cm)"}; + AxisSpec dcaZAxis = {NDcaBins, -dcaZAxisHalfRange, dcaZAxisHalfRange, "DCA_{z} (cm)"}; + // Tiny PID has 255 discrete values from -6.35 to 6.35 in 0.05 steps. + // Extend the histogram limits by half a step so every decoded value is + // located at a bin centre instead of on a floating-point bin boundary. + constexpr int NTinyTPCPidBins = aod::pidtpc_tiny::binning::nbins + 1; + constexpr float HalfTinyTPCPidBinWidth = 0.5f * aod::pidtpc_tiny::binning::bin_width; + constexpr int NTinyTOFPidBins = aod::pidtof_tiny::binning::nbins + 1; + constexpr float HalfTinyTOFPidBinWidth = 0.5f * aod::pidtof_tiny::binning::bin_width; + AxisSpec nSigmaTPCAxis = {NTinyTPCPidBins, + aod::pidtpc_tiny::binning::binned_min - HalfTinyTPCPidBinWidth, + aod::pidtpc_tiny::binning::binned_max + HalfTinyTPCPidBinWidth, + "TPC N#sigma"}; + AxisSpec nSigmaTOFAxis = {NTinyTOFPidBins, + aod::pidtof_tiny::binning::binned_min - HalfTinyTOFPidBinWidth, + aod::pidtof_tiny::binning::binned_max + HalfTinyTOFPidBinWidth, + "TOF N#sigma"}; + + if (processTrackDataEnabled || processTrackMCEnabled) { + qaRegistry.add("QA/hGoodTrackIndices", "hGoodTrackIndices", kTH1D, {idxAxis}); + if (processTrackMCEnabled) { + qaRegistry.add("QA/hGoodMCTrackIndices", "hGoodMCTrackIndices", kTH1D, {idxAxis}); + } + if (cfgDetailTrackQA) { + if (!FilterForDerivedTables.cfgBypassTrackFill) { + qaRegistry.add("QA/h4TrackPtEtaPhi", "ResoTracks pT, eta, phi", kTH3F, {ptAxis, etaAxis, phiAxis}); + qaRegistry.add("QA/h2TrackDCAxyVsPt", "ResoTracks DCAxy vs pT", kTH2F, {dcaPtAxis, dcaXYAxis}); + qaRegistry.add("QA/h2TrackDCAzVsPt", "ResoTracks DCAz vs pT", kTH2F, {dcaPtAxis, dcaZAxis}); + qaRegistry.add("QA/h4TrackTPCnSigma", "ResoTracks TPC nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTPCAxis, nSigmaTPCAxis, nSigmaTPCAxis}); + qaRegistry.add("QA/h4TrackTOFnSigma", "ResoTracks TOF nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTOFAxis, nSigmaTOFAxis, nSigmaTOFAxis}); + } + if (FilterForDerivedTables.cfgFillMicroTracks) { + qaRegistry.add("QA/h4MicroTrackPtEtaPhi", "ResoMicroTracks pT, eta, phi", kTH3F, {ptAxis, etaAxis, phiAxis}); + qaRegistry.add("QA/h2MicroTrackDCAxyVsPt", "ResoMicroTracks DCAxy vs pT", kTH2F, {dcaPtAxis, dcaXYAxis}); + qaRegistry.add("QA/h2MicroTrackDCAzVsPt", "ResoMicroTracks DCAz vs pT", kTH2F, {dcaPtAxis, dcaZAxis}); + qaRegistry.add("QA/h4MicroTrackTPCnSigma", "ResoMicroTracks TPC nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTPCAxis, nSigmaTPCAxis, nSigmaTPCAxis}); + qaRegistry.add("QA/h4MicroTrackTOFnSigma", "ResoMicroTracks TOF nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTOFAxis, nSigmaTOFAxis, nSigmaTOFAxis}); + } + if (FilterForDerivedTables.cfgFillUltraMicroTracks) { + qaRegistry.add("QA/h4UltraMicroTrackPtEtaPhi", "ResoUltraMicroTracks pT, eta, phi", kTH3F, {ptAxis, etaAxis, phiAxis}); + qaRegistry.add("QA/h2UltraMicroTrackDCAxyVsPt", "ResoUltraMicroTracks DCAxy vs pT", kTH2F, {dcaPtAxis, dcaXYAxis}); + qaRegistry.add("QA/h2UltraMicroTrackDCAzVsPt", "ResoUltraMicroTracks DCAz vs pT", kTH2F, {dcaPtAxis, dcaZAxis}); + qaRegistry.add("QA/h4UltraMicroTrackTPCnSigma", "ResoUltraMicroTracks TPC nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTPCAxis, nSigmaTPCAxis, nSigmaTPCAxis}); + qaRegistry.add("QA/h4UltraMicroTrackTOFnSigma", "ResoUltraMicroTracks TOF nSigma Pi, Ka, Pr as pT", kTHnSparseD, {ptAxis, nSigmaTOFAxis, nSigmaTOFAxis, nSigmaTOFAxis}); + } + } } - qaRegistry.add("QA/hTrackPt", "Track pT", kTH1F, {ptAxis}); - qaRegistry.add("QA/hTrackEta", "Track eta", kTH1F, {etaAxis}); - qaRegistry.add("QA/hTrackPhi", "Track phi", kTH1F, {phiAxis}); if (doprocessV0Data || doprocessV0MC) { qaRegistry.add("QA/hGoodV0Indices", "hGoodV0Indices", kTH1D, {idxAxis}); - if (doprocessMC) { + if (doprocessV0MC) { qaRegistry.add("QA/hGoodMCV0Indices", "hGoodMCV0Indices", kTH1D, {idxAxis}); } AxisSpec radiusAxis = {100, 0.0, 200.0, "V0 Radius"}; @@ -738,14 +875,14 @@ struct ResonanceDaughterInitializer { AxisSpec radiusAxis = {100, 0.0, 200.0, "Cascade Radius"}; AxisSpec cosPAAxis = {100, 0.97, 1.0, "Cascade CosPA"}; qaRegistry.add("QA/hGoodCascIndices", "hGoodCascIndices", kTH1D, {idxAxis}); - if (doprocessMC) { + if (doprocessCascMC) { qaRegistry.add("QA/hGoodMCCascIndices", "hGoodMCCascIndices", kTH1D, {idxAxis}); } qaRegistry.add("QA/hCascRadius", "Cascade Radius", kTH1F, {radiusAxis}); qaRegistry.add("QA/hCascCosPA", "Cascade CosPA", kTH1F, {cosPAAxis}); } } - if (doprocessData || doprocessMC) { + if (processTrackDataEnabled || processTrackMCEnabled) { LOGF(info, "ResonanceDaughterInitializer initialized with tracks"); } if (doprocessV0Data || doprocessV0MC) { @@ -756,11 +893,11 @@ struct ResonanceDaughterInitializer { } // Check if the module is initialized with both data and MC - if ((doprocessData && doprocessMC) || (doprocessV0Data && doprocessV0MC) || (doprocessCascData && doprocessCascMC)) { + if ((processTrackDataEnabled && processTrackMCEnabled) || (doprocessV0Data && doprocessV0MC) || (doprocessCascData && doprocessCascMC)) { LOGF(fatal, "ResonanceDaughterInitializer initialized with both data and MC"); } // Check if none of the processes are enabled - if (!doprocessDummy && !doprocessData && !doprocessMC && !doprocessV0Data && !doprocessV0MC && !doprocessCascData && !doprocessCascMC) { + if (!doprocessDummy && !processTrackDataEnabled && !processTrackMCEnabled && !doprocessV0Data && !doprocessV0MC && !doprocessCascData && !doprocessCascMC) { LOGF(fatal, "ResonanceDaughterInitializer not initialized, enable at least one process"); } } @@ -768,19 +905,37 @@ struct ResonanceDaughterInitializer { bool filterMicroTrack(T const& track) { // if no selection is requested, return true - if (!FilterForDerivedTables.cfgFillPionMicroTracks && !FilterForDerivedTables.cfgFillKaonMicroTracks && !FilterForDerivedTables.cfgFillProtonMicroTracks) + if (!FilterForDerivedTables.cfgFillPionMicroTracks && !FilterForDerivedTables.cfgFillKaonMicroTracks && !FilterForDerivedTables.cfgFillProtonMicroTracks) { return true; + } if (FilterForDerivedTables.cfgFillPionMicroTracks) { - if (std::abs(track.tpcNSigmaPi()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } } if (FilterForDerivedTables.cfgFillKaonMicroTracks) { - if (std::abs(track.tpcNSigmaKa()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaKa()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } } if (FilterForDerivedTables.cfgFillProtonMicroTracks) { - if (std::abs(track.tpcNSigmaPr()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaPr()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } + } + return false; + } + + template + bool filterUltraMicroTrack(T const& track) + { + switch (ultraMicroPidSpecies) { + case UltraMicroPidSpecies::Pion: + return std::abs(track.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut; + case UltraMicroPidSpecies::Kaon: + return std::abs(track.tpcNSigmaKa()) < TrackCuts.pidnSigmaPreSelectionCut; + case UltraMicroPidSpecies::Proton: + return std::abs(track.tpcNSigmaPr()) < TrackCuts.pidnSigmaPreSelectionCut; } return false; } @@ -789,23 +944,359 @@ struct ResonanceDaughterInitializer { bool filterTrack(T const& track) { // if no selection is requested, return true - if (!FilterForDerivedTables.cfgFillPionTracks && !FilterForDerivedTables.cfgFillKaonTracks && !FilterForDerivedTables.cfgFillProtonTracks) + if (!FilterForDerivedTables.cfgFillPionTracks && !FilterForDerivedTables.cfgFillKaonTracks && !FilterForDerivedTables.cfgFillProtonTracks) { return true; + } if (FilterForDerivedTables.cfgFillPionTracks) { - if (std::abs(track.tpcNSigmaPi()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } } if (FilterForDerivedTables.cfgFillKaonTracks) { - if (std::abs(track.tpcNSigmaKa()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaKa()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } } if (FilterForDerivedTables.cfgFillProtonTracks) { - if (std::abs(track.tpcNSigmaPr()) < pidnSigmaPreSelectionCut) + if (std::abs(track.tpcNSigmaPr()) < TrackCuts.pidnSigmaPreSelectionCut) { return true; + } } return false; } + template + bool filterV0(CollisionType const& collision, V0Type const& v0, TrackType const&) + { + if (!FilterForDerivedTables.cfgFillK0s && !FilterForDerivedTables.cfgFillLambda0) { + return true; + } + if (!SecondaryCuts.cfgSecondaryRequire) { + return true; + } + if (v0.dcaV0daughters() > SecondaryCuts.cfgSecondaryDauDCAMax || + std::abs(v0.dcapostopv()) < SecondaryCuts.cfgSecondaryDauPosDCAtoPVMin || + std::abs(v0.dcanegtopv()) < SecondaryCuts.cfgSecondaryDauNegDCAtoPVMin || + v0.pt() < SecondaryCuts.cfgSecondaryPtMin || + v0.v0radius() < SecondaryCuts.cfgSecondaryRadiusMin || + v0.v0radius() > SecondaryCuts.cfgSecondaryRadiusMax || + v0.dcav0topv() > SecondaryCuts.cfgSecondaryDCAtoPVMax || + v0.v0cosPA() < SecondaryCuts.cfgSecondaryCosPAMin) { + return false; + } + if (SecondaryCuts.cfgSecondaryArmenterosCut && + v0.qtarm() < SecondaryCuts.cfgSecondaryparamArmenterosCut * std::abs(v0.alpha())) { + return false; + } + + const float decayLengthOverMomentum = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()); + const auto posTrack = v0.template posTrack_as(); + const auto negTrack = v0.template negTrack_as(); + const bool bypassDaughterPID = SecondaryCuts.cfgByPassDauPIDSelection; + bool selected = false; + if (FilterForDerivedTables.cfgFillK0s) { + const bool passesK0DaughterPID = bypassDaughterPID || + (std::abs(posTrack.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut && + std::abs(negTrack.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut); + const bool passesK0 = std::fabs(v0.yK0Short()) <= SecondaryCuts.cfgSecondaryRapidityMax && + decayLengthOverMomentum * MassK0Short <= SecondaryCuts.cfgSecondaryProperLifetimeMax && + std::fabs(v0.mK0Short() - MassK0Short) <= SecondaryCuts.cfgSecondaryMassWindow && + passesK0DaughterPID && + (!SecondaryCuts.cfgSecondaryCrossMassHypothesisCut || + (std::fabs(v0.mLambda() - MassLambda0) >= SecondaryCuts.cfgSecondaryCrossMassCutWindow && + std::fabs(v0.mAntiLambda() - MassLambda0Bar) >= SecondaryCuts.cfgSecondaryCrossMassCutWindow)); + selected = passesK0; + } + if (FilterForDerivedTables.cfgFillLambda0) { + const bool passesLambdaPID = bypassDaughterPID || + (std::abs(posTrack.tpcNSigmaPr()) < TrackCuts.pidnSigmaPreSelectionCut && + std::abs(negTrack.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut); + const bool passesAntiLambdaPID = bypassDaughterPID || + (std::abs(posTrack.tpcNSigmaPi()) < TrackCuts.pidnSigmaPreSelectionCut && + std::abs(negTrack.tpcNSigmaPr()) < TrackCuts.pidnSigmaPreSelectionCut); + const bool passesLambdaMassAndPID = + (std::fabs(v0.mLambda() - MassLambda0) <= SecondaryCuts.cfgSecondaryMassWindow && passesLambdaPID) || + (std::fabs(v0.mAntiLambda() - MassLambda0Bar) <= SecondaryCuts.cfgSecondaryMassWindow && passesAntiLambdaPID); + const bool passesLambda = std::fabs(v0.yLambda()) <= SecondaryCuts.cfgSecondaryRapidityMax && + decayLengthOverMomentum * MassLambda0 <= SecondaryCuts.cfgSecondaryProperLifetimeMax && + passesLambdaMassAndPID && + (!SecondaryCuts.cfgSecondaryCrossMassHypothesisCut || + std::fabs(v0.mK0Short() - MassK0Short) >= SecondaryCuts.cfgSecondaryCrossMassCutWindow); + selected = selected || passesLambda; + } + return selected; + } + + template + bool isMicroTrackSelected(CollisionType const&, TrackType const& track) + { + if (!std::isfinite(track.dcaXY()) || !std::isfinite(track.dcaZ())) { + return false; + } + if (std::fabs(track.dcaXY()) > TrackCuts.cMaxDCArToPVcut) { + return false; + } + if (std::fabs(track.dcaZ()) > TrackCuts.cMaxDCAzToPVcut || std::fabs(track.dcaZ()) < TrackCuts.cMinDCAzToPVcut) { + return false; + } + + return true; + } + + template + bool isTrackSelected(CollisionType const&, TrackType const& track) + { + if (cfgFillQA) { + qaRegistry.fill(HIST("QA/hGoodTrackIndices"), 0.5); + if constexpr (isMC) { + qaRegistry.fill(HIST("QA/hGoodMCTrackIndices"), 0.5); + } + } + if (!std::isfinite(track.dcaXY()) || !std::isfinite(track.dcaZ())) { + return false; + } + if (std::fabs(track.dcaXY()) > TrackCuts.cMaxDCArToPVcut) { + return false; + } + if (cfgFillQA) { + qaRegistry.fill(HIST("QA/hGoodTrackIndices"), 1.5); + } + if (std::fabs(track.dcaZ()) > TrackCuts.cMaxDCAzToPVcut || std::fabs(track.dcaZ()) < TrackCuts.cMinDCAzToPVcut) { + return false; + } + if (cfgFillQA) { + qaRegistry.fill(HIST("QA/hGoodTrackIndices"), 2.5); + qaRegistry.fill(HIST("QA/hGoodTrackIndices"), 7.5); + } + return true; + } + + template + bool isV0Selected(CollisionType const&, V0Type const& v0, TrackType const&, bool fillSelectionQA = true) + { + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodV0Indices"), 0.5); + } + + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + if (posTrack.tpcNClsCrossedRows() < TrackCuts.mincrossedrows || negTrack.tpcNClsCrossedRows() < TrackCuts.mincrossedrows) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodV0Indices"), 1.5); + } + if (std::fabs(posTrack.dcaXY()) < V0Cuts.cMinV0PosDCArToPVcut || + std::fabs(negTrack.dcaXY()) < V0Cuts.cMinV0NegDCArToPVcut) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodV0Indices"), 2.5); + } + if (v0.v0radius() > V0Cuts.cMaxV0Radius || v0.v0radius() < V0Cuts.cMinV0Radius) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodV0Indices"), 3.5); + } + if (v0.v0cosPA() < V0Cuts.cMinV0CosPA) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodV0Indices"), 4.5); + if constexpr (isMC) { + qaRegistry.fill(HIST("QA/hGoodMCV0Indices"), 0.5); + } + } + return true; + } + + template + bool isCascSelected(CollisionType const& collision, CascType const& casc, TrackType const&, bool fillSelectionQA = true) + { + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 0.5); + } + + auto bachelor = casc.template bachelor_as(); + if (bachelor.tpcNClsCrossedRows() < CascadeCuts.cfgMinCrossedRowsCascBach) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 1.5); + } + if (std::fabs(bachelor.dcaXY()) < CascadeCuts.cMinCascBachDCArToPVcut || + std::fabs(bachelor.dcaXY()) > CascadeCuts.cMaxCascBachDCArToPVcut) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 2.5); + } + if (casc.dcaV0daughters() > CascadeCuts.cMaxCascDCAV0Daughters || + casc.dcacascdaughters() > CascadeCuts.cMaxCascDCACascDaughters) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 3.5); + } + if (casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()) < CascadeCuts.cMinCascCosPA || + casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < CascadeCuts.cMinCascV0CosPA) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 4.5); + } + auto v0Radius = casc.v0radius(); + if (v0Radius > CascadeCuts.cMaxCascV0Radius || v0Radius < CascadeCuts.cMinCascV0Radius) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 5.5); + } + auto cascRadius = casc.cascradius(); + if (cascRadius > CascadeCuts.cMaxCascRadius || cascRadius < CascadeCuts.cMinCascRadius) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 6.5); + } + if (std::abs(casc.mXi() - MassXiMinus) > CascadeCuts.cCascMassResol) { + return false; + } + if (cfgFillQA && fillSelectionQA) { + qaRegistry.fill(HIST("QA/hGoodCascIndices"), 7.5); + if constexpr (isMC) { + qaRegistry.fill(HIST("QA/hGoodMCCascIndices"), 0.5); + } + } + return true; + } + + /// @brief Check whether a collision has at least one V0 that would be written + template + bool hasSelectedV0(CollisionType const& collision, V0Type const& v0s, TrackType const& tracks) + { + for (auto const& v0 : v0s) { + if (isV0Selected(collision, v0, tracks, false) && filterV0(collision, v0, tracks)) { + return true; + } + } + return false; + } + + /// @brief Check whether a collision has at least one cascade that would be written + template + bool hasSelectedCascade(CollisionType const& collision, CascType const& cascades, TrackType const& tracks) + { + for (auto const& casc : cascades) { + if (isCascSelected(collision, casc, tracks, false)) { + return true; + } + } + return false; + } + + static bool quantizeP(float p, int16_t& quantized) + { + if (!std::isfinite(p)) { + return false; + } + const double rounded = std::round(static_cast(p) * MomentumQuantizationScale); + if (rounded < static_cast(std::numeric_limits::min()) || + rounded > static_cast(std::numeric_limits::max())) { + return false; + } + quantized = static_cast(rounded); + return true; + } + + float tightDCAThreshold(float pt) const + { + if (!std::isfinite(pt) || pt <= 0.f) { + return -1.f; + } + if (TrackCuts.cfgTightDCAPtCoefficient.value == 0.f) { + return TrackCuts.cfgTightDCAOffset.value; + } + const float threshold = TrackCuts.cfgTightDCAOffset.value + + TrackCuts.cfgTightDCAPtCoefficient.value / + std::pow(pt, TrackCuts.cfgTightDCAPtPower.value); + return std::isfinite(threshold) ? threshold : -1.f; + } + + template + void fillUltraMicroTracks(CollisionType const& collision, TrackType const& tracks) + { + // Loop over tracks + for (auto const& track : tracks) { + if (!isMicroTrackSelected(collision, track)) { + continue; + } + if (!filterUltraMicroTrack(track)) { + continue; + } + if (!o2::aod::resoultramicrodaughter::DCAEncoding::isValid(track.dcaXY()) || + !o2::aod::resoultramicrodaughter::DCAEncoding::isValid(track.dcaZ())) { + continue; + } + o2::aod::resoultramicrodaughter::DCAEncoding dcaEncoding(track.dcaXY(), track.dcaZ()); + int16_t px1000 = 0; + int16_t py1000 = 0; + int16_t pz1000 = 0; + if (!quantizeP(track.px(), px1000) || + !quantizeP(track.py(), py1000) || + !quantizeP(track.pz(), pz1000)) { + if (!warnedUltraMicroMomentumRange) { + LOGF(warn, "Skipping ultra-micro tracks with non-finite or out-of-range momentum components"); + warnedUltraMicroMomentumRange = true; + } + continue; + } + uint8_t trackFlags = (track.passedITSRefit() << 0) | + (track.passedTPCRefit() << 1) | + (track.isGlobalTrackWoDCA() << 2) | + (track.isGlobalTrack() << 3) | + (track.isPrimaryTrack() << 4) | + (track.isPVContributor() << 5) | + (track.hasTOF() << 6) | + ((track.sign() > 0) << 7); // sign +1: 1, -1: 0 + uint8_t pidFlag = 0; + switch (ultraMicroPidSpecies) { + case UltraMicroPidSpecies::Pion: + pidFlag = static_cast(o2::aod::resoultramicrodaughter::PidNSigma( + track.tpcNSigmaPi(), track.tofNSigmaPi(), track.hasTOF())); + break; + case UltraMicroPidSpecies::Kaon: + pidFlag = static_cast(o2::aod::resoultramicrodaughter::PidNSigma( + track.tpcNSigmaKa(), track.tofNSigmaKa(), track.hasTOF())); + break; + case UltraMicroPidSpecies::Proton: + pidFlag = static_cast(o2::aod::resoultramicrodaughter::PidNSigma( + track.tpcNSigmaPr(), track.tofNSigmaPr(), track.hasTOF())); + break; + } + if (cfgFillQA && cfgDetailTrackQA) { + qaRegistry.fill(HIST("QA/h4UltraMicroTrackPtEtaPhi"), track.pt(), track.eta(), track.phi()); + qaRegistry.fill(HIST("QA/h2UltraMicroTrackDCAxyVsPt"), track.pt(), track.dcaXY()); + qaRegistry.fill(HIST("QA/h2UltraMicroTrackDCAzVsPt"), track.pt(), track.dcaZ()); + qaRegistry.fill(HIST("QA/h4UltraMicroTrackTPCnSigma"), track.pt(), track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()); + if (track.hasTOF()) { + qaRegistry.fill(HIST("QA/h4UltraMicroTrackTOFnSigma"), track.pt(), track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()); + } + } + reso2ultramicrotrks(collision.globalIndex(), + px1000, + py1000, + pz1000, + pidFlag, + static_cast(dcaEncoding), + trackFlags); + if (!FilterForDerivedTables.cfgBypassTrackIndexFill) { + resoUltraMicroTrackTracks(track.globalIndex()); + } + } + } + /** * @brief Fills track data * @@ -820,14 +1311,21 @@ struct ResonanceDaughterInitializer { { // Loop over tracks for (auto const& track : tracks) { - if (!filterMicroTrack(track)) + if (!isMicroTrackSelected(collision, track)) { continue; - o2::aod::resomicrodaughter::ResoMicroTrackSelFlag trackSelFlag(track.dcaXY(), track.dcaZ()); - if (std::abs(track.dcaXY()) < (0.004 + (0.013 / track.pt()))) { - trackSelFlag.setDCAxy0(); } - if (std::abs(track.dcaZ()) < (0.004 + (0.013 / track.pt()))) { // TODO: check this - trackSelFlag.setDCAz0(); + if (!filterMicroTrack(track)) { + continue; + } + o2::aod::resomicrodaughter::ResoMicroTrackSelFlag trackSelFlag(track.dcaXY(), track.dcaZ()); + if (TrackCuts.cfgApplyTightDCAPtDepSelection) { + const float dcaThreshold = tightDCAThreshold(track.pt()); + if (dcaThreshold >= 0.f && std::abs(track.dcaXY()) < dcaThreshold) { + trackSelFlag.setDCAxy0(); + } + if (dcaThreshold >= 0.f && std::abs(track.dcaZ()) < dcaThreshold) { + trackSelFlag.setDCAz0(); + } } uint8_t trackFlags = (track.passedITSRefit() << 0) | (track.passedTPCRefit() << 1) | @@ -837,6 +1335,15 @@ struct ResonanceDaughterInitializer { (track.isPVContributor() << 5) | (track.hasTOF() << 6) | ((track.sign() > 0) << 7); // sign +1: 1, -1: 0 + if (cfgFillQA && cfgDetailTrackQA) { + qaRegistry.fill(HIST("QA/h4MicroTrackPtEtaPhi"), track.pt(), track.eta(), track.phi()); + qaRegistry.fill(HIST("QA/h2MicroTrackDCAxyVsPt"), track.pt(), track.dcaXY()); + qaRegistry.fill(HIST("QA/h2MicroTrackDCAzVsPt"), track.pt(), track.dcaZ()); + qaRegistry.fill(HIST("QA/h4MicroTrackTPCnSigma"), track.pt(), track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()); + if (track.hasTOF()) { + qaRegistry.fill(HIST("QA/h4MicroTrackTOFnSigma"), track.pt(), track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()); + } + } reso2microtrks(collision.globalIndex(), track.px(), track.py(), @@ -846,7 +1353,7 @@ struct ResonanceDaughterInitializer { static_cast(o2::aod::resomicrodaughter::PidNSigma(std::abs(track.tpcNSigmaPr()), std::abs(track.tofNSigmaPr()), track.hasTOF())), static_cast(trackSelFlag), trackFlags); - if (!cfgBypassTrackIndexFill) { + if (!FilterForDerivedTables.cfgBypassTrackIndexFill) { resoMicroTrackTracks(track.globalIndex()); } } @@ -864,18 +1371,16 @@ struct ResonanceDaughterInitializer { template void fillTracks(CollisionType const& collision, TrackType const& tracks) { - if (cfgBypassTrackFill) { + if (FilterForDerivedTables.cfgBypassTrackFill) { return; } // Loop over tracks for (auto const& track : tracks) { - if (!filterTrack(track)) + if (!isTrackSelected(collision, track)) { + continue; + } + if (!filterTrack(track)) { continue; - if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodTrackIndices"), 0); - qaRegistry.fill(HIST("QA/hTrackPt"), track.pt()); - qaRegistry.fill(HIST("QA/hTrackEta"), track.eta()); - qaRegistry.fill(HIST("QA/hTrackPhi"), track.phi()); } uint8_t trackFlags = (track.passedITSRefit() << 0) | (track.passedTPCRefit() << 1) | @@ -885,24 +1390,33 @@ struct ResonanceDaughterInitializer { (track.isPVContributor() << 5) | (track.hasTOF() << 6) | ((track.sign() > 0) << 7); // sign +1: 1, -1: 0 + if (cfgFillQA && cfgDetailTrackQA) { + qaRegistry.fill(HIST("QA/h4TrackPtEtaPhi"), track.pt(), track.eta(), track.phi()); + qaRegistry.fill(HIST("QA/h2TrackDCAxyVsPt"), track.pt(), track.dcaXY()); + qaRegistry.fill(HIST("QA/h2TrackDCAzVsPt"), track.pt(), track.dcaZ()); + qaRegistry.fill(HIST("QA/h4TrackTPCnSigma"), track.pt(), track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()); + if (track.hasTOF()) { + qaRegistry.fill(HIST("QA/h4TrackTOFnSigma"), track.pt(), track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()); + } + } reso2trks(collision.globalIndex(), track.pt(), track.px(), track.py(), track.pz(), - (uint8_t)track.tpcNClsCrossedRows(), - (uint8_t)track.tpcNClsFound(), - static_cast(track.dcaXY() * 10000), - static_cast(track.dcaZ() * 10000), - (int8_t)(track.tpcNSigmaPi() * 10), - (int8_t)(track.tpcNSigmaKa() * 10), - (int8_t)(track.tpcNSigmaPr() * 10), - (int8_t)(track.tofNSigmaPi() * 10), - (int8_t)(track.tofNSigmaKa() * 10), - (int8_t)(track.tofNSigmaPr() * 10), - (int8_t)(track.tpcSignal() * 10), + static_cast(track.tpcNClsCrossedRows()), + static_cast(track.tpcNClsFound()), + static_cast(std::round(track.dcaXY() * 10000)), + static_cast(std::round(track.dcaZ() * 10000)), + static_cast(std::round(track.tpcNSigmaPi() * 10)), + static_cast(std::round(track.tpcNSigmaKa() * 10)), + static_cast(std::round(track.tpcNSigmaPr() * 10)), + static_cast(std::round(track.tofNSigmaPi() * 10)), + static_cast(std::round(track.tofNSigmaKa() * 10)), + static_cast(std::round(track.tofNSigmaPr() * 10)), + static_cast(std::round(track.tpcSignal() * 100)), trackFlags); - if (!cfgBypassTrackIndexFill) { + if (!FilterForDerivedTables.cfgBypassTrackIndexFill) { resoTrackTracks(track.globalIndex()); } if constexpr (isMC) { @@ -943,7 +1457,7 @@ struct ResonanceDaughterInitializer { LOGF(debug, " mother index lMother: %d", lMother.globalIndex()); for (auto const& lDaughter : lMother.template daughters_as()) { LOGF(debug, " daughter index lDaughter: %d", lDaughter.globalIndex()); - if (lDaughter.globalIndex() != 0 && lDaughter.globalIndex() != theMcParticle.globalIndex()) { + if (lDaughter.globalIndex() != theMcParticle.globalIndex()) { lSiblingsIndeces.push_back(lDaughter.globalIndex()); } } @@ -953,12 +1467,9 @@ struct ResonanceDaughterInitializer { // ------ std::vector mothers = {-1, -1}; std::vector motherPDGs = {-1, -1}; - int siblings[2] = {0, 0}; - std::vector siblingsTemp = {-1, -1}; + std::array siblings{-1, -1}; + std::vector siblingsTemp{}; if (track.has_mcParticle()) { - if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodMCTrackIndices"), 0); - } // Get the MC particle const auto& particle = track.mcParticle(); if (particle.has_mothers()) { @@ -966,18 +1477,18 @@ struct ResonanceDaughterInitializer { motherPDGs = getMothersPDGCodes(particle); siblingsTemp = getSiblingsIndeces(particle); } - while (mothers.size() > 2) { - mothers.pop_back(); - motherPDGs.pop_back(); - } - if (siblingsTemp.size() > 0) + mothers.resize(StoredMCRelationCount, -1); + motherPDGs.resize(StoredMCRelationCount, -1); + if (!siblingsTemp.empty()) { siblings[0] = siblingsTemp[0]; - if (siblingsTemp.size() > 1) + } + if (siblingsTemp.size() > 1) { siblings[1] = siblingsTemp[1]; + } reso2mctracks(particle.pdgCode(), mothers[0], motherPDGs[0], - siblings, + siblings.data(), particle.isPhysicalPrimary(), particle.producedByGenerator()); } else { @@ -985,7 +1496,7 @@ struct ResonanceDaughterInitializer { reso2mctracks(0, mothers[0], motherPDGs[0], - siblings, + siblings.data(), 0, 0); } @@ -1003,49 +1514,52 @@ struct ResonanceDaughterInitializer { * @param tracks Track data */ template - void fillV0s(CollisionType const& collision, V0Type const& v0s, TrackType const&) + void fillV0s(CollisionType const& collision, V0Type const& v0s, TrackType const& tracks) { - int childIDs[2] = {0, 0}; // these IDs are necessary to keep track of the children for (auto const& v0 : v0s) { + if (!isV0Selected(collision, v0, tracks)) { + continue; + } + if (!filterV0(collision, v0, tracks)) { + continue; + } if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodV0Indices"), 0); qaRegistry.fill(HIST("QA/hV0Radius"), v0.v0radius()); qaRegistry.fill(HIST("QA/hV0CosPA"), v0.v0cosPA()); } - childIDs[0] = v0.posTrackId(); - childIDs[1] = v0.negTrackId(); + const std::array childIDs{v0.posTrackId(), v0.negTrackId()}; // Original track IDs for downstream pair-level shared-daughter rejection reso2v0s(collision.globalIndex(), v0.pt(), v0.px(), v0.py(), v0.pz(), - childIDs, + childIDs.data(), (int8_t)(v0.template posTrack_as().tpcNSigmaPi() * 10), (int8_t)(v0.template posTrack_as().tpcNSigmaKa() * 10), (int8_t)(v0.template posTrack_as().tpcNSigmaPr() * 10), (int8_t)(v0.template negTrack_as().tpcNSigmaPi() * 10), (int8_t)(v0.template negTrack_as().tpcNSigmaKa() * 10), (int8_t)(v0.template negTrack_as().tpcNSigmaPr() * 10), - (int8_t)(v0.template negTrack_as().tofNSigmaPi() * 10), - (int8_t)(v0.template negTrack_as().tofNSigmaKa() * 10), - (int8_t)(v0.template negTrack_as().tofNSigmaPr() * 10), (int8_t)(v0.template posTrack_as().tofNSigmaPi() * 10), (int8_t)(v0.template posTrack_as().tofNSigmaKa() * 10), (int8_t)(v0.template posTrack_as().tofNSigmaPr() * 10), + (int8_t)(v0.template negTrack_as().tofNSigmaPi() * 10), + (int8_t)(v0.template negTrack_as().tofNSigmaKa() * 10), + (int8_t)(v0.template negTrack_as().tofNSigmaPr() * 10), v0.v0cosPA(), v0.dcaV0daughters(), v0.dcapostopv(), v0.dcanegtopv(), v0.dcav0topv(), - 0, - 0, + static_cast(v0.template posTrack_as().tpcNClsCrossedRows()), + static_cast(v0.template negTrack_as().tpcNClsCrossedRows()), v0.mLambda(), v0.mAntiLambda(), v0.mK0Short(), v0.v0radius(), v0.x(), v0.y(), v0.z(), v0.alpha(), v0.qtarm()); - if (!cfgBypassTrackIndexFill) { - resoV0V0s(v0.globalIndex()); + if (!FilterForDerivedTables.cfgBypassTrackIndexFill) { + resoV0V0s(v0.v0Id()); } if constexpr (isMC) { fillMCV0(v0); @@ -1079,13 +1593,27 @@ struct ResonanceDaughterInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pt lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { LOGF(debug, " daughter index lDaughter: %d", lDaughter.globalIndex()); - if (lDaughter.globalIndex() != 0) { - lDaughtersIndeces.push_back(lDaughter.globalIndex()); - } + lDaughtersIndeces.push_back(lDaughter.globalIndex()); } return lDaughtersIndeces; }; @@ -1093,44 +1621,43 @@ struct ResonanceDaughterInitializer { std::vector lDaughtersPDGs{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { LOGF(debug, " daughter pdgcode lDaughter: %d", lDaughter.pdgCode()); - if (lDaughter.globalIndex() != 0) { - lDaughtersPDGs.push_back(lDaughter.pdgCode()); - } + lDaughtersPDGs.push_back(lDaughter.pdgCode()); } return lDaughtersPDGs; }; // ------ std::vector mothers = {-1, -1}; std::vector motherPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; if (v0.has_mcParticle()) { - if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodMCV0Indices"), 0); - } auto v0mc = v0.mcParticle(); if (v0mc.has_mothers()) { mothers = getMothersIndeces(v0mc); motherPDGs = getMothersPDGCodes(v0mc); + mothersPts = getMothersPt(v0mc); + mothersRaps = getMothersRap(v0mc); } - while (mothers.size() > 2) { - mothers.pop_back(); - motherPDGs.pop_back(); - } + mothers.resize(StoredMCRelationCount, -1); + motherPDGs.resize(StoredMCRelationCount, -1); + mothersPts.resize(StoredMCRelationCount, -1.0f); + mothersRaps.resize(StoredMCRelationCount, -1.0f); if (v0mc.has_daughters()) { daughters = getDaughtersIndeces(v0mc); daughterPDGs = getDaughtersPDGCodes(v0mc); } - while (daughters.size() > 2) { + if (daughters.size() > StoredMCRelationCount) { LOGF(info, "daughters.size() is larger than 2"); - daughters.pop_back(); - daughterPDGs.pop_back(); } + daughters.resize(StoredMCRelationCount, -1); + daughterPDGs.resize(StoredMCRelationCount, -1); reso2mcv0s(v0mc.pdgCode(), mothers[0], motherPDGs[0], - 0, - 0, + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1141,8 +1668,8 @@ struct ResonanceDaughterInitializer { reso2mcv0s(0, mothers[0], motherPDGs[0], - 0, - 0, + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1164,24 +1691,23 @@ struct ResonanceDaughterInitializer { * @param tracks Track data */ template - void fillCascades(CollisionType const& collision, CascType const& cascades, TrackType const&) + void fillCascades(CollisionType const& collision, CascType const& cascades, TrackType const& tracks) { - int childIDs[3] = {0, 0, 0}; // these IDs are necessary to keep track of the children for (auto const& casc : cascades) { + if (!isCascSelected(collision, casc, tracks)) { + continue; + } if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodCascIndices"), 0); qaRegistry.fill(HIST("QA/hCascRadius"), casc.cascradius()); qaRegistry.fill(HIST("QA/hCascCosPA"), casc.casccosPA(collision.posX(), collision.posY(), collision.posZ())); } - childIDs[0] = casc.posTrackId(); - childIDs[1] = casc.negTrackId(); - childIDs[2] = casc.bachelorId(); + const std::array childIDs{casc.posTrackId(), casc.negTrackId(), casc.bachelorId()}; // Original track IDs for downstream pair-level shared-daughter rejection reso2cascades(collision.globalIndex(), casc.pt(), casc.px(), casc.py(), casc.pz(), - childIDs, + childIDs.data(), (int8_t)(casc.template posTrack_as().tpcNSigmaPi() * 10), (int8_t)(casc.template posTrack_as().tpcNSigmaKa() * 10), (int8_t)(casc.template posTrack_as().tpcNSigmaPr() * 10), @@ -1211,14 +1737,14 @@ struct ResonanceDaughterInitializer { casc.dcaXYCascToPV(), casc.dcaZCascToPV(), casc.sign(), - 0, - 0, - 0, + static_cast(casc.template posTrack_as().tpcNClsCrossedRows()), + static_cast(casc.template negTrack_as().tpcNClsCrossedRows()), + static_cast(casc.template bachelor_as().tpcNClsCrossedRows()), casc.mLambda(), casc.mXi(), casc.v0radius(), casc.cascradius(), casc.x(), casc.y(), casc.z()); - if (!cfgBypassTrackIndexFill) { - resoCascadeCascades(casc.globalIndex()); + if (!FilterForDerivedTables.cfgBypassTrackIndexFill) { + resoCascadeCascades(casc.cascadeId()); } if constexpr (isMC) { fillMCCascade(casc); @@ -1252,13 +1778,27 @@ struct ResonanceDaughterInitializer { } return lMothersPDGs; }; + auto getMothersPt = [&](auto const& theMcParticle) { + std::vector lMothersPts{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pt lMother: %f", lMother.pt()); + lMothersPts.push_back(lMother.pt()); + } + return lMothersPts; + }; + auto getMothersRap = [&](auto const& theMcParticle) { + std::vector lMothersRaps{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother rap lMother: %f", lMother.y()); + lMothersRaps.push_back(lMother.y()); + } + return lMothersRaps; + }; auto getDaughtersIndeces = [&](auto const& theMcParticle) { std::vector lDaughtersIndeces{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { LOGF(debug, " daughter index lDaughter: %d", lDaughter.globalIndex()); - if (lDaughter.globalIndex() != 0) { - lDaughtersIndeces.push_back(lDaughter.globalIndex()); - } + lDaughtersIndeces.push_back(lDaughter.globalIndex()); } return lDaughtersIndeces; }; @@ -1266,9 +1806,7 @@ struct ResonanceDaughterInitializer { std::vector lDaughtersPDGs{}; for (auto const& lDaughter : theMcParticle.template daughters_as()) { LOGF(debug, " daughter pdgcode lDaughter: %d", lDaughter.pdgCode()); - if (lDaughter.globalIndex() != 0) { - lDaughtersPDGs.push_back(lDaughter.pdgCode()); - } + lDaughtersPDGs.push_back(lDaughter.pdgCode()); } return lDaughtersPDGs; }; @@ -1277,33 +1815,34 @@ struct ResonanceDaughterInitializer { std::vector motherPDGs = {-1, -1}; std::vector daughters = {-1, -1}; std::vector daughterPDGs = {-1, -1}; + std::vector mothersPts = {-1.0f, -1.0f}; + std::vector mothersRaps = {-1.0f, -1.0f}; if (casc.has_mcParticle()) { - if (cfgFillQA) { - qaRegistry.fill(HIST("QA/hGoodMCCascIndices"), 0); - } auto cascmc = casc.mcParticle(); if (cascmc.has_mothers()) { mothers = getMothersIndeces(cascmc); motherPDGs = getMothersPDGCodes(cascmc); + mothersPts = getMothersPt(cascmc); + mothersRaps = getMothersRap(cascmc); } - while (mothers.size() > 2) { - mothers.pop_back(); - motherPDGs.pop_back(); - } + mothers.resize(StoredMCRelationCount, -1); + motherPDGs.resize(StoredMCRelationCount, -1); + mothersPts.resize(StoredMCRelationCount, -1.0f); + mothersRaps.resize(StoredMCRelationCount, -1.0f); if (cascmc.has_daughters()) { daughters = getDaughtersIndeces(cascmc); daughterPDGs = getDaughtersPDGCodes(cascmc); } - while (daughters.size() > 2) { + if (daughters.size() > StoredMCRelationCount) { LOGF(info, "daughters.size() is larger than 2"); - daughters.pop_back(); - daughterPDGs.pop_back(); } + daughters.resize(StoredMCRelationCount, -1); + daughterPDGs.resize(StoredMCRelationCount, -1); reso2mccascades(cascmc.pdgCode(), mothers[0], motherPDGs[0], - 0, - 0, + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1314,8 +1853,8 @@ struct ResonanceDaughterInitializer { reso2mccascades(0, mothers[0], motherPDGs[0], - 0, - 0, + mothersPts[0], + mothersRaps[0], daughters[0], daughters[1], daughterPDGs[0], @@ -1335,22 +1874,66 @@ struct ResonanceDaughterInitializer { } PROCESS_SWITCH(ResonanceDaughterInitializer, processDummy, "Process dummy", true); + /** + * @brief Fills track tables for one original collision + * + * @tparam isMC Boolean indicating if it's MC + * @tparam CollisionType Type of reduced collision + * @tparam TrackTableType Type of input track table + * @tparam PresliceType Type of the original-collision preslice + * @param collision Reduced collision with the original collision index + * @param tracks Input track table + */ + template + void fillTrackTablesForCollision(CollisionType const& collision, TrackTableType const& tracks, PresliceType const& perCollision) + { + auto tracksThisCollision = tracks.sliceBy(perCollision, collision.collisionId()); + fillTracks(collision, tracksThisCollision); + if (FilterForDerivedTables.cfgFillMicroTracks) { + fillMicroTracks(collision, tracksThisCollision); + } + if (FilterForDerivedTables.cfgFillUltraMicroTracks) { + fillUltraMicroTracks(collision, tracksThisCollision); + } + } + /** * @brief Processes data tracks * * @param collision Collision data * @param tracks Track data */ - void processData(aod::ResoCollision const& collision, + void processData(ResoCollisionWithIndex::iterator const& collision, soa::Filtered const& tracks) { - fillTracks(collision, tracks); - if (cfgFillMicroTracks) { - fillMicroTracks(collision, tracks); - } + fillTrackTablesForCollision(collision, tracks, tracksPerCollision); } PROCESS_SWITCH(ResonanceDaughterInitializer, processData, "Process tracks for data", false); + /** + * @brief Processes data tracks with configurable selected-V0 and selected-cascade gates + */ + void processDataWithPairGate(ResoCollisionWithIndex::iterator const& collision, + soa::Filtered const& tracks, + aod::ResoV0Candidates const& v0s, + aod::ResoCascadesCandidates const& cascades) + { + if (FilterForDerivedTables.cfgBypassNoPairV0s) { + auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, collision.collisionId()); + if (!hasSelectedV0(collision, v0sThisCollision, tracks)) { + return; + } + } + if (FilterForDerivedTables.cfgBypassNoPairCascades) { + auto cascadesThisCollision = cascades.sliceBy(cascadesPerCollision, collision.collisionId()); + if (!hasSelectedCascade(collision, cascadesThisCollision, tracks)) { + return; + } + } + fillTrackTablesForCollision(collision, tracks, tracksPerCollision); + } + PROCESS_SWITCH(ResonanceDaughterInitializer, processDataWithPairGate, "Process data tracks with configurable pair gates", false); + /** * @brief Processes MC tracks * @@ -1358,17 +1941,35 @@ struct ResonanceDaughterInitializer { * @param tracks Track data * @param mcParticles MC particles */ - void processMC(aod::ResoCollision const& collision, + void processMC(ResoCollisionWithIndex::iterator const& collision, soa::Filtered const& tracks, aod::McParticles const&) { - fillTracks(collision, tracks); - if (cfgFillMicroTracks) { - fillMicroTracks(collision, tracks); - } + fillTrackTablesForCollision(collision, tracks, tracksMCPerCollision); } PROCESS_SWITCH(ResonanceDaughterInitializer, processMC, "Process tracks for MC", false); + /** + * @brief Processes MC tracks with configurable V0 and cascade candidate gates + */ + void processMCWithPairGate(ResoCollisionWithIndex::iterator const& collision, + soa::Filtered const& tracks, + aod::ResoV0CandidatesMC const& v0s, + aod::ResoCascadesCandidatesMC const& cascades, + aod::McParticles const&) + { + auto v0sThisCollision = v0s.sliceBy(v0sMCPerCollision, collision.collisionId()); + if (FilterForDerivedTables.cfgBypassNoPairV0s && v0sThisCollision.size() < 1) { + return; + } + auto cascadesThisCollision = cascades.sliceBy(cascadesMCPerCollision, collision.collisionId()); + if (FilterForDerivedTables.cfgBypassNoPairCascades && cascadesThisCollision.size() < 1) { + return; + } + fillTrackTablesForCollision(collision, tracks, tracksMCPerCollision); + } + PROCESS_SWITCH(ResonanceDaughterInitializer, processMCWithPairGate, "Process MC tracks with configurable pair gates", false); + /** * @brief Processes V0 data * @@ -1376,9 +1977,10 @@ struct ResonanceDaughterInitializer { * @param v0s V0 data * @param tracks Track data */ - void processV0Data(aod::ResoCollision const& collision, aod::ResoV0Candidates const& v0s, aod::ResoTrackCandidates const& tracks) + void processV0Data(ResoCollisionWithIndex::iterator const& collision, aod::ResoV0Candidates const& v0s, aod::ResoTrackCandidates const& tracks) { - fillV0s(collision, v0s, tracks); + auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, collision.collisionId()); + fillV0s(collision, v0sThisCollision, tracks); } PROCESS_SWITCH(ResonanceDaughterInitializer, processV0Data, "Process V0s for data", false); @@ -1389,9 +1991,10 @@ struct ResonanceDaughterInitializer { * @param v0s V0 data * @param tracks Track data */ - void processV0MC(aod::ResoCollision const& collision, aod::ResoV0CandidatesMC const& v0s, aod::ResoTrackCandidatesMC const& tracks) + void processV0MC(ResoCollisionWithIndex::iterator const& collision, aod::ResoV0CandidatesMC const& v0s, aod::ResoTrackCandidatesMC const& tracks, aod::McParticles const&) { - fillV0s(collision, v0s, tracks); + auto v0sThisCollision = v0s.sliceBy(v0sMCPerCollision, collision.collisionId()); + fillV0s(collision, v0sThisCollision, tracks); } PROCESS_SWITCH(ResonanceDaughterInitializer, processV0MC, "Process V0s for MC", false); @@ -1402,9 +2005,10 @@ struct ResonanceDaughterInitializer { * @param cascades Cascade data * @param tracks Track data */ - void processCascData(aod::ResoCollision const& collision, aod::ResoCascadesCandidates const& cascades, aod::ResoTrackCandidates const& tracks) + void processCascData(ResoCollisionWithIndex::iterator const& collision, aod::ResoCascadesCandidates const& cascades, aod::ResoTrackCandidates const& tracks) { - fillCascades(collision, cascades, tracks); + auto cascadesThisCollision = cascades.sliceBy(cascadesPerCollision, collision.collisionId()); + fillCascades(collision, cascadesThisCollision, tracks); } PROCESS_SWITCH(ResonanceDaughterInitializer, processCascData, "Process Cascades for data", false); @@ -1415,16 +2019,17 @@ struct ResonanceDaughterInitializer { * @param cascades Cascade data * @param tracks Track data */ - void processCascMC(aod::ResoCollision const& collision, aod::ResoCascadesCandidatesMC const& cascades, aod::ResoTrackCandidatesMC const& tracks) + void processCascMC(ResoCollisionWithIndex::iterator const& collision, aod::ResoCascadesCandidatesMC const& cascades, aod::ResoTrackCandidatesMC const& tracks, aod::McParticles const&) { - fillCascades(collision, cascades, tracks); + auto cascadesThisCollision = cascades.sliceBy(cascadesMCPerCollision, collision.collisionId()); + fillCascades(collision, cascadesThisCollision, tracks); } PROCESS_SWITCH(ResonanceDaughterInitializer, processCascMC, "Process Cascades for MC", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(context), + adaptAnalysisTask(context)}; }