diff --git a/PWGCF/EbyEFluctuations/Tasks/partNumFluc.cxx b/PWGCF/EbyEFluctuations/Tasks/partNumFluc.cxx index 51a2474a8a9..d459822b071 100644 --- a/PWGCF/EbyEFluctuations/Tasks/partNumFluc.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/partNumFluc.cxx @@ -57,11 +57,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include @@ -70,16 +70,16 @@ #include #include -#define C_CS(cs) /* NOLINT(cppcoreguidelines-macro-usage) */ \ - [](std::index_sequence) constexpr -> ConstStr<(cs)[indices]...> { \ - static_assert(std::is_array_v> && \ - std::same_as>, char>); \ - return ConstStr<(cs)[indices]...>{}; \ +#define C_CS(cs) /* NOLINT(cppcoreguidelines-macro-usage) */ \ + [](const std::index_sequence) constexpr -> ConstStr<(cs)[indices]...> { \ + static_assert(std::is_array_v> && \ + std::same_as>, char>); \ + return ConstStr<(cs)[indices]...>{}; \ }(std::make_index_sequence{}) -#define C_SV(sv) /* NOLINT(cppcoreguidelines-macro-usage) */ \ - [](std::index_sequence) constexpr -> ConstStr<(sv)[indices]...> { \ - static_assert(std::same_as, std::string_view>); \ - return ConstStr<(sv)[indices]...>{}; \ +#define C_SV(sv) /* NOLINT(cppcoreguidelines-macro-usage) */ \ + [](const std::index_sequence) constexpr -> ConstStr<(sv)[indices]...> { \ + static_assert(std::same_as, std::string_view>); \ + return ConstStr<(sv)[indices]...>{}; \ }(std::make_index_sequence<(sv).size()>{}) using namespace o2; @@ -88,49 +88,74 @@ using namespace o2::framework::expressions; namespace o2::aod { +using JoinedMcCollisions = soa::Join; using JoinedCollisions = soa::Join; -using JoinedTracks = soa::Join; using JoinedCollisionsWithMc = soa::Join; +using JoinedTracks = soa::Join; using JoinedTracksWithMc = soa::Join; -using JoinedMcCollisions = soa::Join; - -namespace mini_mc_collision -{ -DECLARE_SOA_COLUMN(NMcParticlesP, nMcParticlesP, std::uint16_t); -DECLARE_SOA_COLUMN(NMcParticlesM, nMcParticlesM, std::uint16_t); -} // namespace mini_mc_collision - -DECLARE_SOA_TABLE(MiniMcCollisions, "AOD", "MINIMCCOLLISION", soa::Index<>, mini_mc_collision::NMcParticlesP, mini_mc_collision::NMcParticlesM); -using MiniMcCollision = MiniMcCollisions::iterator; namespace mini_collision { -DECLARE_SOA_COLUMN(Vz, vz, std::int8_t); -DECLARE_SOA_COLUMN(Centrality, centrality, std::uint16_t); -DECLARE_SOA_COLUMN(NTracksP, nTracksP, std::uint16_t); -DECLARE_SOA_COLUMN(NTracksM, nTracksM, std::uint16_t); +DECLARE_SOA_COLUMN(Code, code, std::uint16_t); } // namespace mini_collision -DECLARE_SOA_TABLE(MiniCollisions, "AOD", "MINICOLLISION", soa::Index<>, mini_collision::Vz, mini_collision::Centrality, mini_collision::NTracksP, mini_collision::NTracksM); +DECLARE_SOA_TABLE(MiniCollisions, "AOD", "MINICOLLISION", soa::Index<>, mini_collision::Code); using MiniCollision = MiniCollisions::iterator; namespace mini_mc_particle { -DECLARE_SOA_INDEX_COLUMN(MiniMcCollision, miniMcCollision); -DECLARE_SOA_COLUMN(SignedEfficiency, signedEfficiency, std::int16_t); +DECLARE_SOA_INDEX_COLUMN(MiniCollision, miniCollision); +DECLARE_SOA_COLUMN(Code, code, std::uint16_t); } // namespace mini_mc_particle -DECLARE_SOA_TABLE(MiniMcParticles, "AOD", "MINIMCPARTICLE", soa::Index<>, mini_mc_particle::MiniMcCollisionId, mini_mc_particle::SignedEfficiency); +DECLARE_SOA_TABLE(MiniMcParticles, "AOD", "MINIMCPARTICLE", soa::Index<>, mini_mc_particle::MiniCollisionId, mini_mc_particle::Code); using MiniMcParticle = MiniMcParticles::iterator; namespace mini_track { DECLARE_SOA_INDEX_COLUMN(MiniCollision, miniCollision); -DECLARE_SOA_COLUMN(SignedEfficiency, signedEfficiency, std::int16_t); +DECLARE_SOA_COLUMN(Code, code, std::uint32_t); } // namespace mini_track -DECLARE_SOA_TABLE(MiniTracks, "AOD", "MINITRACK", soa::Index<>, mini_track::MiniCollisionId, mini_track::SignedEfficiency); +DECLARE_SOA_TABLE(MiniTracks, "AOD", "MINITRACK", soa::Index<>, mini_track::MiniCollisionId, mini_track::Code); using MiniTrack = MiniTracks::iterator; + +namespace tiny_mc_collision +{ +DECLARE_SOA_COLUMN(NMcParticlesP, nMcParticlesP, std::uint16_t); +DECLARE_SOA_COLUMN(NMcParticlesM, nMcParticlesM, NMcParticlesP::type); +} // namespace tiny_mc_collision + +DECLARE_SOA_TABLE(TinyMcCollisions, "AOD", "TINYMCCOLLISION", soa::Index<>, tiny_mc_collision::NMcParticlesP, tiny_mc_collision::NMcParticlesM); +using TinyMcCollision = TinyMcCollisions::iterator; + +namespace tiny_collision +{ +DECLARE_SOA_COLUMN(Code, code, mini_collision::Code::type); +DECLARE_SOA_COLUMN(NTracksP, nTracksP, std::uint16_t); +DECLARE_SOA_COLUMN(NTracksM, nTracksM, NTracksP::type); +} // namespace tiny_collision + +DECLARE_SOA_TABLE(TinyCollisions, "AOD", "TINYCOLLISION", soa::Index<>, tiny_collision::Code, tiny_collision::NTracksP, tiny_collision::NTracksM); +using TinyCollision = TinyCollisions::iterator; + +namespace tiny_mc_particle +{ +DECLARE_SOA_INDEX_COLUMN(TinyMcCollision, tinyMcCollision); +DECLARE_SOA_COLUMN(SignedEfficiency, signedEfficiency, std::int16_t); +} // namespace tiny_mc_particle + +DECLARE_SOA_TABLE(TinyMcParticles, "AOD", "TINYMCPARTICLE", soa::Index<>, tiny_mc_particle::TinyMcCollisionId, tiny_mc_particle::SignedEfficiency); +using TinyMcParticle = TinyMcParticles::iterator; + +namespace tiny_track +{ +DECLARE_SOA_INDEX_COLUMN(TinyCollision, tinyCollision); +DECLARE_SOA_COLUMN(SignedEfficiency, signedEfficiency, std::int16_t); +} // namespace tiny_track + +DECLARE_SOA_TABLE(TinyTracks, "AOD", "TINYTRACK", soa::Index<>, tiny_track::TinyCollisionId, tiny_track::SignedEfficiency); +using TinyTrack = TinyTracks::iterator; } // namespace o2::aod namespace @@ -160,10 +185,7 @@ inline constexpr std::int32_t NOrderKeys{[]() consteval noexcept -> std::int32_t return std::accumulate(counts.begin(), counts.end(), 0); }()}; inline constexpr std::array, NOrderKeys> OrderKeys{[]() consteval noexcept -> std::array, NOrderKeys> { - std::array, NOrderKeys> result{}; - std::array current{}; - std::int32_t index{}; - const auto fillOrderKeys{[¤t, &index, &result](const auto& self, const std::int32_t position, const std::int32_t sum, const std::int32_t target) consteval noexcept -> void { + constexpr auto FillOrderKeys{[](const auto& self, std::array& current, std::int32_t& index, std::array, NOrderKeys>& result, const std::int32_t position, const std::int32_t sum, const std::int32_t target) consteval noexcept -> void { if (position == NExponentKeys) { if (sum == target) { result[index++] = current; @@ -172,13 +194,17 @@ inline constexpr std::array, NOrderKeys> } const std::int32_t weight{ExponentKeys[position].first}; - for (std::int32_t power{}; sum + power * weight <= target; ++power) { + for (std::int32_t const& power : std::views::iota(0, (target - sum) / weight + 1)) { current[position] = static_cast(power); - self(self, position + 1, sum + power * weight, target); + self(self, current, index, result, position + 1, sum + power * weight, target); } }}; + + std::array, NOrderKeys> result{}; + std::array current{}; + std::int32_t index{}; for (std::int32_t const& target : std::views::iota(0, MaxOrder + 1)) { - fillOrderKeys(fillOrderKeys, 0, 0, target); + FillOrderKeys(FillOrderKeys, current, index, result, 0, 0, target); } return result; }()}; @@ -205,6 +231,32 @@ inline constexpr std::array bool { + if (OrderKeys[0] != std::array{}) { + return false; + } + for (std::int32_t const& iOrderKey : std::views::iota(1, NOrderKeys)) { + const auto& [iOrderKeyParent, factor]{OrderKeyProductLinks[iOrderKey]}; + const auto& [iExponentKey, power]{factor}; + if (iOrderKeyParent < 0 || iOrderKeyParent >= iOrderKey || iExponentKey < 0 || iExponentKey >= NExponentKeys || power <= 0) { + return false; + } + for (std::int32_t const& iExponentKeyAfter : std::views::iota(iExponentKey + 1, NExponentKeys)) { + if (OrderKeys[iOrderKey][iExponentKeyAfter] > 0) { + return false; + } + } + std::array orderKeyReconstructed{OrderKeys[iOrderKeyParent]}; + if (orderKeyReconstructed[iExponentKey] != 0) { + return false; + } + orderKeyReconstructed[iExponentKey] = power; + if (orderKeyReconstructed != OrderKeys[iOrderKey]) { + return false; + } + } + return true; +}()); } // namespace fluctuation_calculator_base class FluctuationCalculatorTrack @@ -259,23 +311,18 @@ class FluctuationCalculatorTrack template concept IsValidEnum = ((std::is_enum_v> && requires { std::remove_cvref_t::N; }) && ...); template -concept IsValid = IsValidEnum...> && ((EValues != std::remove_cvref_t::N) && ...); +concept IsValidEnumValue = IsValidEnum...> && ((EValues != std::remove_cvref_t::N) && ...); template requires IsValidEnum -constexpr std::int32_t toI(const E e) +constexpr std::int32_t toI(const E e) noexcept { return static_cast(e); } - template requires IsValidEnum constexpr std::int32_t NEs{toI(E::N)}; -template - requires IsValidEnum -struct EnumInfo; - enum class NameKind { Default = 0, Lower, @@ -283,43 +330,44 @@ enum class NameKind { DisplayLower, N }; -template - requires IsValidEnum -constexpr std::string_view getName(const std::int32_t index) -{ - if constexpr (NameKindValue == NameKind::Lower) { - return EnumInfo>::NamesLower.at(index); - } else if constexpr (NameKindValue == NameKind::Display) { - return EnumInfo>::DisplayNames.at(index); - } else if constexpr (NameKindValue == NameKind::DisplayLower) { - return EnumInfo>::DisplayNamesLower.at(index); - } else { - return EnumInfo>::Names.at(index); - } -} -template - requires IsValidEnum -constexpr std::string_view getName(const E e) -{ - return getName(toI(e)); -} -template - requires IsValidEnum -std::vector getDisplayNames() -{ - if constexpr (requires { EnumInfo>::DisplayNames; }) { - return {EnumInfo>::DisplayNames.begin(), EnumInfo>::DisplayNames.end()}; - } else { - return {EnumInfo>::Names.begin(), EnumInfo>::Names.end()}; - } -} - enum class DataMode { McMcParticle = 0, McTrack, RawTrack, N }; +enum class Selection { + Loose = 0, + Default, + Tight, + N +}; +enum class RangeEdge { + Min = 0, + Max, + N +}; +enum class McEventSelection { + All = 0, + Good, + Inel, + Vz, + NRecoCollisions, + N +}; +enum class EventSelection { + All = 0, + Good, + Run, + Rct, + Inel, + Bits, + Centrality, + Vz, + Occupancy, + NPvContributors, + N +}; enum class CentralityDefinition { Ft0a = 0, Ft0c, @@ -366,17 +414,17 @@ enum class ParticleSpeciesAll { Proton, N }; +enum class ChargeSpecies { + Plus = 0, + Minus, + N +}; enum class ParticleNumber { Charge = 0, Kaon, Proton, N }; -enum class ChargeSpecies { - Plus = 0, - Minus, - N -}; enum class ChargeNumber { Plus = 0, Minus, @@ -385,6 +433,17 @@ enum class ChargeNumber { N }; +template + requires IsValidEnum +struct EnumInfo; +template <> +struct EnumInfo { + static constexpr std::array> Names{"Loose", "Default", "Tight"}; +}; +template <> +struct EnumInfo { + static constexpr std::array> Names{"Min", "Max"}; +}; template <> struct EnumInfo { static constexpr std::array> Names{"Mean", "Sigma"}; @@ -428,23 +487,54 @@ struct EnumInfo { static constexpr std::array> ParticleSpeciesValues{ParticleSpecies::N, ParticleSpecies::Pion, ParticleSpecies::Kaon, ParticleSpecies::Proton}; }; template <> +struct EnumInfo { + static constexpr std::array> Names{"P", "M"}; + static constexpr std::array> NamesLower{"p", "m"}; + static constexpr std::array> Titles{"+", "#minus"}; +}; +template <> struct EnumInfo { static constexpr std::array> Names{"Ch", "Ka", "Pr"}; static constexpr std::array> DisplayNames{"Charge", "Kaon", "Proton"}; static constexpr std::array> DisplayNamesLower{"charge", "kaon", "proton"}; static constexpr std::array> Titles{"h", "K", "p"}; -}; -template <> -struct EnumInfo { - static constexpr std::array> Names{"P", "M"}; - static constexpr std::array> NamesLower{"p", "m"}; - static constexpr std::array> Titles{"+", "#minus"}; + static constexpr std::array> ParticleSpeciesAllValues{ParticleSpeciesAll::All, ParticleSpeciesAll::Kaon, ParticleSpeciesAll::Proton}; }; template <> struct EnumInfo { static constexpr std::array> Names{"P", "M", "T", "N"}; }; +template + requires IsValidEnum +constexpr std::string_view getName(const std::int32_t index) +{ + if constexpr (NameKindValue == NameKind::Lower) { + return EnumInfo>::NamesLower.at(index); + } else if constexpr (NameKindValue == NameKind::Display) { + return EnumInfo>::DisplayNames.at(index); + } else if constexpr (NameKindValue == NameKind::DisplayLower) { + return EnumInfo>::DisplayNamesLower.at(index); + } else { + return EnumInfo>::Names.at(index); + } +} +template + requires IsValidEnum +constexpr std::string_view getName(const E e) +{ + return getName(toI(e)); +} +template + requires IsValidEnum +std::vector getDisplayNames() +{ + if constexpr (requires { EnumInfo>::DisplayNames; }) { + return {EnumInfo>::DisplayNames.begin(), EnumInfo>::DisplayNames.end()}; + } else { + return {EnumInfo>::Names.begin(), EnumInfo>::Names.end()}; + } +} template requires IsValidEnum constexpr std::string_view getTitle(const std::int32_t index) @@ -453,7 +543,7 @@ constexpr std::string_view getTitle(const std::int32_t index) } template requires IsValidEnum -constexpr To getValue(const E e) +constexpr To getValue(const E e) noexcept { if constexpr (std::is_same_v, Detector>) { return EnumInfo::DetectorValues[toI(e)]; @@ -467,15 +557,183 @@ constexpr To getValue(const E e) return To::N; } } -constexpr std::int32_t getPdgCode(const ParticleSpecies particleSpecies, const ChargeSpecies chargeSpecies) +constexpr std::int32_t getPdgCode(const ParticleSpecies particleSpecies, const ChargeSpecies chargeSpecies) noexcept { return EnumInfo::PdgCodes[toI(particleSpecies)][toI(chargeSpecies)]; } -constexpr double getMass(const ParticleSpecies particleSpecies) +constexpr double getMass(const ParticleSpecies particleSpecies) noexcept { return EnumInfo::Masses[toI(particleSpecies)]; } +namespace mini_collision_codec +{ +inline constexpr std::int32_t NBinsVz{40}; +inline constexpr std::array> RangeVz{-10., 10.}; +inline constexpr std::int32_t NBinsCentrality{550}; +inline constexpr std::array BinEdgesCentrality{[]() consteval noexcept -> std::array { + constexpr std::array EdgesSelection{0., 0.001, 0.01, 0.1, 1., 10., 100.}; + constexpr std::array NBinsPerSection{100, 90, 90, 90, 90, 90}; + std::array edges{}; + std::int32_t index{}; + edges[index++] = EdgesSelection.front(); + for (std::int32_t const& iSection : std::views::iota(0, static_cast(NBinsPerSection.size()))) { + const double width{(EdgesSelection[iSection + 1] - EdgesSelection[iSection]) / NBinsPerSection[iSection]}; + for (std::int32_t const& iBin : std::views::iota(1, NBinsPerSection[iSection] + 1)) { + edges[index++] = EdgesSelection[iSection] + iBin * width; + } + } + return edges; +}()}; +inline constexpr std::int32_t NStatesGoodNPvContributors{2}; +static_assert(static_cast(NBinsVz) * NBinsCentrality * NStatesGoodNPvContributors <= static_cast(std::numeric_limits::max()) + 1); + +std::optional encode(const double vz, const double centrality, const bool isGoodNPvContributors) noexcept +{ + if (!std::isfinite(vz) || !std::isfinite(centrality) || vz < RangeVz[toI(RangeEdge::Min)] || RangeVz[toI(RangeEdge::Max)] <= vz || centrality < BinEdgesCentrality.front() || BinEdgesCentrality.back() <= centrality) { + return std::nullopt; + } + + const std::int32_t vzBinIndex{static_cast(std::floor((vz - RangeVz[toI(RangeEdge::Min)]) / ((RangeVz[toI(RangeEdge::Max)] - RangeVz[toI(RangeEdge::Min)]) / NBinsVz)))}; + if (vzBinIndex < 0 || NBinsVz <= vzBinIndex) { + return std::nullopt; + } + + aod::mini_collision::Code::type code{static_cast(vzBinIndex)}; + aod::mini_collision::Code::type codeSpan{NBinsVz}; + code += (static_cast(std::ranges::upper_bound(BinEdgesCentrality, centrality) - BinEdgesCentrality.begin()) - 1) * codeSpan; + codeSpan *= NBinsCentrality; + code += static_cast(isGoodNPvContributors) * codeSpan; + return code; +} +} // namespace mini_collision_codec + +namespace mini_track_codec +{ +inline constexpr std::int32_t NStatesPvContributor{2}; +inline constexpr std::int32_t NBinsPt{18}; +inline constexpr std::int32_t NBinsEta{16}; +inline constexpr std::array> RangePt{0.2, 2.}; +inline constexpr std::array> RangeEta{-0.8, 0.8}; +static_assert(static_cast(NStatesPvContributor) * NEs * NEs * NEs * NEs * NEs * NEs * NEs * NEs * NBinsPt * NBinsEta * NEs * NEs * NEs <= static_cast(std::numeric_limits::max()) + 1); + +struct ConfigSelection { + std::array> minItsNCls{}; + std::array> maxItsChi2NCls{}; + std::array> minTpcNCls{}; + std::array> maxTpcChi2NCls{}; + std::array> minTpcNCrossedRows{}; + std::array> maxTpcNClsSharedRatio{}; + std::array>, NEs> maxAbsNSigmaDca{}; + std::array>, NEs> maxAbsNSigmaPid{}; +}; + +std::optional encode(const ConfigSelection& configSelection, const bool isPvContributor, const double itsNCls, const double itsChi2NCls, const double tpcNCls, const double tpcChi2NCls, const double tpcNCrossedRows, const double tpcNClsSharedRatio, const std::array>& absNSigmaDca, const double pt, const double eta, const std::int32_t sign, const ParticleSpeciesAll particleSpeciesAll, const double absNSigmaPid) noexcept +{ + if (!std::isfinite(pt) || !std::isfinite(eta) || pt < RangePt[toI(RangeEdge::Min)] || RangePt[toI(RangeEdge::Max)] <= pt || eta < RangeEta[toI(RangeEdge::Min)] || RangeEta[toI(RangeEdge::Max)] <= eta) { + return std::nullopt; + } + + const std::int32_t ptBinIndex{static_cast(std::floor((pt - RangePt[toI(RangeEdge::Min)]) / ((RangePt[toI(RangeEdge::Max)] - RangePt[toI(RangeEdge::Min)]) / NBinsPt)))}; + if (ptBinIndex < 0 || NBinsPt <= ptBinIndex) { + return std::nullopt; + } + + const std::int32_t etaBinIndex{static_cast(std::floor((eta - RangeEta[toI(RangeEdge::Min)]) / ((RangeEta[toI(RangeEdge::Max)] - RangeEta[toI(RangeEdge::Min)]) / NBinsEta)))}; + if (etaBinIndex < 0 || NBinsEta <= etaBinIndex) { + return std::nullopt; + } + + static constexpr auto GetState{ + [] + requires IsValidEnumValue + (const double value, const std::array>& cuts) noexcept -> std::int32_t { + for (std::int32_t const& iSelection : std::views::iota(0, NEs) | std::views::reverse) { + if constexpr (RangeEdgeValue == RangeEdge::Min) { + if (value > cuts[iSelection]) { + return iSelection; + } + } else { + if (value < cuts[iSelection]) { + return iSelection; + } + } + } + return -1; + }}; + + aod::mini_track::Code::type code{static_cast(isPvContributor)}; + aod::mini_track::Code::type codeSpan{NStatesPvContributor}; + for (std::int32_t const& state : std::to_array( + // cppcheck-suppress internalAstError + {GetState.template operator()(itsNCls, configSelection.minItsNCls), + GetState.template operator()(itsChi2NCls, configSelection.maxItsChi2NCls), + GetState.template operator()(tpcNCls, configSelection.minTpcNCls), + GetState.template operator()(tpcChi2NCls, configSelection.maxTpcChi2NCls), + GetState.template operator()(tpcNCrossedRows, configSelection.minTpcNCrossedRows), + GetState.template operator()(tpcNClsSharedRatio, configSelection.maxTpcNClsSharedRatio), + GetState.template operator()(absNSigmaDca[toI(DcaAxis::Xy)], configSelection.maxAbsNSigmaDca[toI(DcaAxis::Xy)]), + GetState.template operator()(absNSigmaDca[toI(DcaAxis::Z)], configSelection.maxAbsNSigmaDca[toI(DcaAxis::Z)])})) { + if (state < 0) { + return std::nullopt; + } + + code += state * codeSpan; + codeSpan *= NEs; + } + code += ptBinIndex * codeSpan; + codeSpan *= NBinsPt; + code += etaBinIndex * codeSpan; + codeSpan *= NBinsEta; + code += (sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)) * codeSpan; + codeSpan *= NEs; + code += toI(particleSpeciesAll) * codeSpan; + codeSpan *= NEs; + if (particleSpeciesAll != ParticleSpeciesAll::All) { + // clang-format v20.1.3 + // clang-format off + const std::int32_t statePid{GetState.template operator()(absNSigmaPid, configSelection.maxAbsNSigmaPid[toI(getValue(particleSpeciesAll))])}; + // clang-format on + if (statePid < 0) { + return std::nullopt; + } + code += statePid * codeSpan; + } + return code; +} +} // namespace mini_track_codec + +namespace mini_mc_particle_codec +{ +static_assert(static_cast(mini_track_codec::NBinsPt) * mini_track_codec::NBinsEta * NEs * NEs <= static_cast(std::numeric_limits::max()) + 1); + +std::optional encode(const double pt, const double eta, const std::int32_t sign, const ParticleSpecies particleSpecies) noexcept +{ + if (!std::isfinite(pt) || !std::isfinite(eta) || pt < mini_track_codec::RangePt[toI(RangeEdge::Min)] || mini_track_codec::RangePt[toI(RangeEdge::Max)] <= pt || eta < mini_track_codec::RangeEta[toI(RangeEdge::Min)] || mini_track_codec::RangeEta[toI(RangeEdge::Max)] <= eta) { + return std::nullopt; + } + + const std::int32_t ptBinIndex{static_cast(std::floor((pt - mini_track_codec::RangePt[toI(RangeEdge::Min)]) / ((mini_track_codec::RangePt[toI(RangeEdge::Max)] - mini_track_codec::RangePt[toI(RangeEdge::Min)]) / mini_track_codec::NBinsPt)))}; + if (ptBinIndex < 0 || mini_track_codec::NBinsPt <= ptBinIndex) { + return std::nullopt; + } + + const std::int32_t etaBinIndex{static_cast(std::floor((eta - mini_track_codec::RangeEta[toI(RangeEdge::Min)]) / ((mini_track_codec::RangeEta[toI(RangeEdge::Max)] - mini_track_codec::RangeEta[toI(RangeEdge::Min)]) / mini_track_codec::NBinsEta)))}; + if (etaBinIndex < 0 || mini_track_codec::NBinsEta <= etaBinIndex) { + return std::nullopt; + } + + aod::mini_mc_particle::Code::type code{static_cast(ptBinIndex)}; + aod::mini_mc_particle::Code::type codeSpan{mini_track_codec::NBinsPt}; + code += etaBinIndex * codeSpan; + codeSpan *= mini_track_codec::NBinsEta; + code += (sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)) * codeSpan; + codeSpan *= NEs; + code += toI(particleSpecies) * codeSpan; + return code; +} +} // namespace mini_mc_particle_codec + template requires std::is_arithmetic_v constexpr std::int32_t nEnabled(const Configurable>& cfg) @@ -496,36 +754,36 @@ double interpolate(const TH3* const h, const double x, const double y, const dou return 0.; } - static constexpr auto GetBinIndicesWeights{[](const TAxis* const axis, const double position) -> std::array, 2> { + static constexpr auto GetBinIndicesWeights{[](const TAxis* const axis, const double position) -> std::array, NEs> { if (!axis) { return {}; } - const std::int32_t n{axis->GetNbins()}; - if (n == 1 || position <= axis->GetBinCenter(1)) { + const std::int32_t nBins{axis->GetNbins()}; + if (nBins == 1 || position <= axis->GetBinCenter(1)) { return {{{1, 1.}, {1, 0.}}}; } - if (position >= axis->GetBinCenter(n)) { - return {{{n, 1.}, {n, 0.}}}; + if (position >= axis->GetBinCenter(nBins)) { + return {{{nBins, 1.}, {nBins, 0.}}}; } - const std::int32_t bin{axis->FindFixBin(position)}; - const std::int32_t lower{position < axis->GetBinCenter(bin) ? bin - 1 : bin}; - const std::int32_t upper{lower + 1}; - const double fraction{(position - axis->GetBinCenter(lower)) / (axis->GetBinCenter(upper) - axis->GetBinCenter(lower))}; + const std::int32_t binIndex{axis->FindFixBin(position)}; + const std::int32_t binIndexLower{position < axis->GetBinCenter(binIndex) ? binIndex - 1 : binIndex}; + const std::int32_t binIndexUpper{binIndexLower + 1}; + const double fraction{(position - axis->GetBinCenter(binIndexLower)) / (axis->GetBinCenter(binIndexUpper) - axis->GetBinCenter(binIndexLower))}; - return {{{lower, 1. - fraction}, {upper, fraction}}}; + return {{{binIndexLower, 1. - fraction}, {binIndexUpper, fraction}}}; }}; - const std::array, 2> xb{GetBinIndicesWeights(h->GetXaxis(), x)}; - const std::array, 2> yb{GetBinIndicesWeights(h->GetYaxis(), y)}; - const std::array, 2> zb{GetBinIndicesWeights(h->GetZaxis(), z)}; + const std::array, NEs> binIndicesWeightsX{GetBinIndicesWeights(h->GetXaxis(), x)}; + const std::array, NEs> binIndicesWeightsY{GetBinIndicesWeights(h->GetYaxis(), y)}; + const std::array, NEs> binIndicesWeightsZ{GetBinIndicesWeights(h->GetZaxis(), z)}; double result{}; - for (const auto& [ix, wx] : xb) { - for (const auto& [iy, wy] : yb) { - for (const auto& [iz, wz] : zb) { - result += wx * wy * wz * h->GetBinContent(ix, iy, iz); + for (const auto& [iBinX, weightX] : binIndicesWeightsX) { + for (const auto& [iBinY, weightY] : binIndicesWeightsY) { + for (const auto& [iBinZ, weightZ] : binIndicesWeightsZ) { + result += weightX * weightY * weightZ * h->GetBinContent(iBinX, iBinY, iBinZ); } } } @@ -535,14 +793,14 @@ double interpolate(const TH3* const h, const double x, const double y, const dou struct PartNumFluc { struct HolderCcdb { - [[maybe_unused]] static constexpr std::int32_t NDimensionsEfficiency{4}; + static constexpr std::int32_t NDimensionsEfficiency{4}; const TList* lCcdb{}; std::map> runNumbersIndicesGroupIndices; std::int32_t runGroupIndexCurrent{}; - std::array, NEs>, NEs>, NEs> fPtMeasureDca{}; - std::array>, NEs>, NEs> hCentralityPtEtaShiftNSigmaPid{}; - std::array>, NEs>, NEs> hVzCentralityPtEtaEfficiency{}; + std::array, NEs>, NEs>, NEs> calibrationsPtMeasureDca{}; + std::array>, NEs>, NEs> hsCentralityPtEtaShiftNSigmaPid{}; + std::array>, NEs>, NEs> hsVzCentralityPtEtaEfficiency{}; } holderCcdb{}; struct HolderMcEvent { @@ -550,12 +808,12 @@ struct PartNumFluc { std::array>, NEs> numbers{}; std::array>, NEs> numbersEff{}; - void clear() { *this = {}; } + void clear() noexcept { *this = {}; } } holderMcEvent{}; struct HolderEvent { - static constexpr std::pair RangeCentrality{0., 100.}; - static constexpr bool isValidCentrality(const double value) { return RangeCentrality.first <= value && value <= RangeCentrality.second; } + static constexpr std::array> RangeCentrality{0., 100.}; + static constexpr bool isValidCentrality(const double value) noexcept { return RangeCentrality[toI(RangeEdge::Min)] <= value && value < RangeCentrality[toI(RangeEdge::Max)]; } std::int32_t runNumber{}; std::int32_t runIndex{}; @@ -568,14 +826,14 @@ struct PartNumFluc { double centralityCalibration{}; double centrality{}; std::int32_t subgroupIndex{}; - std::array>, NEs> numbers{}; + std::array>, NEs> numbers{}; - void clear() { *this = {}; } - [[nodiscard]] std::int32_t getNGlobalTracks() const { return std::accumulate(nGlobalTracks.begin(), nGlobalTracks.end(), 0); } - [[nodiscard]] std::int32_t getNPvContributors() const { return std::accumulate(nPvContributors.begin(), nPvContributors.end(), 0); } + void clear() noexcept { *this = {}; } + [[nodiscard]] std::int32_t getNGlobalTracks() const noexcept { return std::accumulate(nGlobalTracks.begin(), nGlobalTracks.end(), 0); } + [[nodiscard]] std::int32_t getNPvContributors() const noexcept { return std::accumulate(nPvContributors.begin(), nPvContributors.end(), 0); } template - requires IsValid - [[nodiscard]] double getMeasureDca() const + requires IsValidEnumValue + [[nodiscard]] double getMeasureDca() const noexcept { const std::int32_t sumNGlobalTracks{getNGlobalTracks()}; if (sumNGlobalTracks == 0) { @@ -596,7 +854,7 @@ struct PartNumFluc { return sumDca / sumNGlobalTracks; } } - [[nodiscard]] std::int32_t getNTofBeta() const { return std::accumulate(nTofBeta.begin(), nTofBeta.end(), 0); } + [[nodiscard]] std::int32_t getNTofBeta() const noexcept { return std::accumulate(nTofBeta.begin(), nTofBeta.end(), 0); } } holderEvent{}; struct HolderMcParticle { @@ -605,18 +863,18 @@ struct PartNumFluc { double pt{}; double eta{}; - void clear() { *this = {}; } + void clear() noexcept { *this = {}; } } holderMcParticle{}; struct HolderTrack { static constexpr double TruncationAbsNSigmaPid{999.}; - static constexpr double truncateNSigmaPid(const double value, const double shift = {}) + static constexpr double truncateNSigmaPid(const double value, const double shift = {}) noexcept { const double valueShifted{value - shift}; return std::abs(value) < TruncationAbsNSigmaPid && std::abs(valueShifted) < TruncationAbsNSigmaPid ? valueShifted : -TruncationAbsNSigmaPid; } - [[nodiscard]] double getNSigmaPidCombined(const std::int32_t particleSpeciesIndex) const + [[nodiscard]] double getNSigmaPidCombined(const std::int32_t particleSpeciesIndex) const noexcept { return truncateNSigmaPid(std::copysign(std::hypot(nSigmaPid[toI(Detector::Tpc)][particleSpeciesIndex], nSigmaPid[toI(Detector::Tof)][particleSpeciesIndex]), nSigmaPid[toI(Detector::Tpc)][particleSpeciesIndex] + nSigmaPid[toI(Detector::Tof)][particleSpeciesIndex])); } @@ -635,27 +893,22 @@ struct PartNumFluc { return a; }()}; - void clear() { *this = {}; } + void clear() noexcept { *this = {}; } } holderTrack{}; struct HolderDerivedData { template - static constexpr T convertFloor(const double value) - { - return std::numeric_limits::lowest() <= value && value <= std::numeric_limits::max() ? static_cast(std::floor(value)) : (std::is_signed_v ? std::numeric_limits::lowest() : std::numeric_limits::max()); - } - template - static constexpr T convertRound(const double value) + static constexpr T convert(const double value) noexcept { - return std::numeric_limits::lowest() <= value && value <= std::numeric_limits::max() ? static_cast(std::round(value)) : (std::is_signed_v ? std::numeric_limits::lowest() : std::numeric_limits::max()); + return std::numeric_limits::lowest() <= value && value <= std::numeric_limits::max() ? static_cast(std::rint(value)) : (std::is_signed_v ? std::numeric_limits::lowest() : std::numeric_limits::max()); } - std::array> nMcParticles{}; - std::array> nTracks{}; - std::vector signedEfficienciesMcParticle{[]() constexpr -> std::vector { std::vector v{}; v.reserve(256); return v; }()}; - std::vector signedEfficienciesTrack{[]() constexpr -> std::vector { std::vector v{}; v.reserve(256); return v; }()}; + std::array> nMcParticles{}; + std::array> nTracks{}; + std::vector signedEfficienciesMcParticle{[]() constexpr -> std::vector { std::vector v{}; v.reserve(256); return v; }()}; + std::vector signedEfficienciesTrack{[]() constexpr -> std::vector { std::vector v{}; v.reserve(256); return v; }()}; - void clear() + void clear() noexcept { nMcParticles = {}; nTracks = {}; @@ -664,15 +917,15 @@ struct PartNumFluc { } } holderDerivedData{}; - std::array, NEs>, NEs> fluctuationCalculatorsTrack{}; - struct : ConfigurableGroup { + std::string prefix{"cgCcdb"}; Configurable cfgCcdbUrl{"cfgCcdbUrl", "http://ccdb-test.cern.ch:8080", "Url of CCDB"}; Configurable cfgCcdbPath{"cfgCcdbPath", "Users/f/fasi/test", "Path in CCDB"}; Configurable cfgCcdbTimestampLatest{"cfgCcdbTimestampLatest", std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(), "Latest timestamp in CCDB"}; - } groupCcdb{}; + } cgCcdb{}; struct : ConfigurableGroup { + std::string prefix{"cgAnalysis"}; Configurable cfgFlagQaRun{"cfgFlagQaRun", false, "Run QA flag"}; Configurable cfgFlagQaEvent{"cfgFlagQaEvent", false, "Event QA flag"}; Configurable cfgFlagQaCentrality{"cfgFlagQaCentrality", false, "Centrality QA flag"}; @@ -686,9 +939,11 @@ struct PartNumFluc { Configurable> cfgFlagsCalculationPurity{"cfgFlagsCalculationPurity", {std::array>{}.data(), NEs, getDisplayNames()}, "Purity calculation flags"}; Configurable> cfgFlagsCalculationFractionPrimary{"cfgFlagsCalculationFractionPrimary", {std::array>{}.data(), NEs, getDisplayNames()}, "Primary fraction calculation flags"}; Configurable> cfgFlagsCalculationFluctuation{"cfgFlagsCalculationFluctuation", {std::array>{}.data(), NEs, getDisplayNames()}, "Fluctuation calculation flags"}; - } groupAnalysis{}; + Configurable> cfgFlagsStorageMiniTable{"cfgFlagsStorageMiniTable", {std::array>{}.data(), NEs, getDisplayNames()}, "Mini table storage flags"}; + } cgAnalysis{}; struct : ConfigurableGroup { + std::string prefix{"cgEvent"}; Configurable cfgFlagRejectionRunBad{"cfgFlagRejectionRunBad", false, "Bad run rejection flag"}; Configurable cfgFlagRejectionRunBadMc{"cfgFlagRejectionRunBadMc", false, "MC bad run rejection flag"}; Configurable cfgLabelFlagsRct{"cfgLabelFlagsRct", "CBT_hadronPID", "RCT flags label"}; @@ -697,84 +952,105 @@ struct PartNumFluc { Configurable cfgFlagInelEvent{"cfgFlagInelEvent", true, "Flag of requiring inelastic event"}; Configurable cfgFlagInelEventMc{"cfgFlagInelEventMc", false, "Flag of requiring inelastic MC event"}; Configurable cfgCutMaxAbsVz{"cfgCutMaxAbsVz", 8., "Maximum absolute z-vertex position (cm)"}; - Configurable cfgCutMaxAbsVzMc{"cfgCutMaxAbsVzMc", 999., "Maximum absolute MC z-vertex position (cm)"}; + Configurable cfgCutMaxAbsVzMc{"cfgCutMaxAbsVzMc", 8., "Maximum absolute MC z-vertex position (cm)"}; + Configurable cfgCutMaxOccupancy{"cfgCutMaxOccupancy", 0, "Maximum occupancy"}; Configurable cfgCutMinDeviationNPvContributors{"cfgCutMinDeviationNPvContributors", -4, "Minimum nPvContributors deviation from nGlobalTracks"}; Configurable cfgIndexDefinitionCentrality{"cfgIndexDefinitionCentrality", 2, "Centrality definition index"}; Configurable cfgFlagDefinitionCentralitySameQa{"cfgFlagDefinitionCentralitySameQa", false, "Flag of using the same centrality definition for QA"}; - ConfigurableAxis cfgAxisCentrality{"cfgAxisCentrality", {18, 0., 90.}, "Centrality axis in fluctuation calculation"}; + Configurable cfgNMultiplicityBinsAxis{"cfgNMultiplicityBinsAxis", 200, "Number of multiplicity bins for axis"}; ConfigurableAxis cfgAxisCentralityCalibration{"cfgAxisCentralityCalibration", {VARIABLE_WIDTH, 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100.}, "Centrality axis in calibration"}; + ConfigurableAxis cfgAxisCentrality{"cfgAxisCentrality", {20, 0., 100.}, "Centrality axis in fluctuation calculation"}; Configurable cfgNSubgroups{"cfgNSubgroups", 20, "Number of subgroups in fluctuation calculation"}; + Configurable cfgFactorStorageMiniTable{"cfgFactorStorageMiniTable", 1., "Mini table storage inverse probability factor"}; Configurable cfgFlagSingleCollisionMc{"cfgFlagSingleCollisionMc", false, "Flag of requiring exactly single collision of MC collision"}; - Configurable cfgFlagBestCollisionMc{"cfgFlagBestCollisionMc", false, "Flag of requiring best collision of MC collision"}; Configurable cfgFlagMcCollisionVz{"cfgFlagMcCollisionVz", false, "Flag of using z-vertex position of MC collision"}; - } groupEvent{}; + } cgEvent{}; struct : ConfigurableGroup { + std::string prefix{"cgTrack"}; Configurable cfgFlagPvContributor{"cfgFlagPvContributor", true, "Flag of requiring PV contributor"}; - Configurable cfgCutMinItsNCls{"cfgCutMinItsNCls", 5, "Minimum number of clusters ITS"}; - Configurable cfgCutMaxItsChi2NCls{"cfgCutMaxItsChi2NCls", 25., "Maximum chi2 per cluster ITS"}; - Configurable cfgCutMinTpcNCls{"cfgCutMinTpcNCls", 55, "Minimum number of clusters TPC"}; + Configurable> cfgCutsMinItsNCls{"cfgCutsMinItsNCls", {std::array>{4, 5, 6}.data(), NEs, getDisplayNames()}, "Minimum numbers of clusters ITS"}; + Configurable> cfgCutsMaxItsChi2NCls{"cfgCutsMaxItsChi2NCls", {std::array>{30., 25., 20.}.data(), NEs, getDisplayNames()}, "Maximum chi2 per cluster ITS"}; + Configurable> cfgCutsMinTpcNCls{"cfgCutsMinTpcNCls", {std::array>{45, 55, 65}.data(), NEs, getDisplayNames()}, "Minimum numbers of clusters TPC"}; Configurable cfgCutMinTpcChi2NCls{"cfgCutMinTpcChi2NCls", 0., "Minimum chi2 per cluster TPC"}; - Configurable cfgCutMaxTpcChi2NCls{"cfgCutMaxTpcChi2NCls", 3.5, "Maximum chi2 per cluster TPC"}; - Configurable cfgCutMaxTpcNClsSharedRatio{"cfgCutMaxTpcNClsSharedRatio", 0.4, "Maximum ratio of shared clusters over clusters TPC"}; - Configurable cfgCutMinTpcNCrossedRows{"cfgCutMinTpcNCrossedRows", 80, "Minimum number of crossed rows TPC"}; + Configurable> cfgCutsMaxTpcChi2NCls{"cfgCutsMaxTpcChi2NCls", {std::array>{4., 3.5, 3.}.data(), NEs, getDisplayNames()}, "Maximum chi2 per cluster TPC"}; + Configurable> cfgCutsMaxTpcNClsSharedRatio{"cfgCutsMaxTpcNClsSharedRatio", {std::array>{0.5, 0.4, 0.3}.data(), NEs, getDisplayNames()}, "Maximum ratios of shared clusters over clusters TPC"}; + Configurable> cfgCutsMinTpcNCrossedRows{"cfgCutsMinTpcNCrossedRows", {std::array>{70, 80, 90}.data(), NEs, getDisplayNames()}, "Minimum numbers of crossed rows TPC"}; Configurable cfgCutMinTpcNCrossedRowsRatio{"cfgCutMinTpcNCrossedRowsRatio", 0.8, "Minimum ratio of crossed rows over findable clusters TPC"}; Configurable cfgFlagRecalibrationDca{"cfgFlagRecalibrationDca", false, "DCA recalibration flag"}; - Configurable> cfgCutsMaxAbsNSigmaDca{"cfgCutsMaxAbsNSigmaDca", {std::array>{2.5, 2.5}.data(), NEs, getDisplayNames()}, "Maximum absolute nSigma values of DCA (cm)"}; - Configurable cfgCutMinPt{"cfgCutMinPt", 0.4, "Minimum pT (GeV/c)"}; - Configurable cfgCutMaxPt{"cfgCutMaxPt", 2., "Maximum pT (GeV/c)"}; + Configurable> cfgCutsMaxAbsNSigmaDca{"cfgCutsMaxAbsNSigmaDca", {std::array * NEs>{3., 2.5, 2., 3., 2.5, 2.}.data(), NEs, NEs, getDisplayNames(), getDisplayNames()}, "Maximum absolute nSigma values of DCA (cm)"}; + Configurable> cfgCutsRangePt{"cfgCutsRangePt", {std::array * NEs>{0.2, 2., 0.3, 2., 0.4, 2.}.data(), NEs, NEs, getDisplayNames(), getDisplayNames()}, "pT ranges (GeV/c)"}; Configurable cfgCutMaxAbsEta{"cfgCutMaxAbsEta", 0.8, "Maximum absolute eta"}; - Configurable> cfgThresholdsPtTofPid{"cfgThresholdsPtTofPid", {std::array>{0.5, 0.5, 0.8}.data(), NEs, getDisplayNames()}, "pT (GeV/c) thresholds for TOF PID"}; + Configurable> cfgThresholdsPtTofPid{"cfgThresholdsPtTofPid", {std::array>{0.4, 0.4, 0.8}.data(), NEs, getDisplayNames()}, "pT (GeV/c) thresholds for TOF PID"}; Configurable> cfgFlagsRecalibrationNSigmaPid{"cfgFlagsRecalibrationNSigmaPid", {std::array>{}.data(), NEs, getDisplayNames()}, "nSigma PID recalibration flags"}; Configurable cfgFlagRejectionOthers{"cfgFlagRejectionOthers", false, "Other particle species rejection flag"}; - Configurable> cfgCutsMaxAbsNSigmaPid{"cfgCutsMaxAbsNSigmaPid", {std::array>{2., 2., 2.}.data(), NEs, getDisplayNames()}, "Maximum absolute nSigma values for PID"}; + Configurable> cfgCutsMaxAbsNSigmaPid{"cfgCutsMaxAbsNSigmaPid", {std::array * NEs>{2.5, 2., 1.5, 2.5, 2., 1.5, 2.5, 2., 1.5}.data(), NEs, NEs, getDisplayNames(), getDisplayNames()}, "Maximum absolute nSigma values for PID"}; Configurable cfgFlagMcParticlePhysicalPrimary{"cfgFlagMcParticlePhysicalPrimary", true, "Flag of requiring physical primary MC particle"}; Configurable cfgFlagMcParticleMomentum{"cfgFlagMcParticleMomentum", true, "Flag of using momentum of MC particle"}; - } groupTrack{}; - - bool doQaAcceptance{}; - bool doQaPhi{}; - bool doQaPid{}; - bool doCalculationYield{}; - bool doCalculationPurity{}; - bool doCalculationFractionPrimary{}; - bool doCalculationFluctuation{}; - - HistogramRegistry hrCalculationFluctuation{"hrCalculationFluctuation", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrCalculationFractionPrimary{"hrCalculationFractionPrimary", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrCalculationPurity{"hrCalculationPurity", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrCalculationYield{"hrCalculationYield", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaMc{"hrQaMc", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaPid{"hrQaPid", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaPhi{"hrQaPhi", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaAcceptance{"hrQaAcceptance", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaDca{"hrQaDca", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaTrack{"hrQaTrack", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaCentrality{"hrQaCentrality", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaEvent{"hrQaEvent", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrQaRun{"hrQaRun", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry hrCounter{"hrCounter", {}, OutputObjHandlingPolicy::AnalysisObject}; - - aod::rctsel::RCTFlagsChecker rctFlagsChecker; + } cgTrack{}; Service pdg{}; Service ccdb{}; + aod::rctsel::RCTFlagsChecker rctFlagsChecker; + Filter filterCollision{aod::evsel::sel8 == true}; Filter filterTrack{requireQualityTracksInFilter() && requireTrackCutInFilter(TrackSelectionFlags::kGoldenChi2)}; Filter filterMcCollision{aod::mccollisionprop::numRecoCollision > 0}; Preslice presliceTracksPerCollision{aod::track::collisionId}; + PresliceUnsorted presliceTracksPerMcParticle{aod::mctracklabel::mcParticleId}; - Produces miniMcCollision{}; Produces miniCollision{}; Produces miniMcParticle{}; Produces miniTrack{}; + Produces tinyMcCollision{}; + Produces tinyCollision{}; + Produces tinyMcParticle{}; + Produces tinyTrack{}; + + HistogramRegistry hrCalculationFluctuation{"CalculationFluctuation", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrCalculationFractionPrimary{"CalculationFractionPrimary", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrCalculationPurity{"CalculationPurity", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrCalculationYield{"CalculationYield", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaMc{"QaMc", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaPid{"QaPid", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaPhi{"QaPhi", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaAcceptance{"QaAcceptance", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaDca{"QaDca", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaTrack{"QaTrack", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaCentrality{"QaCentrality", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaEvent{"QaEvent", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrQaRun{"QaRun", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry hrCounter{"Counter", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + bool doQaAcceptance{}; + bool doQaPhi{}; + bool doQaPid{}; + bool doCalculationYield{}; + bool doCalculationPurity{}; + bool doCalculationFractionPrimary{}; + bool doCalculationFluctuation{}; + bool doStorageMiniTable{}; + std::array> rangePtAll{std::numeric_limits::max(), std::numeric_limits::lowest()}; + std::unique_ptr configSelection; + + std::array, NEs>, NEs> fluctuationCalculatorsTrackMcParticle{}; + std::array, NEs>, NEs> fluctuationCalculatorsTrackTrack{}; void init(const InitContext&) { gRandom->SetSeed(0); + doQaAcceptance = isEnabled(cgAnalysis.cfgFlagsQaAcceptance); + doQaPhi = isEnabled(cgAnalysis.cfgFlagsQaPhi); + doQaPid = isEnabled(cgAnalysis.cfgFlagsQaPid); + doCalculationYield = isEnabled(cgAnalysis.cfgFlagsCalculationYield); + doCalculationPurity = isEnabled(cgAnalysis.cfgFlagsCalculationPurity); + doCalculationFractionPrimary = isEnabled(cgAnalysis.cfgFlagsCalculationFractionPrimary); + doCalculationFluctuation = isEnabled(cgAnalysis.cfgFlagsCalculationFluctuation); + doStorageMiniTable = isEnabled(cgAnalysis.cfgFlagsStorageMiniTable); + if (doProcessRaw.value == doProcessMc.value) { LOG(fatal) << "Identical values of doProcessRaw and doProcessMc!"; } @@ -783,30 +1059,98 @@ struct PartNumFluc { } else { LOG(info) << "Enabling raw data process."; } + if (nEnabled(cgAnalysis.cfgFlagsCalculationFluctuation) > 1) { + LOG(fatal) << "Invalid " << cgAnalysis.cfgFlagsCalculationFluctuation.name << "!"; + } + if (nEnabled(cgAnalysis.cfgFlagsStorageMiniTable) > 1) { + LOG(fatal) << "Invalid " << cgAnalysis.cfgFlagsStorageMiniTable.name << "!"; + } + if ((cgAnalysis.cfgFlagQaEvent.value || cgAnalysis.cfgFlagQaCentrality.value || doCalculationFluctuation) && cgEvent.cfgNMultiplicityBinsAxis.value <= 0) { + LOG(fatal) << "Invalid " << cgEvent.cfgNMultiplicityBinsAxis.name << "!"; + } + if (doCalculationFluctuation && cgEvent.cfgNSubgroups.value <= 0) { + LOG(fatal) << "Invalid " << cgEvent.cfgNSubgroups.name << "!"; + } + + if (doStorageMiniTable) { + configSelection = std::make_unique(); + for (std::int32_t const& iSelection : std::views::iota(0, NEs)) { + configSelection->minItsNCls[iSelection] = cgTrack.cfgCutsMinItsNCls.value.get(iSelection); + configSelection->maxItsChi2NCls[iSelection] = cgTrack.cfgCutsMaxItsChi2NCls.value.get(iSelection); + configSelection->minTpcNCls[iSelection] = cgTrack.cfgCutsMinTpcNCls.value.get(iSelection); + configSelection->maxTpcChi2NCls[iSelection] = cgTrack.cfgCutsMaxTpcChi2NCls.value.get(iSelection); + configSelection->minTpcNCrossedRows[iSelection] = cgTrack.cfgCutsMinTpcNCrossedRows.value.get(iSelection); + configSelection->maxTpcNClsSharedRatio[iSelection] = cgTrack.cfgCutsMaxTpcNClsSharedRatio.value.get(iSelection); + } + for (std::int32_t const& iDcaAxis : std::views::iota(0, NEs)) { + for (std::int32_t const& iSelection : std::views::iota(0, NEs)) { + configSelection->maxAbsNSigmaDca[iDcaAxis][iSelection] = cgTrack.cfgCutsMaxAbsNSigmaDca.value.get(getName(iDcaAxis).data(), getName(iSelection).data()); + } + } + for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { + for (std::int32_t const& iSelection : std::views::iota(0, NEs)) { + configSelection->maxAbsNSigmaPid[iParticleSpecies][iSelection] = cgTrack.cfgCutsMaxAbsNSigmaPid.value.get(getName(iParticleSpecies).data(), getName(iSelection).data()); + } + } + + static constexpr auto ValidateCuts{ + [] + requires IsValidEnumValue + (const auto& cuts, const auto& name) constexpr -> void { + for (std::int32_t const& iSelection : std::views::iota(0, NEs)) { + if (!std::isfinite(cuts[iSelection])) { + LOG(fatal) << "Non-finite value in " << name << "!"; + } + if (iSelection == 0) { + continue; + } + if constexpr (RangeEdgeValue == RangeEdge::Min) { + if (!(cuts[iSelection - 1] < cuts[iSelection])) { + LOG(fatal) << "Values in " << name << " must satisfy " << getName(Selection::Loose) << " < " << getName(Selection::Default) << " < " << getName(Selection::Tight) << "!"; + } + } else { + if (!(cuts[iSelection - 1] > cuts[iSelection])) { + LOG(fatal) << "Values in " << name << " must satisfy " << getName(Selection::Loose) << " > " << getName(Selection::Default) << " > " << getName(Selection::Tight) << "!"; + } + } + } + }}; + ValidateCuts.template operator()(configSelection->minItsNCls, cgTrack.cfgCutsMinItsNCls.name); + ValidateCuts.template operator()(configSelection->maxItsChi2NCls, cgTrack.cfgCutsMaxItsChi2NCls.name); + ValidateCuts.template operator()(configSelection->minTpcNCls, cgTrack.cfgCutsMinTpcNCls.name); + ValidateCuts.template operator()(configSelection->maxTpcChi2NCls, cgTrack.cfgCutsMaxTpcChi2NCls.name); + ValidateCuts.template operator()(configSelection->minTpcNCrossedRows, cgTrack.cfgCutsMinTpcNCrossedRows.name); + ValidateCuts.template operator()(configSelection->maxTpcNClsSharedRatio, cgTrack.cfgCutsMaxTpcNClsSharedRatio.name); + for (std::int32_t const& iDcaAxis : std::views::iota(0, NEs)) { + ValidateCuts.template operator()(configSelection->maxAbsNSigmaDca[iDcaAxis], cgTrack.cfgCutsMaxAbsNSigmaDca.name); + } + for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { + ValidateCuts.template operator()(configSelection->maxAbsNSigmaPid[iParticleSpecies], cgTrack.cfgCutsMaxAbsNSigmaPid.name); + } + } - doQaAcceptance = isEnabled(groupAnalysis.cfgFlagsQaAcceptance); - doQaPhi = isEnabled(groupAnalysis.cfgFlagsQaPhi); - doQaPid = isEnabled(groupAnalysis.cfgFlagsQaPid); - doCalculationYield = isEnabled(groupAnalysis.cfgFlagsCalculationYield); - doCalculationPurity = isEnabled(groupAnalysis.cfgFlagsCalculationPurity); - doCalculationFractionPrimary = isEnabled(groupAnalysis.cfgFlagsCalculationFractionPrimary); - doCalculationFluctuation = isEnabled(groupAnalysis.cfgFlagsCalculationFluctuation); + if (doStorageMiniTable || static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Charge)))) { + for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { + rangePtAll[toI(RangeEdge::Min)] = std::min(rangePtAll[toI(RangeEdge::Min)], cgTrack.cfgCutsRangePt.value.get(iParticleSpecies, toI(RangeEdge::Min))); + rangePtAll[toI(RangeEdge::Max)] = std::max(rangePtAll[toI(RangeEdge::Max)], cgTrack.cfgCutsRangePt.value.get(iParticleSpecies, toI(RangeEdge::Max))); + } + } - rctFlagsChecker.init(groupEvent.cfgLabelFlagsRct.value, static_cast(groupEvent.cfgFlagsRct.value.get("ZDC")), static_cast(groupEvent.cfgFlagsRct.value.get("Acceptance")), static_cast(groupEvent.cfgFlagsRct.value.get("Table"))); + rctFlagsChecker.init(cgEvent.cfgLabelFlagsRct.value, static_cast(cgEvent.cfgFlagsRct.value.get("ZDC")), static_cast(cgEvent.cfgFlagsRct.value.get("Acceptance")), static_cast(cgEvent.cfgFlagsRct.value.get("Table"))); - ccdb->setURL(groupCcdb.cfgCcdbUrl.value); + ccdb->setURL(cgCcdb.cfgCcdbUrl.value); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(true); - if (groupCcdb.cfgCcdbTimestampLatest.value >= 0) { - ccdb->setCreatedNotAfter(groupCcdb.cfgCcdbTimestampLatest.value); + if (cgCcdb.cfgCcdbTimestampLatest.value >= 0) { + ccdb->setCreatedNotAfter(cgCcdb.cfgCcdbTimestampLatest.value); } readCcdb(); std::int32_t nRunsBad{}; - for (const auto& runNumberIndexGroupIndex : holderCcdb.runNumbersIndicesGroupIndices) { - const std::int32_t runGroupIndex{runNumberIndexGroupIndex.second.second}; - if (runGroupIndex == 0 || (groupEvent.cfgFlagRejectionRunBad.value && runGroupIndex < 0)) { + for (const auto& [runNumber, runIndexGroupIndex] : holderCcdb.runNumbersIndicesGroupIndices) { + const std::int32_t runGroupIndex{runIndexGroupIndex.second}; + if (runGroupIndex == 0 || (cgEvent.cfgFlagRejectionRunBad.value && runGroupIndex < 0)) { ++nRunsBad; } } @@ -821,7 +1165,7 @@ struct PartNumFluc { LOG(info) << "Number of bad runs: " << nRunsBad; } for (const auto& [runNumber, runIndexGroupIndex] : holderCcdb.runNumbersIndicesGroupIndices) { - if (runIndexGroupIndex.second == 0 || (groupEvent.cfgFlagRejectionRunBad.value && runIndexGroupIndex.second < 0)) { + if (runIndexGroupIndex.second == 0 || (cgEvent.cfgFlagRejectionRunBad.value && runIndexGroupIndex.second < 0)) { LOG(info) << "Enabling rejecting run: " << runNumber << " (" << runIndexGroupIndex.second << ")"; } else { LOG(info) << "Enabling processing run: " << runNumber << " (" << std::abs(runIndexGroupIndex.second) << ")"; @@ -829,32 +1173,32 @@ struct PartNumFluc { } } - if (groupEvent.cfgLabelFlagsRct.value.empty()) { + if (cgEvent.cfgLabelFlagsRct.value.empty()) { LOG(info) << "No RCT flags label enabled."; } else { - LOG(info) << "Enabling RCT flags label: " << groupEvent.cfgLabelFlagsRct.value; + LOG(info) << "Enabling RCT flags label: " << cgEvent.cfgLabelFlagsRct.value; } - if (static_cast(groupEvent.cfgFlagsRct.value.get("ZDC"))) { + if (static_cast(cgEvent.cfgFlagsRct.value.get("ZDC"))) { LOG(info) << "Enabling RCT flag: ZDC"; } - if (static_cast(groupEvent.cfgFlagsRct.value.get("Acceptance"))) { + if (static_cast(cgEvent.cfgFlagsRct.value.get("Acceptance"))) { LOG(info) << "Enabling RCT flag: acceptance"; } - if (static_cast(groupEvent.cfgFlagsRct.value.get("Table"))) { + if (static_cast(cgEvent.cfgFlagsRct.value.get("Table"))) { LOG(info) << "Enabling RCT flag: table"; } - if ((groupEvent.cfgBitsSelectionEvent.value & ((std::uint64_t{1} << aod::evsel::EventSelectionFlags::kNsel) - 1)) == 0) { + if ((cgEvent.cfgBitsSelectionEvent.value & ((std::uint64_t{1} << aod::evsel::EventSelectionFlags::kNsel) - 1)) == 0) { LOG(info) << "No event selection bit enabled."; } else { - for (std::int32_t const& iEvSel : std::views::iota(0, aod::evsel::EventSelectionFlags::kNsel)) { - if (static_cast((groupEvent.cfgBitsSelectionEvent.value >> iEvSel) & 1)) { - LOG(info) << "Enabling event selection bit: " << aod::evsel::selectionLabels[iEvSel]; + for (std::int32_t const& iBit : std::views::iota(0, aod::evsel::EventSelectionFlags::kNsel)) { + if (static_cast((cgEvent.cfgBitsSelectionEvent.value >> iBit) & 1)) { + LOG(info) << "Enabling event selection bit: " << aod::evsel::selectionLabels[iBit]; } } } - switch (groupEvent.cfgIndexDefinitionCentrality.value) { + switch (cgEvent.cfgIndexDefinitionCentrality.value) { case toI(CentralityDefinition::Ft0a): LOG(info) << "Enabling centrality definition: FT0A"; break; @@ -866,12 +1210,12 @@ struct PartNumFluc { break; } - hrCounter.add("hNEvents", ";;No. of Events", {HistType::kTH1D, {{10 + aod::evsel::EventSelectionFlags::kNsel, -0.5, 9.5 + static_cast(aod::evsel::EventSelectionFlags::kNsel), "Selection"}}}); + hrCounter.add("hNEvents", ";;No. of Events", {HistType::kTH1D, {{NEs + aod::evsel::EventSelectionFlags::kNsel, -0.5, static_cast(NEs + aod::evsel::EventSelectionFlags::kNsel) - 0.5, "Selection"}}}); if (doProcessMc.value) { - hrCounter.add("hNMcEvents", ";;No. of MC Events", {HistType::kTH1D, {{10, -0.5, 9.5, "Selection"}}}); + hrCounter.add("hNMcEvents", ";;No. of MC Events", {HistType::kTH1D, {{NEs, -0.5, static_cast(NEs) - 0.5, "Selection"}}}); } - if (groupAnalysis.cfgFlagQaRun.value) { + if (cgAnalysis.cfgFlagQaRun.value) { LOG(info) << "Enabling run QA."; const HistogramConfigSpec hcsQaRun{HistType::kTProfile, {{static_cast(holderCcdb.runNumbersIndicesGroupIndices.size()), -0.5, holderCcdb.runNumbersIndicesGroupIndices.size() - 0.5, "Run Index"}}}; @@ -912,72 +1256,71 @@ struct PartNumFluc { {std::format("{}NSigma{}", getName(Detector::Tof), getName(ParticleSpecies::Kaon)), std::format("{} #LT#it{{n}}#it{{#sigma}}_{{K}}#GT", getName(Detector::Tof)), true}, {std::format("{}NSigma{}", getName(Detector::Tof), getName(ParticleSpecies::Proton)), std::format("{} #LT#it{{n}}#it{{#sigma}}_{{p}}#GT", getName(Detector::Tof)), true}})) { if (!isChargeSpeciesSeparated) { - hrQaRun.add(std::format("QaRun/pRunIndex{}", name).c_str(), std::format(";;{}", title).c_str(), hcsQaRun); + hrQaRun.add(std::format("pRunIndex{}", name).c_str(), std::format(";;{}", title).c_str(), hcsQaRun); } else { - hrQaRun.add(std::format("QaRun/pRunIndex{}_{}", name, getName(ChargeSpecies::Plus)).c_str(), std::format(";;{} (#it{{q}}>0)", title).c_str(), hcsQaRun); - hrQaRun.add(std::format("QaRun/pRunIndex{}_{}", name, getName(ChargeSpecies::Minus)).c_str(), std::format(";;{} (#it{{q}}<0)", title).c_str(), hcsQaRun); + hrQaRun.add(std::format("pRunIndex{}_{}", name, getName(ChargeSpecies::Plus)).c_str(), std::format(";;{} (#it{{q}}>0)", title).c_str(), hcsQaRun); + hrQaRun.add(std::format("pRunIndex{}_{}", name, getName(ChargeSpecies::Minus)).c_str(), std::format(";;{} (#it{{q}}<0)", title).c_str(), hcsQaRun); } } } - if (groupAnalysis.cfgFlagQaEvent.value) { + if (cgAnalysis.cfgFlagQaEvent.value) { LOG(info) << "Enabling event QA."; - const AxisSpec asNTracks{200, -0.5, 199.5}; - const HistogramConfigSpec hcsQaEvent{HistType::kTHnSparseD, {asNTracks, asNTracks}}; + const AxisSpec asMultiplicity{cgEvent.cfgNMultiplicityBinsAxis.value, -0.5, cgEvent.cfgNMultiplicityBinsAxis.value - 0.5}; + const HistogramConfigSpec hcsQaEvent{HistType::kTHnSparseD, {asMultiplicity, asMultiplicity}}; - hrQaEvent.add("QaEvent/hVxVy", "", {HistType::kTHnSparseD, {{150, -0.15, 0.15, "#it{V}_{#it{x}} (cm)"}, {150, -0.15, 0.15, "#it{V}_{#it{y}} (cm)"}}}); - hrQaEvent.add("QaEvent/hVz", "", {HistType::kTH1D, {{300, -15., 15., "#it{V}_{#it{z}} (cm)"}}}); - hrQaEvent.add("QaEvent/hNPvContributorsNGlobalTracks", ";nPvContributors;nGlobalTracks;", hcsQaEvent); - hrQaEvent.add("QaEvent/hNGlobalTracksMeanDcaXy", ";nGlobalTracks;", {HistType::kTHnSparseD, {asNTracks, {250, -0.25, 0.25, "#LTDCA_{#it{xy}}#GT_{event} (cm)"}}}); - hrQaEvent.add("QaEvent/hNGlobalTracksMeanDcaXy_nPvContributorsCut", ";nGlobalTracks;", {HistType::kTHnSparseD, {asNTracks, {250, -0.25, 0.25, "#LTDCA_{#it{xy}}#GT_{event} (cm)"}}}); - hrQaEvent.add("QaEvent/hNGlobalTracksMeanDcaZ", ";nGlobalTracks;", {HistType::kTHnSparseD, {asNTracks, {200, -2., 2., "#LTDCA_{#it{z}}#GT_{event} (cm)"}}}); - hrQaEvent.add("QaEvent/hNGlobalTracksMeanDcaZ_nPvContributorsCut", ";nGlobalTracks;", {HistType::kTHnSparseD, {asNTracks, {200, -2., 2., "#LTDCA_{#it{z}}#GT_{event} (cm)"}}}); - hrQaEvent.add("QaEvent/hNTofBetaNGlobalTracks", ";nTofBeta;nGlobalTracks;", hcsQaEvent); - hrQaEvent.add("QaEvent/hNTofBetaNGlobalTracks_nPvContributorsCut", ";nTofBeta;nGlobalTracks;", hcsQaEvent); + hrQaEvent.add("hVxVy", "", {HistType::kTHnSparseD, {{150, -0.15, 0.15, "#it{V}_{#it{x}} (cm)"}, {150, -0.15, 0.15, "#it{V}_{#it{y}} (cm)"}}}); + hrQaEvent.add("hVz", "", {HistType::kTH1D, {{300, -15., 15., "#it{V}_{#it{z}} (cm)"}}}); + hrQaEvent.add("hNPvContributorsNGlobalTracks", ";nPvContributors;nGlobalTracks;", hcsQaEvent); + hrQaEvent.add("hNGlobalTracksMeanDcaXy", ";nGlobalTracks;", {HistType::kTHnSparseD, {asMultiplicity, {250, -0.25, 0.25, "#LTDCA_{#it{xy}}#GT_{event} (cm)"}}}); + hrQaEvent.add("hNGlobalTracksMeanDcaXy_nPvContributorsCut", ";nGlobalTracks;", {HistType::kTHnSparseD, {asMultiplicity, {250, -0.25, 0.25, "#LTDCA_{#it{xy}}#GT_{event} (cm)"}}}); + hrQaEvent.add("hNGlobalTracksMeanDcaZ", ";nGlobalTracks;", {HistType::kTHnSparseD, {asMultiplicity, {200, -2., 2., "#LTDCA_{#it{z}}#GT_{event} (cm)"}}}); + hrQaEvent.add("hNGlobalTracksMeanDcaZ_nPvContributorsCut", ";nGlobalTracks;", {HistType::kTHnSparseD, {asMultiplicity, {200, -2., 2., "#LTDCA_{#it{z}}#GT_{event} (cm)"}}}); + hrQaEvent.add("hNTofBetaNGlobalTracks", ";nTofBeta;nGlobalTracks;", hcsQaEvent); + hrQaEvent.add("hNTofBetaNGlobalTracks_nPvContributorsCut", ";nTofBeta;nGlobalTracks;", hcsQaEvent); } - if (groupAnalysis.cfgFlagQaCentrality.value) { + if (cgAnalysis.cfgFlagQaCentrality.value) { LOG(info) << "Enabling centrality QA."; - hrQaCentrality.add("QaCentrality/hCentralitySelection", "", {HistType::kTHnSparseD, {{100, 0., 100., "Centrality (%)"}, {10 + aod::evsel::EventSelectionFlags::kNsel, -0.5, 9.5 + static_cast(aod::evsel::EventSelectionFlags::kNsel), "Selection"}}}); - hrQaCentrality.add("QaCentrality/hCentralityMultiplicity", "", {HistType::kTHnSparseD, {{100, 0., 100., "Centrality (%)"}, {200, -0.5, 199.5, "Multiplicity"}}}); + hrQaCentrality.add("hCentralitySelection", "", {HistType::kTHnSparseD, {{100, 0., 100., "Centrality (%)"}, {NEs + aod::evsel::EventSelectionFlags::kNsel, -0.5, static_cast(NEs + aod::evsel::EventSelectionFlags::kNsel) - 0.5, "Selection"}}}); + hrQaCentrality.add("hCentralityMultiplicity", "", {HistType::kTHnSparseD, {{100, 0., 100., "Centrality (%)"}, {cgEvent.cfgNMultiplicityBinsAxis.value, -0.5, cgEvent.cfgNMultiplicityBinsAxis.value - 0.5, "Multiplicity"}}}); } - if (groupAnalysis.cfgFlagQaTrack.value) { + if (cgAnalysis.cfgFlagQaTrack.value) { LOG(info) << "Enabling track QA."; - for (const auto& [name, configSpec] : std::to_array>( - {{"ItsNCls", {HistType::kTH1D, {{10, -0.5, 9.5, "ITS nClusters"}}}}, - {"ItsChi2NCls", {HistType::kTH1D, {{80, 0., 40., "ITS #it{#chi}^{2}/nClusters"}}}}, + for (const auto& [name, hcs] : std::to_array>( + {{"ItsNClsChi2NCls", {HistType::kTHnSparseD, {{10, -0.5, 9.5, "ITS nClusters"}, {80, 0., 40., "ITS #it{#chi}^{2}/nClusters"}}}}, + {"TpcNClsChi2NCls", {HistType::kTHnSparseD, {{180, -0.5, 179.5, "TPC nClusters"}, {100, 0., 5., "TPC #it{#chi}^{2}/nClusters"}}}}, {"TpcNClsNClsShared", {HistType::kTHnSparseD, {{180, -0.5, 179.5, "TPC nClusters"}, {180, -0.5, 179.5, "TPC nSharedClusters"}}}}, - {"TpcChi2NCls", {HistType::kTH1D, {{100, 0., 5., "TPC #it{#chi}^{2}/nClusters"}}}}, {"TpcNClsFindableNCrossedRows", {HistType::kTHnSparseD, {{180, -0.5, 179.5, "TPC nFindableClusters"}, {180, -0.5, 179.5, "TPC nCrossedRows"}}}}})) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaTrack.add(std::format("QaTrack/h{}_{}", name, getName(iChargeSpecies)).c_str(), "", configSpec); + hrQaTrack.add(std::format("h{}_{}", name, getName(iChargeSpecies)).c_str(), "", hcs); } } } - if (groupAnalysis.cfgFlagQaDca.value) { + if (cgAnalysis.cfgFlagQaDca.value) { LOG(info) << "Enabling DCA QA."; const AxisSpec asPt{40, 0., 2., "#it{p}_{T} (GeV/#it{c})"}; - const HistogramConfigSpec hcsQaDcaProfile{HistType::kTProfile3D, {{groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, asPt, {24, -1.2, 1.2, "#it{#eta}"}}}; + const HistogramConfigSpec hcsQaDcaProfile{HistType::kTProfile3D, {{cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, asPt, {24, -1.2, 1.2, "#it{#eta}"}}}; - for (const auto& [name, title, configSpec] : std::to_array>( + for (const auto& [name, title, hcs] : std::to_array>( {{"hPtDcaXy", "", {HistType::kTHnSparseD, {asPt, {250, -0.25, 0.25, "DCA_{#it{xy}} (cm)"}}}}, {"pCentralityPtEtaDcaXy", ";;#LTDCA_{#it{xy}}#GT (cm)", hcsQaDcaProfile}, {"hPtDcaZ", "", {HistType::kTHnSparseD, {asPt, {250, -0.5, 0.5, "DCA_{#it{z}} (cm)"}}}}, {"pCentralityPtEtaDcaZ", ";;#LTDCA_{#it{z}}#GT (cm)", hcsQaDcaProfile}})) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaDca.add(std::format("QaDca/{}_{}", name, getName(iChargeSpecies)).c_str(), title.data(), configSpec); + hrQaDca.add(std::format("{}_{}", name, getName(iChargeSpecies)).c_str(), title.data(), hcs); } } } for (std::int32_t const& iParticleSpeciesAll : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsQaAcceptance.value.get(iParticleSpeciesAll))) { + if (!static_cast(cgAnalysis.cfgFlagsQaAcceptance.value.get(iParticleSpeciesAll))) { continue; } @@ -985,107 +1328,107 @@ struct PartNumFluc { for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaAcceptance.add(std::format("QaAcceptance/h{}Pt_{}Edge{}{}", iParticleSpeciesAll == toI(ParticleSpeciesAll::All) ? "Eta" : "Rapidity", getName(iPidStrategy), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), "", {HistType::kTHnSparseD, {{300, -1.5, 1.5, iParticleSpeciesAll == toI(ParticleSpeciesAll::All) ? "#it{#eta}" : "#it{y}"}, {250, 0., 2.5, "#it{p}_{T} (GeV/#it{c})"}}}); + hrQaAcceptance.add(std::format("h{}Pt_{}Edge{}{}", iParticleSpeciesAll == toI(ParticleSpeciesAll::All) ? "Eta" : "Rapidity", getName(iPidStrategy), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), "", {HistType::kTHnSparseD, {{300, -1.5, 1.5, iParticleSpeciesAll == toI(ParticleSpeciesAll::All) ? "#it{#eta}" : "#it{y}"}, {250, 0., 2.5, "#it{p}_{T} (GeV/#it{c})"}}}); } } } for (std::int32_t const& iParticleSpeciesAll : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsQaPhi.value.get(iParticleSpeciesAll))) { + if (!static_cast(cgAnalysis.cfgFlagsQaPhi.value.get(iParticleSpeciesAll))) { continue; } LOG(info) << "Enabling " << getName(iParticleSpeciesAll) << " phi QA."; - const HistogramConfigSpec hcsQaPhi{HistType::kTHnSparseF, {{groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {24, -1.2, 1.2, "#it{#eta}"}, {360, 0., constants::math::TwoPI, "#it{#varphi} (rad)"}}}; + const HistogramConfigSpec hcsQaPhi{HistType::kTHnSparseF, {{cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {24, -1.2, 1.2, "#it{#eta}"}, {360, 0., constants::math::TwoPI, "#it{#varphi} (rad)"}}}; for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPhi.add(std::format("QaPhi/hCentralityPtEtaPhi_{}{}{}{}", doProcessMc.value ? "mc" : "", doProcessMc.value ? getName(iPidStrategy) : getName(iPidStrategy), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), "", hcsQaPhi); + hrQaPhi.add(std::format("hCentralityPtEtaPhi_{}{}{}{}", doProcessMc.value ? "mc" : "", doProcessMc.value ? getName(iPidStrategy) : getName(iPidStrategy), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), "", hcsQaPhi); } } } for (std::int32_t const& iParticleSpeciesAll : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsQaPid.value.get(iParticleSpeciesAll))) { + if (!static_cast(cgAnalysis.cfgFlagsQaPid.value.get(iParticleSpeciesAll))) { continue; } LOG(info) << "Enabling " << getName(iParticleSpeciesAll) << " PID QA."; - const AxisSpec asCentrality{groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}; + const AxisSpec asCentrality{cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}; if (iParticleSpeciesAll == toI(ParticleSpeciesAll::All)) { const AxisSpec asPOverQ{350, -3.5, 3.5, "#it{p}/#it{q} (GeV/#it{c})"}; const AxisSpec asEta{48, -1.2, 1.2, "#it{#eta}"}; - hrQaPid.add("QaPid/hCentralityPOverQEtaTpcLnDeDx", "", {HistType::kTHnSparseF, {asCentrality, asPOverQ, asEta, {240, 3., 9., "TPC ln(d#it{E}/d#it{x} (a.u.))"}}}); - hrQaPid.add("QaPid/hCentralityPOverQEtaTofInverseBeta", "", {HistType::kTHnSparseF, {asCentrality, asPOverQ, asEta, {120, 0.5, 3.5, "TOF 1/#it{#beta}"}}}); + hrQaPid.add("hCentralityPOverQEtaTpcLnDeDx", "", {HistType::kTHnSparseF, {asCentrality, asPOverQ, asEta, {240, 3., 9., "TPC ln(d#it{E}/d#it{x} (a.u.))"}}}); + hrQaPid.add("hCentralityPOverQEtaTofInverseBeta", "", {HistType::kTHnSparseF, {asCentrality, asPOverQ, asEta, {120, 0.5, 3.5, "TOF 1/#it{#beta}"}}}); } else { const HistogramConfigSpec hcsQaPid{HistType::kTHnSparseF, {asCentrality, {40, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {32, -0.8, 0.8, "#it{#eta}"}, {300, -30., 30.}}}; if (doProcessMc.value) { for (std::int32_t const& iDetector : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPid.add(std::format("QaPid/hCentralityPtEta{}NSigma{}_mc{}{}", getName(iDetector), getName(iParticleSpeciesAll), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(iDetector), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); + hrQaPid.add(std::format("hCentralityPtEta{}NSigma{}_mc{}{}", getName(iDetector), getName(iParticleSpeciesAll), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(iDetector), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); } } } else { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPid.add(std::format("QaPid/hCentralityPtEta{}NSigma{}_{}", getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tpc), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); + hrQaPid.add(std::format("hCentralityPtEta{}NSigma{}_{}", getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tpc), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); } for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPid.add(std::format("QaPid/hCentralityPtEta{}NSigma{}_{}{}{}", getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(Detector::Tof), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tpc), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); + hrQaPid.add(std::format("hCentralityPtEta{}NSigma{}_{}{}{}", getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(Detector::Tof), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tpc), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); } for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPid.add(std::format("QaPid/hCentralityPtEta{}NSigma{}_{}{}{}", getName(Detector::Tof), getName(iParticleSpeciesAll), getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tof), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); + hrQaPid.add(std::format("hCentralityPtEta{}NSigma{}_{}{}{}", getName(Detector::Tof), getName(iParticleSpeciesAll), getName(Detector::Tpc), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(Detector::Tof), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); } for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrQaPid.add(std::format("QaPid/hCentralityPtEta{}NSigma{}_{}", getName(PidStrategy::TpcTof), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(PidStrategy::TpcTof), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); + hrQaPid.add(std::format("hCentralityPtEta{}NSigma{}_{}", getName(PidStrategy::TpcTof), getName(iParticleSpeciesAll), getName(iChargeSpecies)).c_str(), std::format(";;;;{} #it{{n}}#it{{#sigma}}_{{{}}};", getName(PidStrategy::TpcTof), getTitle(iParticleSpeciesAll)).c_str(), hcsQaPid); } } } } if (doProcessMc.value) { - if (groupAnalysis.cfgFlagQaMc.value) { + if (cgAnalysis.cfgFlagQaMc.value) { LOG(info) << "Enabling MC QA."; - const double maxAbsVz{std::ceil(groupEvent.cfgFlagMcCollisionVz.value ? groupEvent.cfgCutMaxAbsVzMc.value : groupEvent.cfgCutMaxAbsVz.value)}; + const double maxAbsVz{std::ceil(cgEvent.cfgFlagMcCollisionVz.value ? cgEvent.cfgCutMaxAbsVzMc.value : cgEvent.cfgCutMaxAbsVz.value)}; const AxisSpec asCentrality{20, 0., 100., "Centrality (%)"}; - hrQaMc.add("QaMc/hCentralityVzMcDeltaVz", "", {HistType::kTHnSparseF, {asCentrality, {static_cast(maxAbsVz) * 20, -maxAbsVz, maxAbsVz, "#it{V}_{#it{z}}^{Gen} (cm)"}, {200, -0.2, 0.2, "#it{V}_{#it{z}}^{Rec}#minus#it{V}_{#it{z}}^{Gen} (cm)"}}}); - hrQaMc.add("QaMc/hCentralityPtMcEtaMcDeltaPt", "", {HistType::kTHnSparseF, {asCentrality, {200, 0., 2., "#it{p}_{T}^{Gen} (GeV/#it{c})"}, {24, -1.2, 1.2, "#it{#eta}_{Gen}"}, {320, -0.8, 0.8, "#it{p}_{T}^{Rec}#minus#it{p}_{T}^{Gen} (GeV/#it{c})"}}}); - hrQaMc.add("QaMc/hCentralityPtMcEtaMcDeltaEta", "", {HistType::kTHnSparseF, {asCentrality, {20, 0., 2., "#it{p}_{T}^{Gen} (GeV/#it{c})"}, {240, -1.2, 1.2, "#it{#eta}_{Gen}"}, {160, -0.4, 0.4, "#it{#eta}_{Rec}#minus#it{#eta}_{Gen}"}}}); + hrQaMc.add("hCentralityVzMcDeltaVz", "", {HistType::kTHnSparseF, {asCentrality, {static_cast(maxAbsVz) * 20, -maxAbsVz, maxAbsVz, "#it{V}_{#it{z}}^{Gen} (cm)"}, {200, -0.2, 0.2, "#it{V}_{#it{z}}^{Rec}#minus#it{V}_{#it{z}}^{Gen} (cm)"}}}); + hrQaMc.add("hCentralityPtMcEtaMcDeltaPt", "", {HistType::kTHnSparseF, {asCentrality, {200, 0., 2., "#it{p}_{T}^{Gen} (GeV/#it{c})"}, {24, -1.2, 1.2, "#it{#eta}_{Gen}"}, {320, -0.8, 0.8, "#it{p}_{T}^{Rec}#minus#it{p}_{T}^{Gen} (GeV/#it{c})"}}}); + hrQaMc.add("hCentralityPtMcEtaMcDeltaEta", "", {HistType::kTHnSparseF, {asCentrality, {20, 0., 2., "#it{p}_{T}^{Gen} (GeV/#it{c})"}, {240, -1.2, 1.2, "#it{#eta}_{Gen}"}, {160, -0.4, 0.4, "#it{#eta}_{Rec}#minus#it{#eta}_{Gen}"}}}); } } for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsCalculationYield.value.get(iParticleSpecies))) { + if (!static_cast(cgAnalysis.cfgFlagsCalculationYield.value.get(iParticleSpecies))) { continue; } LOG(info) << "Enabling " << getName(iParticleSpecies) << " yield calculation."; - const double maxAbsVz{std::ceil(doProcessMc.value && groupEvent.cfgFlagMcCollisionVz.value ? groupEvent.cfgCutMaxAbsVzMc.value : groupEvent.cfgCutMaxAbsVz.value)}; - const HistogramConfigSpec hcsCalculationYield{HistType::kTHnSparseF, {{static_cast(maxAbsVz) * 2, -maxAbsVz, maxAbsVz, "#it{V}_{#it{z}} (cm)"}, {groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {40, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {32, -0.8, 0.8, "#it{#eta}"}}}; + const double maxAbsVz{std::ceil(doProcessMc.value && cgEvent.cfgFlagMcCollisionVz.value ? cgEvent.cfgCutMaxAbsVzMc.value : cgEvent.cfgCutMaxAbsVz.value)}; + const HistogramConfigSpec hcsCalculationYield{HistType::kTHnSparseF, {{static_cast(maxAbsVz) * 2, -maxAbsVz, maxAbsVz, "#it{V}_{#it{z}} (cm)"}, {cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {40, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {32, -0.8, 0.8, "#it{#eta}"}}}; if (doProcessMc.value) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrCalculationYield.add(std::format("CalculationYield/hVzCentralityPtMcEtaMc_mc{}{}", getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); + hrCalculationYield.add(std::format("hVzCentralityPtMcEtaMc_mc{}{}", getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); } for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - if (groupTrack.cfgFlagMcParticleMomentum.value) { - hrCalculationYield.add(std::format("CalculationYield/hVzCentralityPtMcEtaMc_mc{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); + if (cgTrack.cfgFlagMcParticleMomentum.value) { + hrCalculationYield.add(std::format("hVzCentralityPtMcEtaMc_mc{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); } else { - hrCalculationYield.add(std::format("CalculationYield/hVzCentralityPtEta_mc{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); + hrCalculationYield.add(std::format("hVzCentralityPtEta_mc{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); } } } } else { for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrCalculationYield.add(std::format("CalculationYield/hVzCentralityPtEta_{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); + hrCalculationYield.add(std::format("hVzCentralityPtEta_{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationYield); } } } @@ -1093,17 +1436,17 @@ struct PartNumFluc { if (doProcessMc.value) { for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsCalculationPurity.value.get(iParticleSpecies))) { + if (!static_cast(cgAnalysis.cfgFlagsCalculationPurity.value.get(iParticleSpecies))) { continue; } LOG(info) << "Enabling " << getName(iParticleSpecies) << " purity calculation."; - const HistogramConfigSpec hcsCalculationPurity{HistType::kTProfile3D, {{groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}}}; + const HistogramConfigSpec hcsCalculationPurity{HistType::kTProfile3D, {{cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}}}; for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrCalculationPurity.add(std::format("CalculationPurity/pCentralityPtEtaPurity{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationPurity); + hrCalculationPurity.add(std::format("pCentralityPtEtaPurity{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationPurity); } } } @@ -1111,53 +1454,51 @@ struct PartNumFluc { if (doProcessMc.value) { for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsCalculationFractionPrimary.value.get(iParticleSpecies))) { + if (!static_cast(cgAnalysis.cfgFlagsCalculationFractionPrimary.value.get(iParticleSpecies))) { continue; } LOG(info) << "Enabling " << getName(iParticleSpecies) << " primary fraction calculation."; - const HistogramConfigSpec hcsCalculationFractionPrimary{HistType::kTProfile3D, {{groupEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}}}; + const HistogramConfigSpec hcsCalculationFractionPrimary{HistType::kTProfile3D, {{cgEvent.cfgAxisCentralityCalibration, "Centrality (%)"}, {20, 0., 2., "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}}}; for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - hrCalculationFractionPrimary.add(std::format("CalculationFractionPrimary/pCentralityPtEtaFractionPrimary{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationFractionPrimary); + hrCalculationFractionPrimary.add(std::format("pCentralityPtEtaFractionPrimary{}{}{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies)).c_str(), "", hcsCalculationFractionPrimary); } } } } - if (nEnabled(groupAnalysis.cfgFlagsCalculationFluctuation) > 1) { - LOG(fatal) << "Invalid " << groupAnalysis.cfgFlagsCalculationFluctuation.name << "!"; - } - if (doCalculationFluctuation && groupEvent.cfgNSubgroups.value <= 0) { - LOG(fatal) << "Invalid " << groupEvent.cfgNSubgroups.name << "!"; - } for (std::int32_t const& iParticleNumber : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsCalculationFluctuation.value.get(iParticleNumber))) { + if (!static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(iParticleNumber))) { continue; } LOG(info) << "Enabling " << getName(iParticleNumber) << " number fluctuation calculation."; - const AxisSpec asCentrality{groupEvent.cfgAxisCentrality, "Centrality (%)"}; - const HistogramConfigSpec hcsDistribution{HistType::kTHnSparseD, {asCentrality, {200, -0.5, 199.5}, {200, -0.5, 199.5}}}; - const HistogramConfigSpec hcsFluctuationCalculator{HistType::kTH3D, {asCentrality, {groupEvent.cfgNSubgroups.value, -0.5, groupEvent.cfgNSubgroups.value - 0.5, "Subgroup Index"}, {fluctuation_calculator_base::NOrderKeys, -0.5, fluctuation_calculator_base::NOrderKeys - 0.5, "Order Key Index"}}}; + const AxisSpec asCentrality{cgEvent.cfgAxisCentrality, "Centrality (%)"}; + const AxisSpec asMultiplicity{cgEvent.cfgNMultiplicityBinsAxis.value, -0.5, cgEvent.cfgNMultiplicityBinsAxis.value - 0.5}; + const HistogramConfigSpec hcsDistribution{HistType::kTHnSparseD, {asCentrality, asMultiplicity, asMultiplicity}}; + const HistogramConfigSpec hcsFluctuationCalculator{HistType::kTH3D, {asCentrality, {cgEvent.cfgNSubgroups.value, -0.5, cgEvent.cfgNSubgroups.value - 0.5, "Subgroup Index"}, {fluctuation_calculator_base::NOrderKeys, -0.5, fluctuation_calculator_base::NOrderKeys - 0.5, "Order Key Index"}}}; for (std::int32_t const& iChargeNumber : std::views::iota(0, NEs)) { - fluctuationCalculatorsTrack[iParticleNumber][iChargeNumber] = std::make_unique(); + if (doProcessMc.value) { + fluctuationCalculatorsTrackMcParticle[iParticleNumber][iChargeNumber] = std::make_unique(); + } + fluctuationCalculatorsTrackTrack[iParticleNumber][iChargeNumber] = std::make_unique(); } if (doProcessMc.value) { - hrCalculationFluctuation.add(std::format("CalculationFluctuation/hCentralityN{}{}N{}{}_mc", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); - hrCalculationFluctuation.add(std::format("CalculationFluctuation/hCentralityN{}{}N{}{}_mcEff", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); + hrCalculationFluctuation.add(std::format("hCentralityN{}{}N{}{}_mc", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); + hrCalculationFluctuation.add(std::format("hCentralityN{}{}N{}{}_mcEff", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); for (std::int32_t const& iChargeNumber : std::views::iota(0, NEs)) { - hrCalculationFluctuation.add(std::format("CalculationFluctuation/hFluctuationCalculator{}{}_mc", getName(iParticleNumber), getName(iChargeNumber)).c_str(), "", hcsFluctuationCalculator); + hrCalculationFluctuation.add(std::format("hFluctuationCalculator{}{}_mc", getName(iParticleNumber), getName(iChargeNumber)).c_str(), "", hcsFluctuationCalculator); } } - hrCalculationFluctuation.add(std::format("CalculationFluctuation/hCentralityN{}{}N{}{}", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); + hrCalculationFluctuation.add(std::format("hCentralityN{}{}N{}{}", getName(iParticleNumber), getName(ChargeSpecies::Plus), getName(iParticleNumber), getName(ChargeSpecies::Minus)).c_str(), std::format(";;#it{{N}}({}^{{{}}});#it{{N}}({}^{{{}}});", getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Plus)), getTitle(iParticleNumber), getTitle(toI(ChargeSpecies::Minus))).c_str(), hcsDistribution); for (std::int32_t const& iChargeNumber : std::views::iota(0, NEs)) { - hrCalculationFluctuation.add(std::format("CalculationFluctuation/hFluctuationCalculator{}{}", getName(iParticleNumber), getName(iChargeNumber)).c_str(), "", hcsFluctuationCalculator); + hrCalculationFluctuation.add(std::format("hFluctuationCalculator{}{}", getName(iParticleNumber), getName(iChargeNumber)).c_str(), "", hcsFluctuationCalculator); } } } @@ -1166,7 +1507,7 @@ struct PartNumFluc { void readCcdb() { if constexpr (DoInit) { - holderCcdb.lCcdb = ccdb->get(groupCcdb.cfgCcdbPath.value); + holderCcdb.lCcdb = ccdb->get(cgCcdb.cfgCcdbPath.value); if (!holderCcdb.lCcdb || holderCcdb.lCcdb->IsA() != TList::Class()) { LOG(fatal) << "Invalid CCDB object!"; } @@ -1176,17 +1517,17 @@ struct PartNumFluc { LOG(fatal) << "Invalid gRunNumberGroupIndex!"; } for (std::int32_t const& iRun : std::views::iota(0, gRunNumberGroupIndex->GetN())) { - holderCcdb.runNumbersIndicesGroupIndices[std::llrint(gRunNumberGroupIndex->GetX()[iRun])] = {iRun, static_cast(std::llrint(gRunNumberGroupIndex->GetY()[iRun]))}; + holderCcdb.runNumbersIndicesGroupIndices[static_cast(std::rint(gRunNumberGroupIndex->GetX()[iRun]))] = {iRun, static_cast(std::rint(gRunNumberGroupIndex->GetY()[iRun]))}; } - if (groupEvent.cfgFlagRejectionRunBadMc.value) { + if (cgEvent.cfgFlagRejectionRunBadMc.value) { const TGraph* const gRunNumberGroupIndexMc{dynamic_cast(holderCcdb.lCcdb->FindObject("gRunNumberGroupIndex_mc"))}; if (!gRunNumberGroupIndexMc || gRunNumberGroupIndexMc->IsA() != TGraph::Class()) { LOG(fatal) << "Invalid gRunNumberGroupIndex_mc!"; } for (std::int32_t const& iRun : std::views::iota(0, gRunNumberGroupIndexMc->GetN())) { - if (std::llrint(gRunNumberGroupIndexMc->GetY()[iRun]) <= 0) { - if (const auto iter{holderCcdb.runNumbersIndicesGroupIndices.find(std::llrint(gRunNumberGroupIndexMc->GetX()[iRun]))}; iter != holderCcdb.runNumbersIndicesGroupIndices.end() && iter->second.second > 0) { + if (static_cast(std::rint(gRunNumberGroupIndexMc->GetY()[iRun])) <= 0) { + if (const auto iter{holderCcdb.runNumbersIndicesGroupIndices.find(static_cast(std::rint(gRunNumberGroupIndexMc->GetX()[iRun])))}; iter != holderCcdb.runNumbersIndicesGroupIndices.end() && iter->second.second > 0) { iter->second.second = -iter->second.second; } } @@ -1199,10 +1540,10 @@ struct PartNumFluc { } holderCcdb.runGroupIndexCurrent = runGroupIndex; - holderCcdb.fPtMeasureDca = {}; - holderCcdb.hCentralityPtEtaShiftNSigmaPid = {}; - holderCcdb.hVzCentralityPtEtaEfficiency = {}; - if (!groupTrack.cfgFlagRecalibrationDca.value && !isEnabled(groupTrack.cfgFlagsRecalibrationNSigmaPid) && !doCalculationFluctuation) { + holderCcdb.calibrationsPtMeasureDca = {}; + holderCcdb.hsCentralityPtEtaShiftNSigmaPid = {}; + holderCcdb.hsVzCentralityPtEtaEfficiency = {}; + if (!cgTrack.cfgFlagRecalibrationDca.value && !isEnabled(cgTrack.cfgFlagsRecalibrationNSigmaPid) && !doCalculationFluctuation) { return; } @@ -1212,11 +1553,11 @@ struct PartNumFluc { LOG(fatal) << "Invalid " << nameList << "!"; } - if (groupTrack.cfgFlagRecalibrationDca.value) { + if (cgTrack.cfgFlagRecalibrationDca.value) { for (std::int32_t const& iDcaMeasure : std::views::iota(0, NEs)) { for (std::int32_t const& iDcaAxis : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - std::pair& calibration{holderCcdb.fPtMeasureDca[iDcaMeasure][iDcaAxis][iChargeSpecies]}; + std::pair& calibration{holderCcdb.calibrationsPtMeasureDca[iDcaMeasure][iDcaAxis][iChargeSpecies]}; const std::string nameFormula{std::format("fPt{}Dca{}{}{}_runGroup{}", getName(iDcaMeasure), getName(iDcaAxis), getName(iChargeSpecies), doProcessMc.value ? "_mc" : "", runGroupIndex)}; calibration.first = dynamic_cast(lRunGroup->FindObject(nameFormula.c_str())); if (!calibration.first || calibration.first->GetNdim() != 1 || calibration.first->GetNpar() <= 0) { @@ -1237,14 +1578,14 @@ struct PartNumFluc { } for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { - if (!static_cast(groupTrack.cfgFlagsRecalibrationNSigmaPid.value.get(iParticleSpecies))) { + if (!static_cast(cgTrack.cfgFlagsRecalibrationNSigmaPid.value.get(iParticleSpecies))) { continue; } for (std::int32_t const& iDetector : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { const std::string name{std::format("hCentralityPtEtaShift{}NSigma{}{}{}_runGroup{}", getName(iDetector), getName(iParticleSpecies), getName(iChargeSpecies), doProcessMc.value ? "_mc" : "", runGroupIndex)}; - const TH3*& histogram{holderCcdb.hCentralityPtEtaShiftNSigmaPid[iDetector][iParticleSpecies][iChargeSpecies]}; + const TH3*& histogram{holderCcdb.hsCentralityPtEtaShiftNSigmaPid[iDetector][iParticleSpecies][iChargeSpecies]}; histogram = dynamic_cast(lRunGroup->FindObject(name.c_str())); if (!histogram) { LOG(fatal) << "Invalid " << name << "!"; @@ -1255,14 +1596,14 @@ struct PartNumFluc { } for (std::int32_t const& iParticleSpecies : std::views::iota(0, NEs)) { - if (!static_cast(groupAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Charge))) && (iParticleSpecies != toI(ParticleSpecies::Kaon) || !static_cast(groupAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Kaon)))) && (iParticleSpecies != toI(ParticleSpecies::Proton) || !static_cast(groupAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Proton))))) { + if (!static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Charge))) && (iParticleSpecies != toI(ParticleSpecies::Kaon) || !static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Kaon)))) && (iParticleSpecies != toI(ParticleSpecies::Proton) || !static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumber::Proton))))) { continue; } for (std::int32_t const& iPidStrategy : std::views::iota(0, NEs)) { for (std::int32_t const& iChargeSpecies : std::views::iota(0, NEs)) { - const std::string name{std::format("hVzCentralityPtEtaEfficiency{}{}{}_runGroup{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies), runGroupIndex)}; - const THnBase*& histogram{holderCcdb.hVzCentralityPtEtaEfficiency[iPidStrategy][iParticleSpecies][iChargeSpecies]}; + const std::string name{std::format("hsVzCentralityPtEtaEfficiency{}{}{}_runGroup{}", getName(iPidStrategy), getName(iParticleSpecies), getName(iChargeSpecies), runGroupIndex)}; + const THnBase*& histogram{holderCcdb.hsVzCentralityPtEtaEfficiency[iPidStrategy][iParticleSpecies][iChargeSpecies]}; histogram = dynamic_cast(lRunGroup->FindObject(name.c_str())); if (!histogram || histogram->GetNdimensions() != HolderCcdb::NDimensionsEfficiency) { LOG(fatal) << "Invalid " << name << "!"; @@ -1275,40 +1616,25 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue double getEfficiency(const bool doUseMcParticleMomentum) const { - const THnBase* const hVzCentralityPtEtaEfficiency{holderCcdb.hVzCentralityPtEtaEfficiency[toI(PidStrategyValue)][toI(ParticleSpeciesValue)][toI(ChargeSpeciesValue)]}; - return hVzCentralityPtEtaEfficiency ? hVzCentralityPtEtaEfficiency->GetBinContent(hVzCentralityPtEtaEfficiency->GetBin(std::array{doProcessMc.value && groupEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, doUseMcParticleMomentum ? holderMcParticle.pt : holderTrack.pt, doUseMcParticleMomentum ? holderMcParticle.eta : holderTrack.eta}.data())) : 0.; + const THnBase* const hsVzCentralityPtEtaEfficiency{holderCcdb.hsVzCentralityPtEtaEfficiency[toI(PidStrategyValue)][toI(ParticleSpeciesValue)][toI(ChargeSpeciesValue)]}; + return hsVzCentralityPtEtaEfficiency ? hsVzCentralityPtEtaEfficiency->GetBinContent(hsVzCentralityPtEtaEfficiency->GetBin(std::array{doProcessMc.value && cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, doUseMcParticleMomentum ? holderMcParticle.pt : holderTrack.pt, doUseMcParticleMomentum ? holderMcParticle.eta : holderTrack.eta}.data())) : 0.; } template - requires IsValid + requires IsValidEnumValue double getShiftNSigmaPid() const { if constexpr (DoRecalibrate) { - if (groupTrack.cfgFlagsRecalibrationNSigmaPid.value.get(toI(ParticleSpeciesValue))) { - return interpolate(holderCcdb.hCentralityPtEtaShiftNSigmaPid[toI(DetectorValue)][toI(ParticleSpeciesValue)][holderTrack.sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)], holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta); + if (cgTrack.cfgFlagsRecalibrationNSigmaPid.value.get(toI(ParticleSpeciesValue))) { + return interpolate(holderCcdb.hsCentralityPtEtaShiftNSigmaPid[toI(DetectorValue)][toI(ParticleSpeciesValue)][holderTrack.sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)], holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta); } } return 0.; } - template - void setNSigmaPid(const T& track) - { - if (holderTrack.hasPid[toI(Detector::Tpc)]) { - holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Pion)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaPi(), getShiftNSigmaPid()); - holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Kaon)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaKa(), getShiftNSigmaPid()); - holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Proton)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaPr(), getShiftNSigmaPid()); - } - if (holderTrack.hasPid[toI(Detector::Tof)]) { - holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Pion)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaPi(), getShiftNSigmaPid()); - holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Kaon)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaKa(), getShiftNSigmaPid()); - holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Proton)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaPr(), getShiftNSigmaPid()); - } - } - double getMeasureDca(const std::pair& calibration) const { static thread_local std::vector parametersPtMeasureDcaScratch{}; @@ -1316,7 +1642,7 @@ struct PartNumFluc { const TFormula* const fPtMeasureDca{calibration.first}; const TH3* const hCentralityEtaParameterPtMeasureDca{calibration.second}; const std::int32_t nParameters{fPtMeasureDca->GetNpar()}; - if (static_cast(parametersPtMeasureDcaScratch.size()) < nParameters) { + if (std::cmp_less(parametersPtMeasureDcaScratch.size(), nParameters)) { parametersPtMeasureDcaScratch.resize(nParameters); } const std::int32_t centralityBinIndex{std::clamp(hCentralityEtaParameterPtMeasureDca->GetXaxis()->FindFixBin(holderEvent.centralityCalibration), 1, hCentralityEtaParameterPtMeasureDca->GetNbinsX())}; @@ -1327,8 +1653,21 @@ struct PartNumFluc { return fPtMeasureDca->EvalPar(&holderTrack.pt, parametersPtMeasureDcaScratch.data()); } - template - requires IsValid + template + requires IsValidEnumValue + double getAbsNSigmaDca() const + { + if (!cgTrack.cfgFlagRecalibrationDca.value) { + return std::abs(holderTrack.dca[toI(DcaAxisValue)]); + } + + const std::int32_t chargeSpeciesIndex{holderTrack.sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)}; + const double sigma{getMeasureDca(holderCcdb.calibrationsPtMeasureDca[toI(DcaMeasure::Sigma)][toI(DcaAxisValue)][chargeSpeciesIndex])}; + return sigma > 0. ? std::abs((holderTrack.dca[toI(DcaAxisValue)] - getMeasureDca(holderCcdb.calibrationsPtMeasureDca[toI(DcaMeasure::Mean)][toI(DcaAxisValue)][chargeSpeciesIndex])) / sigma) : std::numeric_limits::infinity(); + } + + template + requires IsValidEnumValue bool isPid(const bool doRejectOthers) const { if constexpr (ParticleSpeciesAllValue == ParticleSpeciesAll::All) { @@ -1348,10 +1687,10 @@ struct PartNumFluc { } else { constexpr std::int32_t ParticleSpeciesIndex{toI(getValue(ParticleSpeciesAllValue))}; if constexpr (PidStrategyAllValue == PidStrategyAll::TpcTofSeparated) { - if (!(std::abs(holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]) < groupTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex))) { + if (!(std::abs(holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]) < cgTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex, toI(SelectionValue)))) { return false; } - if (!(std::abs(holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]) < groupTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex))) { + if (!(std::abs(holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]) < cgTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex, toI(SelectionValue)))) { return false; } if (doRejectOthers && !(std::abs(holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]) < std::min(std::abs(holderTrack.nSigmaPid[toI(Detector::Tof)][(ParticleSpeciesIndex + 1) % NEs]), std::abs(holderTrack.nSigmaPid[toI(Detector::Tof)][(ParticleSpeciesIndex + 2) % NEs])))) { @@ -1359,7 +1698,7 @@ struct PartNumFluc { } } else if constexpr (PidStrategyAllValue == PidStrategyAll::TpcTofCombined) { const double absNSigmaPidCombined{std::abs(holderTrack.getNSigmaPidCombined(ParticleSpeciesIndex))}; - if (!(absNSigmaPidCombined < groupTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex))) { + if (!(absNSigmaPidCombined < cgTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex, toI(SelectionValue)))) { return false; } if (doRejectOthers && !(absNSigmaPidCombined < std::min(std::abs(holderTrack.getNSigmaPidCombined((ParticleSpeciesIndex + 1) % NEs)), std::abs(holderTrack.getNSigmaPidCombined((ParticleSpeciesIndex + 2) % NEs))))) { @@ -1367,7 +1706,7 @@ struct PartNumFluc { } } else { constexpr std::int32_t DetectorIndex{toI(getValue(PidStrategyAllValue))}; - if (!(std::abs(holderTrack.nSigmaPid[DetectorIndex][ParticleSpeciesIndex]) < groupTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex))) { + if (!(std::abs(holderTrack.nSigmaPid[DetectorIndex][ParticleSpeciesIndex]) < cgTrack.cfgCutsMaxAbsNSigmaPid.value.get(ParticleSpeciesIndex, toI(SelectionValue)))) { return false; } if (doRejectOthers && !(std::abs(holderTrack.nSigmaPid[DetectorIndex][ParticleSpeciesIndex]) < std::min(std::abs(holderTrack.nSigmaPid[DetectorIndex][(ParticleSpeciesIndex + 1) % NEs]), std::abs(holderTrack.nSigmaPid[DetectorIndex][(ParticleSpeciesIndex + 2) % NEs])))) { @@ -1379,7 +1718,7 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue bool isPid() const { if constexpr (ParticleSpeciesAllValue == ParticleSpeciesAll::All) { @@ -1389,92 +1728,69 @@ struct PartNumFluc { } } + template + requires IsValidEnumValue bool isGoodMomentum(const bool doUseMcParticleMomentum) const { const double pt{doUseMcParticleMomentum ? holderMcParticle.pt : holderTrack.pt}; - const double eta{doUseMcParticleMomentum ? holderMcParticle.eta : holderTrack.eta}; - if (!(groupTrack.cfgCutMinPt.value < pt) || !(pt < groupTrack.cfgCutMaxPt.value)) { - return false; - } - if (!(std::abs(eta) < groupTrack.cfgCutMaxAbsEta.value)) { - return false; + if constexpr (ParticleSpeciesAllValue == ParticleSpeciesAll::All) { + if (!(rangePtAll[toI(RangeEdge::Min)] < pt) || !(pt < rangePtAll[toI(RangeEdge::Max)])) { + return false; + } + } else { + const std::int32_t particleSpeciesIndex{toI(getValue(ParticleSpeciesAllValue))}; + if (!(cgTrack.cfgCutsRangePt.value.get(particleSpeciesIndex, toI(RangeEdge::Min)) < pt) || !(pt < cgTrack.cfgCutsRangePt.value.get(particleSpeciesIndex, toI(RangeEdge::Max)))) { + return false; + } } - return true; + return std::abs(doUseMcParticleMomentum ? holderMcParticle.eta : holderTrack.eta) < cgTrack.cfgCutMaxAbsEta.value; } bool isGoodDca() const { - if (!groupTrack.cfgFlagRecalibrationDca.value) { - for (std::int32_t const& iDcaAxis : std::views::iota(0, NEs)) { - if (!(std::abs(holderTrack.dca[iDcaAxis]) < groupTrack.cfgCutsMaxAbsNSigmaDca.value.get(iDcaAxis))) { - return false; - } - } - } else { - const std::int32_t chargeSpeciesIndex{holderTrack.sign > 0 ? toI(ChargeSpecies::Plus) : toI(ChargeSpecies::Minus)}; - const std::array, NEs>, NEs>, NEs>& fPtMeasureDca{holderCcdb.fPtMeasureDca}; - for (std::int32_t const& iDcaAxis : std::views::iota(0, NEs)) { - const double mean{getMeasureDca(fPtMeasureDca[toI(DcaMeasure::Mean)][iDcaAxis][chargeSpeciesIndex])}; - const double sigma{getMeasureDca(fPtMeasureDca[toI(DcaMeasure::Sigma)][iDcaAxis][chargeSpeciesIndex])}; - if (!(std::abs(holderTrack.dca[iDcaAxis] - mean) < groupTrack.cfgCutsMaxAbsNSigmaDca.value.get(iDcaAxis) * sigma)) { - return false; - } - } - } - return true; + return getAbsNSigmaDca() < cgTrack.cfgCutsMaxAbsNSigmaDca.value.get(toI(DcaAxis::Xy), toI(Selection::Default)) && getAbsNSigmaDca() < cgTrack.cfgCutsMaxAbsNSigmaDca.value.get(toI(DcaAxis::Z), toI(Selection::Default)); } template bool isGoodTrack(const T& track) const { - if (groupTrack.cfgFlagPvContributor.value && !track.isPVContributor()) { + if (cgTrack.cfgFlagPvContributor.value && !track.isPVContributor()) { return false; } - if (!(track.itsNCls() > groupTrack.cfgCutMinItsNCls.value)) { + if (!(track.itsNCls() > cgTrack.cfgCutsMinItsNCls.value.get(toI(Selection::Default)))) { return false; } - if (!(track.itsChi2NCl() < groupTrack.cfgCutMaxItsChi2NCls.value)) { + if (!(track.itsChi2NCl() < cgTrack.cfgCutsMaxItsChi2NCls.value.get(toI(Selection::Default)))) { return false; } - if (!(track.tpcNClsFound() > groupTrack.cfgCutMinTpcNCls.value)) { + if (!(track.tpcNClsFound() > cgTrack.cfgCutsMinTpcNCls.value.get(toI(Selection::Default)))) { return false; } - if (!(groupTrack.cfgCutMinTpcChi2NCls.value < track.tpcChi2NCl()) || !(track.tpcChi2NCl() < groupTrack.cfgCutMaxTpcChi2NCls.value)) { + if (!(cgTrack.cfgCutMinTpcChi2NCls.value < track.tpcChi2NCl()) || !(track.tpcChi2NCl() < cgTrack.cfgCutsMaxTpcChi2NCls.value.get(toI(Selection::Default)))) { return false; } - if (!(track.tpcFractionSharedCls() < groupTrack.cfgCutMaxTpcNClsSharedRatio.value)) { + if (!(track.tpcFractionSharedCls() < cgTrack.cfgCutsMaxTpcNClsSharedRatio.value.get(toI(Selection::Default)))) { return false; } - if (!(track.tpcNClsCrossedRows() > groupTrack.cfgCutMinTpcNCrossedRows.value)) { + if (!(track.tpcNClsCrossedRows() > cgTrack.cfgCutsMinTpcNCrossedRows.value.get(toI(Selection::Default)))) { return false; } - if (!(track.tpcCrossedRowsOverFindableCls() > groupTrack.cfgCutMinTpcNCrossedRowsRatio.value)) { + if (!(track.tpcCrossedRowsOverFindableCls() > cgTrack.cfgCutMinTpcNCrossedRowsRatio.value)) { return false; } return true; } - template - bool isGoodMcParticle(const MP& mcParticle) const - { - if constexpr (IsMc) { - if (!mcParticle.isPhysicalPrimary()) { - return false; - } - } - return true; - } - template - requires IsValid + requires IsValidEnumValue void fillQaRunByTrackByChargeSpecies(const T& track) { const auto fill{[this](const auto& name, const auto value) -> void { - hrQaRun.fill(C_CS("QaRun/pRunIndex") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.runIndex, value); + hrQaRun.fill(C_CS("pRunIndex") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.runIndex, value); }}; const auto fillNSigmaPidByDetectorParticleSpecies{ [this, &fill] - requires IsValid + requires IsValidEnumValue () -> void { const double nSigmaPid{holderTrack.nSigmaPid[toI(DetectorValue)][toI(ParticleSpeciesValue)]}; if (std::abs(nSigmaPid) < HolderTrack::TruncationAbsNSigmaPid) { @@ -1507,11 +1823,11 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue void fillQaRunByEventByChargeSpecies() { const auto fill{[this](const auto& name, const auto value) -> void { - hrQaRun.fill(C_CS("QaRun/pRunIndex") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.runIndex, value); + hrQaRun.fill(C_CS("pRunIndex") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.runIndex, value); }}; fill(C_CS("NGlobalTracks"), holderEvent.nGlobalTracks[toI(ChargeSpeciesValue)]); @@ -1526,30 +1842,29 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue void fillQaTrackByChargeSpecies(const T& track) { const auto fill{[this](const auto& name, const auto... positionAndWeight) -> void { - hrQaTrack.fill(C_CS("QaTrack/h") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), positionAndWeight...); + hrQaTrack.fill(C_CS("h") + name + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), positionAndWeight...); }}; - fill(C_CS("ItsNCls"), track.itsNCls()); - fill(C_CS("ItsChi2NCls"), track.itsChi2NCl()); + fill(C_CS("ItsNClsChi2NCls"), track.itsNCls(), track.itsChi2NCl()); + fill(C_CS("TpcNClsChi2NCls"), track.tpcNClsFound(), track.tpcChi2NCl()); fill(C_CS("TpcNClsNClsShared"), track.tpcNClsFound(), track.tpcNClsShared()); - fill(C_CS("TpcChi2NCls"), track.tpcChi2NCl()); fill(C_CS("TpcNClsFindableNCrossedRows"), track.tpcNClsFindable(), track.tpcNClsCrossedRows()); } template - requires IsValid + requires IsValidEnumValue void fillQaDcaByChargeSpecies() { const auto fillByDcaAxis{ [this] - requires IsValid + requires IsValidEnumValue () -> void { - hrQaDca.fill(C_CS("QaDca/hPtDca") + C_SV(getName(DcaAxisValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderTrack.pt, holderTrack.dca[toI(DcaAxisValue)]); - hrQaDca.fill(C_CS("QaDca/pCentralityPtEtaDca") + C_SV(getName(DcaAxisValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.dca[toI(DcaAxisValue)]); + hrQaDca.fill(C_CS("hPtDca") + C_SV(getName(DcaAxisValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderTrack.pt, holderTrack.dca[toI(DcaAxisValue)]); + hrQaDca.fill(C_CS("pCentralityPtEtaDca") + C_SV(getName(DcaAxisValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.dca[toI(DcaAxisValue)]); }}; fillByDcaAxis.template operator()(); @@ -1557,23 +1872,23 @@ struct PartNumFluc { } template - requires IsValid - void fillQaAcceptancebyParticleSpeciesAll(const T& track) + requires IsValidEnumValue + void fillQaAcceptanceByParticleSpeciesAll(const T& track) { - if (!groupAnalysis.cfgFlagsQaAcceptance.value.get(toI(ParticleSpeciesAllValue))) { + if (!cgAnalysis.cfgFlagsQaAcceptance.value.get(toI(ParticleSpeciesAllValue))) { return; } const auto fillByChargeSpecies{ [this] - requires IsValid + requires IsValidEnumValue (const auto& name, const double value) -> void { const auto fillByPidStrategy{ [this, &name, value] - requires IsValid + requires IsValidEnumValue () -> void { if (isPid(PidStrategyValue), ParticleSpeciesAllValue>(false)) { - hrQaAcceptance.fill(C_CS("QaAcceptance/h") + name + C_CS("Pt_") + C_SV(getName(PidStrategyValue)) + C_CS("Edge") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), value, holderTrack.pt); // NOLINT(clang-analyzer-core.NonNullParamChecker) + hrQaAcceptance.fill(C_CS("h") + name + C_CS("Pt_") + C_SV(getName(PidStrategyValue)) + C_CS("Edge") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), value, holderTrack.pt); // NOLINT(clang-analyzer-core.NonNullParamChecker) } }}; @@ -1597,28 +1912,28 @@ struct PartNumFluc { } template - requires IsValid && (DataModeValue != DataMode::McMcParticle) + requires IsValidEnumValue && (DataModeValue != DataMode::McMcParticle) void fillQaPhiByParticleSpeciesAll() { - if (!groupAnalysis.cfgFlagsQaPhi.value.get(toI(ParticleSpeciesAllValue))) { + if (!cgAnalysis.cfgFlagsQaPhi.value.get(toI(ParticleSpeciesAllValue))) { return; } const auto fillByChargeSpecies{ [this] - requires IsValid + requires IsValidEnumValue () -> void { const auto fillByPidStrategy{ [this] - requires IsValid + requires IsValidEnumValue () -> void { if constexpr (DataModeValue == DataMode::McTrack) { if (isPid() && isPid(PidStrategyValue), ParticleSpeciesAllValue>(false)) { - hrQaPhi.fill(C_CS("QaPhi/hCentralityPtEtaPhi_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.phi); + hrQaPhi.fill(C_CS("hCentralityPtEtaPhi_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.phi); } } else { // DataModeValue == DataMode::RawTrack if (isPid(PidStrategyValue), ParticleSpeciesAllValue>(false)) { - hrQaPhi.fill(C_CS("QaPhi/hCentralityPtEtaPhi_") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.phi); + hrQaPhi.fill(C_CS("hCentralityPtEtaPhi_") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.phi); } } }}; @@ -1635,40 +1950,40 @@ struct PartNumFluc { } template - requires IsValid && (DataModeValue != DataMode::McMcParticle) + requires IsValidEnumValue && (DataModeValue != DataMode::McMcParticle) void fillQaPidByParticleSpeciesAll(const T& track) { - if (!groupAnalysis.cfgFlagsQaPid.value.get(toI(ParticleSpeciesAllValue))) { + if (!cgAnalysis.cfgFlagsQaPid.value.get(toI(ParticleSpeciesAllValue))) { return; } if constexpr (ParticleSpeciesAllValue == ParticleSpeciesAll::All) { if (isPid(PidStrategy::Tpc), ParticleSpeciesAll::All>(false)) { - hrQaPid.fill(C_CS("QaPid/hCentralityPOverQEtaTpcLnDeDx"), holderEvent.centralityCalibration, track.p() / holderTrack.sign, holderTrack.eta, std::log(track.tpcSignal())); + hrQaPid.fill(C_CS("hCentralityPOverQEtaTpcLnDeDx"), holderEvent.centralityCalibration, track.p() / holderTrack.sign, holderTrack.eta, std::log(track.tpcSignal())); } if (isPid(PidStrategy::TpcTof), ParticleSpeciesAll::All>(false)) { - hrQaPid.fill(C_CS("QaPid/hCentralityPOverQEtaTofInverseBeta"), holderEvent.centralityCalibration, track.p() / holderTrack.sign, holderTrack.eta, 1. / track.beta()); + hrQaPid.fill(C_CS("hCentralityPOverQEtaTofInverseBeta"), holderEvent.centralityCalibration, track.p() / holderTrack.sign, holderTrack.eta, 1. / track.beta()); } } else { constexpr std::int32_t ParticleSpeciesIndex{toI(getValue(ParticleSpeciesAllValue))}; const auto fillByChargeSpecies{ [this] - requires IsValid + requires IsValidEnumValue () -> void { if constexpr (DataModeValue == DataMode::McTrack) { if (isPid()) { - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_mc") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(Detector::Tof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_mc") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_mc") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(Detector::Tof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_mc") + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]); } } else { // DataModeValue == DataMode::RawTrack - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); if (isPid(false)) { - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(Detector::Tof)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(Detector::Tpc)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(Detector::Tof)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]); } if (isPid(false)) { - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(Detector::Tof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(Detector::Tpc)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(Detector::Tof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(Detector::Tpc)) + C_SV(getName(ParticleSpeciesAllValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.nSigmaPid[toI(Detector::Tof)][ParticleSpeciesIndex]); } - hrQaPid.fill(C_CS("QaPid/hCentralityPtEta") + C_SV(getName(PidStrategy::TpcTof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.getNSigmaPidCombined(ParticleSpeciesIndex)); + hrQaPid.fill(C_CS("hCentralityPtEta") + C_SV(getName(PidStrategy::TpcTof)) + C_CS("NSigma") + C_SV(getName(ParticleSpeciesAllValue)) + C_CS("_") + C_SV(getName(ChargeSpeciesValue)), holderEvent.centralityCalibration, holderTrack.pt, holderTrack.eta, holderTrack.getNSigmaPidCombined(ParticleSpeciesIndex)); } }}; @@ -1681,10 +1996,10 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue void fillCalculationYieldByParticleSpecies() { - if (!groupAnalysis.cfgFlagsCalculationYield.value.get(toI(ParticleSpeciesValue))) { + if (!cgAnalysis.cfgFlagsCalculationYield.value.get(toI(ParticleSpeciesValue))) { return; } @@ -1703,28 +2018,28 @@ struct PartNumFluc { const auto fillByChargeSpecies{ [this] - requires IsValid + requires IsValidEnumValue () -> void { if constexpr (DataModeValue == DataMode::McMcParticle) { if (isPid(ParticleSpeciesValue), ChargeSpeciesValue>()) { - hrCalculationYield.fill(C_CS("CalculationYield/hVzCentralityPtMcEtaMc_mc") + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), groupEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderMcParticle.pt, holderMcParticle.eta); + hrCalculationYield.fill(C_CS("hVzCentralityPtMcEtaMc_mc") + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderMcParticle.pt, holderMcParticle.eta); } } else { const auto fillByPidStrategy{ [this] - requires IsValid + requires IsValidEnumValue () -> void { if constexpr (DataModeValue == DataMode::McTrack) { - if (isPid(ParticleSpeciesValue), ChargeSpeciesValue>() && isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value)) { - if (groupTrack.cfgFlagMcParticleMomentum.value) { - hrCalculationYield.fill(C_CS("CalculationYield/hVzCentralityPtMcEtaMc_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), groupEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderMcParticle.pt, holderMcParticle.eta); + if (isPid(ParticleSpeciesValue), ChargeSpeciesValue>() && isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value)) { + if (cgTrack.cfgFlagMcParticleMomentum.value) { + hrCalculationYield.fill(C_CS("hVzCentralityPtMcEtaMc_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderMcParticle.pt, holderMcParticle.eta); } else { - hrCalculationYield.fill(C_CS("CalculationYield/hVzCentralityPtEta_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), groupEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderTrack.pt, holderTrack.eta); + hrCalculationYield.fill(C_CS("hVzCentralityPtEta_mc") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, holderTrack.pt, holderTrack.eta); } } } else { // DataModeValue == DataMode::RawTrack - if (isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value)) { - hrCalculationYield.fill(C_CS("CalculationYield/hVzCentralityPtEta_") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.vz, holderEvent.centrality, holderTrack.pt, holderTrack.eta); + if (isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value)) { + hrCalculationYield.fill(C_CS("hVzCentralityPtEta_") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.vz, holderEvent.centrality, holderTrack.pt, holderTrack.eta); } } }}; @@ -1742,23 +2057,23 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue void fillCalculationPurityByParticleSpecies() { - if (!groupAnalysis.cfgFlagsCalculationPurity.value.get(toI(ParticleSpeciesValue))) { + if (!cgAnalysis.cfgFlagsCalculationPurity.value.get(toI(ParticleSpeciesValue))) { return; } const auto fillByChargeSpecies{ [this] - requires IsValid + requires IsValidEnumValue () -> void { const auto fillByPidStrategy{ [this] - requires IsValid + requires IsValidEnumValue () -> void { - if (isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value)) { - hrCalculationPurity.fill(C_CS("CalculationPurity/pCentralityPtEtaPurity") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centrality, holderTrack.pt, holderTrack.eta, isPid(ParticleSpeciesValue), ChargeSpeciesValue>() ? 1. : 0.); + if (isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value)) { + hrCalculationPurity.fill(C_CS("pCentralityPtEtaPurity") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centrality, holderTrack.pt, holderTrack.eta, isPid(ParticleSpeciesValue), ChargeSpeciesValue>() ? 1. : 0.); } }}; @@ -1774,23 +2089,23 @@ struct PartNumFluc { } template - requires IsValid + requires IsValidEnumValue void fillCalculationFractionPrimaryByParticleSpecies(const MP& mcParticle) { - if (!groupAnalysis.cfgFlagsCalculationFractionPrimary.value.get(toI(ParticleSpeciesValue))) { + if (!cgAnalysis.cfgFlagsCalculationFractionPrimary.value.get(toI(ParticleSpeciesValue))) { return; } const auto fillByChargeSpecies{ [this, &mcParticle] - requires IsValid + requires IsValidEnumValue () -> void { const auto fillByPidStrategy{ [this, &mcParticle] - requires IsValid + requires IsValidEnumValue () -> void { - if (isPid(ParticleSpeciesValue), ChargeSpeciesValue>() && isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value)) { - hrCalculationFractionPrimary.fill(C_CS("CalculationFractionPrimary/pCentralityPtEtaFractionPrimary") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centrality, holderTrack.pt, holderTrack.eta, mcParticle.isPhysicalPrimary() ? 1. : 0.); + if (isPid(ParticleSpeciesValue), ChargeSpeciesValue>() && isPid(PidStrategyValue), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value)) { + hrCalculationFractionPrimary.fill(C_CS("pCentralityPtEtaFractionPrimary") + C_SV(getName(PidStrategyValue)) + C_SV(getName(ParticleSpeciesValue)) + C_SV(getName(ChargeSpeciesValue)), holderEvent.centrality, holderTrack.pt, holderTrack.eta, mcParticle.isPhysicalPrimary() ? 1. : 0.); } }}; @@ -1808,19 +2123,22 @@ struct PartNumFluc { void initCalculationFluctuation() { for (std::int32_t const& iParticleNumber : std::views::iota(0, NEs)) { - if (static_cast(groupAnalysis.cfgFlagsCalculationFluctuation.value.get(iParticleNumber))) { + if (static_cast(cgAnalysis.cfgFlagsCalculationFluctuation.value.get(iParticleNumber))) { for (std::int32_t const& iChargeNumber : std::views::iota(0, NEs)) { - fluctuationCalculatorsTrack[iParticleNumber][iChargeNumber]->clear(); + if (doProcessMc.value) { + fluctuationCalculatorsTrackMcParticle[iParticleNumber][iChargeNumber]->clear(); + } + fluctuationCalculatorsTrackTrack[iParticleNumber][iChargeNumber]->clear(); } } } } template - requires IsValid + requires IsValidEnumValue void calculateFluctuationByParticleNumber() { - if (!groupAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumberValue))) { + if (!cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumberValue))) { return; } @@ -1845,49 +2163,47 @@ struct PartNumFluc { } else { return false; } - }(groupTrack.cfgFlagMcParticleMomentum.value)}; - if (!isGoodMomentum(doUseMcParticleMomentum)) { - return; - } - - if constexpr (DataModeValue == DataMode::McMcParticle) { - ++holderDerivedData.nMcParticles[toI(chargeSign > 0 ? ChargeSpecies::Plus : ChargeSpecies::Minus)]; - } else { - ++holderDerivedData.nTracks[toI(chargeSign > 0 ? ChargeSpecies::Plus : ChargeSpecies::Minus)]; + }(cgTrack.cfgFlagMcParticleMomentum.value)}; + const ChargeSpecies chargeSpecies{chargeSign > 0 ? ChargeSpecies::Plus : ChargeSpecies::Minus}; + if (isGoodMomentum(ParticleNumberValue)>(doUseMcParticleMomentum)) { + if constexpr (DataModeValue == DataMode::McMcParticle) { + ++holderDerivedData.nMcParticles[toI(chargeSpecies)]; + } else { + ++holderDerivedData.nTracks[toI(chargeSpecies)]; + } } const auto calculateByParticleSpecies{ [this, chargeSign, doUseMcParticleMomentum] - requires IsValid && (ParticleNumberValue == ParticleNumber::Charge || (ParticleNumberValue == ParticleNumber::Kaon && ParticleSpeciesValue == ParticleSpecies::Kaon) || (ParticleNumberValue == ParticleNumber::Proton && ParticleSpeciesValue == ParticleSpecies::Proton)) - () -> void { - const auto calculateByChargeSpecies{ - [this, chargeSign, doUseMcParticleMomentum] - requires IsValid - () -> void { - if constexpr (DataModeValue != DataMode::RawTrack) { - if (!isPid(ParticleSpeciesValue), ChargeSpeciesValue>()) { - return; - } - } + requires IsValidEnumValue && (getValue(ParticleNumberValue) == ParticleSpeciesAll::All || getValue(ParticleNumberValue) == getValue(ParticleSpeciesValue)) + () -> bool { + if (!isGoodMomentum(ParticleSpeciesValue)>(doUseMcParticleMomentum) || (DataModeValue != DataMode::RawTrack && (chargeSign > 0 ? !isPid(ParticleSpeciesValue), ChargeSpecies::Plus>() : !isPid(ParticleSpeciesValue), ChargeSpecies::Minus>()))) { + return false; + } - const bool doUseTofPid{[](const bool doUseMcParticleMomentumValue, const double ptMcParticle, const double ptTrack, const double thresholdPtTofPid) constexpr -> bool { - if constexpr (DataModeValue == DataMode::McMcParticle) { - return ptMcParticle >= thresholdPtTofPid; - } else if constexpr (DataModeValue == DataMode::McTrack) { - return (doUseMcParticleMomentumValue ? ptMcParticle : ptTrack) >= thresholdPtTofPid; - } else { - return ptTrack >= thresholdPtTofPid; - } - }(doUseMcParticleMomentum, holderMcParticle.pt, holderTrack.pt, groupTrack.cfgThresholdsPtTofPid.value.get(toI(ParticleSpeciesValue)))}; - if constexpr (DataModeValue != DataMode::McMcParticle) { - if (!(doUseTofPid ? isPid(PidStrategy::TpcTof), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value) : isPid(PidStrategy::Tpc), getValue(ParticleSpeciesValue)>(groupTrack.cfgFlagRejectionOthers.value))) { - return; - } - } + const bool doUseTofPid{[](const bool doUseMcParticleMomentumValue, const double ptMcParticle, const double ptTrack, const double thresholdPtTofPid) constexpr -> bool { + if constexpr (DataModeValue == DataMode::McMcParticle) { + return ptMcParticle >= thresholdPtTofPid; + } else if constexpr (DataModeValue == DataMode::McTrack) { + return (doUseMcParticleMomentumValue ? ptMcParticle : ptTrack) >= thresholdPtTofPid; + } else { + return ptTrack >= thresholdPtTofPid; + } + }(doUseMcParticleMomentum, holderMcParticle.pt, holderTrack.pt, cgTrack.cfgThresholdsPtTofPid.value.get(toI(ParticleSpeciesValue)))}; + if constexpr (DataModeValue != DataMode::McMcParticle) { + if (!(doUseTofPid ? isPid(PidStrategy::TpcTof), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value) : isPid(PidStrategy::Tpc), getValue(ParticleSpeciesValue)>(cgTrack.cfgFlagRejectionOthers.value))) { + return false; + } + } - const double efficiency{doUseTofPid ? getEfficiency(doUseMcParticleMomentum) : getEfficiency(doUseMcParticleMomentum)}; // NOLINT(clang-analyzer-core.NullDereference) + const auto calculateByChargeSpecies{ + [this, chargeSign, doUseMcParticleMomentum, doUseTofPid] + requires IsValidEnumValue + () -> void { + const double efficiency{doUseTofPid ? getEfficiency(doUseMcParticleMomentum) : getEfficiency(doUseMcParticleMomentum)}; const auto fill{ [this, efficiency]() -> void { + std::array, NEs>, NEs>& fluctuationCalculatorsTrack{DataModeValue == DataMode::McMcParticle ? fluctuationCalculatorsTrackMcParticle : fluctuationCalculatorsTrackTrack}; if constexpr (ChargeSpeciesValue == ChargeSpecies::Plus) { fluctuationCalculatorsTrack[toI(ParticleNumberValue)][toI(ChargeNumber::Plus)]->fill(1., efficiency); fluctuationCalculatorsTrack[toI(ParticleNumberValue)][toI(ChargeNumber::Net)]->fill(1., efficiency); @@ -1903,59 +2219,60 @@ struct PartNumFluc { ++holderMcEvent.numbersEff[toI(ParticleNumberValue)][toI(ChargeSpeciesValue)]; fill(); } - holderDerivedData.signedEfficienciesMcParticle.push_back(HolderDerivedData::convertRound(std::copysign(std::numeric_limits::max(), chargeSign) * efficiency)); + holderDerivedData.signedEfficienciesMcParticle.push_back(HolderDerivedData::convert(std::copysign(std::numeric_limits::max(), chargeSign) * efficiency)); } else { ++holderEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpeciesValue)]; fill(); - holderDerivedData.signedEfficienciesTrack.push_back(HolderDerivedData::convertRound(std::copysign(std::numeric_limits::max(), chargeSign) * efficiency)); + holderDerivedData.signedEfficienciesTrack.push_back(HolderDerivedData::convert(std::copysign(std::numeric_limits::max(), chargeSign) * efficiency)); } }}; - if (chargeSign > 0) { // NOLINT(clang-analyzer-core.NullDereference) + if (chargeSign > 0) { calculateByChargeSpecies.template operator()(); } else { calculateByChargeSpecies.template operator()(); } + return true; }}; - if constexpr (ParticleNumberValue == ParticleNumber::Kaon) { - calculateByParticleSpecies.template operator()(); - } else if constexpr (ParticleNumberValue == ParticleNumber::Proton) { - calculateByParticleSpecies.template operator()(); - } else { // ParticleNumberValue == ParticleNumber::Charge - calculateByParticleSpecies.template operator()(); - calculateByParticleSpecies.template operator()(); - calculateByParticleSpecies.template operator()(); + if constexpr (getValue(ParticleNumberValue) == ParticleSpeciesAll::All) { + if (!calculateByParticleSpecies.template operator()() && !calculateByParticleSpecies.template operator()()) { + calculateByParticleSpecies.template operator()(); + } + } else { + calculateByParticleSpecies.template operator()(getValue(ParticleNumberValue))>(); } } - template - requires IsValid + template + requires IsValidEnumValue void fillCalculationFluctuationByParticleNumber() { - if (!groupAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumberValue))) { + if (!cgAnalysis.cfgFlagsCalculationFluctuation.value.get(toI(ParticleNumberValue))) { return; } - if constexpr (DataModeValue == DataMode::McMcParticle) { - hrCalculationFluctuation.fill(C_CS("CalculationFluctuation/hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)) + C_CS("_mc"), holderEvent.centrality, holderMcEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderMcEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); - hrCalculationFluctuation.fill(C_CS("CalculationFluctuation/hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)) + C_CS("_mcEff"), holderEvent.centrality, holderMcEvent.numbersEff[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderMcEvent.numbersEff[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); - } else { - hrCalculationFluctuation.fill(C_CS("CalculationFluctuation/hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)), holderEvent.centrality, holderEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); + if (doProcessMc.value) { + hrCalculationFluctuation.fill(C_CS("hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)) + C_CS("_mc"), holderEvent.centrality, holderMcEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderMcEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); + hrCalculationFluctuation.fill(C_CS("hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)) + C_CS("_mcEff"), holderEvent.centrality, holderMcEvent.numbersEff[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderMcEvent.numbersEff[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); } + hrCalculationFluctuation.fill(C_CS("hCentralityN") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Plus)) + C_CS("N") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeSpecies::Minus)), holderEvent.centrality, holderEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Plus)], holderEvent.numbers[toI(ParticleNumberValue)][toI(ChargeSpecies::Minus)]); const auto fillByChargeNumber{ [this] - requires IsValid + requires IsValidEnumValue () -> void { - const std::array products{fluctuationCalculatorsTrack[toI(ParticleNumberValue)][toI(ChargeNumberValue)]->getProducts()}; - for (std::int32_t const& iOrderKey : std::views::iota(0, fluctuation_calculator_base::NOrderKeys)) { - if constexpr (DataModeValue == DataMode::McMcParticle) { - hrCalculationFluctuation.fill(C_CS("CalculationFluctuation/hFluctuationCalculator") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeNumberValue)) + C_CS("_mc"), holderEvent.centrality, holderEvent.subgroupIndex, iOrderKey, products[iOrderKey]); - } else { - hrCalculationFluctuation.fill(C_CS("CalculationFluctuation/hFluctuationCalculator") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeNumberValue)), holderEvent.centrality, holderEvent.subgroupIndex, iOrderKey, products[iOrderKey]); + if (doProcessMc.value) { + const std::array products{fluctuationCalculatorsTrackMcParticle[toI(ParticleNumberValue)][toI(ChargeNumberValue)]->getProducts()}; + for (std::int32_t const& iOrderKey : std::views::iota(0, fluctuation_calculator_base::NOrderKeys)) { + hrCalculationFluctuation.fill(C_CS("hFluctuationCalculator") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeNumberValue)) + C_CS("_mc"), holderEvent.centrality, holderEvent.subgroupIndex, iOrderKey, products[iOrderKey]); } } + + const std::array products{fluctuationCalculatorsTrackTrack[toI(ParticleNumberValue)][toI(ChargeNumberValue)]->getProducts()}; + for (std::int32_t const& iOrderKey : std::views::iota(0, fluctuation_calculator_base::NOrderKeys)) { + hrCalculationFluctuation.fill(C_CS("hFluctuationCalculator") + C_SV(getName(ParticleNumberValue)) + C_SV(getName(ChargeNumberValue)), holderEvent.centrality, holderEvent.subgroupIndex, iOrderKey, products[iOrderKey]); + } }}; fillByChargeNumber.template operator()(); @@ -1964,8 +2281,8 @@ struct PartNumFluc { fillByChargeNumber.template operator()(); } - template - bool initTrack(const T& track) + template + bool setTrack(const T& track) { if (std::abs(track.sign()) != 1) { return false; @@ -1980,7 +2297,197 @@ struct PartNumFluc { holderTrack.phi = track.phi(); holderTrack.hasPid[toI(Detector::Tpc)] = (track.hasTPC() && track.tpcSignal() > 0.); holderTrack.hasPid[toI(Detector::Tof)] = (track.hasTOF() && track.beta() > 0.); - setNSigmaPid(track); + if (holderTrack.hasPid[toI(Detector::Tpc)]) { + holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Pion)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaPi(), getShiftNSigmaPid()); + holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Kaon)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaKa(), getShiftNSigmaPid()); + holderTrack.nSigmaPid[toI(Detector::Tpc)][toI(ParticleSpecies::Proton)] = HolderTrack::truncateNSigmaPid(track.tpcNSigmaPr(), getShiftNSigmaPid()); + } + if (holderTrack.hasPid[toI(Detector::Tof)]) { + holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Pion)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaPi(), getShiftNSigmaPid()); + holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Kaon)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaKa(), getShiftNSigmaPid()); + holderTrack.nSigmaPid[toI(Detector::Tof)][toI(ParticleSpecies::Proton)] = HolderTrack::truncateNSigmaPid(track.tofNSigmaPr(), getShiftNSigmaPid()); + } + + return true; + } + + template + requires IsValidEnumValue + void fillMiniTableMc(const MPs& mcParticles, const Ts& tracks, const bool isGoodNPvContributors) + { + if (!static_cast(cgAnalysis.cfgFlagsStorageMiniTable.value.get(toI(ParticleNumberValue)))) { + return; + } + + const std::optional codeCollision{mini_collision_codec::encode(cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, isGoodNPvContributors)}; + if (!codeCollision.has_value()) { + return; + } + + miniCollision(*codeCollision); + + for (const auto& mcParticle : mcParticles) { + initMcParticle(mcParticle); + + if (mcParticle.isPhysicalPrimary()) { + const auto fillByParticleSpecies{ + [this] + requires IsValidEnumValue && (getValue(ParticleNumberValue) == ParticleSpeciesAll::All || getValue(ParticleNumberValue) == getValue(ParticleSpeciesValue)) + () -> bool { + if (!isGoodMomentum(ParticleSpeciesValue)>(true) || (holderMcParticle.charge > 0 ? !isPid(ParticleSpeciesValue), ChargeSpecies::Plus>() : !isPid(ParticleSpeciesValue), ChargeSpecies::Minus>())) { + return false; + } + + const std::optional codeMcParticle{mini_mc_particle_codec::encode(holderMcParticle.pt, holderMcParticle.eta, holderMcParticle.charge, ParticleSpeciesValue)}; + if (codeMcParticle.has_value()) { + miniMcParticle(miniCollision.lastIndex(), *codeMcParticle); + } + return true; + }}; + + if constexpr (ParticleNumberValue == ParticleNumber::Charge) { + if (!fillByParticleSpecies.template operator()() && !fillByParticleSpecies.template operator()()) { + fillByParticleSpecies.template operator()(); + } + } else { + fillByParticleSpecies.template operator()(getValue(ParticleNumberValue))>(); + } + } + + if (!cgTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary()) { + const auto& tracksMatched{tracks.sliceBy(presliceTracksPerMcParticle, mcParticle.globalIndex())}; + for (const auto& track : tracksMatched) { + if (!setTrack(track) || !(cgTrack.cfgCutMinTpcChi2NCls.value < track.tpcChi2NCl()) || !(track.tpcCrossedRowsOverFindableCls() > cgTrack.cfgCutMinTpcNCrossedRowsRatio.value)) { + continue; + } + + const auto fillByParticleSpeciesAll{ + [this, &track] + requires IsValidEnumValue && (ParticleSpeciesAllValue != ParticleSpeciesAll::All) + () -> bool { + const bool doUseMcParticleMomentum{cgTrack.cfgFlagMcParticleMomentum.value}; + if (!isGoodMomentum(doUseMcParticleMomentum) || (holderTrack.sign > 0 ? !isPid() : !isPid())) { + return false; + } + + constexpr std::int32_t ParticleSpeciesIndex{toI(getValue(ParticleSpeciesAllValue))}; + const double pt{doUseMcParticleMomentum ? holderMcParticle.pt : holderTrack.pt}; + const bool doUseTofPid{pt >= cgTrack.cfgThresholdsPtTofPid.value.get(ParticleSpeciesIndex)}; + if (!(doUseTofPid ? isPid(PidStrategy::TpcTof), ParticleSpeciesAllValue, Selection::Loose>(cgTrack.cfgFlagRejectionOthers.value) : isPid(PidStrategy::Tpc), ParticleSpeciesAllValue, Selection::Loose>(cgTrack.cfgFlagRejectionOthers.value))) { + return false; + } + + const std::optional codeTrack{mini_track_codec::encode(*configSelection, track.isPVContributor(), track.itsNCls(), track.itsChi2NCl(), track.tpcNClsFound(), track.tpcChi2NCl(), track.tpcNClsCrossedRows(), track.tpcFractionSharedCls(), {getAbsNSigmaDca(), getAbsNSigmaDca()}, pt, doUseMcParticleMomentum ? holderMcParticle.eta : holderTrack.eta, holderTrack.sign, ParticleSpeciesAllValue, doUseTofPid ? std::abs(holderTrack.getNSigmaPidCombined(ParticleSpeciesIndex)) : std::abs(holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex]))}; + if (codeTrack.has_value()) { + miniTrack(miniCollision.lastIndex(), *codeTrack); + } + return true; + }}; + + if constexpr (ParticleNumberValue == ParticleNumber::Charge) { + if (!fillByParticleSpeciesAll.template operator()() && !fillByParticleSpeciesAll.template operator()()) { + fillByParticleSpeciesAll.template operator()(); + } + } else { + fillByParticleSpeciesAll.template operator()(ParticleNumberValue)>(); + } + } + } + } + } + + template + requires IsValidEnumValue + void fillMiniTableRaw(const T& tracks, const bool isGoodNPvContributors) + { + if (!static_cast(cgAnalysis.cfgFlagsStorageMiniTable.value.get(toI(ParticleNumberValue)))) { + return; + } + + const std::optional codeCollision{mini_collision_codec::encode(holderEvent.vz, holderEvent.centrality, isGoodNPvContributors)}; + if (!codeCollision.has_value()) { + return; + } + + miniCollision(*codeCollision); + + for (const auto& track : tracks) { + if (!setTrack(track) || !(cgTrack.cfgCutMinTpcChi2NCls.value < track.tpcChi2NCl()) || !(track.tpcCrossedRowsOverFindableCls() > cgTrack.cfgCutMinTpcNCrossedRowsRatio.value)) { + continue; + } + + const auto fill{ + [this, &track](const ParticleSpeciesAll particleSpeciesAll, const double absNSigmaPid) -> void { + const std::optional codeTrack{mini_track_codec::encode(*configSelection, track.isPVContributor(), track.itsNCls(), track.itsChi2NCl(), track.tpcNClsFound(), track.tpcChi2NCl(), track.tpcNClsCrossedRows(), track.tpcFractionSharedCls(), {getAbsNSigmaDca(), getAbsNSigmaDca()}, holderTrack.pt, holderTrack.eta, holderTrack.sign, particleSpeciesAll, absNSigmaPid)}; + if (codeTrack.has_value()) { + miniTrack(miniCollision.lastIndex(), *codeTrack); + } + }}; + const auto fillByParticleSpeciesAll{ + [this, &fill] + requires IsValidEnumValue + () -> bool { + if constexpr (ParticleSpeciesAllValue == ParticleSpeciesAll::All) { + if (!isGoodMomentum(ParticleNumberValue)>(false)) { + return false; + } + + fill(ParticleSpeciesAll::All, 0.); + } else { + if (!isGoodMomentum(false)) { + return false; + } + + constexpr std::int32_t ParticleSpeciesIndex{toI(getValue(ParticleSpeciesAllValue))}; + const bool doUseTofPid{holderTrack.pt >= cgTrack.cfgThresholdsPtTofPid.value.get(ParticleSpeciesIndex)}; + if (!(doUseTofPid ? isPid(PidStrategy::TpcTof), ParticleSpeciesAllValue, Selection::Loose>(cgTrack.cfgFlagRejectionOthers.value) : isPid(PidStrategy::Tpc), ParticleSpeciesAllValue, Selection::Loose>(cgTrack.cfgFlagRejectionOthers.value))) { + return false; + } + + fill(ParticleSpeciesAllValue, doUseTofPid ? std::abs(holderTrack.getNSigmaPidCombined(ParticleSpeciesIndex)) : std::abs(holderTrack.nSigmaPid[toI(Detector::Tpc)][ParticleSpeciesIndex])); + } + return true; + }}; + + if constexpr (ParticleNumberValue == ParticleNumber::Charge) { + if (!fillByParticleSpeciesAll.template operator()() && !fillByParticleSpeciesAll.template operator()() && !fillByParticleSpeciesAll.template operator()()) { + fillByParticleSpeciesAll.template operator()(); + } + } else { + if (!fillByParticleSpeciesAll.template operator()(ParticleNumberValue)>()) { + fillByParticleSpeciesAll.template operator()(); + } + } + } + } + + template + void initMcParticle(const MP& mcParticle) + { + holderMcParticle.clear(); + holderMcParticle.pdgCode = mcParticle.pdgCode(); + const TParticlePDG* const particlePdg{pdg->GetParticle(mcParticle.pdgCode())}; + if (particlePdg) { + holderMcParticle.charge = static_cast(std::rint(particlePdg->Charge())); + } else { + switch (std::abs(holderMcParticle.pdgCode) / 100000000) { + case 10: + holderMcParticle.charge = holderMcParticle.pdgCode / 10000 % 1000; + break; + default: + break; + } + } + holderMcParticle.pt = mcParticle.pt(); + holderMcParticle.eta = mcParticle.eta(); + } + + template + bool initTrack(const T& track) + { + if (!setTrack(track)) { + return false; + } if constexpr (DoInitEvent) { if (track.isPrimaryTrack()) { @@ -1998,17 +2505,15 @@ struct PartNumFluc { } } - if (groupAnalysis.cfgFlagQaRun.value && track.isPrimaryTrack()) { + if (cgAnalysis.cfgFlagQaRun.value && track.isPrimaryTrack()) { if (holderTrack.sign > 0) { fillQaRunByTrackByChargeSpecies(track); } else { fillQaRunByTrackByChargeSpecies(track); } } - - return true; } else { - if (groupAnalysis.cfgFlagQaTrack.value && track.isPrimaryTrack()) { + if (cgAnalysis.cfgFlagQaTrack.value && track.isPrimaryTrack()) { if (holderTrack.sign > 0) { fillQaTrackByChargeSpecies(track); } else { @@ -2020,7 +2525,7 @@ struct PartNumFluc { return false; } - if (groupAnalysis.cfgFlagQaDca.value) { + if (cgAnalysis.cfgFlagQaDca.value) { if (holderTrack.sign > 0) { fillQaDcaByChargeSpecies(); } else { @@ -2032,49 +2537,55 @@ struct PartNumFluc { return false; } - { - const double vz{doProcessMc.value && groupEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz}; - if (doQaAcceptance && (holderTrack.eta * vz > 0. && std::abs(vz) > (doProcessMc.value && groupEvent.cfgFlagMcCollisionVz.value ? groupEvent.cfgCutMaxAbsVzMc.value : groupEvent.cfgCutMaxAbsVz.value) - 1.)) { - fillQaAcceptancebyParticleSpeciesAll(track); - fillQaAcceptancebyParticleSpeciesAll(track); - fillQaAcceptancebyParticleSpeciesAll(track); - fillQaAcceptancebyParticleSpeciesAll(track); + if (doQaAcceptance) { + const double vz{doProcessMc.value && cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz}; + if (holderTrack.eta * vz > 0. && std::abs(vz) > (doProcessMc.value && cgEvent.cfgFlagMcCollisionVz.value ? cgEvent.cfgCutMaxAbsVzMc.value : cgEvent.cfgCutMaxAbsVz.value) - 1.) { + fillQaAcceptanceByParticleSpeciesAll(track); + fillQaAcceptanceByParticleSpeciesAll(track); + fillQaAcceptanceByParticleSpeciesAll(track); + fillQaAcceptanceByParticleSpeciesAll(track); } } - - return true; } + + return true; } - template - bool initMcParticle(const MP& mcParticle) + template + McEventSelection initMcEvent(const MC& mcCollision) { - holderMcParticle.clear(); - holderMcParticle.pdgCode = mcParticle.pdgCode(); - const TParticlePDG* const particlePdg{pdg->GetParticle(mcParticle.pdgCode())}; - if (particlePdg) { - holderMcParticle.charge = std::llrint(particlePdg->Charge()); - } else { - switch (std::abs(holderMcParticle.pdgCode) / 100000000) { - case 10: - holderMcParticle.charge = holderMcParticle.pdgCode / 10000 % 1000; - break; - default: - break; - } + holderMcEvent.clear(); + holderMcEvent.vz = mcCollision.posZ(); + + const auto fillMcEventSelection{ + [this](const McEventSelection selection) -> McEventSelection { + hrCounter.fill(C_CS("hNMcEvents"), toI(selection)); + return selection; + }}; + + fillMcEventSelection(McEventSelection::All); + + if (cgEvent.cfgFlagInelEventMc.value && !mcCollision.isInelGt0()) { + return fillMcEventSelection(McEventSelection::Inel); + } + + if (!(std::abs(holderMcEvent.vz) < cgEvent.cfgCutMaxAbsVzMc.value)) { + return fillMcEventSelection(McEventSelection::Vz); + } + + if (cgEvent.cfgFlagSingleCollisionMc.value && mcCollision.numRecoCollision() != 1) { + return fillMcEventSelection(McEventSelection::NRecoCollisions); } - holderMcParticle.pt = mcParticle.pt(); - holderMcParticle.eta = mcParticle.eta(); - return isGoodMcParticle(mcParticle); + return fillMcEventSelection(McEventSelection::Good); } template - bool initEvent(const C& collision, const Ts& tracks) + EventSelection initEvent(const C& collision, const Ts& tracks) { holderEvent.clear(); holderEvent.vz = collision.posZ(); - switch (groupEvent.cfgIndexDefinitionCentrality.value) { + switch (cgEvent.cfgIndexDefinitionCentrality.value) { case toI(CentralityDefinition::Ft0a): holderEvent.centrality = collision.centFT0A(); break; @@ -2085,20 +2596,25 @@ struct PartNumFluc { holderEvent.centrality = collision.centFT0M(); break; } - holderEvent.centralityCalibration = (groupEvent.cfgFlagDefinitionCentralitySameQa.value ? holderEvent.centrality : collision.centNTPV()); + holderEvent.centralityCalibration = (cgEvent.cfgFlagDefinitionCentralitySameQa.value ? holderEvent.centrality : collision.centNTPV()); - const auto fillEventSelection{[this](const auto... selections) -> void { - (hrCounter.fill(C_CS("hNEvents"), selections), ...); - if (groupAnalysis.cfgFlagQaCentrality.value) { - (hrQaCentrality.fill(C_CS("QaCentrality/hCentralitySelection"), holderEvent.centrality, selections), ...); - } - }}; + const auto fillEventSelection{ + [this](const EventSelection selection, const auto... selectionsAdditional) -> EventSelection + requires((std::is_arithmetic_v> && ...)) + { + hrCounter.fill(C_CS("hNEvents"), toI(selection)); + (hrCounter.fill(C_CS("hNEvents"), selectionsAdditional), ...); + if (cgAnalysis.cfgFlagQaCentrality.value) { + hrQaCentrality.fill(C_CS("hCentralitySelection"), holderEvent.centrality, toI(selection)); + (hrQaCentrality.fill(C_CS("hCentralitySelection"), holderEvent.centrality, selectionsAdditional), ...); + } + return selection; + }}; - fillEventSelection(0.); + fillEventSelection(EventSelection::All); if (!collision.has_foundBC()) { - fillEventSelection(2.); - return false; + return fillEventSelection(EventSelection::Run); } const auto& foundBc{collision.template foundBC_as()}; @@ -2107,66 +2623,66 @@ struct PartNumFluc { { const auto iter{holderCcdb.runNumbersIndicesGroupIndices.find(holderEvent.runNumber)}; if (iter == holderCcdb.runNumbersIndicesGroupIndices.end()) { - fillEventSelection(2.); - return false; + return fillEventSelection(EventSelection::Run); } std::tie(holderEvent.runIndex, holderEvent.runGroupIndex) = iter->second; } - if (holderEvent.runGroupIndex == 0 || (groupEvent.cfgFlagRejectionRunBad.value && holderEvent.runGroupIndex < 0)) { - fillEventSelection(2.); - return false; + if (holderEvent.runGroupIndex == 0 || (cgEvent.cfgFlagRejectionRunBad.value && holderEvent.runGroupIndex < 0)) { + return fillEventSelection(EventSelection::Run); } if (rctFlagsChecker.any() && !rctFlagsChecker.checkTable(collision)) { - fillEventSelection(3.); - return false; + return fillEventSelection(EventSelection::Rct); } - if (groupEvent.cfgFlagInelEvent.value && !collision.isInelGt0()) { - fillEventSelection(4.); - return false; + if (cgEvent.cfgFlagInelEvent.value && !collision.isInelGt0()) { + return fillEventSelection(EventSelection::Inel); } - for (std::int32_t const& iEvSel : std::views::iota(0, aod::evsel::EventSelectionFlags::kNsel)) { - if (((groupEvent.cfgBitsSelectionEvent.value >> iEvSel) & 1) && !collision.selection_bit(iEvSel)) { - fillEventSelection(5., 10. + iEvSel); - return false; + for (std::int32_t const& iBit : std::views::iota(0, aod::evsel::EventSelectionFlags::kNsel)) { + if (((cgEvent.cfgBitsSelectionEvent.value >> iBit) & 1) && !collision.selection_bit(iBit)) { + return fillEventSelection(EventSelection::Bits, NEs + iBit); } } if (!HolderEvent::isValidCentrality(holderEvent.centrality) || !HolderEvent::isValidCentrality(holderEvent.centralityCalibration)) { - fillEventSelection(6.); - return false; + return fillEventSelection(EventSelection::Centrality); } - if (groupAnalysis.cfgFlagQaEvent.value) { - hrQaEvent.fill(C_CS("QaEvent/hVxVy"), collision.posX(), collision.posY()); - hrQaEvent.fill(C_CS("QaEvent/hVz"), holderEvent.vz); + if (cgAnalysis.cfgFlagQaEvent.value) { + hrQaEvent.fill(C_CS("hVxVy"), collision.posX(), collision.posY()); + hrQaEvent.fill(C_CS("hVz"), holderEvent.vz); } - if (!(std::abs(holderEvent.vz) < groupEvent.cfgCutMaxAbsVz.value)) { - fillEventSelection(7.); - return false; + if (!(std::abs(holderEvent.vz) < cgEvent.cfgCutMaxAbsVz.value)) { + return fillEventSelection(EventSelection::Vz); + } + + if (cgEvent.cfgCutMaxOccupancy.value > 0) { + const std::int32_t occupancy{collision.trackOccupancyInTimeRange()}; + if (occupancy < 0 || occupancy >= cgEvent.cfgCutMaxOccupancy.value) { + return fillEventSelection(EventSelection::Occupancy); + } } - if (groupAnalysis.cfgFlagQaRun.value) { - hrQaRun.fill(C_CS("QaRun/pRunIndexVx"), holderEvent.runIndex, collision.posX()); - hrQaRun.fill(C_CS("QaRun/pRunIndexVy"), holderEvent.runIndex, collision.posY()); - hrQaRun.fill(C_CS("QaRun/pRunIndexVz"), holderEvent.runIndex, holderEvent.vz); - hrQaRun.fill(C_CS("QaRun/pRunIndexMultiplicityFt0a"), holderEvent.runIndex, collision.multZeqFT0A()); - hrQaRun.fill(C_CS("QaRun/pRunIndexMultiplicityFt0c"), holderEvent.runIndex, collision.multZeqFT0C()); + if (cgAnalysis.cfgFlagQaRun.value) { + hrQaRun.fill(C_CS("pRunIndexVx"), holderEvent.runIndex, collision.posX()); + hrQaRun.fill(C_CS("pRunIndexVy"), holderEvent.runIndex, collision.posY()); + hrQaRun.fill(C_CS("pRunIndexVz"), holderEvent.runIndex, holderEvent.vz); + hrQaRun.fill(C_CS("pRunIndexMultiplicityFt0a"), holderEvent.runIndex, collision.multZeqFT0A()); + hrQaRun.fill(C_CS("pRunIndexMultiplicityFt0c"), holderEvent.runIndex, collision.multZeqFT0C()); if (const double centrality{collision.centNTPV()}; HolderEvent::isValidCentrality(centrality)) { - hrQaRun.fill(C_CS("QaRun/pRunIndexCentralityNtpv"), holderEvent.runIndex, centrality); + hrQaRun.fill(C_CS("pRunIndexCentralityNtpv"), holderEvent.runIndex, centrality); } if (const double centrality{collision.centFT0A()}; HolderEvent::isValidCentrality(centrality)) { - hrQaRun.fill(C_CS("QaRun/pRunIndexCentralityFt0a"), holderEvent.runIndex, centrality); + hrQaRun.fill(C_CS("pRunIndexCentralityFt0a"), holderEvent.runIndex, centrality); } if (const double centrality{collision.centFT0C()}; HolderEvent::isValidCentrality(centrality)) { - hrQaRun.fill(C_CS("QaRun/pRunIndexCentralityFt0c"), holderEvent.runIndex, centrality); + hrQaRun.fill(C_CS("pRunIndexCentralityFt0c"), holderEvent.runIndex, centrality); } if (const double centrality{collision.centFT0M()}; HolderEvent::isValidCentrality(centrality)) { - hrQaRun.fill(C_CS("QaRun/pRunIndexCentralityFt0m"), holderEvent.runIndex, centrality); + hrQaRun.fill(C_CS("pRunIndexCentralityFt0m"), holderEvent.runIndex, centrality); } } @@ -2182,79 +2698,185 @@ struct PartNumFluc { } } - if (groupAnalysis.cfgFlagQaRun.value) { + if (cgAnalysis.cfgFlagQaRun.value) { fillQaRunByEventByChargeSpecies(); fillQaRunByEventByChargeSpecies(); } - if (groupAnalysis.cfgFlagQaEvent.value) { - hrQaEvent.fill(C_CS("QaEvent/hNPvContributorsNGlobalTracks"), holderEvent.getNPvContributors(), holderEvent.getNGlobalTracks()); + if (cgAnalysis.cfgFlagQaEvent.value) { + hrQaEvent.fill(C_CS("hNPvContributorsNGlobalTracks"), holderEvent.getNPvContributors(), holderEvent.getNGlobalTracks()); if (holderEvent.getNGlobalTracks() > 0) { - hrQaEvent.fill(C_CS("QaEvent/hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Xy)), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); - hrQaEvent.fill(C_CS("QaEvent/hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Z)), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); + hrQaEvent.fill(C_CS("hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Xy)), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); + hrQaEvent.fill(C_CS("hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Z)), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); } - hrQaEvent.fill(C_CS("QaEvent/hNTofBetaNGlobalTracks"), holderEvent.getNTofBeta(), holderEvent.getNGlobalTracks()); + hrQaEvent.fill(C_CS("hNTofBetaNGlobalTracks"), holderEvent.getNTofBeta(), holderEvent.getNGlobalTracks()); } - if (!(holderEvent.getNPvContributors() - holderEvent.getNGlobalTracks() > groupEvent.cfgCutMinDeviationNPvContributors.value)) { - fillEventSelection(8.); - return false; + if (!(holderEvent.getNPvContributors() - holderEvent.getNGlobalTracks() > cgEvent.cfgCutMinDeviationNPvContributors.value)) { + return fillEventSelection(EventSelection::NPvContributors); } - fillEventSelection(1.); readCcdb(); - if (groupAnalysis.cfgFlagQaEvent.value) { + if (cgAnalysis.cfgFlagQaEvent.value) { if (holderEvent.getNGlobalTracks() > 0) { - hrQaEvent.fill(C_CS("QaEvent/hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Xy)) + C_CS("_nPvContributorsCut"), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); - hrQaEvent.fill(C_CS("QaEvent/hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Z)) + C_CS("_nPvContributorsCut"), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); + hrQaEvent.fill(C_CS("hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Xy)) + C_CS("_nPvContributorsCut"), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); + hrQaEvent.fill(C_CS("hNGlobalTracks") + C_SV(getName(DcaMeasure::Mean)) + C_CS("Dca") + C_SV(getName(DcaAxis::Z)) + C_CS("_nPvContributorsCut"), holderEvent.getNGlobalTracks(), holderEvent.getMeasureDca()); } - hrQaEvent.fill(C_CS("QaEvent/hNTofBetaNGlobalTracks_nPvContributorsCut"), holderEvent.getNTofBeta(), holderEvent.getNGlobalTracks()); + hrQaEvent.fill(C_CS("hNTofBetaNGlobalTracks_nPvContributorsCut"), holderEvent.getNTofBeta(), holderEvent.getNGlobalTracks()); } - if (groupAnalysis.cfgFlagQaCentrality.value) { - hrQaCentrality.fill(C_CS("QaCentrality/hCentralityMultiplicity"), holderEvent.centrality, collision.multNTracksPVeta1()); + if (cgAnalysis.cfgFlagQaCentrality.value) { + hrQaCentrality.fill(C_CS("hCentralityMultiplicity"), holderEvent.centrality, collision.multNTracksPVeta1()); } - return true; + return fillEventSelection(EventSelection::Good); } - template - bool initMcEvent(const MC& mcCollision) + void processMc(const soa::Filtered::iterator& mcCollision, const aod::McParticles& mcParticles, const soa::SmallGroups& collisions, const soa::Filtered& tracksUngrouped, const aod::BCsWithTimestamps&) { - holderMcEvent.clear(); - holderMcEvent.vz = mcCollision.posZ(); + const McEventSelection mcEventSelection{initMcEvent(mcCollision)}; + if (mcEventSelection != McEventSelection::Good) { + return; + } - hrCounter.fill(C_CS("hNMcEvents"), 0.); + for (const auto& collision : collisions) { + if (collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } - if (groupEvent.cfgFlagInelEventMc.value && !mcCollision.isInelGt0()) { - hrCounter.fill(C_CS("hNMcEvents"), 2.); - return false; - } + const auto& tracks{tracksUngrouped.sliceBy(presliceTracksPerCollision, collision.globalIndex())}; - if (!(std::abs(holderMcEvent.vz) < groupEvent.cfgCutMaxAbsVzMc.value)) { - hrCounter.fill(C_CS("hNMcEvents"), 3.); - return false; - } + const EventSelection eventSelection{initEvent(collision, tracks)}; + if ((eventSelection == EventSelection::Good || toI(eventSelection) >= toI(EventSelection::NPvContributors)) && doStorageMiniTable && gRandom->Rndm() * cgEvent.cfgFactorStorageMiniTable.value < 1.) { + const bool isGoodNPvContributors{eventSelection != EventSelection::NPvContributors}; + readCcdb(); + fillMiniTableMc(mcParticles, tracks, isGoodNPvContributors); + fillMiniTableMc(mcParticles, tracks, isGoodNPvContributors); + fillMiniTableMc(mcParticles, tracks, isGoodNPvContributors); + } + if (eventSelection != EventSelection::Good) { + continue; + } - if (groupEvent.cfgFlagSingleCollisionMc.value && mcCollision.numRecoCollision() != 1) { - hrCounter.fill(C_CS("hNMcEvents"), 4.); - return false; - } + if (cgAnalysis.cfgFlagQaMc.value) { + hrQaMc.fill(C_CS("hCentralityVzMcDeltaVz"), holderEvent.centralityCalibration, holderMcEvent.vz, holderEvent.vz - holderMcEvent.vz); + } + + if (cgAnalysis.cfgFlagQaTrack.value || cgAnalysis.cfgFlagQaDca.value || doQaAcceptance || doQaPhi || doQaPid || cgAnalysis.cfgFlagQaMc.value || doCalculationYield || doCalculationPurity || doCalculationFractionPrimary || doCalculationFluctuation) { + if (doCalculationFluctuation) { + holderEvent.subgroupIndex = gRandom->Integer(cgEvent.cfgNSubgroups.value); + initCalculationFluctuation(); + holderDerivedData.clear(); + } - hrCounter.fill(C_CS("hNMcEvents"), 1.); + for (const auto& mcParticle : mcParticles) { + initMcParticle(mcParticle); - return true; + if (mcParticle.isPhysicalPrimary()) { + if (doCalculationYield) { + fillCalculationYieldByParticleSpecies(); + fillCalculationYieldByParticleSpecies(); + fillCalculationYieldByParticleSpecies(); + } + + if (doCalculationFluctuation) { + calculateFluctuationByParticleNumber(); + calculateFluctuationByParticleNumber(); + calculateFluctuationByParticleNumber(); + } + } + + const auto& tracksMatched{tracks.sliceBy(presliceTracksPerMcParticle, mcParticle.globalIndex())}; + for (const auto& track : tracksMatched) { + if (!initTrack(track)) { + continue; + } + + if (doQaPhi) { + fillQaPhiByParticleSpeciesAll(); + fillQaPhiByParticleSpeciesAll(); + fillQaPhiByParticleSpeciesAll(); + fillQaPhiByParticleSpeciesAll(); + } + + if (doQaPid) { + fillQaPidByParticleSpeciesAll(track); + fillQaPidByParticleSpeciesAll(track); + fillQaPidByParticleSpeciesAll(track); + fillQaPidByParticleSpeciesAll(track); + } + + if (cgAnalysis.cfgFlagQaMc.value && (!cgTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { + hrQaMc.fill(C_CS("hCentralityPtMcEtaMcDeltaPt"), holderEvent.centralityCalibration, holderMcParticle.pt, holderMcParticle.eta, holderTrack.pt - holderMcParticle.pt); + hrQaMc.fill(C_CS("hCentralityPtMcEtaMcDeltaEta"), holderEvent.centralityCalibration, holderMcParticle.pt, holderMcParticle.eta, holderTrack.eta - holderMcParticle.eta); + } + + if (doCalculationYield && (!cgTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { + fillCalculationYieldByParticleSpecies(); + fillCalculationYieldByParticleSpecies(); + fillCalculationYieldByParticleSpecies(); + } + + if (doCalculationPurity && (!cgTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { + fillCalculationPurityByParticleSpecies(); + fillCalculationPurityByParticleSpecies(); + fillCalculationPurityByParticleSpecies(); + } + + if (doCalculationFractionPrimary) { + fillCalculationFractionPrimaryByParticleSpecies(mcParticle); + fillCalculationFractionPrimaryByParticleSpecies(mcParticle); + fillCalculationFractionPrimaryByParticleSpecies(mcParticle); + } + + if (doCalculationFluctuation && (!cgTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { + calculateFluctuationByParticleNumber(); + calculateFluctuationByParticleNumber(); + calculateFluctuationByParticleNumber(); + } + } + } + + if (doCalculationFluctuation) { + fillCalculationFluctuationByParticleNumber(); + fillCalculationFluctuationByParticleNumber(); + fillCalculationFluctuationByParticleNumber(); + if (const std::optional codeCollision{mini_collision_codec::encode(cgEvent.cfgFlagMcCollisionVz.value ? holderMcEvent.vz : holderEvent.vz, holderEvent.centrality, true)}; codeCollision.has_value()) { + tinyMcCollision(holderDerivedData.nMcParticles[toI(ChargeSpecies::Plus)], holderDerivedData.nMcParticles[toI(ChargeSpecies::Minus)]); + tinyCollision(*codeCollision, holderDerivedData.nTracks[toI(ChargeSpecies::Plus)], holderDerivedData.nTracks[toI(ChargeSpecies::Minus)]); + for (aod::tiny_mc_particle::SignedEfficiency::type const& signedEfficiency : holderDerivedData.signedEfficienciesMcParticle) { + tinyMcParticle(tinyMcCollision.lastIndex(), signedEfficiency); + } + for (aod::tiny_track::SignedEfficiency::type const& signedEfficiency : holderDerivedData.signedEfficienciesTrack) { + tinyTrack(tinyCollision.lastIndex(), signedEfficiency); + } + } + } + } + } } void processRaw(const soa::Filtered::iterator& collision, const soa::Filtered& tracks, const aod::BCsWithTimestamps&) { - if (!initEvent(collision, tracks) || (!groupAnalysis.cfgFlagQaTrack.value && !groupAnalysis.cfgFlagQaDca.value && !doQaAcceptance && !doQaPhi && !doQaPid && !doCalculationYield && !doCalculationFluctuation)) { + const EventSelection eventSelection{initEvent(collision, tracks)}; + if ((eventSelection == EventSelection::Good || toI(eventSelection) >= toI(EventSelection::NPvContributors)) && doStorageMiniTable && gRandom->Rndm() * cgEvent.cfgFactorStorageMiniTable.value < 1.) { + const bool isGoodNPvContributors{eventSelection != EventSelection::NPvContributors}; + readCcdb(); + fillMiniTableRaw(tracks, isGoodNPvContributors); + fillMiniTableRaw(tracks, isGoodNPvContributors); + fillMiniTableRaw(tracks, isGoodNPvContributors); + } + if (eventSelection != EventSelection::Good) { + return; + } + + if (!cgAnalysis.cfgFlagQaTrack.value && !cgAnalysis.cfgFlagQaDca.value && !doQaAcceptance && !doQaPhi && !doQaPid && !doCalculationYield && !doCalculationFluctuation) { return; } if (doCalculationFluctuation) { - holderEvent.subgroupIndex = gRandom->Integer(groupEvent.cfgNSubgroups.value); + holderEvent.subgroupIndex = gRandom->Integer(cgEvent.cfgNSubgroups.value); initCalculationFluctuation(); holderDerivedData.clear(); } @@ -2292,153 +2914,20 @@ struct PartNumFluc { } if (doCalculationFluctuation) { - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); - miniCollision(HolderDerivedData::convertFloor(holderEvent.vz * 10.), HolderDerivedData::convertFloor(holderEvent.centrality * 500.), holderDerivedData.nTracks[toI(ChargeSpecies::Plus)], holderDerivedData.nTracks[toI(ChargeSpecies::Minus)]); - for (std::int16_t const& signedEfficiency : holderDerivedData.signedEfficienciesTrack) { - miniTrack(miniCollision.lastIndex(), signedEfficiency); - } - } - } - - void processMc(const soa::Filtered::iterator& mcCollision, const aod::McParticles& mcParticles, const soa::SmallGroups& collisions, const soa::Filtered& tracksUngrouped, const aod::BCsWithTimestamps&) - { - if (!initMcEvent(mcCollision)) { - return; - } - - for (const auto& collision : collisions) { - if (groupEvent.cfgFlagBestCollisionMc.value && collision.globalIndex() != mcCollision.bestCollisionIndex()) { - continue; - } - - const auto& tracks{tracksUngrouped.sliceBy(presliceTracksPerCollision, collision.globalIndex())}; - - if (!initEvent(collision, tracks)) { - continue; - } - - if (groupAnalysis.cfgFlagQaMc.value) { - hrQaMc.fill(C_CS("QaMc/hCentralityVzMcDeltaVz"), holderEvent.centralityCalibration, holderMcEvent.vz, holderEvent.vz - holderMcEvent.vz); - } - - if (doCalculationYield || doCalculationFluctuation) { - if (doCalculationFluctuation) { - holderEvent.subgroupIndex = gRandom->Integer(groupEvent.cfgNSubgroups.value); - initCalculationFluctuation(); - holderDerivedData.clear(); - } - - for (const auto& mcParticle : mcParticles) { - if (!initMcParticle(mcParticle)) { - continue; - } - - if (doCalculationYield) { - fillCalculationYieldByParticleSpecies(); - fillCalculationYieldByParticleSpecies(); - fillCalculationYieldByParticleSpecies(); - } - - if (doCalculationFluctuation) { - calculateFluctuationByParticleNumber(); - calculateFluctuationByParticleNumber(); - calculateFluctuationByParticleNumber(); - } - } - - if (doCalculationFluctuation) { - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); + fillCalculationFluctuationByParticleNumber(); + fillCalculationFluctuationByParticleNumber(); + fillCalculationFluctuationByParticleNumber(); + if (const std::optional codeCollision{mini_collision_codec::encode(holderEvent.vz, holderEvent.centrality, true)}; codeCollision.has_value()) { + tinyCollision(*codeCollision, holderDerivedData.nTracks[toI(ChargeSpecies::Plus)], holderDerivedData.nTracks[toI(ChargeSpecies::Minus)]); + for (aod::tiny_track::SignedEfficiency::type const& signedEfficiency : holderDerivedData.signedEfficienciesTrack) { + tinyTrack(tinyCollision.lastIndex(), signedEfficiency); } } - - if (groupAnalysis.cfgFlagQaTrack.value || groupAnalysis.cfgFlagQaDca.value || doQaAcceptance || doQaPhi || doQaPid || groupAnalysis.cfgFlagQaMc.value || doCalculationYield || doCalculationPurity || doCalculationFractionPrimary || doCalculationFluctuation) { - if (doCalculationFluctuation) { - initCalculationFluctuation(); - } - - for (const auto& track : tracks) { - if (!track.has_mcParticle()) { - continue; - } - - const auto& mcParticle{track.template mcParticle_as()}; - if (!mcParticle.has_mcCollision() || mcParticle.mcCollisionId() != mcCollision.globalIndex()) { - continue; - } - - if (!initTrack(track) || !initMcParticle(mcParticle)) { - continue; - } - - if (doQaPhi) { - fillQaPhiByParticleSpeciesAll(); - fillQaPhiByParticleSpeciesAll(); - fillQaPhiByParticleSpeciesAll(); - fillQaPhiByParticleSpeciesAll(); - } - - if (doQaPid) { - fillQaPidByParticleSpeciesAll(track); - fillQaPidByParticleSpeciesAll(track); - fillQaPidByParticleSpeciesAll(track); - fillQaPidByParticleSpeciesAll(track); - } - - if (groupAnalysis.cfgFlagQaMc.value && (!groupTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { - hrQaMc.fill(C_CS("QaMc/hCentralityPtMcEtaMcDeltaPt"), holderEvent.centralityCalibration, holderMcParticle.pt, holderMcParticle.eta, holderTrack.pt - holderMcParticle.pt); - hrQaMc.fill(C_CS("QaMc/hCentralityPtMcEtaMcDeltaEta"), holderEvent.centralityCalibration, holderMcParticle.pt, holderMcParticle.eta, holderTrack.eta - holderMcParticle.eta); - } - - if (doCalculationYield && (!groupTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { - fillCalculationYieldByParticleSpecies(); - fillCalculationYieldByParticleSpecies(); - fillCalculationYieldByParticleSpecies(); - } - - if (doCalculationPurity && (!groupTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { - fillCalculationPurityByParticleSpecies(); - fillCalculationPurityByParticleSpecies(); - fillCalculationPurityByParticleSpecies(); - } - - if (doCalculationFractionPrimary) { - fillCalculationFractionPrimaryByParticleSpecies(mcParticle); - fillCalculationFractionPrimaryByParticleSpecies(mcParticle); - fillCalculationFractionPrimaryByParticleSpecies(mcParticle); - } - - if (doCalculationFluctuation && (!groupTrack.cfgFlagMcParticlePhysicalPrimary.value || mcParticle.isPhysicalPrimary())) { - calculateFluctuationByParticleNumber(); - calculateFluctuationByParticleNumber(); - calculateFluctuationByParticleNumber(); - } - } - - if (doCalculationFluctuation) { - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); - fillCalculationFluctuationByParticleNumber(); - miniMcCollision(holderDerivedData.nMcParticles[toI(ChargeSpecies::Plus)], holderDerivedData.nMcParticles[toI(ChargeSpecies::Minus)]); - miniCollision(HolderDerivedData::convertFloor(holderEvent.vz * 10.), HolderDerivedData::convertFloor(holderEvent.centrality * 500.), holderDerivedData.nTracks[toI(ChargeSpecies::Plus)], holderDerivedData.nTracks[toI(ChargeSpecies::Minus)]); - for (std::int16_t const& signedEfficiency : holderDerivedData.signedEfficienciesMcParticle) { - miniMcParticle(miniMcCollision.lastIndex(), signedEfficiency); - } - for (std::int16_t const& signedEfficiency : holderDerivedData.signedEfficienciesTrack) { - miniTrack(miniCollision.lastIndex(), signedEfficiency); - } - } - } - - break; } } - PROCESS_SWITCH_FULL(PartNumFluc, processRaw, ProcessRaw, "Flag of processing raw data", true); PROCESS_SWITCH_FULL(PartNumFluc, processMc, ProcessMc, "Flag of processing MC data", false); + PROCESS_SWITCH_FULL(PartNumFluc, processRaw, ProcessRaw, "Flag of processing raw data", true); }; WorkflowSpec defineDataProcessing(const ConfigContext& configContext)