diff --git a/PWGLF/Core/Omega2012AnalysisCore.h b/PWGLF/Core/Omega2012AnalysisCore.h new file mode 100644 index 00000000000..d0d3600d7b4 --- /dev/null +++ b/PWGLF/Core/Omega2012AnalysisCore.h @@ -0,0 +1,1269 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file Omega2012AnalysisCore.h +/// \brief Omega(2012) selection, candidate enumeration and truth classification, separately for the two decay modes +/// \author Bong-Hwi Lim +/// +/// Mode A (DecayMode::XiK0s): Omega(2012)- -> Xi- K0S. +/// Mode B (DecayMode::Xi1530K): Omega(2012)- -> Xi(1530)0 K- -> Xi- pi+ K- (charged kaon track). +/// The two modes are never merged: each has its own enumeration, hooks, cut flow, pass bits and truth classification. +/// The Xi selection is common to both modes. The core owns the selection and its cut-flow instrumentation +/// (CutFlow/*, ML/*); the tasks own their output histograms and fill them from the hooks. + +#ifndef PWGLF_CORE_OMEGA2012ANALYSISCORE_H_ +#define PWGLF_CORE_OMEGA2012ANALYSISCORE_H_ + +#include "PWGLF/Core/Omega2012MlFeatures.h" +#include "PWGLF/Core/ResoAnalysisSelectionCore.h" + +#include +#include +#include +#include +#include +#include + +#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::omega2012 +{ + +using o2::analysis::resonance::EventCuts; +using o2::analysis::resonance::PIDCutConfig; +using o2::analysis::resonance::ResoAnalysisSelectionCore; +using o2::analysis::resonance::TrackCuts; +using o2::analysis::resonance::TrackStage; + +enum class DecayMode : uint8_t { + XiK0s = 0, // Omega(2012)- -> Xi- K0S + Xi1530K = 1 // Omega(2012)- -> Xi(1530)0 K- -> Xi- pi+ K- +}; + +inline constexpr int PdgXi1530Zero = 3324; // Xi(1530)0, not in o2::constants::physics::Pdg +inline constexpr float SmallNumber = 1e-10f; // avoids division by zero (value of the original task) +inline constexpr float MaxDCAV0ToPV = 1.0f; // maximum K0S DCA to PV (value of the original task) +inline constexpr float CascDCAzPtExponent = -1.1f; +inline constexpr const char* TrackCutsPrefix = "trk."; // JSON prefix of the mode-B track selection +inline constexpr int NXiK0sCandidateStages = 4; +inline constexpr int NXi1530KCandidateStages = 4; +inline constexpr int NGeneratedChannels = 3; + +// Last stage passed by a cascade (Xi selection, common to both modes). +enum XiStage : int { + kXiInput = 0, // failed |eta| or pT + kXiKinematics = 1, // passed |eta| and pT + kXiDCA = 2, // passed the cascade DCA to PV + kXiV0Topology = 3, // passed the Lambda topology and mass window + kXiCascTopology = 4, // passed the cascade topology + kXiMass = 5, // passed the Xi mass window + kXiSelected = kXiMass, + kXiNStages = 6 +}; + +// Last stage passed by a V0 (K0S selection of mode A). +enum K0sStage : int { + kK0sInput = 0, // failed |eta| or pT + kK0sKinematics = 1, // passed |eta| and pT + kK0sTopology = 2, // passed cosPA, daughter DCAs, radius and DCA to PV + kK0sLifetime = 3, // passed the proper lifetime and the minimum qT + kK0sMass = 4, // passed the K0S mass window and the (anti)Lambda rejection + kK0sDaughters = 5, // passed the daughter pion TPC nSigma and crossed rows + kK0sArmenteros = 6, // passed the Armenteros qT > coefficient * |alpha| cut + kK0sSelected = kK0sArmenteros, + kK0sNStages = 7 +}; + +// The value is the species index of the PID configuration in ResoAnalysisSelectionCore. +enum class Species : int { + Pion = 0, + Kaon = 1 +}; + +// Cumulative pass bits of mode A (see LFOmega2012MlTables.h) +enum XiK0sPassBit : uint16_t { + kXiK0sPassLoose = 1, + kXiK0sPassXi = 2, + kXiK0sPassK0s = 4, + kXiK0sPassKinematic = 8 +}; +inline constexpr uint16_t XiK0sPassBitsSelected = kXiK0sPassLoose | kXiK0sPassXi | kXiK0sPassK0s | kXiK0sPassKinematic; + +// Cumulative pass bits of mode B (see LFOmega2012MlTables.h) +enum Xi1530KPassBit : uint16_t { + kXi1530KPassLoose = 1, + kXi1530KPassXi = 2, + kXi1530KPassQuality = 4, + kXi1530KPassPID = 8, + kXi1530KPassWindow = 16 +}; +inline constexpr uint16_t Xi1530KPassBitsSelected = kXi1530KPassLoose | kXi1530KPassXi | kXi1530KPassQuality | kXi1530KPassPID | kXi1530KPassWindow; + +enum class XiK0sTruth : uint8_t { + None = 0, + Matched = 1 +}; + +enum class Xi1530KTruth : uint8_t { + None = 0, + Matched = 1 +}; + +// Immediate decay channel of a generated Omega(2012) +enum class GeneratedChannel : uint8_t { + Other = 0, + XiK0s = 1, + Xi1530K = 2 +}; + +// Configurable groups without prefix: the JSON keys are the plain configurable names of the original task. + +/// Xi (cascade) selection, common to both modes. cMinPtcut is also the K0S minimum pT (as in the original task). +struct XiCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cMinPtcut{"cMinPtcut", 0.15, "Minimum pT for candidates"}; + o2::framework::Configurable cMaxEtaCut{"cMaxEtaCut", 0.8, "Maximum |eta|"}; + o2::framework::Configurable cDCAxyToPVByPtCascP0{"cDCAxyToPVByPtCascP0", 999., "Cascade DCAxy p0"}; + o2::framework::Configurable cDCAxyToPVByPtCascExp{"cDCAxyToPVByPtCascExp", 1., "Cascade DCAxy exp"}; + o2::framework::Configurable cDCAxyToPVAsPtForCasc{"cDCAxyToPVAsPtForCasc", true, "Use pt-dep DCAxy cut (casc)"}; + o2::framework::Configurable cDCAzToPVAsPtForCasc{"cDCAzToPVAsPtForCasc", true, "Use pt-dep DCAz cut (casc)"}; + // V0 topology inside cascade (Lambda) + o2::framework::Configurable cDCALambdaDaugtherscut{"cDCALambdaDaugtherscut", 0.7, "Λ daughters DCA cut"}; + o2::framework::Configurable cDCALambdaToPVcut{"cDCALambdaToPVcut", 0.02, "Λ DCA to PV min"}; + o2::framework::Configurable cDCAPionToPVcut{"cDCAPionToPVcut", 0.06, "π DCA to PV min"}; + o2::framework::Configurable cDCAProtonToPVcut{"cDCAProtonToPVcut", 0.07, "p DCA to PV min"}; + o2::framework::Configurable cV0CosPACutPtDepP0{"cV0CosPACutPtDepP0", 0.25, "V0 CosPA p0"}; + o2::framework::Configurable cV0CosPACutPtDepP1{"cV0CosPACutPtDepP1", 0.022, "V0 CosPA p1"}; + o2::framework::Configurable cMaxV0radiuscut{"cMaxV0radiuscut", 200., "V0 radius max"}; + o2::framework::Configurable cMinV0radiuscut{"cMinV0radiuscut", 2.5, "V0 radius min"}; + o2::framework::Configurable cMasswindowV0cut{"cMasswindowV0cut", 0.005, "Λ mass window for cascade V0"}; + // Cascade topology + o2::framework::Configurable cDCABachlorToPVcut{"cDCABachlorToPVcut", 0.06, "Bachelor DCA to PV min"}; + o2::framework::Configurable cDCAXiDaugthersCutPtRangeLower{"cDCAXiDaugthersCutPtRangeLower", 1., "Xi pt low boundary"}; + o2::framework::Configurable cDCAXiDaugthersCutPtRangeUpper{"cDCAXiDaugthersCutPtRangeUpper", 4., "Xi pt high boundary"}; + o2::framework::Configurable cDCAXiDaugthersCutPtDepLower{"cDCAXiDaugthersCutPtDepLower", 0.8, "Xi daugh DCA (pt cDCAXiDaugthersCutPtDepMiddle{"cDCAXiDaugthersCutPtDepMiddle", 0.5, "Xi daugh DCA (low<=pt cDCAXiDaugthersCutPtDepUpper{"cDCAXiDaugthersCutPtDepUpper", 0.2, "Xi daugh DCA (pt>=high)"}; + o2::framework::Configurable cCosPACascCutPtDepP0{"cCosPACascCutPtDepP0", 0.2, "Cascade CosPA p0"}; + o2::framework::Configurable cCosPACascCutPtDepP1{"cCosPACascCutPtDepP1", 0.022, "Cascade CosPA p1"}; + o2::framework::Configurable cMaxCascradiuscut{"cMaxCascradiuscut", 200., "Cascade radius max"}; + o2::framework::Configurable cMinCascradiuscut{"cMinCascradiuscut", 1.1, "Cascade radius min"}; + o2::framework::Configurable cMasswindowCasccut{"cMasswindowCasccut", 0.008, "Xi mass window"}; + o2::framework::Configurable cMassXiminus{"cMassXiminus", 1.32171, "Xi mass (GeV/c^2)"}; // PDG +}; + +/// K0S selection and Xi-K0S kinematic (opening-angle) cut of mode A +struct K0sCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cKinCuts{"cKinCuts", false, "Kinematic cuts for Xi-K0s opening angle"}; + o2::framework::Configurable> cKinCutsPt{"cKinCutsPt", {0.0, 0.4, 0.6, 0.8, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0, 4.0, 5.0, 6.0, 1e10}, "Omega(2012) pT bins for kinematic cuts"}; + o2::framework::Configurable> cKinLowerCutsAlpha{"cKinLowerCutsAlpha", {1.5, 1.0, 0.5, 0.3, 0.2, 0.15, 0.1, 0.08, 0.07, 0.06, 0.04, 0.02, 0.02}, "Lower cut on Xi-K0s opening angle"}; + o2::framework::Configurable> cKinUpperCutsAlpha{"cKinUpperCutsAlpha", {3.0, 2.0, 1.5, 1.4, 1.0, 0.8, 0.6, 0.5, 0.45, 0.35, 0.3, 0.25, 0.2}, "Upper cut on Xi-K0s opening angle"}; + o2::framework::Configurable cK0sMinCosPA{"cK0sMinCosPA", 0.98, "K0s minimum pointing angle cosine"}; + o2::framework::Configurable cK0sMaxDaughDCA{"cK0sMaxDaughDCA", 0.5, "K0s daughter DCA Maximum"}; + o2::framework::Configurable cK0sMassWindow{"cK0sMassWindow", 0.025, "Mass window for K0s selection (GeV/c^2)"}; + o2::framework::Configurable cMaxV0Etacut{"cMaxV0Etacut", 0.8, "V0 maximum eta cut"}; + o2::framework::Configurable cK0sProperLifetimeMax{"cK0sProperLifetimeMax", 20.0, "K0s proper lifetime max (cm/c)"}; + o2::framework::Configurable cK0sArmenterosQtMin{"cK0sArmenterosQtMin", 0.0, "K0s Armenteros qt min"}; + o2::framework::Configurable cK0sArmenterosAlphaCoeff{"cK0sArmenterosAlphaCoeff", 0.2, "K0s Armenteros alpha coefficient"}; + o2::framework::Configurable cK0sDauPosDCAtoPVMin{"cK0sDauPosDCAtoPVMin", 0.05, "K0s positive daughter DCA to PV min"}; + o2::framework::Configurable cK0sDauNegDCAtoPVMin{"cK0sDauNegDCAtoPVMin", 0.05, "K0s negative daughter DCA to PV min"}; + o2::framework::Configurable cK0sRadiusMin{"cK0sRadiusMin", 0.5, "K0s decay radius min"}; + o2::framework::Configurable cK0sRadiusMax{"cK0sRadiusMax", 200.0, "K0s decay radius max"}; + o2::framework::Configurable cK0sCrossMassRejection{"cK0sCrossMassRejection", true, "Enable Lambda mass rejection for K0s"}; + o2::framework::Configurable cK0sCrossMassRejectionWindow{"cK0sCrossMassRejectionWindow", 0.01, "Lambda mass rejection window for K0s (GeV/c^2)"}; + o2::framework::Configurable cK0sDaughterPiTPCNSigmaMax{"cK0sDaughterPiTPCNSigmaMax", 5.0, "Maximum TPC NSigma for K0s daughter pions"}; + o2::framework::Configurable cK0sPosDaughterMinCrossedRows{"cK0sPosDaughterMinCrossedRows", 50, "Minimum TPC crossed rows for K0s positive daughter"}; + o2::framework::Configurable cK0sNegDaughterMinCrossedRows{"cK0sNegDaughterMinCrossedRows", 50, "Minimum TPC crossed rows for K0s negative daughter"}; +}; + +/// Xi(1530)0 K- selection of mode B +struct Xi1530KCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cXi1530Mass{"cXi1530Mass", 1.53, "Xi(1530) mass (GeV/c^2)"}; + o2::framework::Configurable cXi1530MassWindow{"cXi1530MassWindow", 0.01, "Xi(1530) mass window (GeV/c^2)"}; + o2::framework::Configurable cXi1530UseMassWindow{"cXi1530UseMassWindow", true, "Require m(Xi pi) inside the Xi(1530) mass window"}; + o2::framework::Configurable cXi1530KFillWrongSign{"cXi1530KFillWrongSign", true, "Fill the wrong-sign control histograms (xi1530K_wrongSign/)"}; + o2::framework::Configurable cByPassTOF{"cByPassTOF", false, "Bypass the TOF nSigma selection of the pion and the kaon"}; +}; + +/// Pion PID of mode B (keys of the original three-body pion selection) +struct PionPidCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cPionTPCNSigmaMax{"cPionTPCNSigmaMax", 3.0, "Maximum TPC NSigma for pion"}; + o2::framework::Configurable cPionTOFNSigmaMax{"cPionTOFNSigmaMax", 3.0, "Maximum TOF NSigma for pion"}; + o2::framework::Configurable cPionUsePtDepPID{"cPionUsePtDepPID", false, "Use pT-dependent PID cuts for pion"}; + o2::framework::Configurable> cPionPIDPtBins{"cPionPIDPtBins", {0.0f, 0.5f, 0.8f, 2.0f, 999.0f}, "pT bin edges for pion PID cuts"}; + o2::framework::Configurable> cPionTPCNSigmaCuts{"cPionTPCNSigmaCuts", {3.0f, 3.0f, 2.0f, 2.0f}, "TPC NSigma cuts per pT bin (pion)"}; + o2::framework::Configurable> cPionTOFNSigmaCuts{"cPionTOFNSigmaCuts", {3.0f, 3.0f, 3.0f, 3.0f}, "TOF NSigma cuts per pT bin (pion)"}; + o2::framework::Configurable> cPionTOFRequired{"cPionTOFRequired", {0, 0, 1, 1}, "Require TOF per pT bin (pion)"}; +}; + +/// Kaon PID of mode B (same shape as the pion keys) +struct KaonPidCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cKaonTPCNSigmaMax{"cKaonTPCNSigmaMax", 3.0, "Maximum TPC NSigma for kaon"}; + o2::framework::Configurable cKaonTOFNSigmaMax{"cKaonTOFNSigmaMax", 3.0, "Maximum TOF NSigma for kaon"}; + o2::framework::Configurable cKaonUsePtDepPID{"cKaonUsePtDepPID", false, "Use pT-dependent PID cuts for kaon"}; + o2::framework::Configurable> cKaonPIDPtBins{"cKaonPIDPtBins", {0.0f, 0.5f, 0.8f, 2.0f, 999.0f}, "pT bin edges for kaon PID cuts"}; + o2::framework::Configurable> cKaonTPCNSigmaCuts{"cKaonTPCNSigmaCuts", {3.0f, 3.0f, 2.0f, 2.0f}, "TPC NSigma cuts per pT bin (kaon)"}; + o2::framework::Configurable> cKaonTOFNSigmaCuts{"cKaonTOFNSigmaCuts", {3.0f, 3.0f, 3.0f, 3.0f}, "TOF NSigma cuts per pT bin (kaon)"}; + o2::framework::Configurable> cKaonTOFRequired{"cKaonTOFRequired", {0, 0, 1, 1}, "Require TOF per pT bin (kaon)"}; +}; + +/// Common candidate selection +struct CandidateCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cfgRapidityCut{"cfgRapidityCut", 0.5, "Rapidity cut"}; +}; + +/// Track selection of the mode-B pion and kaon: the common resonance TrackCuts with the JSON prefix "trk." +/// (the Xi selection owns cMinPtcut). The framework prefixes a group only through a `prefix` data member, which +/// TrackCuts does not have and a derived group cannot add (members must belong to one class for the structured +/// binding of the option registration); the names are therefore prefixed here in the same way ("trk.cMinPtcut"). +/// Defaults that differ from TrackCuts follow the original pion selection: |eta| < 0.8 (cPionEtaMax), +/// 70 found TPC clusters for ResoTracks (cPionTPCNClusMin) and DCAz < 0.15 cm (cPionDCAzMax was 0.2 cm, which is +/// beyond the 0.15 cm range of the micro001 DCA grid). +inline TrackCuts makeXi1530KTrackCuts() +{ + TrackCuts cuts; + cuts.cMaxEtacut.value = 0.8f; + cuts.cMaxDCAzToPVcut.value = 0.15; + cuts.cfgTPCcluster.value = 70; + o2::framework::homogeneous_apply_refs( + [](auto& option) { + if constexpr (requires { option.name; }) { + option.name.insert(0, TrackCutsPrefix); + } + return true; + }, + cuts); + return cuts; +} + +/// Process functions enabled in a task; they decide the configuration checks and the registered histograms. +struct ProcessModes { + bool xiK0s = false; // any mode-A process function + bool xi1530K = false; // any mode-B process function + bool microTracks = false; // any mode-B process function reading micro tracks (quantised DCA and nSigma) + bool mcReco = false; // reconstructed MC + bool mcGen = false; // generated Omega(2012) parents (ResoMCParents_001) + bool mixing = false; // event mixing +}; + +/// Loose-stage traversal for the export hook (as the K1 LooseStageOptions). +/// With the defaults and without an export hook, the enumeration applies only the conventional selection. +struct LooseStageOptions { + bool audit = false; // fill ML//looseCutflow and ML//looseMassPtActivity + bool exportSelected = false; // hand candidates to the export hook at the selected stage instead of the loose stage +}; + +/// Mode-A candidate handed to the candidate and export hooks. +/// alpha is computed when cKinCuts is on or the values were requested; passesKinCut is true when cKinCuts is off. +struct XiK0sCandidateValues { + ROOT::Math::PxPyPzEVector omega; // xi + k0s + ROOT::Math::PxPyPzEVector xi; // (pt, eta, phi, mXi) + ROOT::Math::PxPyPzEVector k0s; // (pt, eta, phi, MassK0Short) + float alpha = 0.f; // Xi-K0S opening angle + bool passesKinCut = true; + bool inRapidity = false; // |y| < cfgRapidityCut +}; + +/// Mode-B candidate handed to the candidate and export hooks. +/// massXiK, massPiK, openingAngleXi1530K and cosThetaStar are computed only when the values were requested. +struct Xi1530KCandidateValues { + ROOT::Math::PxPyPzEVector omega; // xi + pion + kaon + ROOT::Math::PxPyPzEVector xi; // (pt, eta, phi, mXi) + ROOT::Math::PxPyPzEVector pion; // (pt, eta, phi, MassPionCharged) + ROOT::Math::PxPyPzEVector kaon; // (pt, eta, phi, MassKaonCharged) + ROOT::Math::PxPyPzEVector xi1530; // xi + pion + float massXiPi = 0.f; + float massXiK = 0.f; + float massPiK = 0.f; + float openingAngleXi1530K = 0.f; + float cosThetaStar = 0.f; + uint8_t chargePattern = o2::analysis::omega2012ml::kSignalPattern; + bool inXi1530Window = false; + bool inRapidity = false; // |y| < cfgRapidityCut +}; + +// Truth classification, separately per mode + +/// Mode A: reconstructed Xi and K0S from the same Omega(2012) (criteria of the original task). +/// A K0S from an intermediate K0bar has motherPDG 311 and is not matched; the truth table keeps motherPDG for audit. +template +XiK0sTruth classifyXiK0sTruth(const Xi& xi, const V0& v0) +{ + if (std::abs(xi.pdgCode()) != kXiMinus || std::abs(v0.pdgCode()) != kK0Short) { + return XiK0sTruth::None; + } + if (xi.motherId() < 0 || xi.motherId() != v0.motherId()) { + return XiK0sTruth::None; + } + if (std::abs(xi.motherPDG()) != o2::constants::physics::Pdg::kOmega2012Minus || xi.motherPDG() != v0.motherPDG()) { + return XiK0sTruth::None; + } + return XiK0sTruth::Matched; +} + +template +bool hasSibling(const Track& track, int particleId) +{ + if (particleId < 0) { + return false; + } + const auto siblings = track.siblingIds(); + return siblings[0] == particleId || siblings[1] == particleId; +} + +/// Mode B: Xi and pion from the same Xi(1530)0, kaon from the Omega(2012) whose other daughter is that Xi(1530)0. +/// Cascades carry no sibling IDs, so the kaon carries the link (siblingIds contains the Xi(1530)0 ID). +template +Xi1530KTruth classifyXi1530KTruth(const Xi& xi, const Track& pion, const Track& kaon) +{ + if (std::abs(xi.pdgCode()) != kXiMinus) { + return Xi1530KTruth::None; + } + const int xiCharge = xi.pdgCode() > 0 ? -1 : 1; // Xi- has PDG code +3312 + if (pion.pdgCode() != -xiCharge * kPiPlus || kaon.pdgCode() != xiCharge * kKPlus) { + return Xi1530KTruth::None; + } + if (xi.motherId() < 0 || xi.motherId() != pion.motherId() || std::abs(xi.motherPDG()) != PdgXi1530Zero || pion.motherPDG() != xi.motherPDG()) { + return Xi1530KTruth::None; + } + // Omega(2012)- and Xi- both have positive PDG codes + if (kaon.motherId() < 0 || kaon.motherPDG() != -xiCharge * o2::constants::physics::Pdg::kOmega2012Minus) { + return Xi1530KTruth::None; + } + return hasSibling(kaon, xi.motherId()) ? Xi1530KTruth::Matched : Xi1530KTruth::None; +} + +/// Immediate channel of a generated Omega(2012) from the PDG codes of its two daughters (charge conjugates included). +inline GeneratedChannel classifyGeneratedOmega2012(int pdg, int daughter1, int daughter2) +{ + const int sign = pdg > 0 ? 1 : -1; + auto isPair = [&](int a, int b) { + if (a == sign * PdgXi1530Zero && b == -sign * kKPlus) { + return GeneratedChannel::Xi1530K; + } + if (a == sign * kXiMinus && (b == kK0Short || b == -sign * kK0)) { + return GeneratedChannel::XiK0s; + } + return GeneratedChannel::Other; + }; + const auto channel = isPair(daughter1, daughter2); + return channel != GeneratedChannel::Other ? channel : isPair(daughter2, daughter1); +} + +/// Omega(2012) selection and candidate enumeration of the Omega(2012) tasks. +/// The task owns the configurable groups and the histogram registry and passes them in init(). +class Omega2012AnalysisCore +{ + public: + void init(o2::framework::HistogramRegistry& histos, + EventCuts const& eventCuts, XiCuts const& xiCuts, K0sCuts const& k0sCuts, + Xi1530KCuts const& xi1530KCuts, TrackCuts const& trackCuts, + PionPidCuts const& pionPidCuts, KaonPidCuts const& kaonPidCuts, + CandidateCuts const& candidateCuts, ProcessModes const& modes, LooseStageOptions const& looseOptions = {}) + { + mXiCuts = xiCuts; + mK0sCuts = k0sCuts; + mXi1530KCuts = xi1530KCuts; + mCandidateCuts = candidateCuts; + mLooseOptions = looseOptions; + + // The order follows Species + std::vector pid(2); + auto& pion = pid[static_cast(Species::Pion)]; + pion.species = "Pion"; + pion.maxTPCnSigma = pionPidCuts.cPionTPCNSigmaMax.value; + pion.maxTOFnSigma = pionPidCuts.cPionTOFNSigmaMax.value; + pion.usePtDependent = pionPidCuts.cPionUsePtDepPID.value; + pion.ptBins = pionPidCuts.cPionPIDPtBins.value; + pion.tpcNSigmaCuts = pionPidCuts.cPionTPCNSigmaCuts.value; + pion.tofNSigmaCuts = pionPidCuts.cPionTOFNSigmaCuts.value; + pion.tofRequired = pionPidCuts.cPionTOFRequired.value; + pion.maxTPCName = "cPionTPCNSigmaMax"; + pion.maxTOFName = "cPionTOFNSigmaMax"; + pion.tpcCutsName = "cPionTPCNSigmaCuts"; + pion.tofCutsName = "cPionTOFNSigmaCuts"; + auto& kaon = pid[static_cast(Species::Kaon)]; + kaon.species = "Kaon"; + kaon.maxTPCnSigma = kaonPidCuts.cKaonTPCNSigmaMax.value; + kaon.maxTOFnSigma = kaonPidCuts.cKaonTOFNSigmaMax.value; + kaon.usePtDependent = kaonPidCuts.cKaonUsePtDepPID.value; + kaon.ptBins = kaonPidCuts.cKaonPIDPtBins.value; + kaon.tpcNSigmaCuts = kaonPidCuts.cKaonTPCNSigmaCuts.value; + kaon.tofNSigmaCuts = kaonPidCuts.cKaonTOFNSigmaCuts.value; + kaon.tofRequired = kaonPidCuts.cKaonTOFRequired.value; + kaon.maxTPCName = "cKaonTPCNSigmaMax"; + kaon.maxTOFName = "cKaonTOFNSigmaMax"; + kaon.tpcCutsName = "cKaonTPCNSigmaCuts"; + kaon.tofCutsName = "cKaonTOFNSigmaCuts"; + mSelection.init(eventCuts, trackCuts, std::move(pid), mXi1530KCuts.cByPassTOF.value, modes.microTracks); + + registerHistograms(histos, modes); + } + + template + bool passesEventCuts(const CollisionType& collision) + { + return mSelection.passesEventCuts(collision); + } + + template + bool passesMCEventCuts(const CollisionType& collision) + { + return mSelection.passesMCEventCuts(collision); + } + + [[nodiscard]] bool kinCutsEnabled() const { return mK0sCuts.cKinCuts.value; } + [[nodiscard]] bool fillWrongSign() const { return mXi1530KCuts.cXi1530KFillWrongSign.value; } + + // Xi selection (common to both modes); equivalent to cascprimaryTrackCut && casctopCut of the original task. + template + [[nodiscard]] int xiSelectionStage(const CascadeType& c) const + { + const auto& cut = mXiCuts; + if (std::abs(c.eta()) > cut.cMaxEtaCut.value || std::abs(c.pt()) < cut.cMinPtcut.value) { + return kXiInput; + } + if (cut.cDCAxyToPVAsPtForCasc.value && std::abs(c.dcaXYCascToPV()) > (cut.cDCAxyToPVByPtCascP0.value + cut.cDCAxyToPVByPtCascExp.value * c.pt())) { + return kXiKinematics; + } + // The DCAz cut uses the DCAxy parameters, as in the original task + if (cut.cDCAzToPVAsPtForCasc.value && std::abs(c.dcaZCascToPV()) > (cut.cDCAxyToPVByPtCascP0.value + cut.cDCAxyToPVByPtCascExp.value * std::pow(c.pt(), CascDCAzPtExponent))) { + return kXiKinematics; + } + + // V0 (Lambda) topology inside the cascade + if (std::abs(c.daughDCA()) > cut.cDCALambdaDaugtherscut.value || std::abs(c.dcav0topv()) < cut.cDCALambdaToPVcut.value) { + return kXiDCA; + } + if (c.sign() < 0) { // Xi- + if (std::abs(c.dcanegtopv()) < cut.cDCAPionToPVcut.value || std::abs(c.dcapostopv()) < cut.cDCAProtonToPVcut.value) { + return kXiDCA; + } + } else { // Anti-Xi + if (std::abs(c.dcanegtopv()) < cut.cDCAProtonToPVcut.value || std::abs(c.dcapostopv()) < cut.cDCAPionToPVcut.value) { + return kXiDCA; + } + } + if (c.v0CosPA() < std::cos(cut.cV0CosPACutPtDepP0.value - cut.cV0CosPACutPtDepP1.value * c.pt())) { + return kXiDCA; + } + if (c.transRadius() > cut.cMaxV0radiuscut.value || c.transRadius() < cut.cMinV0radiuscut.value) { + return kXiDCA; + } + if (std::abs(c.mLambda() - o2::constants::physics::MassLambda) > cut.cMasswindowV0cut.value) { + return kXiDCA; + } + + // Cascade topology + if (std::abs(c.dcabachtopv()) < cut.cDCABachlorToPVcut.value) { + return kXiV0Topology; + } + if (c.pt() < cut.cDCAXiDaugthersCutPtRangeLower.value) { + if (c.cascDaughDCA() > cut.cDCAXiDaugthersCutPtDepLower.value) { + return kXiV0Topology; + } + } else if (c.pt() < cut.cDCAXiDaugthersCutPtRangeUpper.value) { + if (c.cascDaughDCA() > cut.cDCAXiDaugthersCutPtDepMiddle.value) { + return kXiV0Topology; + } + } else { + if (c.cascDaughDCA() > cut.cDCAXiDaugthersCutPtDepUpper.value) { + return kXiV0Topology; + } + } + if (c.cascCosPA() < std::cos(cut.cCosPACascCutPtDepP0.value - cut.cCosPACascCutPtDepP1.value * c.pt())) { + return kXiV0Topology; + } + if (c.cascTransRadius() > cut.cMaxCascradiuscut.value || c.cascTransRadius() < cut.cMinCascradiuscut.value) { + return kXiV0Topology; + } + if (std::abs(c.mXi() - cut.cMassXiminus.value) > cut.cMasswindowCasccut.value) { + return kXiCascTopology; + } + return kXiMass; + } + + // K0S proper lifetime (cm/c), the expression of the original task + template + static double k0sProperLifetime(const CollisionType& collision, const V0Type& v0) + { + float dx = v0.decayVtxX() - collision.posX(); + float dy = v0.decayVtxY() - collision.posY(); + float dz = v0.decayVtxZ() - collision.posZ(); + float l = std::sqrt(dx * dx + dy * dy + dz * dz); + float p = std::sqrt(v0.px() * v0.px() + v0.py() * v0.py() + v0.pz() * v0.pz()); + return (l / (p + SmallNumber)) * o2::constants::physics::MassK0Short; + } + + // K0S selection of mode A; equivalent to v0CutEnhanced of the original task. + // The collision is the one of the V0 (its primary vertex defines the proper lifetime). + template + [[nodiscard]] int k0sSelectionStage(const CollisionType& collision, const V0Type& v0) const + { + const auto& cut = mK0sCuts; + if (std::abs(v0.eta()) > cut.cMaxV0Etacut.value || v0.pt() < mXiCuts.cMinPtcut.value) { + return kK0sInput; + } + if (v0.v0CosPA() < cut.cK0sMinCosPA.value || v0.daughDCA() > cut.cK0sMaxDaughDCA.value) { + return kK0sKinematics; + } + if (std::abs(v0.dcapostopv()) < cut.cK0sDauPosDCAtoPVMin.value || std::abs(v0.dcanegtopv()) < cut.cK0sDauNegDCAtoPVMin.value) { + return kK0sKinematics; + } + auto radius = v0.transRadius(); + if (radius < cut.cK0sRadiusMin.value || radius > cut.cK0sRadiusMax.value) { + return kK0sKinematics; + } + if (std::abs(v0.dcav0topv()) > MaxDCAV0ToPV) { + return kK0sKinematics; + } + if (k0sProperLifetime(collision, v0) > cut.cK0sProperLifetimeMax.value) { + return kK0sTopology; + } + if (v0.qtarm() < cut.cK0sArmenterosQtMin.value) { + return kK0sTopology; + } + if (std::abs(v0.mK0Short() - o2::constants::physics::MassK0Short) > cut.cK0sMassWindow.value) { + return kK0sLifetime; + } + if (cut.cK0sCrossMassRejection.value) { + if (std::abs(v0.mLambda() - o2::constants::physics::MassLambda) < cut.cK0sCrossMassRejectionWindow.value || + std::abs(v0.mAntiLambda() - o2::constants::physics::MassLambda) < cut.cK0sCrossMassRejectionWindow.value) { + return kK0sLifetime; + } + } + if (std::abs(v0.daughterTPCNSigmaPosPi()) >= cut.cK0sDaughterPiTPCNSigmaMax.value || + std::abs(v0.daughterTPCNSigmaNegPi()) >= cut.cK0sDaughterPiTPCNSigmaMax.value) { + return kK0sMass; + } + if (v0.nCrossedRowsPos() <= cut.cK0sPosDaughterMinCrossedRows.value || v0.nCrossedRowsNeg() <= cut.cK0sNegDaughterMinCrossedRows.value) { + return kK0sMass; + } + if (v0.qtarm() < cut.cK0sArmenterosAlphaCoeff.value * std::fabs(v0.alpha())) { + return kK0sDaughters; + } + return kK0sArmenteros; + } + + // Xi-K0S opening angle and the pT-dependent opening-angle window of the original task (kinCuts) + template + [[nodiscard]] bool passesKinematicCut(const FirstVecT& firstDaughter, const SecondVecT& secondDaughter, const MotherVecT& mother, float& alpha) const + { + auto firstP = std::sqrt(firstDaughter.Px() * firstDaughter.Px() + firstDaughter.Py() * firstDaughter.Py() + firstDaughter.Pz() * firstDaughter.Pz()); + auto secondP = std::sqrt(secondDaughter.Px() * secondDaughter.Px() + secondDaughter.Py() * secondDaughter.Py() + secondDaughter.Pz() * secondDaughter.Pz()); + if (firstP < SmallNumber || secondP < SmallNumber) { + alpha = 0.f; + return false; + } + + auto cosAlpha = (firstDaughter.Px() * secondDaughter.Px() + firstDaughter.Py() * secondDaughter.Py() + firstDaughter.Pz() * secondDaughter.Pz()) / (firstP * secondP); + if (cosAlpha > 1.) { + cosAlpha = 1.; + } else if (cosAlpha < -1.) { + cosAlpha = -1.; + } + alpha = std::acos(cosAlpha); + + const auto& kinCutsPt = mK0sCuts.cKinCutsPt.value; + const auto& kinLowerCutsAlpha = mK0sCuts.cKinLowerCutsAlpha.value; + const auto& kinUpperCutsAlpha = mK0sCuts.cKinUpperCutsAlpha.value; + + int kinCutsSize = static_cast(kinUpperCutsAlpha.size()); + if (kinCutsSize > static_cast(kinLowerCutsAlpha.size())) { + kinCutsSize = static_cast(kinLowerCutsAlpha.size()); + } + if (kinCutsSize > static_cast(kinCutsPt.size()) - 1) { + kinCutsSize = static_cast(kinCutsPt.size()) - 1; + } + + for (int i = 0; i < kinCutsSize; ++i) { + if ((mother.Pt() > kinCutsPt[i] && mother.Pt() <= kinCutsPt[i + 1]) && (alpha < kinLowerCutsAlpha[i] || alpha > kinUpperCutsAlpha[i])) { + return false; + } + } + return true; + } + + // Full selection stage of a mode-B track (quality, TOF requirement, PID) + template + int trackSelectionStage(const TrackType& track) + { + return speciesStage(track, mSelection.trackQualityStage(track)); + } + + template + int pionSelectionStage(const TrackType& track) + { + return trackSelectionStage(track); + } + + template + int kaonSelectionStage(const TrackType& track) + { + return trackSelectionStage(track); + } + + // Mode A: (Xi, K0S) enumeration of one collision, or of one mixed pair (Xi from collision, K0S from v0Collision). + // The order of the original task is kept: K0S selection, Xi selection, then Xi (outer) x K0S (inner). + // Hooks (nullptr to skip): + // - onK0s(v0, properLifetime, selected): every V0, in table order, + // - onXi(xi, selected): every cascade, in table order, + // - onCandidate(xi, v0, XiK0sCandidateValues): every pair of selected objects without shared daughters + // (the daughter-ID check is skipped in mixed events); kinematic and rapidity flags are in the values, + // - onExport(collision, xi, v0, values, passBits): same-event candidates at the loose or selected stage. + template + void forEachXiK0sCandidate(o2::framework::HistogramRegistry& histos, const CollisionType& collision, const V0CollisionType& v0Collision, + const CascadesType& cascades, const V0sType& v0s, bool computeValues, + XiHook onXi = nullptr, K0sHook onK0s = nullptr, CandidateHook onCandidate = nullptr, ExportHook onExport = nullptr) + { + constexpr bool HasXiHook = !std::is_same_v; + constexpr bool HasK0sHook = !std::is_same_v; + constexpr bool HasCandidateHook = !std::is_same_v; + constexpr bool HasExportHook = !std::is_same_v; + constexpr bool FillCutFlow = !IsMix; + bool visitLoose = false; + bool checkValidity = false; + if constexpr (!IsMix) { + visitLoose = mLooseOptions.audit || (!mLooseOptions.exportSelected && HasExportHook); + checkValidity = visitLoose || (mLooseOptions.exportSelected && HasExportHook); + } + + // K0S selection cache + using V0Row = std::decay_t; + std::vector> k0sList; + k0sList.reserve(v0s.size()); + for (const auto& v0 : v0s) { + const int stage = k0sSelectionStage(v0Collision, v0); + const bool selected = stage == kK0sSelected; + if constexpr (HasK0sHook) { + onK0s(v0, k0sProperLifetime(v0Collision, v0), selected); + } + if constexpr (FillCutFlow) { + for (int i = 0; i <= stage; ++i) { + histos.fill(HIST("CutFlow/xiK0s/k0s"), i); + } + } + if (selected || visitLoose) { + const auto indices = v0.indices(); + k0sList.push_back({v0, {indices[0], indices[1]}, selected}); + } + } + + // Xi selection cache + using XiRow = std::decay_t; + std::vector> xiList; + xiList.reserve(cascades.size()); + for (const auto& xi : cascades) { + const int stage = xiSelectionStage(xi); + const bool selected = stage == kXiSelected; + if constexpr (HasXiHook) { + onXi(xi, selected); + } + if constexpr (FillCutFlow) { + for (int i = 0; i <= stage; ++i) { + histos.fill(HIST("CutFlow/xiK0s/xi"), i); + } + } + if (selected || visitLoose) { + const auto indices = xi.cascadeIndices(); + xiList.push_back({xi, {indices[0], indices[1], indices[2]}, selected}); + } + } + + const bool kinCutsOn = mK0sCuts.cKinCuts.value; + const float rapidityMax = mCandidateCuts.cfgRapidityCut.value; + for (const auto& cachedXi : xiList) { + const auto& xi = cachedXi.row; + for (const auto& cachedK0s : k0sList) { + const auto& v0 = cachedK0s.row; + const bool objectsSelected = cachedXi.selected && cachedK0s.selected; + bool truthMatched = false; + if constexpr (IsMC && !IsMix) { + truthMatched = classifyXiK0sTruth(xi, v0) == XiK0sTruth::Matched; + } + auto countCandidate = [&](int stage) { + if constexpr (FillCutFlow) { + if (objectsSelected) { + histos.fill(HIST("CutFlow/xiK0s/candidates"), stage, 0); + if (truthMatched) { + histos.fill(HIST("CutFlow/xiK0s/candidates"), stage, 1); + } + } + } + }; + countCandidate(0); + if constexpr (!IsMix) { + if (sharesAnyDaughterId(cachedXi.daughterIds, cachedK0s.daughterIds)) { + continue; + } + } + countCandidate(1); + + // 4-vectors (construction of the original task) + XiK0sCandidateValues values; + values.xi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); + values.k0s = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(v0.pt(), v0.eta(), v0.phi(), o2::constants::physics::MassK0Short)); + values.omega = values.xi + values.k0s; + if (kinCutsOn || computeValues || checkValidity) { + float alpha = 0.f; + const bool kinCutFlag = passesKinematicCut(values.xi, values.k0s, values.omega, alpha); + values.alpha = alpha; + if (kinCutsOn) { + values.passesKinCut = kinCutFlag; + } + } + values.inRapidity = !(std::abs(values.omega.Rapidity()) >= rapidityMax); + if (values.passesKinCut) { + countCandidate(2); + if (values.inRapidity) { + countCandidate(3); + } + } + + if constexpr (!IsMix) { + if (checkValidity) { + exportXiK0s(histos, collision, xi, v0, values, cachedXi.selected, cachedK0s.selected, truthMatched, visitLoose, onExport); + } + } + if constexpr (HasCandidateHook) { + if (objectsSelected) { + onCandidate(xi, v0, values); + } + } + } + } + } + + // Mode B: (Xi, pion, kaon) enumeration of one collision, or of one mixed pair (Xi from collision, tracks from the other). + // trackIds: ResoTrackTracks for full ResoTracks (positional source track IDs), unused for ResoMicroTracks_001 (trackId()). + // Hooks (nullptr to skip): + // - onXi(xi, selected): every cascade, in table order, + // - onTrack(track, pionStage, kaonStage): every track, in table order (TrackStage values), + // - onCandidate(xi, pion, kaon, Xi1530KCandidateValues): every triplet of selected objects with five distinct + // source track IDs (no Xi-daughter check in mixed events) inside the Xi(1530) window; all charge patterns, + // - onExport(collision, xi, pion, kaon, values, passBits): same-event micro candidates at the loose or selected stage. + template + void forEachXi1530KCandidate(o2::framework::HistogramRegistry& histos, const CollisionType& collision, + const CascadesType& cascades, const TracksType& tracks, const TrackIdsType& trackIds, bool computeValues, + XiHook onXi = nullptr, TrackHook onTrack = nullptr, CandidateHook onCandidate = nullptr, ExportHook onExport = nullptr) + { + constexpr bool HasXiHook = !std::is_same_v; + constexpr bool HasTrackHook = !std::is_same_v; + constexpr bool HasCandidateHook = !std::is_same_v; + constexpr bool HasExportHook = !std::is_same_v; + constexpr bool FillCutFlow = !IsMix; + bool visitLoose = false; + bool checkValidity = false; + if constexpr (IsResoMicrotrack && !IsMix) { + visitLoose = mLooseOptions.audit || (!mLooseOptions.exportSelected && HasExportHook); + checkValidity = visitLoose || (mLooseOptions.exportSelected && HasExportHook); + } + + // Xi selection cache + using XiRow = std::decay_t; + std::vector> xiList; + xiList.reserve(cascades.size()); + for (const auto& xi : cascades) { + const int stage = xiSelectionStage(xi); + const bool selected = stage == kXiSelected; + if constexpr (HasXiHook) { + onXi(xi, selected); + } + if constexpr (FillCutFlow) { + for (int i = 0; i <= stage; ++i) { + histos.fill(HIST("CutFlow/xi1530K/xi"), i); + } + } + if (selected || visitLoose) { + const auto indices = xi.cascadeIndices(); + xiList.push_back({xi, {indices[0], indices[1], indices[2]}, selected}); + } + } + + // Track selection cache: every track is selected once per collision, for both species + const std::size_t nTracks = tracks.size(); + if (nTracks == 0) { + return; + } + const int64_t firstIndex = tracks.begin().index(); + std::vector trackCache(nTracks); + for (const auto& track : tracks) { + auto& entry = trackCache[getCacheIndex(track, firstIndex, nTracks)]; + const int qualityStage = mSelection.trackQualityStage(track); + const int pionStage = speciesStage(track, qualityStage); + const int kaonStage = speciesStage(track, qualityStage); + entry.sourceId = sourceTrackId(track, trackIds); + entry.quality = qualityStage == TrackStage::kTrkClusters; + entry.pionSelected = pionStage == TrackStage::kTrkPID; + entry.kaonSelected = kaonStage == TrackStage::kTrkPID; + if constexpr (HasTrackHook) { + onTrack(track, pionStage, kaonStage); + } + if constexpr (FillCutFlow) { + for (int i = 0; i <= pionStage; ++i) { + histos.fill(HIST("CutFlow/xi1530K/tracks"), i, static_cast(Species::Pion)); + } + for (int i = 0; i <= kaonStage; ++i) { + histos.fill(HIST("CutFlow/xi1530K/tracks"), i, static_cast(Species::Kaon)); + } + } + } + + const float rapidityMax = mCandidateCuts.cfgRapidityCut.value; + const bool useWindow = mXi1530KCuts.cXi1530UseMassWindow.value; + const float windowCenter = mXi1530KCuts.cXi1530Mass.value; + const float windowWidth = mXi1530KCuts.cXi1530MassWindow.value; + for (const auto& cachedXi : xiList) { + const auto& xi = cachedXi.row; + const int xiSign = xi.sign() < 0 ? -1 : 1; + const ROOT::Math::PxPyPzEVector pXi(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); + for (const auto& pion : tracks) { + const auto& pionEntry = trackCache[getCacheIndex(pion, firstIndex, nTracks)]; + if (!pionEntry.pionSelected && !visitLoose) { + continue; + } + if constexpr (!IsMix) { + if (sharesDaughterId(cachedXi.daughterIds, pionEntry.sourceId)) { + continue; + } + } + const bool pairSelected = cachedXi.selected && pionEntry.pionSelected; + const bool pionSignal = pion.sign() != xiSign; + const ROOT::Math::PxPyPzEVector pPion(ROOT::Math::PtEtaPhiMVector(pion.pt(), pion.eta(), pion.phi(), o2::constants::physics::MassPionCharged)); + const ROOT::Math::PxPyPzEVector pXi1530 = pXi + pPion; + const float massXiPi = pXi1530.M(); + const bool inWindow = !useWindow || std::abs(massXiPi - windowCenter) < windowWidth; + if constexpr (FillCutFlow) { + if (pairSelected) { + countXi1530K(histos, 0, pionSignal, false); + if (inWindow) { + countXi1530K(histos, 1, pionSignal, false); + } + } + } + if (!(pairSelected && inWindow) && !visitLoose) { + continue; + } + for (const auto& kaon : tracks) { + if (kaon.index() == pion.index()) { + continue; + } + const auto& kaonEntry = trackCache[getCacheIndex(kaon, firstIndex, nTracks)]; + if (!kaonEntry.kaonSelected && !visitLoose) { + continue; + } + if (kaonEntry.sourceId == pionEntry.sourceId) { + continue; + } + if constexpr (!IsMix) { + if (sharesDaughterId(cachedXi.daughterIds, kaonEntry.sourceId)) { + continue; + } + } + const bool tripletSelected = pairSelected && inWindow && kaonEntry.kaonSelected; + Xi1530KCandidateValues values; + values.xi = pXi; + values.pion = pPion; + values.kaon = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(kaon.pt(), kaon.eta(), kaon.phi(), o2::constants::physics::MassKaonCharged)); + values.xi1530 = pXi1530; + values.omega = pXi1530 + values.kaon; + values.massXiPi = massXiPi; + values.chargePattern = o2::analysis::omega2012ml::chargePattern(xiSign, pion.sign(), kaon.sign()); + values.inXi1530Window = inWindow; + values.inRapidity = !(std::abs(values.omega.Rapidity()) >= rapidityMax); + if (computeValues || checkValidity) { + values.massXiK = (pXi + values.kaon).M(); + values.massPiK = (pPion + values.kaon).M(); + values.openingAngleXi1530K = ROOT::Math::VectorUtil::Angle(pXi1530, values.kaon); + values.cosThetaStar = static_cast(o2::analysis::omega2012ml::detail::cosThetaStar( + o2::analysis::omega2012ml::detail::LorentzVector(pXi1530), o2::analysis::omega2012ml::detail::LorentzVector(values.omega))); + } + bool truthMatched = false; + if constexpr (IsMC && !IsMix) { + truthMatched = classifyXi1530KTruth(xi, pion, kaon) == Xi1530KTruth::Matched; + } + if constexpr (FillCutFlow) { + if (tripletSelected) { + const bool signal = values.chargePattern == o2::analysis::omega2012ml::kSignalPattern; + countXi1530K(histos, 2, signal, truthMatched); + if (values.inRapidity) { + countXi1530K(histos, 3, signal, truthMatched); + } + } + } + if constexpr (IsResoMicrotrack && !IsMix) { + if (checkValidity) { + exportXi1530K(histos, collision, xi, pion, kaon, values, cachedXi.selected, + pionEntry.quality && kaonEntry.quality, pionEntry.pionSelected && kaonEntry.kaonSelected, + truthMatched, visitLoose, onExport); + } + } + if constexpr (HasCandidateHook) { + if (tripletSelected) { + onCandidate(xi, pion, kaon, values); + } + } + } + } + } + } + + // Generated Omega(2012) parents of a selected reconstructed MC collision (ResoMCParents_001). The optional callback + // receives (parent, immediate channel) for the parents inside the rapidity window. Split reconstructed collisions + // repeat parent sets: this is not an unconditional generated denominator. + template + void forEachGeneratedOmega2012(o2::framework::HistogramRegistry& histos, const ParentsType& resoParents, Callback callback = nullptr) + { + for (const auto& part : resoParents) { + if (std::abs(part.pdgCode()) != o2::constants::physics::Pdg::kOmega2012Minus) { + continue; + } + const GeneratedChannel channel = classifyGeneratedOmega2012(part.pdgCode(), part.daughterPDG1(), part.daughterPDG2()); + histos.fill(HIST("CutFlow/generated"), 0, static_cast(channel)); + if (!(std::abs(part.y()) < mCandidateCuts.cfgRapidityCut.value)) { + continue; + } + histos.fill(HIST("CutFlow/generated"), 1, static_cast(channel)); + if constexpr (!std::is_same_v) { + callback(part, channel); + } + } + } + + template + static bool sharesAnyDaughterId(const std::array& first, const std::array& second) + { + for (const auto& firstId : first) { + for (const auto& secondId : second) { + if (firstId == secondId) { + return true; + } + } + } + return false; + } + + template + static bool sharesDaughterId(const std::array& daughters, int64_t trackId) + { + for (const auto& daughterId : daughters) { + if (static_cast(daughterId) == trackId) { + return true; + } + } + return false; + } + + // Source track ID: trackId() of ResoMicroTracks_001, or the positional ResoTrackTracks row of a ResoTracks row. + template + static int64_t sourceTrackId(const TrackType& track, const TrackIdsType& trackIds) + { + if constexpr (requires { track.trackId(); }) { + return static_cast(track.trackId()); + } else { + const auto rowIndex = track.globalIndex(); + if (rowIndex >= 0 && rowIndex < static_cast(trackIds.size())) { + return static_cast(trackIds.rawIteratorAt(rowIndex).trackId()); + } + return static_cast(rowIndex); + } + } + + private: + template + struct CachedObject { + Row row; + std::array daughterIds{}; + bool selected = false; + }; + + struct TrackCacheEntry { + int64_t sourceId = -1; + bool quality = false; + bool pionSelected = false; + bool kaonSelected = false; + }; + + template + int speciesStage(const TrackType& track, int qualityStage) + { + if (qualityStage < TrackStage::kTrkClusters) { + return qualityStage; + } + constexpr int SpeciesIndex = static_cast(S); + if (!mSelection.passesTOFRequired(SpeciesIndex, track)) { + return TrackStage::kTrkClusters; + } + const bool hasTOF = track.hasTOF(); + const double tpcNSigma = (S == Species::Pion) ? track.tpcNSigmaPi() : track.tpcNSigmaKa(); + double tofNSigma = std::numeric_limits::quiet_NaN(); // TOF value is only valid with hasTOF + if (hasTOF) { + tofNSigma = (S == Species::Pion) ? track.tofNSigmaPi() : track.tofNSigmaKa(); + } + if (!mSelection.passesPID(SpeciesIndex, track.pt(), hasTOF, tpcNSigma, tofNSigma)) { + return TrackStage::kTrkTOFRequired; + } + return TrackStage::kTrkPID; + } + + template + static std::size_t getCacheIndex(const TrackType& track, int64_t firstIndex, std::size_t size) + { + const int64_t index = static_cast(track.index()) - firstIndex; + if (index < 0 || index >= static_cast(size)) { + LOG(fatal) << "Track index " << track.index() << " is outside the selection cache [" << firstIndex << ", " << firstIndex + static_cast(size) << ")"; + } + return static_cast(index); + } + + static void countXi1530K(o2::framework::HistogramRegistry& histos, int stage, bool signalPattern, bool truthMatched) + { + histos.fill(HIST("CutFlow/xi1530K/candidates"), stage, 0); + if (signalPattern) { + histos.fill(HIST("CutFlow/xi1530K/candidates"), stage, 1); + } + if (truthMatched) { + histos.fill(HIST("CutFlow/xi1530K/candidates"), stage, 2); + } + } + + // Loose / selected stage of a mode-A same-event candidate + template + void exportXiK0s(o2::framework::HistogramRegistry& histos, const CollisionType& collision, const XiType& xi, const V0Type& v0, + XiK0sCandidateValues const& values, bool xiSelected, bool k0sSelected, bool truthMatched, + bool visitLoose, ExportHook onExport) + { + constexpr bool HasExportHook = !std::is_same_v; + namespace ml = o2::analysis::omega2012ml; + const auto canonical = ml::canonicalizeXiK0s(ml::makeCascadeSnapshot(collision, xi), ml::makeV0Snapshot(collision, v0)); + const bool valid = canonical.status == ml::BuildStatus::Ok && + std::isfinite(values.omega.M()) && std::isfinite(values.omega.Pt()) && std::isfinite(values.omega.Rapidity()) && + ml::buildXiK0sFeatures(canonical.candidate).status == ml::BuildStatus::Ok; + if (!valid) { + return; + } + auto countMl = [&](int stage) { + if (mLooseOptions.audit) { + histos.fill(HIST("ML/xiK0s/looseCutflow"), stage, 0); + if (truthMatched) { + histos.fill(HIST("ML/xiK0s/looseCutflow"), stage, 1); + } + } + }; + countMl(0); + if (!values.inRapidity) { + return; + } + const bool selected = xiSelected && k0sSelected && values.passesKinCut; + if (visitLoose) { + countMl(1); + if (mLooseOptions.audit) { + histos.fill(HIST("ML/xiK0s/looseMassPtActivity"), values.omega.M(), values.omega.Pt(), collision.cent()); + } + uint16_t passBits = kXiK0sPassLoose; + if (xiSelected) { + passBits |= kXiK0sPassXi; + countMl(2); + if (k0sSelected) { + passBits |= kXiK0sPassK0s; + countMl(3); + if (values.passesKinCut) { + passBits |= kXiK0sPassKinematic; + countMl(4); + } + } + } + if constexpr (HasExportHook) { + if (!mLooseOptions.exportSelected) { + onExport(collision, xi, v0, values, passBits); + } + } + } + if constexpr (HasExportHook) { + if (mLooseOptions.exportSelected && selected) { + onExport(collision, xi, v0, values, XiK0sPassBitsSelected); + } + } + } + + // Loose / selected stage of a mode-B same-event micro candidate + template + void exportXi1530K(o2::framework::HistogramRegistry& histos, const CollisionType& collision, const XiType& xi, + const TrackType& pion, const TrackType& kaon, Xi1530KCandidateValues const& values, + bool xiSelected, bool qualityBoth, bool pidBoth, bool truthMatched, bool visitLoose, ExportHook onExport) + { + constexpr bool HasExportHook = !std::is_same_v; + namespace ml = o2::analysis::omega2012ml; + const auto canonical = ml::canonicalizeXi1530K(ml::makeCascadeSnapshot(collision, xi), ml::makeTrackSnapshot(pion), ml::makeTrackSnapshot(kaon)); + const bool valid = canonical.status == ml::BuildStatus::Ok && + std::isfinite(values.omega.M()) && std::isfinite(values.omega.Pt()) && std::isfinite(values.omega.Rapidity()) && + ml::buildXi1530KFeatures(canonical.candidate).status == ml::BuildStatus::Ok; + if (!valid) { + return; + } + const int category = values.chargePattern == ml::kSignalPattern ? 0 : 1; + auto countMl = [&](int stage) { + if (mLooseOptions.audit) { + histos.fill(HIST("ML/xi1530K/looseCutflow"), stage, category); + if (truthMatched) { + histos.fill(HIST("ML/xi1530K/looseCutflow"), stage, 2); + } + } + }; + countMl(0); + if (!values.inRapidity) { + return; + } + const bool selected = xiSelected && qualityBoth && pidBoth && values.inXi1530Window; + if (visitLoose) { + countMl(1); + if (mLooseOptions.audit && category == 0) { + histos.fill(HIST("ML/xi1530K/looseMassPtActivity"), values.omega.M(), values.omega.Pt(), collision.cent()); + } + uint16_t passBits = kXi1530KPassLoose; + if (xiSelected) { + passBits |= kXi1530KPassXi; + countMl(2); + if (qualityBoth) { + passBits |= kXi1530KPassQuality; + countMl(3); + if (pidBoth) { + passBits |= kXi1530KPassPID; + countMl(4); + if (values.inXi1530Window) { + passBits |= kXi1530KPassWindow; + countMl(5); + } + } + } + } + if constexpr (HasExportHook) { + if (!mLooseOptions.exportSelected) { + onExport(collision, xi, pion, kaon, values, passBits); + } + } + } + if constexpr (HasExportHook) { + if (mLooseOptions.exportSelected && selected) { + onExport(collision, xi, pion, kaon, values, Xi1530KPassBitsSelected); + } + } + } + + // Cut-flow instrumentation of the selection, per decay mode; the output histograms belong to the tasks. + void registerHistograms(o2::framework::HistogramRegistry& histos, ProcessModes const& modes) + { + using o2::framework::HistType; + const std::array xiLabels{"input", "|#eta|, p_{T}", "DCA to PV", "#Lambda topology", "cascade topology", "#Xi mass"}; + if (modes.xiK0s) { + auto xiFlow = histos.add("CutFlow/xiK0s/xi", "XiK0s mode: cascades, once per selected collision;stage;cascades", HistType::kTH1D, {{kXiNStages, -0.5, static_cast(kXiNStages) - 0.5}}); + for (std::size_t i = 0; i < xiLabels.size(); ++i) { + xiFlow->GetXaxis()->SetBinLabel(i + 1, xiLabels[i]); + } + auto k0sFlow = histos.add("CutFlow/xiK0s/k0s", "XiK0s mode: V0s, once per selected collision;stage;V0s", HistType::kTH1D, {{kK0sNStages, -0.5, static_cast(kK0sNStages) - 0.5}}); + const std::array k0sLabels{"input", "|#eta|, p_{T}", "topology", "lifetime, min q_{T}", "mass, #Lambda rejection", "daughters", "Armenteros"}; + for (std::size_t i = 0; i < k0sLabels.size(); ++i) { + k0sFlow->GetXaxis()->SetBinLabel(i + 1, k0sLabels[i]); + } + auto candidateFlow = histos.add("CutFlow/xiK0s/candidates", "XiK0s mode: selected Xi x selected K0S, same event;stage;category", HistType::kTH2D, + {{NXiK0sCandidateStages, -0.5, NXiK0sCandidateStages - 0.5}, {2, -0.5, 1.5}}); + const std::array candidateLabels{"selected pairs", "distinct daughters", "kinematic cut", "rapidity"}; + for (std::size_t i = 0; i < candidateLabels.size(); ++i) { + candidateFlow->GetXaxis()->SetBinLabel(i + 1, candidateLabels[i]); + } + candidateFlow->GetYaxis()->SetBinLabel(1, "all"); + candidateFlow->GetYaxis()->SetBinLabel(2, "XiK0s truth"); + if (mLooseOptions.audit) { + auto flow = histos.add("ML/xiK0s/looseCutflow", "XiK0s mode: valid canonical candidates;stage;category", HistType::kTH2D, {{5, -0.5, 4.5}, {2, -0.5, 1.5}}); + const std::array labels{"valid canonical", "loose acceptance", "Xi selection", "K0S selection", "kinematic cut"}; + for (std::size_t i = 0; i < labels.size(); ++i) { + flow->GetXaxis()->SetBinLabel(i + 1, labels[i]); + } + flow->GetYaxis()->SetBinLabel(1, "all"); + flow->GetYaxis()->SetBinLabel(2, "XiK0s truth"); + histos.add("ML/xiK0s/looseMassPtActivity", "XiK0s mode loose candidates;mass (GeV/c^{2});pT (GeV/c);centrality", HistType::kTH3D, + {{300, 1.6, 2.8}, {{0., 0.5, 1., 2., 3., 5., 8., 15., 30., 100.}, "pT"}, {{0., 10., 30., 50., 70., 100., 110.}, "centrality"}}); + } + } + if (modes.xi1530K) { + auto xiFlow = histos.add("CutFlow/xi1530K/xi", "Xi1530K mode: cascades, once per selected collision;stage;cascades", HistType::kTH1D, {{kXiNStages, -0.5, static_cast(kXiNStages) - 0.5}}); + for (std::size_t i = 0; i < xiLabels.size(); ++i) { + xiFlow->GetXaxis()->SetBinLabel(i + 1, xiLabels[i]); + } + constexpr int NTrackStages = TrackStage::kTrkNStages; + auto trackFlow = histos.add("CutFlow/xi1530K/tracks", "Xi1530K mode: tracks, once per selected collision;stage;species", HistType::kTH2D, + {{NTrackStages, -0.5, NTrackStages - 0.5}, {2, -0.5, 1.5}}); + const std::array trackLabels{"input", "pT", "eta", "DCAxy", "DCAz", "track flags", "clusters / crossed rows", "TOF required", "PID"}; + for (std::size_t i = 0; i < trackLabels.size(); ++i) { + trackFlow->GetXaxis()->SetBinLabel(i + 1, trackLabels[i]); + } + trackFlow->GetYaxis()->SetBinLabel(1, "pion"); + trackFlow->GetYaxis()->SetBinLabel(2, "kaon"); + auto candidateFlow = histos.add("CutFlow/xi1530K/candidates", "Xi1530K mode: same event;stage;category", HistType::kTH2D, + {{NXi1530KCandidateStages, -0.5, NXi1530KCandidateStages - 0.5}, {3, -0.5, 2.5}}); + const std::array candidateLabels{"selected Xi #pi pairs", "Xi(1530) window", "selected Xi #pi K triplets", "rapidity"}; + for (std::size_t i = 0; i < candidateLabels.size(); ++i) { + candidateFlow->GetXaxis()->SetBinLabel(i + 1, candidateLabels[i]); + } + candidateFlow->GetYaxis()->SetBinLabel(1, "all"); + candidateFlow->GetYaxis()->SetBinLabel(2, "signal charge pattern"); + candidateFlow->GetYaxis()->SetBinLabel(3, "Xi1530K truth"); + if (mLooseOptions.audit) { + auto flow = histos.add("ML/xi1530K/looseCutflow", "Xi1530K mode: valid canonical candidates;stage;category", HistType::kTH2D, {{6, -0.5, 5.5}, {3, -0.5, 2.5}}); + const std::array labels{"valid canonical", "loose acceptance", "Xi selection", "track quality", "TOF + PID", "Xi(1530) window"}; + for (std::size_t i = 0; i < labels.size(); ++i) { + flow->GetXaxis()->SetBinLabel(i + 1, labels[i]); + } + flow->GetYaxis()->SetBinLabel(1, "signal charge pattern"); + flow->GetYaxis()->SetBinLabel(2, "wrong-sign control"); + flow->GetYaxis()->SetBinLabel(3, "Xi1530K truth"); + histos.add("ML/xi1530K/looseMassPtActivity", "Xi1530K mode loose signal-pattern candidates;mass (GeV/c^{2});pT (GeV/c);centrality", HistType::kTH3D, + {{300, 1.8, 3.0}, {{0., 0.5, 1., 2., 3., 5., 8., 15., 30., 100.}, "pT"}, {{0., 10., 30., 50., 70., 100., 110.}, "centrality"}}); + } + } + if (modes.mcGen) { + auto generated = histos.add("CutFlow/generated", "Generated Omega(2012) parent rows of selected reconstructed events;stage;immediate channel", HistType::kTH2D, + {{2, -0.5, 1.5}, {NGeneratedChannels, -0.5, NGeneratedChannels - 0.5}}); + generated->GetXaxis()->SetBinLabel(1, "all parent rows"); + generated->GetXaxis()->SetBinLabel(2, "rapidity window"); + generated->GetYaxis()->SetBinLabel(1, "other / unresolved"); + generated->GetYaxis()->SetBinLabel(2, "Xi K0S"); + generated->GetYaxis()->SetBinLabel(3, "Xi(1530) K"); + } + } + + ResoAnalysisSelectionCore mSelection; + XiCuts mXiCuts; + K0sCuts mK0sCuts; + Xi1530KCuts mXi1530KCuts; + CandidateCuts mCandidateCuts; + LooseStageOptions mLooseOptions; +}; + +} // namespace o2::analysis::omega2012 + +#endif // PWGLF_CORE_OMEGA2012ANALYSISCORE_H_ diff --git a/PWGLF/Core/Omega2012MlFeatures.h b/PWGLF/Core/Omega2012MlFeatures.h new file mode 100644 index 00000000000..345a6c70c38 --- /dev/null +++ b/PWGLF/Core/Omega2012MlFeatures.h @@ -0,0 +1,1387 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file Omega2012MlFeatures.h +/// \brief Snapshots, canonical candidates and the two separate feature contracts of the Omega(2012) decay modes +/// \author Bong-Hwi Lim +/// +/// Mode A (XiK0s): Omega(2012)- -> Xi- K0S. Mode B (Xi1530K): Omega(2012)- -> Xi(1530)0 K- -> Xi- pi+ K-. +/// The two modes have independent, separately versioned contracts (feature_contract_omega2012_v1.json); +/// no feature vector, status or projection is shared between them. + +#ifndef PWGLF_CORE_OMEGA2012MLFEATURES_H_ +#define PWGLF_CORE_OMEGA2012MLFEATURES_H_ + +#include "PWGLF/DataModel/LFResonanceTables.h" + +#include +#include + +#include +#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::omega2012ml +{ +inline constexpr std::size_t NItsLayers = 7; +inline constexpr std::size_t NCascadeDaughters = 3; +inline constexpr std::size_t NV0Daughters = 2; +inline constexpr std::size_t NPidSpecies = 3; +inline constexpr std::size_t XiK0sNMasterFeatures = 63; +inline constexpr std::size_t Xi1530KNMasterFeatures = 126; +inline constexpr int8_t SaturatedLow = std::numeric_limits::lowest(); // x10 nSigma code of -inf / missing TOF +inline constexpr int8_t SaturatedHigh = std::numeric_limits::max(); // x10 nSigma code of NaN / +inf +inline constexpr float NSigmaScale = 10.f; +inline constexpr float MicroPidOverflow = 3.5f; + +enum class Profile : uint8_t { DetectorV1, + RelationalV1, + StudyV1 }; +enum class BuildStatus : uint8_t { + Ok, + InvalidChargePattern, + ReusedDaughter, + InvalidMomentum, + InvalidKinematics, + InvalidContract +}; + +/// Charge pattern of a mode-B candidate relative to the Xi charge (metadata, never a feature). +enum ChargePattern : uint8_t { + kSignalPattern = 0, // pion opposite to the Xi charge, kaon equal to it (Xi- pi+ K- and conjugate) + kWrongSignPion = 1, // pion with the Xi charge + kWrongSignKaon = 2, // kaon opposite to the Xi charge + kWrongSignBoth = 3 // both wrong +}; + +/// All ResoCascades scalars of one cascade plus derived decay lengths. +/// Daughter arrays: index 0 positive, 1 negative, 2 bachelor; nSigma: [daughter * 3 + species], species pi, ka, pr. +struct CascadeSnapshot { + int64_t sourceRow = -1; // row of the cascade in the input ResoCascades table + int8_t sign = 0; + float px = 0.f, py = 0.f, pz = 0.f; + float mXi = 0.f, mLambda = 0.f; + float v0CosPA = 0.f, cascCosPA = 0.f, v0DaughDCA = 0.f, cascDaughDCA = 0.f; + float dcaPosToPV = 0.f, dcaNegToPV = 0.f, dcaBachToPV = 0.f, dcaV0ToPV = 0.f, dcaXYToPV = 0.f, dcaZToPV = 0.f; + float v0Radius = 0.f, cascRadius = 0.f; + float decayVtxX = 0.f, decayVtxY = 0.f, decayVtxZ = 0.f; + std::array tpcNSigma10{}; + std::array tofNSigma10{}; + std::array crossedRows{}; + std::array daughterIds{}; + float decayLength = 0.f; // 3D distance from the primary vertex to the decay vertex + float properLength = 0.f; // decayLength * MassXiMinus / p +}; + +/// All ResoV0s scalars of one V0 plus derived decay lengths. Daughter arrays: index 0 positive, 1 negative. +struct V0Snapshot { + int64_t sourceRow = -1; // row of the V0 in the input ResoV0s table + float px = 0.f, py = 0.f, pz = 0.f; + float mK0Short = 0.f, mLambda = 0.f, mAntiLambda = 0.f; + float cosPA = 0.f, daughDCA = 0.f, dcaPosToPV = 0.f, dcaNegToPV = 0.f, dcaV0ToPV = 0.f, radius = 0.f; + float decayVtxX = 0.f, decayVtxY = 0.f, decayVtxZ = 0.f; + float alpha = 0.f, qtarm = 0.f; + std::array tpcNSigma10{}; + std::array tofNSigma10{}; + std::array crossedRows{}; + std::array daughterIds{}; + float decayLength = 0.f; // 3D distance from the primary vertex to the decay vertex + float properLength = 0.f; // decayLength * MassK0Short / p (proper lifetime c*tau) +}; + +/// ResoMicroTracks_001 row: raw packed bytes and the decoded public accessors. +/// The encoding rules are those of the K1 micro001 contract; the type is separate so the K1 contract stays immutable. +struct TrackSnapshot { + int64_t sourceTrackId = -1; + int8_t charge = 0; + float px = 0.f, py = 0.f, pz = 0.f; + uint8_t pidNSigmaPiFlag = 0, pidNSigmaKaFlag = 0, pidNSigmaPrFlag = 0; + uint8_t trackSelectionFlags = 0, trackFlags = 0, tpcNClsCrossedRows = 0, itsClusterMap = 0; + std::array tpcNSigma{}; // pi, ka, pr + std::array tofNSigma{}; // pi, ka, pr + float dcaXY = std::numeric_limits::quiet_NaN(); + float dcaZ = std::numeric_limits::quiet_NaN(); + bool passedPtDependentDCAxy = false, passedPtDependentDCAz = false; + bool hasTOF = false, isPVContributor = false; + std::array itsHit{}; +}; + +struct XiK0sCandidateSnapshot { + CascadeSnapshot xi{}; + V0Snapshot k0s{}; +}; + +struct Xi1530KCandidateSnapshot { + CascadeSnapshot xi{}; + TrackSnapshot pion{}; + TrackSnapshot kaon{}; + uint8_t chargePattern = kSignalPattern; // metadata +}; + +template +struct CanonicalizationResult { + Snapshot candidate{}; + BuildStatus status = BuildStatus::InvalidChargePattern; + explicit operator bool() const { return status == BuildStatus::Ok; } +}; + +/// Candidate kinematics from the stored 3-momenta and PDG masses (audit only, never features) +struct KinematicAudit { + float mass = 0.f; + float pt = 0.f; + float eta = 0.f; + float phi = 0.f; + float scalarSumPt = 0.f; +}; + +template +struct FeaturePack { + std::array master{}; + BuildStatus status = BuildStatus::InvalidContract; + KinematicAudit kinematics{}; +}; +using XiK0sFeaturePack = FeaturePack; +using Xi1530KFeaturePack = FeaturePack; + +namespace detail +{ +template +float decayLength(Collision const& collision, Row const& row) +{ + const float dx = row.decayVtxX() - collision.posX(); + const float dy = row.decayVtxY() - collision.posY(); + const float dz = row.decayVtxZ() - collision.posZ(); + return std::sqrt(dx * dx + dy * dy + dz * dz); +} + +inline float properLength(float length, float px, float py, float pz, double mass) +{ + const float p = std::sqrt(px * px + py * py + pz * pz); + if (!(p > 0.f)) { + return std::numeric_limits::quiet_NaN(); + } + return static_cast(length * mass / p); +} +} // namespace detail + +template +CascadeSnapshot makeCascadeSnapshot(Collision const& collision, Cascade const& row) +{ + CascadeSnapshot out; + out.sourceRow = static_cast(row.globalIndex()); + out.sign = static_cast(row.sign()); + out.px = row.px(); + out.py = row.py(); + out.pz = row.pz(); + out.mXi = row.mXi(); + out.mLambda = row.mLambda(); + out.v0CosPA = row.v0CosPA(); + out.cascCosPA = row.cascCosPA(); + out.v0DaughDCA = row.daughDCA(); + out.cascDaughDCA = row.cascDaughDCA(); + out.dcaPosToPV = row.dcapostopv(); + out.dcaNegToPV = row.dcanegtopv(); + out.dcaBachToPV = row.dcabachtopv(); + out.dcaV0ToPV = row.dcav0topv(); + out.dcaXYToPV = row.dcaXYCascToPV(); + out.dcaZToPV = row.dcaZCascToPV(); + out.v0Radius = row.transRadius(); + out.cascRadius = row.cascTransRadius(); + out.decayVtxX = row.decayVtxX(); + out.decayVtxY = row.decayVtxY(); + out.decayVtxZ = row.decayVtxZ(); + out.tpcNSigma10 = {row.daughterTPCNSigmaPosPi10(), row.daughterTPCNSigmaPosKa10(), row.daughterTPCNSigmaPosPr10(), + row.daughterTPCNSigmaNegPi10(), row.daughterTPCNSigmaNegKa10(), row.daughterTPCNSigmaNegPr10(), + row.daughterTPCNSigmaBachPi10(), row.daughterTPCNSigmaBachKa10(), row.daughterTPCNSigmaBachPr10()}; + out.tofNSigma10 = {row.daughterTOFNSigmaPosPi10(), row.daughterTOFNSigmaPosKa10(), row.daughterTOFNSigmaPosPr10(), + row.daughterTOFNSigmaNegPi10(), row.daughterTOFNSigmaNegKa10(), row.daughterTOFNSigmaNegPr10(), + row.daughterTOFNSigmaBachPi10(), row.daughterTOFNSigmaBachKa10(), row.daughterTOFNSigmaBachPr10()}; + out.crossedRows = {row.nCrossedRowsPos(), row.nCrossedRowsNeg(), row.nCrossedRowsBach()}; + const auto indices = row.cascadeIndices(); + out.daughterIds = {indices[0], indices[1], indices[2]}; + out.decayLength = detail::decayLength(collision, row); + out.properLength = detail::properLength(out.decayLength, out.px, out.py, out.pz, o2::constants::physics::MassXiMinus); + return out; +} + +template +V0Snapshot makeV0Snapshot(Collision const& collision, V0 const& row) +{ + V0Snapshot out; + out.sourceRow = static_cast(row.globalIndex()); + out.px = row.px(); + out.py = row.py(); + out.pz = row.pz(); + out.mK0Short = row.mK0Short(); + out.mLambda = row.mLambda(); + out.mAntiLambda = row.mAntiLambda(); + out.cosPA = row.v0CosPA(); + out.daughDCA = row.daughDCA(); + out.dcaPosToPV = row.dcapostopv(); + out.dcaNegToPV = row.dcanegtopv(); + out.dcaV0ToPV = row.dcav0topv(); + out.radius = row.transRadius(); + out.decayVtxX = row.decayVtxX(); + out.decayVtxY = row.decayVtxY(); + out.decayVtxZ = row.decayVtxZ(); + out.alpha = row.alpha(); + out.qtarm = row.qtarm(); + out.tpcNSigma10 = {row.daughterTPCNSigmaPosPi10(), row.daughterTPCNSigmaPosKa10(), row.daughterTPCNSigmaPosPr10(), + row.daughterTPCNSigmaNegPi10(), row.daughterTPCNSigmaNegKa10(), row.daughterTPCNSigmaNegPr10()}; + out.tofNSigma10 = {row.daughterTOFNSigmaPosPi10(), row.daughterTOFNSigmaPosKa10(), row.daughterTOFNSigmaPosPr10(), + row.daughterTOFNSigmaNegPi10(), row.daughterTOFNSigmaNegKa10(), row.daughterTOFNSigmaNegPr10()}; + out.crossedRows = {row.nCrossedRowsPos(), row.nCrossedRowsNeg()}; + const auto indices = row.indices(); + out.daughterIds = {indices[0], indices[1]}; + out.decayLength = detail::decayLength(collision, row); + out.properLength = detail::properLength(out.decayLength, out.px, out.py, out.pz, o2::constants::physics::MassK0Short); + return out; +} + +template +TrackSnapshot makeTrackSnapshot(MicroTrack const& row) +{ + TrackSnapshot out; + out.sourceTrackId = static_cast(row.trackId()); + out.charge = static_cast(row.sign()); + out.px = static_cast(row.px()); + out.py = static_cast(row.py()); + out.pz = static_cast(row.pz()); + out.pidNSigmaPiFlag = static_cast(row.pidNSigmaPiFlag()); + out.pidNSigmaKaFlag = static_cast(row.pidNSigmaKaFlag()); + out.pidNSigmaPrFlag = static_cast(row.pidNSigmaPrFlag()); + out.trackSelectionFlags = static_cast(row.trackSelectionFlags()); + out.trackFlags = static_cast(row.trackFlags()); + out.tpcNClsCrossedRows = static_cast(row.tpcNClsCrossedRows()); + out.itsClusterMap = static_cast(row.itsClusterMap()); + out.tpcNSigma = {static_cast(row.tpcNSigmaPi()), static_cast(row.tpcNSigmaKa()), static_cast(row.tpcNSigmaPr())}; + out.tofNSigma = {static_cast(row.tofNSigmaPi()), static_cast(row.tofNSigmaKa()), static_cast(row.tofNSigmaPr())}; + out.dcaXY = static_cast(row.dcaXY()); + out.dcaZ = static_cast(row.dcaZ()); + out.passedPtDependentDCAxy = row.passedPtDependentDCAxy(); + out.passedPtDependentDCAz = row.passedPtDependentDCAz(); + out.hasTOF = row.hasTOF(); + out.isPVContributor = row.isPVContributor(); + for (std::size_t layer = 0; layer < NItsLayers; ++layer) { + out.itsHit[layer] = row.hasITSHitInLayer(static_cast(layer)); + } + return out; +} + +inline bool hasFiniteMomentum(float px, float py, float pz) +{ + return std::isfinite(px) && std::isfinite(py) && std::isfinite(pz) && std::hypot(px, py) > 0.f; +} + +inline bool sharesDaughter(CascadeSnapshot const& xi, int64_t trackId) +{ + return std::any_of(xi.daughterIds.begin(), xi.daughterIds.end(), [trackId](int id) { return static_cast(id) == trackId; }); +} + +/// Mode A canonical candidate [Xi, K0s]; the roles are fixed by the object types. +inline CanonicalizationResult canonicalizeXiK0s(CascadeSnapshot const& xi, V0Snapshot const& k0s) +{ + CanonicalizationResult result; + for (const auto& id : k0s.daughterIds) { + if (sharesDaughter(xi, id)) { + result.status = BuildStatus::ReusedDaughter; + return result; + } + } + if (xi.sign != 1 && xi.sign != -1) { + result.status = BuildStatus::InvalidChargePattern; + return result; + } + if (!hasFiniteMomentum(xi.px, xi.py, xi.pz) || !hasFiniteMomentum(k0s.px, k0s.py, k0s.pz)) { + result.status = BuildStatus::InvalidMomentum; + return result; + } + result.candidate.xi = xi; + result.candidate.k0s = k0s; + result.status = BuildStatus::Ok; + return result; +} + +/// Charge pattern of (Xi, pion, kaon); 0 is the Omega(2012) signal pattern. +inline uint8_t chargePattern(int xiSign, int pionSign, int kaonSign) +{ + uint8_t pattern = kSignalPattern; + if (pionSign == xiSign) { + pattern |= kWrongSignPion; + } + if (kaonSign != xiSign) { + pattern |= kWrongSignKaon; + } + return pattern; +} + +/// Mode B canonical candidate [Xi, pion, kaon]; the roles are fixed by the selected species, the charge pattern is metadata. +inline CanonicalizationResult canonicalizeXi1530K(CascadeSnapshot const& xi, TrackSnapshot const& pion, TrackSnapshot const& kaon) +{ + CanonicalizationResult result; + if (pion.sourceTrackId == kaon.sourceTrackId || sharesDaughter(xi, pion.sourceTrackId) || sharesDaughter(xi, kaon.sourceTrackId)) { + result.status = BuildStatus::ReusedDaughter; + return result; + } + if ((xi.sign != 1 && xi.sign != -1) || (pion.charge != 1 && pion.charge != -1) || (kaon.charge != 1 && kaon.charge != -1)) { + result.status = BuildStatus::InvalidChargePattern; + return result; + } + if (!hasFiniteMomentum(xi.px, xi.py, xi.pz) || !hasFiniteMomentum(pion.px, pion.py, pion.pz) || !hasFiniteMomentum(kaon.px, kaon.py, kaon.pz)) { + result.status = BuildStatus::InvalidMomentum; + return result; + } + result.candidate.xi = xi; + result.candidate.pion = pion; + result.candidate.kaon = kaon; + result.candidate.chargePattern = chargePattern(xi.sign, pion.charge, kaon.charge); + result.status = BuildStatus::Ok; + return result; +} + +namespace detail +{ +struct EncodedValue { + float value; + float valid; + float overflow; +}; + +// int8 x10 nSigma of ResoCascades/ResoV0s: saturated codes carry no usable value +inline EncodedValue encodeNSigma10(int8_t code) +{ + if (code == SaturatedLow || code == SaturatedHigh) { + return {.value = 0.f, .valid = 0.f, .overflow = 0.f}; + } + return {.value = static_cast(code) / NSigmaScale, .valid = 1.f, .overflow = 0.f}; +} + +// micro001 nSigma, as the K1 contract +inline EncodedValue encodePID(float decoded) +{ + if (std::isnan(decoded)) { + return {.value = 0.f, .valid = 0.f, .overflow = 0.f}; + } + if (std::isinf(decoded)) { + return {.value = std::signbit(decoded) ? -MicroPidOverflow : MicroPidOverflow, .valid = 1.f, .overflow = 1.f}; + } + return {.value = decoded, .valid = 1.f, .overflow = 0.f}; +} + +// micro001 DCA, as the K1 contract +inline EncodedValue encodeDCA(float decoded) +{ + if (!std::isfinite(decoded)) { + return {.value = 0.f, .valid = 0.f, .overflow = 0.f}; + } + const bool overflow = decoded == o2::aod::resomicrodaughter001::DCAEncoding::MaxDCA; + return {.value = decoded, .valid = 1.f, .overflow = overflow ? 1.f : 0.f}; +} + +template +class Writer +{ + public: + explicit Writer(std::array& output) : out(output) {} + std::size_t index = 0; + // A count beyond N is detected by finishPack(); nothing is written past the end. + void add(float value) + { + if (index < N) { + out[index] = value; + } + ++index; + } + void add(double value) { add(static_cast(value)); } + void add(bool value) { add(value ? 1.f : 0.f); } + void addValid(EncodedValue value) + { + add(value.value); + add(value.valid); + } + void addFull(EncodedValue value) + { + add(value.value); + add(value.valid); + add(value.overflow); + } + + private: + std::array& out; +}; + +using LorentzVector = ROOT::Math::PxPyPzMVector; + +inline LorentzVector fourMomentum(float px, float py, float pz, double mass) +{ + return {px, py, pz, mass}; +} + +template +void appendCascade(Writer& w, CascadeSnapshot const& xi) +{ + w.add(xi.v0CosPA); + w.add(xi.cascCosPA); + w.add(xi.v0DaughDCA); + w.add(xi.cascDaughDCA); + w.add(std::abs(xi.dcaPosToPV)); + w.add(std::abs(xi.dcaNegToPV)); + w.add(std::abs(xi.dcaBachToPV)); + w.add(std::abs(xi.dcaV0ToPV)); + w.add(std::abs(xi.dcaXYToPV)); + w.add(std::abs(xi.dcaZToPV)); + w.add(xi.v0Radius); + w.add(xi.cascRadius); + w.add(static_cast(xi.mLambda) - o2::constants::physics::MassLambda0); + w.add(static_cast(xi.mXi) - o2::constants::physics::MassXiMinus); + w.add(xi.decayLength); + w.add(xi.properLength); + // Roles: the (anti)proton is the positive daughter of a Xi-, the negative daughter of a Xi+ + const bool isXiMinus = xi.sign < 0; + constexpr std::size_t Pi = 0, Pr = 2, Pos = 0, Neg = 1, Bach = 2; + const std::size_t baryon = isXiMinus ? Pos : Neg; + const std::size_t meson = isXiMinus ? Neg : Pos; + w.addValid(encodeNSigma10(xi.tpcNSigma10[baryon * NPidSpecies + Pr])); + w.addValid(encodeNSigma10(xi.tpcNSigma10[meson * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(xi.tpcNSigma10[Bach * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(xi.tofNSigma10[baryon * NPidSpecies + Pr])); + w.addValid(encodeNSigma10(xi.tofNSigma10[meson * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(xi.tofNSigma10[Bach * NPidSpecies + Pi])); + w.add(static_cast(xi.crossedRows[baryon])); + w.add(static_cast(xi.crossedRows[meson])); + w.add(static_cast(xi.crossedRows[Bach])); +} + +template +void appendV0(Writer& w, V0Snapshot const& k0s) +{ + w.add(k0s.cosPA); + w.add(k0s.daughDCA); + w.add(std::abs(k0s.dcaPosToPV)); + w.add(std::abs(k0s.dcaNegToPV)); + w.add(std::abs(k0s.dcaV0ToPV)); + w.add(k0s.radius); + w.add(static_cast(k0s.mK0Short) - o2::constants::physics::MassK0Short); + w.add(static_cast(k0s.mLambda) - o2::constants::physics::MassLambda0); + w.add(static_cast(k0s.mAntiLambda) - o2::constants::physics::MassLambda0); + w.add(k0s.alpha); + w.add(k0s.qtarm); + w.add(k0s.decayLength); + w.add(k0s.properLength); + constexpr std::size_t Pi = 0, Pos = 0, Neg = 1; + w.addValid(encodeNSigma10(k0s.tpcNSigma10[Pos * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(k0s.tpcNSigma10[Neg * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(k0s.tofNSigma10[Pos * NPidSpecies + Pi])); + w.addValid(encodeNSigma10(k0s.tofNSigma10[Neg * NPidSpecies + Pi])); + w.add(static_cast(k0s.crossedRows[Pos])); + w.add(static_cast(k0s.crossedRows[Neg])); +} + +template +void appendTrack(Writer& w, TrackSnapshot const& track) +{ + for (const float& decoded : track.tpcNSigma) { + w.addFull(encodePID(decoded)); + } + for (const float& decoded : track.tofNSigma) { + w.addFull(track.hasTOF ? encodePID(decoded) : EncodedValue{.value = 0.f, .valid = 0.f, .overflow = 0.f}); + } + w.addFull(encodeDCA(track.dcaXY)); + w.addFull(encodeDCA(track.dcaZ)); + w.add(track.passedPtDependentDCAxy); + w.add(track.passedPtDependentDCAz); + w.add(track.hasTOF); + w.add(static_cast(track.tpcNClsCrossedRows)); + for (const bool& hit : track.itsHit) { + w.add(hit); + } + w.add(track.isPVContributor); +} + +// delta eta, sin/cos delta phi, delta R, z_pT of two objects +template +bool appendPair(Writer& w, LorentzVector const& a, LorentzVector const& b) +{ + const double deta = a.Eta() - b.Eta(); + const double dphi = std::remainder(a.Phi() - b.Phi(), o2::constants::math::TwoPI); + const double sumPt = a.Pt() + b.Pt(); + if (!std::isfinite(deta) || !std::isfinite(dphi) || !(sumPt > 0.)) { + return false; + } + w.add(deta); + w.add(std::sin(dphi)); + w.add(std::cos(dphi)); + w.add(std::hypot(deta, dphi)); + w.add(std::min(a.Pt(), b.Pt()) / sumPt); + return true; +} + +// Cosine of the daughter direction in the mother rest frame w.r.t. the mother direction +inline double cosThetaStar(LorentzVector const& daughter, LorentzVector const& mother) +{ + const double p = mother.P(); + if (!(p > 0.) || !(mother.E() > p)) { + return std::numeric_limits::quiet_NaN(); + } + const ROOT::Math::Boost boost{mother.BoostToCM()}; + const auto inRest = boost(daughter); + const double pStar = inRest.P(); + if (!(pStar > 0.)) { + return std::numeric_limits::quiet_NaN(); + } + return inRest.Vect().Dot(mother.Vect()) / (pStar * p); +} + +template +bool finishPack(FeaturePack& pack, std::size_t index, LorentzVector const& mother, double scalarSumPt) +{ + const bool allFinite = std::all_of(pack.master.begin(), pack.master.end(), [](float x) { return std::isfinite(x); }); + if (index != N || !allFinite) { + pack.status = BuildStatus::InvalidContract; + return false; + } + pack.kinematics.mass = static_cast(mother.M()); + pack.kinematics.pt = static_cast(mother.Pt()); + pack.kinematics.eta = static_cast(mother.Eta()); + pack.kinematics.phi = static_cast(mother.Phi()); + pack.kinematics.scalarSumPt = static_cast(scalarSumPt); + pack.status = BuildStatus::Ok; + return true; +} + +template +constexpr bool isStrictlyIncreasingBelow(std::array const& indices, std::size_t bound) +{ + for (std::size_t i = 0; i < N; ++i) { + if (indices[i] >= bound || (i > 0 && indices[i - 1] >= indices[i])) { + return false; + } + } + return true; +} +} // namespace detail + +// Names, projections and SHA-256 are generated from feature_contract_omega2012_v1.json (one section per mode). +inline constexpr std::array XiK0sMasterFeatureNames{ + "xi.v0_cos_pa", + "xi.casc_cos_pa", + "xi.v0_daughter_dca", + "xi.casc_daughter_dca", + "xi.abs_dca_pos_to_pv", + "xi.abs_dca_neg_to_pv", + "xi.abs_dca_bach_to_pv", + "xi.abs_dca_v0_to_pv", + "xi.abs_dca_xy_to_pv", + "xi.abs_dca_z_to_pv", + "xi.v0_radius", + "xi.casc_radius", + "xi.delta_mass_lambda", + "xi.delta_mass_xi", + "xi.decay_length", + "xi.proper_length", + "xi.baryon_tpc_nsigma_pr", + "xi.baryon_tpc_nsigma_pr_valid", + "xi.meson_tpc_nsigma_pi", + "xi.meson_tpc_nsigma_pi_valid", + "xi.bach_tpc_nsigma_pi", + "xi.bach_tpc_nsigma_pi_valid", + "xi.baryon_tof_nsigma_pr", + "xi.baryon_tof_nsigma_pr_valid", + "xi.meson_tof_nsigma_pi", + "xi.meson_tof_nsigma_pi_valid", + "xi.bach_tof_nsigma_pi", + "xi.bach_tof_nsigma_pi_valid", + "xi.baryon_crossed_rows", + "xi.meson_crossed_rows", + "xi.bach_crossed_rows", + "k0s.cos_pa", + "k0s.daughter_dca", + "k0s.abs_dca_pos_to_pv", + "k0s.abs_dca_neg_to_pv", + "k0s.abs_dca_v0_to_pv", + "k0s.radius", + "k0s.delta_mass_k0s", + "k0s.delta_mass_lambda", + "k0s.delta_mass_antilambda", + "k0s.arm_alpha", + "k0s.arm_qt", + "k0s.decay_length", + "k0s.proper_length", + "k0s.pos_tpc_nsigma_pi", + "k0s.pos_tpc_nsigma_pi_valid", + "k0s.neg_tpc_nsigma_pi", + "k0s.neg_tpc_nsigma_pi_valid", + "k0s.pos_tof_nsigma_pi", + "k0s.pos_tof_nsigma_pi_valid", + "k0s.neg_tof_nsigma_pi", + "k0s.neg_tof_nsigma_pi_valid", + "k0s.pos_crossed_rows", + "k0s.neg_crossed_rows", + "xi.pt_fraction", + "k0s.pt_fraction", + "xi__k0s.delta_eta", + "xi__k0s.sin_delta_phi", + "xi__k0s.cos_delta_phi", + "xi__k0s.delta_r", + "xi__k0s.z_pt", + "xi__k0s.opening_angle", + "xi__k0s.cos_theta_star"}; +inline constexpr std::array XiK0sDetectorV1Projection{ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 37, + 38, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 52, + 53}; +inline constexpr std::array XiK0sRelationalV1Projection{ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 37, + 38, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 52, + 53, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + 61, + 62}; +inline constexpr std::string_view XiK0sFeatureContractSha256 = "b767a5643da563246295f543ec7296c9e21b00fd921fef88acc7b19f3b24942c"; + +inline constexpr std::array Xi1530KMasterFeatureNames{ + "xi.v0_cos_pa", + "xi.casc_cos_pa", + "xi.v0_daughter_dca", + "xi.casc_daughter_dca", + "xi.abs_dca_pos_to_pv", + "xi.abs_dca_neg_to_pv", + "xi.abs_dca_bach_to_pv", + "xi.abs_dca_v0_to_pv", + "xi.abs_dca_xy_to_pv", + "xi.abs_dca_z_to_pv", + "xi.v0_radius", + "xi.casc_radius", + "xi.delta_mass_lambda", + "xi.delta_mass_xi", + "xi.decay_length", + "xi.proper_length", + "xi.baryon_tpc_nsigma_pr", + "xi.baryon_tpc_nsigma_pr_valid", + "xi.meson_tpc_nsigma_pi", + "xi.meson_tpc_nsigma_pi_valid", + "xi.bach_tpc_nsigma_pi", + "xi.bach_tpc_nsigma_pi_valid", + "xi.baryon_tof_nsigma_pr", + "xi.baryon_tof_nsigma_pr_valid", + "xi.meson_tof_nsigma_pi", + "xi.meson_tof_nsigma_pi_valid", + "xi.bach_tof_nsigma_pi", + "xi.bach_tof_nsigma_pi_valid", + "xi.baryon_crossed_rows", + "xi.meson_crossed_rows", + "xi.bach_crossed_rows", + "pion.tpc_nsigma_pi", + "pion.tpc_nsigma_pi_valid", + "pion.tpc_nsigma_pi_overflow", + "pion.tpc_nsigma_ka", + "pion.tpc_nsigma_ka_valid", + "pion.tpc_nsigma_ka_overflow", + "pion.tpc_nsigma_pr", + "pion.tpc_nsigma_pr_valid", + "pion.tpc_nsigma_pr_overflow", + "pion.tof_nsigma_pi", + "pion.tof_nsigma_pi_valid", + "pion.tof_nsigma_pi_overflow", + "pion.tof_nsigma_ka", + "pion.tof_nsigma_ka_valid", + "pion.tof_nsigma_ka_overflow", + "pion.tof_nsigma_pr", + "pion.tof_nsigma_pr_valid", + "pion.tof_nsigma_pr_overflow", + "pion.abs_dca_xy", + "pion.abs_dca_xy_valid", + "pion.abs_dca_xy_overflow", + "pion.abs_dca_z", + "pion.abs_dca_z_valid", + "pion.abs_dca_z_overflow", + "pion.passed_ptdep_dca_xy", + "pion.passed_ptdep_dca_z", + "pion.has_tof", + "pion.tpc_crossed_rows", + "pion.its_hit_l0", + "pion.its_hit_l1", + "pion.its_hit_l2", + "pion.its_hit_l3", + "pion.its_hit_l4", + "pion.its_hit_l5", + "pion.its_hit_l6", + "pion.is_pv_contributor", + "kaon.tpc_nsigma_pi", + "kaon.tpc_nsigma_pi_valid", + "kaon.tpc_nsigma_pi_overflow", + "kaon.tpc_nsigma_ka", + "kaon.tpc_nsigma_ka_valid", + "kaon.tpc_nsigma_ka_overflow", + "kaon.tpc_nsigma_pr", + "kaon.tpc_nsigma_pr_valid", + "kaon.tpc_nsigma_pr_overflow", + "kaon.tof_nsigma_pi", + "kaon.tof_nsigma_pi_valid", + "kaon.tof_nsigma_pi_overflow", + "kaon.tof_nsigma_ka", + "kaon.tof_nsigma_ka_valid", + "kaon.tof_nsigma_ka_overflow", + "kaon.tof_nsigma_pr", + "kaon.tof_nsigma_pr_valid", + "kaon.tof_nsigma_pr_overflow", + "kaon.abs_dca_xy", + "kaon.abs_dca_xy_valid", + "kaon.abs_dca_xy_overflow", + "kaon.abs_dca_z", + "kaon.abs_dca_z_valid", + "kaon.abs_dca_z_overflow", + "kaon.passed_ptdep_dca_xy", + "kaon.passed_ptdep_dca_z", + "kaon.has_tof", + "kaon.tpc_crossed_rows", + "kaon.its_hit_l0", + "kaon.its_hit_l1", + "kaon.its_hit_l2", + "kaon.its_hit_l3", + "kaon.its_hit_l4", + "kaon.its_hit_l5", + "kaon.its_hit_l6", + "kaon.is_pv_contributor", + "xi.pt_fraction", + "pion.pt_fraction", + "kaon.pt_fraction", + "xi__pion.delta_eta", + "xi__pion.sin_delta_phi", + "xi__pion.cos_delta_phi", + "xi__pion.delta_r", + "xi__pion.z_pt", + "xi__kaon.delta_eta", + "xi__kaon.sin_delta_phi", + "xi__kaon.cos_delta_phi", + "xi__kaon.delta_r", + "xi__kaon.z_pt", + "pion__kaon.delta_eta", + "pion__kaon.sin_delta_phi", + "pion__kaon.cos_delta_phi", + "pion__kaon.delta_r", + "pion__kaon.z_pt", + "xi1530__kaon.opening_angle", + "xi1530__kaon.cos_theta_star", + "mass_xi_pi", + "mass_xi_k", + "mass_pi_k"}; +inline constexpr std::array Xi1530KDetectorV1Projection{ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 37, + 38, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 52, + 53, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + 61, + 62, + 63, + 64, + 65, + 66, + 67, + 68, + 69, + 70, + 71, + 72, + 73, + 74, + 75, + 76, + 77, + 78, + 79, + 80, + 81, + 82, + 83, + 84, + 85, + 86, + 87, + 88, + 89, + 90, + 91, + 92, + 93, + 94, + 95, + 96, + 97, + 98, + 99, + 100, + 101, + 102}; +inline constexpr std::array Xi1530KRelationalV1Projection{ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 37, + 38, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 52, + 53, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + 61, + 62, + 63, + 64, + 65, + 66, + 67, + 68, + 69, + 70, + 71, + 72, + 73, + 74, + 75, + 76, + 77, + 78, + 79, + 80, + 81, + 82, + 83, + 84, + 85, + 86, + 87, + 88, + 89, + 90, + 91, + 92, + 93, + 94, + 95, + 96, + 97, + 98, + 99, + 100, + 101, + 102, + 103, + 104, + 105, + 106, + 107, + 108, + 109, + 110, + 111, + 112, + 113, + 114, + 115, + 116, + 117, + 118, + 119, + 120, + 121, + 122, + 123}; +inline constexpr std::array Xi1530KStudyV1Projection{ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34, + 35, + 36, + 37, + 38, + 39, + 40, + 41, + 42, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 52, + 53, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + 61, + 62, + 63, + 64, + 65, + 66, + 67, + 68, + 69, + 70, + 71, + 72, + 73, + 74, + 75, + 76, + 77, + 78, + 79, + 80, + 81, + 82, + 83, + 84, + 85, + 86, + 87, + 88, + 89, + 90, + 91, + 92, + 93, + 94, + 95, + 96, + 97, + 98, + 99, + 100, + 101, + 102, + 103, + 104, + 105, + 106, + 107, + 108, + 109, + 110, + 111, + 112, + 113, + 114, + 115, + 116, + 117, + 118, + 119, + 120, + 121, + 122, + 123, + 124, + 125}; +inline constexpr std::string_view Xi1530KFeatureContractSha256 = "ca244fe167a599f9bc60d3909630edca5c797bb7b84d75d80944e6ad1a44aec5"; +static_assert(XiK0sMasterFeatureNames.size() == XiK0sNMasterFeatures, "XiK0s feature name count must match XiK0sNMasterFeatures"); +static_assert(Xi1530KMasterFeatureNames.size() == Xi1530KNMasterFeatures, "Xi1530K feature name count must match Xi1530KNMasterFeatures"); +static_assert(detail::isStrictlyIncreasingBelow(XiK0sDetectorV1Projection, XiK0sNMasterFeatures), "XiK0s DetectorV1 projection must be strictly increasing"); +static_assert(detail::isStrictlyIncreasingBelow(XiK0sRelationalV1Projection, XiK0sNMasterFeatures), "XiK0s RelationalV1 projection must be strictly increasing"); +static_assert(detail::isStrictlyIncreasingBelow(Xi1530KDetectorV1Projection, Xi1530KNMasterFeatures), "Xi1530K DetectorV1 projection must be strictly increasing"); +static_assert(detail::isStrictlyIncreasingBelow(Xi1530KRelationalV1Projection, Xi1530KNMasterFeatures), "Xi1530K RelationalV1 projection must be strictly increasing"); +static_assert(detail::isStrictlyIncreasingBelow(Xi1530KStudyV1Projection, Xi1530KNMasterFeatures), "Xi1530K StudyV1 projection must be strictly increasing"); + +/// Mode A master features (XiK0s contract) +inline XiK0sFeaturePack buildXiK0sFeatures(XiK0sCandidateSnapshot const& candidate) +{ + using o2::constants::physics::MassK0Short; + using o2::constants::physics::MassXiMinus; + XiK0sFeaturePack pack; + const auto& xi = candidate.xi; + const auto& k0s = candidate.k0s; + if (!hasFiniteMomentum(xi.px, xi.py, xi.pz) || !hasFiniteMomentum(k0s.px, k0s.py, k0s.pz)) { + pack.status = BuildStatus::InvalidMomentum; + return pack; + } + const auto pXi = detail::fourMomentum(xi.px, xi.py, xi.pz, MassXiMinus); + const auto pK0s = detail::fourMomentum(k0s.px, k0s.py, k0s.pz, MassK0Short); + const auto mother = pXi + pK0s; + const double sumPt = pXi.Pt() + pK0s.Pt(); + if (!(sumPt > 0.)) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + detail::Writer w{pack.master}; + detail::appendCascade(w, xi); + detail::appendV0(w, k0s); + w.add(pXi.Pt() / sumPt); + w.add(pK0s.Pt() / sumPt); + if (!detail::appendPair(w, pXi, pK0s)) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + w.add(ROOT::Math::VectorUtil::Angle(pXi, pK0s)); + w.add(detail::cosThetaStar(pXi, mother)); + detail::finishPack(pack, w.index, mother, sumPt); + return pack; +} + +/// Mode B master features (Xi1530K contract) +inline Xi1530KFeaturePack buildXi1530KFeatures(Xi1530KCandidateSnapshot const& candidate) +{ + using o2::constants::physics::MassKaonCharged; + using o2::constants::physics::MassPionCharged; + using o2::constants::physics::MassXiMinus; + Xi1530KFeaturePack pack; + const auto& xi = candidate.xi; + const auto& pion = candidate.pion; + const auto& kaon = candidate.kaon; + if (!hasFiniteMomentum(xi.px, xi.py, xi.pz) || !hasFiniteMomentum(pion.px, pion.py, pion.pz) || !hasFiniteMomentum(kaon.px, kaon.py, kaon.pz)) { + pack.status = BuildStatus::InvalidMomentum; + return pack; + } + if (pion.sourceTrackId == kaon.sourceTrackId || sharesDaughter(xi, pion.sourceTrackId) || sharesDaughter(xi, kaon.sourceTrackId)) { + pack.status = BuildStatus::ReusedDaughter; + return pack; + } + const auto pXi = detail::fourMomentum(xi.px, xi.py, xi.pz, MassXiMinus); + const auto pPion = detail::fourMomentum(pion.px, pion.py, pion.pz, MassPionCharged); + const auto pKaon = detail::fourMomentum(kaon.px, kaon.py, kaon.pz, MassKaonCharged); + const auto xi1530 = pXi + pPion; + const auto mother = xi1530 + pKaon; + const double sumPt = pXi.Pt() + pPion.Pt() + pKaon.Pt(); + if (!(sumPt > 0.)) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + detail::Writer w{pack.master}; + detail::appendCascade(w, xi); + detail::appendTrack(w, pion); + detail::appendTrack(w, kaon); + w.add(pXi.Pt() / sumPt); + w.add(pPion.Pt() / sumPt); + w.add(pKaon.Pt() / sumPt); + if (!detail::appendPair(w, pXi, pPion) || !detail::appendPair(w, pXi, pKaon) || !detail::appendPair(w, pPion, pKaon)) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + w.add(ROOT::Math::VectorUtil::Angle(xi1530, pKaon)); + w.add(detail::cosThetaStar(xi1530, mother)); + w.add(xi1530.M()); + w.add((pXi + pKaon).M()); + w.add((pPion + pKaon).M()); + detail::finishPack(pack, w.index, mother, sumPt); + return pack; +} + +template +std::vector project(FeaturePack const& pack, std::array const& indices) +{ + std::vector projected; + projected.reserve(M); + for (const auto& i : indices) { + projected.push_back(pack.master[i]); + } + return projected; +} + +inline std::vector projectFeatures(XiK0sFeaturePack const& pack, Profile profile) +{ + if (pack.status != BuildStatus::Ok) { + return {}; + } + switch (profile) { + case Profile::DetectorV1: + return project(pack, XiK0sDetectorV1Projection); + case Profile::RelationalV1: + return project(pack, XiK0sRelationalV1Projection); + default: + return {}; // mode A has no study-only features + } +} + +inline std::vector projectFeatures(Xi1530KFeaturePack const& pack, Profile profile) +{ + if (pack.status != BuildStatus::Ok) { + return {}; + } + switch (profile) { + case Profile::DetectorV1: + return project(pack, Xi1530KDetectorV1Projection); + case Profile::RelationalV1: + return project(pack, Xi1530KRelationalV1Projection); + case Profile::StudyV1: + return project(pack, Xi1530KStudyV1Projection); + default: + return {}; + } +} +} // namespace o2::analysis::omega2012ml + +#endif // PWGLF_CORE_OMEGA2012MLFEATURES_H_ diff --git a/PWGLF/DataModel/LFOmega2012MlTables.h b/PWGLF/DataModel/LFOmega2012MlTables.h new file mode 100644 index 00000000000..5cca99cbd46 --- /dev/null +++ b/PWGLF/DataModel/LFOmega2012MlTables.h @@ -0,0 +1,239 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file LFOmega2012MlTables.h +/// \brief Derived Omega(2012) training and audit tables, separate for the XiK0s and Xi1530K decay modes +/// \author Bong-Hwi Lim +/// +/// Mode A (XiK0s): Omega(2012)- -> Xi- K0S. Mode B (Xi1530K): Omega(2012)- -> Xi(1530)0 K- -> Xi- pi+ K-. +/// Only the input-object tables (events, cascades, V0s, tracks) are common; every event row belongs to one +/// decay mode (column omDecayMode), so the object rows of the two modes never mix either. +/// +/// Cumulative pass bits (each bit requires all lower bits): +/// mode A: 1 loose (valid canonical candidate inside the rapidity window), 2 Xi selection, 4 K0s selection, +/// 8 Xi-K0s kinematic (opening-angle) cut; selected = 15. +/// mode B: 1 loose (valid canonical candidate inside the rapidity window), 2 Xi selection, 4 track quality of the +/// pion and the kaon, 8 TOF requirement and PID of the pion and the kaon, 16 Xi(1530) mass window; +/// selected = 31. +#ifndef PWGLF_DATAMODEL_LFOMEGA2012MLTABLES_H_ +#define PWGLF_DATAMODEL_LFOMEGA2012MLTABLES_H_ + +#include +#include + +#include + +namespace o2::aod +{ +namespace omega2012ml +{ +// Persisted relations target only these derived tables; source AO2D row numbers are scalar audit values. +// Events +DECLARE_SOA_COLUMN(OmDecayMode, omDecayMode, uint8_t); //! 0 XiK0s, 1 Xi1530K: the mode whose candidates reference this event row +DECLARE_SOA_COLUMN(OmRecoCollisionId, omRecoCollisionId, int64_t); //! row of the reduced ResoCollisions_001 collision in the input DF +DECLARE_SOA_COLUMN(OmPosZ, omPosZ, float); +DECLARE_SOA_COLUMN(OmBField, omBField, float); +DECLARE_SOA_COLUMN(OmCentrality, omCentrality, float); +DECLARE_SOA_COLUMN(OmMultiplicity, omMultiplicity, float); +DECLARE_SOA_COLUMN(OmRecINELgt0, omRecINELgt0, bool); +// Common object columns +DECLARE_SOA_COLUMN(OmSourceRow, omSourceRow, int64_t); //! row of the object in the input ResoCascades / ResoV0s table +DECLARE_SOA_COLUMN(OmPx, omPx, float); +DECLARE_SOA_COLUMN(OmPy, omPy, float); +DECLARE_SOA_COLUMN(OmPz, omPz, float); +DECLARE_SOA_COLUMN(OmV0CosPA, omV0CosPA, float); +DECLARE_SOA_COLUMN(OmV0DaughDCA, omV0DaughDCA, float); +DECLARE_SOA_COLUMN(OmDcaPosToPV, omDcaPosToPV, float); +DECLARE_SOA_COLUMN(OmDcaNegToPV, omDcaNegToPV, float); +DECLARE_SOA_COLUMN(OmDcaV0ToPV, omDcaV0ToPV, float); +DECLARE_SOA_COLUMN(OmV0Radius, omV0Radius, float); +DECLARE_SOA_COLUMN(OmMassLambda, omMassLambda, float); +DECLARE_SOA_COLUMN(OmDecayVtxX, omDecayVtxX, float); +DECLARE_SOA_COLUMN(OmDecayVtxY, omDecayVtxY, float); +DECLARE_SOA_COLUMN(OmDecayVtxZ, omDecayVtxZ, float); +DECLARE_SOA_COLUMN(OmTpcPosPi10, omTpcPosPi10, int8_t); //! TPC nSigma x10 of the positive daughter (pion hypothesis), as in the reduced tables +DECLARE_SOA_COLUMN(OmTpcPosKa10, omTpcPosKa10, int8_t); +DECLARE_SOA_COLUMN(OmTpcPosPr10, omTpcPosPr10, int8_t); +DECLARE_SOA_COLUMN(OmTpcNegPi10, omTpcNegPi10, int8_t); +DECLARE_SOA_COLUMN(OmTpcNegKa10, omTpcNegKa10, int8_t); +DECLARE_SOA_COLUMN(OmTpcNegPr10, omTpcNegPr10, int8_t); +DECLARE_SOA_COLUMN(OmTofPosPi10, omTofPosPi10, int8_t); //! TOF nSigma x10 of the positive daughter (pion hypothesis), as in the reduced tables +DECLARE_SOA_COLUMN(OmTofPosKa10, omTofPosKa10, int8_t); +DECLARE_SOA_COLUMN(OmTofPosPr10, omTofPosPr10, int8_t); +DECLARE_SOA_COLUMN(OmTofNegPi10, omTofNegPi10, int8_t); +DECLARE_SOA_COLUMN(OmTofNegKa10, omTofNegKa10, int8_t); +DECLARE_SOA_COLUMN(OmTofNegPr10, omTofNegPr10, int8_t); +DECLARE_SOA_COLUMN(OmCrossedRowsPos, omCrossedRowsPos, uint8_t); +DECLARE_SOA_COLUMN(OmCrossedRowsNeg, omCrossedRowsNeg, uint8_t); +// Cascade-only columns +DECLARE_SOA_COLUMN(OmSign, omSign, int8_t); +DECLARE_SOA_COLUMN(OmMassXi, omMassXi, float); +DECLARE_SOA_COLUMN(OmCascCosPA, omCascCosPA, float); +DECLARE_SOA_COLUMN(OmCascDaughDCA, omCascDaughDCA, float); +DECLARE_SOA_COLUMN(OmDcaBachToPV, omDcaBachToPV, float); +DECLARE_SOA_COLUMN(OmDcaXYCascToPV, omDcaXYCascToPV, float); +DECLARE_SOA_COLUMN(OmDcaZCascToPV, omDcaZCascToPV, float); +DECLARE_SOA_COLUMN(OmCascRadius, omCascRadius, float); +DECLARE_SOA_COLUMN(OmTpcBachPi10, omTpcBachPi10, int8_t); +DECLARE_SOA_COLUMN(OmTpcBachKa10, omTpcBachKa10, int8_t); +DECLARE_SOA_COLUMN(OmTpcBachPr10, omTpcBachPr10, int8_t); +DECLARE_SOA_COLUMN(OmTofBachPi10, omTofBachPi10, int8_t); +DECLARE_SOA_COLUMN(OmTofBachKa10, omTofBachKa10, int8_t); +DECLARE_SOA_COLUMN(OmTofBachPr10, omTofBachPr10, int8_t); +DECLARE_SOA_COLUMN(OmCrossedRowsBach, omCrossedRowsBach, uint8_t); +// The daughter-ID arrays mirror the persistent ResoCascades / ResoV0s columns. +DECLARE_SOA_COLUMN(OmCascadeIndices, omCascadeIndices, int[3]); //! source track IDs of the cascade daughters (positive, negative, bachelor) +// V0-only columns +DECLARE_SOA_COLUMN(OmMassK0Short, omMassK0Short, float); +DECLARE_SOA_COLUMN(OmMassAntiLambda, omMassAntiLambda, float); +DECLARE_SOA_COLUMN(OmArmAlpha, omArmAlpha, float); +DECLARE_SOA_COLUMN(OmArmQt, omArmQt, float); +DECLARE_SOA_COLUMN(OmV0Indices, omV0Indices, int[2]); //! source track IDs of the V0 daughters (positive, negative) +// Track columns +DECLARE_SOA_COLUMN(OmSourceTrackId, omSourceTrackId, int64_t); //! trackId of the reduced micro track +DECLARE_SOA_COLUMN(OmPidPi, omPidPi, uint8_t); +DECLARE_SOA_COLUMN(OmPidKa, omPidKa, uint8_t); +DECLARE_SOA_COLUMN(OmPidPr, omPidPr, uint8_t); +DECLARE_SOA_COLUMN(OmSelectionFlags, omSelectionFlags, uint8_t); +DECLARE_SOA_COLUMN(OmTrackFlags, omTrackFlags, uint8_t); +DECLARE_SOA_COLUMN(OmCrossedRows, omCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(OmItsClusterMap, omItsClusterMap, uint8_t); +// Candidate columns +DECLARE_SOA_COLUMN(OmMass, omMass, float); +DECLARE_SOA_COLUMN(OmPt, omPt, float); +DECLARE_SOA_COLUMN(OmY, omY, float); +DECLARE_SOA_COLUMN(OmEta, omEta, float); +DECLARE_SOA_COLUMN(OmPhi, omPhi, float); +DECLARE_SOA_COLUMN(OmOpeningAngle, omOpeningAngle, float); //! Xi-K0s opening angle alpha_oa of the kinematic cut +DECLARE_SOA_COLUMN(OmCharge, omCharge, int8_t); //! sign of the Xi +DECLARE_SOA_COLUMN(OmMassXiPi, omMassXiPi, float); +DECLARE_SOA_COLUMN(OmMassXiK, omMassXiK, float); +DECLARE_SOA_COLUMN(OmMassPiK, omMassPiK, float); +DECLARE_SOA_COLUMN(OmChargePattern, omChargePattern, uint8_t); //! 0 signal (pion opposite to the Xi, kaon equal), 1 wrong-sign pion, 2 wrong-sign kaon, 3 both +DECLARE_SOA_COLUMN(OmPassBits, omPassBits, uint16_t); //! cumulative pass bits of the mode, see the file header +DECLARE_SOA_COLUMN(OmXiK0sFeatures, omXiK0sFeatures, float[63]); //! master features of the XiK0s contract (PWGLF/Core/Omega2012MlFeatures.h) +DECLARE_SOA_COLUMN(OmXi1530KFeatures, omXi1530KFeatures, float[126]); //! master features of the Xi1530K contract (PWGLF/Core/Omega2012MlFeatures.h) +DECLARE_SOA_COLUMN(OmFeatureStatus, omFeatureStatus, uint8_t); //! o2::analysis::omega2012ml::BuildStatus; only Ok (0) rows are written +// Truth and generated audit +DECLARE_SOA_COLUMN(OmTruthStatus, omTruthStatus, uint8_t); //! 0 data, 1 matched, 2 unmatched +DECLARE_SOA_COLUMN(OmMotherPdg, omMotherPdg, int32_t); +DECLARE_SOA_COLUMN(OmMotherId, omMotherId, int64_t); +DECLARE_SOA_COLUMN(OmXiMotherPdg, omXiMotherPdg, int32_t); //! immediate-mother PDG of the Xi (any candidate) +DECLARE_SOA_COLUMN(OmV0MotherPdg, omV0MotherPdg, int32_t); //! immediate-mother PDG of the V0 (any candidate; 311 shows a K0bar intermediate) +DECLARE_SOA_COLUMN(OmXi1530Id, omXi1530Id, int64_t); //! immediate-mother ID of the Xi (the Xi(1530)0 candidate) +DECLARE_SOA_COLUMN(OmOriginalMcParticleId, omOriginalMcParticleId, int64_t); +DECLARE_SOA_COLUMN(OmPdg, omPdg, int32_t); +DECLARE_SOA_COLUMN(OmDaughterPdg1, omDaughterPdg1, int32_t); +DECLARE_SOA_COLUMN(OmDaughterPdg2, omDaughterPdg2, int32_t); +DECLARE_SOA_COLUMN(OmGenChannel, omGenChannel, uint8_t); //! immediate channel of a generated Omega(2012): 0 other, 1 XiK0s, 2 Xi1530K +DECLARE_SOA_COLUMN(OmGenPt, omGenPt, float); +DECLARE_SOA_COLUMN(OmGenY, omGenY, float); +} // namespace omega2012ml + +DECLARE_SOA_TABLE(Omega2012MlEvents, "AOD", "OMMLEVENT", + o2::soa::Index<>, omega2012ml::OmDecayMode, omega2012ml::OmRecoCollisionId, + omega2012ml::OmPosZ, omega2012ml::OmBField, omega2012ml::OmCentrality, + omega2012ml::OmMultiplicity, omega2012ml::OmRecINELgt0); +namespace omega2012ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlEvent, omega2012MlEvent, int, Omega2012MlEvents, ""); +} // namespace omega2012ml + +DECLARE_SOA_TABLE(Omega2012MlCascades, "AOD", "OMMLCASCADE", + o2::soa::Index<>, omega2012ml::Omega2012MlEventId, omega2012ml::OmSourceRow, + omega2012ml::OmPx, omega2012ml::OmPy, omega2012ml::OmPz, omega2012ml::OmSign, + omega2012ml::OmMassXi, omega2012ml::OmMassLambda, + omega2012ml::OmV0CosPA, omega2012ml::OmCascCosPA, omega2012ml::OmV0DaughDCA, omega2012ml::OmCascDaughDCA, + omega2012ml::OmDcaPosToPV, omega2012ml::OmDcaNegToPV, omega2012ml::OmDcaBachToPV, omega2012ml::OmDcaV0ToPV, + omega2012ml::OmDcaXYCascToPV, omega2012ml::OmDcaZCascToPV, + omega2012ml::OmV0Radius, omega2012ml::OmCascRadius, + omega2012ml::OmDecayVtxX, omega2012ml::OmDecayVtxY, omega2012ml::OmDecayVtxZ, + omega2012ml::OmTpcPosPi10, omega2012ml::OmTpcPosKa10, omega2012ml::OmTpcPosPr10, + omega2012ml::OmTpcNegPi10, omega2012ml::OmTpcNegKa10, omega2012ml::OmTpcNegPr10, + omega2012ml::OmTpcBachPi10, omega2012ml::OmTpcBachKa10, omega2012ml::OmTpcBachPr10, + omega2012ml::OmTofPosPi10, omega2012ml::OmTofPosKa10, omega2012ml::OmTofPosPr10, + omega2012ml::OmTofNegPi10, omega2012ml::OmTofNegKa10, omega2012ml::OmTofNegPr10, + omega2012ml::OmTofBachPi10, omega2012ml::OmTofBachKa10, omega2012ml::OmTofBachPr10, + omega2012ml::OmCrossedRowsPos, omega2012ml::OmCrossedRowsNeg, omega2012ml::OmCrossedRowsBach, + omega2012ml::OmCascadeIndices); +DECLARE_SOA_TABLE(Omega2012MlV0s, "AOD", "OMMLV0", + o2::soa::Index<>, omega2012ml::Omega2012MlEventId, omega2012ml::OmSourceRow, + omega2012ml::OmPx, omega2012ml::OmPy, omega2012ml::OmPz, + omega2012ml::OmMassK0Short, omega2012ml::OmMassLambda, omega2012ml::OmMassAntiLambda, + omega2012ml::OmV0CosPA, omega2012ml::OmV0DaughDCA, + omega2012ml::OmDcaPosToPV, omega2012ml::OmDcaNegToPV, omega2012ml::OmDcaV0ToPV, omega2012ml::OmV0Radius, + omega2012ml::OmDecayVtxX, omega2012ml::OmDecayVtxY, omega2012ml::OmDecayVtxZ, + omega2012ml::OmArmAlpha, omega2012ml::OmArmQt, + omega2012ml::OmTpcPosPi10, omega2012ml::OmTpcPosKa10, omega2012ml::OmTpcPosPr10, + omega2012ml::OmTpcNegPi10, omega2012ml::OmTpcNegKa10, omega2012ml::OmTpcNegPr10, + omega2012ml::OmTofPosPi10, omega2012ml::OmTofPosKa10, omega2012ml::OmTofPosPr10, + omega2012ml::OmTofNegPi10, omega2012ml::OmTofNegKa10, omega2012ml::OmTofNegPr10, + omega2012ml::OmCrossedRowsPos, omega2012ml::OmCrossedRowsNeg, + omega2012ml::OmV0Indices); +DECLARE_SOA_TABLE(Omega2012MlTracks, "AOD", "OMMLTRACK", + o2::soa::Index<>, omega2012ml::Omega2012MlEventId, omega2012ml::OmSourceTrackId, + omega2012ml::OmPx, omega2012ml::OmPy, omega2012ml::OmPz, + omega2012ml::OmPidPi, omega2012ml::OmPidKa, omega2012ml::OmPidPr, + omega2012ml::OmSelectionFlags, omega2012ml::OmTrackFlags, + omega2012ml::OmCrossedRows, omega2012ml::OmItsClusterMap); +namespace omega2012ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlCascade, omega2012MlCascade, int, Omega2012MlCascades, ""); +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlV0, omega2012MlV0, int, Omega2012MlV0s, ""); +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlPionTrack, omega2012MlPionTrack, int, Omega2012MlTracks, "_Pion"); +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlKaonTrack, omega2012MlKaonTrack, int, Omega2012MlTracks, "_Kaon"); +} // namespace omega2012ml + +// Mode A: Xi K0S +DECLARE_SOA_TABLE(Omega2012MlXiK0sCandidates, "AOD", "OMMLXIK0SCAND", + o2::soa::Index<>, omega2012ml::Omega2012MlEventId, + omega2012ml::Omega2012MlCascadeId, omega2012ml::Omega2012MlV0Id, + omega2012ml::OmMass, omega2012ml::OmPt, omega2012ml::OmY, omega2012ml::OmEta, omega2012ml::OmPhi, + omega2012ml::OmOpeningAngle, omega2012ml::OmCharge, omega2012ml::OmPassBits); +namespace omega2012ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlXiK0sCandidate, omega2012MlXiK0sCandidate, int, Omega2012MlXiK0sCandidates, ""); +} // namespace omega2012ml +DECLARE_SOA_TABLE(Omega2012MlXiK0sInputs, "AOD", "OMMLXIK0SINPUT", + o2::soa::Index<>, omega2012ml::Omega2012MlXiK0sCandidateId, + omega2012ml::OmXiK0sFeatures, omega2012ml::OmFeatureStatus); +DECLARE_SOA_TABLE(Omega2012MlXiK0sTruth, "AOD", "OMMLXIK0STRUTH", + o2::soa::Index<>, omega2012ml::Omega2012MlXiK0sCandidateId, + omega2012ml::OmTruthStatus, omega2012ml::OmMotherPdg, omega2012ml::OmMotherId, + omega2012ml::OmV0MotherPdg, omega2012ml::OmXiMotherPdg); + +// Mode B: Xi(1530)0 K- -> Xi- pi+ K- +DECLARE_SOA_TABLE(Omega2012MlXi1530KCandidates, "AOD", "OMMLXI1530KCAND", + o2::soa::Index<>, omega2012ml::Omega2012MlEventId, + omega2012ml::Omega2012MlCascadeId, omega2012ml::Omega2012MlPionTrackId, omega2012ml::Omega2012MlKaonTrackId, + omega2012ml::OmMass, omega2012ml::OmMassXiPi, omega2012ml::OmMassXiK, omega2012ml::OmMassPiK, + omega2012ml::OmPt, omega2012ml::OmY, omega2012ml::OmEta, omega2012ml::OmPhi, + omega2012ml::OmChargePattern, omega2012ml::OmPassBits); +namespace omega2012ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Omega2012MlXi1530KCandidate, omega2012MlXi1530KCandidate, int, Omega2012MlXi1530KCandidates, ""); +} // namespace omega2012ml +DECLARE_SOA_TABLE(Omega2012MlXi1530KInputs, "AOD", "OMMLXI1530KINP", + o2::soa::Index<>, omega2012ml::Omega2012MlXi1530KCandidateId, + omega2012ml::OmXi1530KFeatures, omega2012ml::OmFeatureStatus); +DECLARE_SOA_TABLE(Omega2012MlXi1530KTruth, "AOD", "OMMLXI1530KTRU", + o2::soa::Index<>, omega2012ml::Omega2012MlXi1530KCandidateId, + omega2012ml::OmTruthStatus, omega2012ml::OmMotherPdg, omega2012ml::OmMotherId, + omega2012ml::OmXi1530Id); + +// Generated Omega(2012) parents of selected reconstructed MC events (ResoMCParents_001) +DECLARE_SOA_TABLE(Omega2012MlGenAudit, "AOD", "OMMLGENAUDIT", + o2::soa::Index<>, omega2012ml::OmRecoCollisionId, omega2012ml::OmOriginalMcParticleId, + omega2012ml::OmPdg, omega2012ml::OmDaughterPdg1, omega2012ml::OmDaughterPdg2, + omega2012ml::OmGenChannel, omega2012ml::OmGenPt, omega2012ml::OmGenY); +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFOMEGA2012MLTABLES_H_ diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index b1f09db8d9a..9c9fd4bf987 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -189,6 +189,11 @@ o2physics_add_dpl_workflow(omega2012analysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(omega2012-training-table + SOURCES omega2012TrainingTable.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(xi1820analysis SOURCES xi1820Analysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Resonances/omega2012Analysis.cxx b/PWGLF/Tasks/Resonances/omega2012Analysis.cxx index cfbbab3d7d6..74146abee60 100644 --- a/PWGLF/Tasks/Resonances/omega2012Analysis.cxx +++ b/PWGLF/Tasks/Resonances/omega2012Analysis.cxx @@ -12,7 +12,27 @@ /// \file omega2012Analysis.cxx /// \brief Invariant Mass Reconstruction of Omega(2012) Resonance /// \author Bong-Hwi Lim - +/// +/// Two decay modes, analysed separately (never merged in any histogram or process function): +/// - mode A (XiK0s): Omega(2012)- -> Xi- K0S. processData, processMixedEvent, processMC, processMCGenerated. +/// Histograms: Event/, QAbefore/, QAafter/, omega2012/, MC/ (unchanged). +/// - mode B (Xi1530K): Omega(2012)- -> Xi(1530)0 K- -> Xi- pi+ K- with a charged kaon track. +/// processXi1530KMicro, processXi1530KTracks, processXi1530KMCMicro, processXi1530KMixedMicro. +/// Histograms: xi1530K/ (signal charge pattern) and xi1530K_wrongSign/ (charge-pattern controls). +/// Selection, candidate enumeration and truth classification: PWGLF/Core/Omega2012AnalysisCore.h. +/// +/// Configuration keys: all mode-A keys are unchanged. The former three-body (Xi pi K0S) process functions +/// processThreeBodyWithTracks / processThreeBodyWithMicroTracks are removed (that final state cannot be an Omega(2012)-). +/// Mode-B pion/kaon track selection: the common resonance TrackCuts with the prefix "trk.": +/// cPionPtMin -> trk.cMinPtcut, cPionEtaMax -> trk.cMaxEtacut, cPionDCAxyMax -> trk.cMaxDCArToPVcut, +/// cPionDCAzMax -> trk.cMaxDCAzToPVcut (default 0.15 cm), cPionTPCNClusMin -> trk.cfgTPCcluster; +/// further trk.* keys as in PWGLF/Core/ResoAnalysisSelectionCore.h. +/// Pion PID keys cPion* are unchanged; the kaon PID keys cKaon* have the same shape. +/// New mode-B keys: cXi1530UseMassWindow, cXi1530KFillWrongSign, cByPassTOF. + +#include "PWGLF/Core/Omega2012AnalysisCore.h" +#include "PWGLF/Core/Omega2012MlFeatures.h" +#include "PWGLF/Core/ResoAnalysisSelectionCore.h" #include "PWGLF/DataModel/LFResonanceTables.h" #include @@ -25,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -33,28 +54,24 @@ #include #include -#include #include -#include -#include -#include +#include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::constants::physics; +using namespace o2::analysis::omega2012; struct Omega2012Analysis { // Constants - static constexpr float kSmallNumber = 1e-10f; // Small number to avoid division by zero - static constexpr float kMaxDCAV0ToPV = 1.0f; // Maximum DCA of V0 to PV - static constexpr int kNumExpectedDaughters = 2; // Expected number of daughters for 2-body decay + static constexpr int NumExpectedDaughters = 2; // Expected number of daughters for 2-body decay SliceCache cache; Preslice perResoCollisionCasc = aod::resodaughter::resoCollisionId; Preslice perResoCollisionV0 = aod::resodaughter::resoCollisionId; Preslice perResoCollisionTrack = aod::resodaughter::resoCollisionId; - Preslice perResoCollisionMicroTrack = aod::resodaughter::resoCollisionId; + Preslice perResoCollisionMicroTrack = aod::resodaughter::resoCollisionId; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Axes @@ -62,109 +79,39 @@ struct Omega2012Analysis { ConfigurableAxis binsPtQA{"binsPtQA", {VARIABLE_WIDTH, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0}, "pT (QA)"}; ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 1., 5., 10., 30., 50., 70., 100., 110.}, "Centrality"}; - // Invariant mass range for Omega(2012) → Xi + K0s + // Invariant mass range for Omega(2012) Configurable cInvMassStart{"cInvMassStart", 1.6, "Invariant mass start (GeV/c^2)"}; Configurable cInvMassEnd{"cInvMassEnd", 2.2, "Invariant mass end (GeV/c^2)"}; Configurable cInvMassBins{"cInvMassBins", 600, "Invariant mass bins"}; - // Basic pre-selections (mirroring refs) - Configurable cMinPtcut{"cMinPtcut", 0.15, "Minimum pT for candidates"}; - Configurable cMaxEtaCut{"cMaxEtaCut", 0.8, "Maximum |eta|"}; - Configurable cfgRapidityCut{"cfgRapidityCut", 0.5, "Rapidity cut"}; - Configurable cRecoINELgt0{"cRecoINELgt0", false, "Apply Reco INEL>0 event selection"}; - Configurable cMCINELgt0{"cMCINELgt0", false, "Require generator INEL>0 in MC processes"}; - Configurable cMCVtxIn10{"cMCVtxIn10", false, "Require generator |vertex z| < 10 cm in MC processes"}; - Configurable cKinCuts{"cKinCuts", false, "Kinematic cuts for Xi-K0s opening angle"}; - Configurable> cKinCutsPt{"cKinCutsPt", {0.0, 0.4, 0.6, 0.8, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0, 4.0, 5.0, 6.0, 1e10}, "Omega(2012) pT bins for kinematic cuts"}; - Configurable> cKinLowerCutsAlpha{"cKinLowerCutsAlpha", {1.5, 1.0, 0.5, 0.3, 0.2, 0.15, 0.1, 0.08, 0.07, 0.06, 0.04, 0.02, 0.02}, "Lower cut on Xi-K0s opening angle"}; - Configurable> cKinUpperCutsAlpha{"cKinUpperCutsAlpha", {3.0, 2.0, 1.5, 1.4, 1.0, 0.8, 0.6, 0.5, 0.45, 0.35, 0.3, 0.25, 0.2}, "Upper cut on Xi-K0s opening angle"}; - // V0 selections (K0s) - Configurable cK0sMinCosPA{"cK0sMinCosPA", 0.98, "K0s minimum pointing angle cosine"}; - Configurable cK0sMaxDaughDCA{"cK0sMaxDaughDCA", 0.5, "K0s daughter DCA Maximum"}; - Configurable cK0sMassWindow{"cK0sMassWindow", 0.025, "Mass window for K0s selection (GeV/c^2)"}; - Configurable cMaxV0Etacut{"cMaxV0Etacut", 0.8, "V0 maximum eta cut"}; - - // Xi (cascade) selections from xi1530Analysisqa.cxx - Configurable cDCAxyToPVByPtCascP0{"cDCAxyToPVByPtCascP0", 999., "Cascade DCAxy p0"}; - Configurable cDCAxyToPVByPtCascExp{"cDCAxyToPVByPtCascExp", 1., "Cascade DCAxy exp"}; - Configurable cDCAxyToPVAsPtForCasc{"cDCAxyToPVAsPtForCasc", true, "Use pt-dep DCAxy cut (casc)"}; - - Configurable cDCAzToPVAsPtForCasc{"cDCAzToPVAsPtForCasc", true, "Use pt-dep DCAz cut (casc)"}; - - // V0 topology inside cascade (Λ) - Configurable cDCALambdaDaugtherscut{"cDCALambdaDaugtherscut", 0.7, "Λ daughters DCA cut"}; - Configurable cDCALambdaToPVcut{"cDCALambdaToPVcut", 0.02, "Λ DCA to PV min"}; - Configurable cDCAPionToPVcut{"cDCAPionToPVcut", 0.06, "π DCA to PV min"}; - Configurable cDCAProtonToPVcut{"cDCAProtonToPVcut", 0.07, "p DCA to PV min"}; - Configurable cV0CosPACutPtDepP0{"cV0CosPACutPtDepP0", 0.25, "V0 CosPA p0"}; - Configurable cV0CosPACutPtDepP1{"cV0CosPACutPtDepP1", 0.022, "V0 CosPA p1"}; - Configurable cMaxV0radiuscut{"cMaxV0radiuscut", 200., "V0 radius max"}; - Configurable cMinV0radiuscut{"cMinV0radiuscut", 2.5, "V0 radius min"}; - Configurable cMasswindowV0cut{"cMasswindowV0cut", 0.005, "Λ mass window for cascade V0"}; - - // Cascade topology - Configurable cDCABachlorToPVcut{"cDCABachlorToPVcut", 0.06, "Bachelor DCA to PV min"}; - Configurable cDCAXiDaugthersCutPtRangeLower{"cDCAXiDaugthersCutPtRangeLower", 1., "Xi pt low boundary"}; - Configurable cDCAXiDaugthersCutPtRangeUpper{"cDCAXiDaugthersCutPtRangeUpper", 4., "Xi pt high boundary"}; - Configurable cDCAXiDaugthersCutPtDepLower{"cDCAXiDaugthersCutPtDepLower", 0.8, "Xi daugh DCA (pt cDCAXiDaugthersCutPtDepMiddle{"cDCAXiDaugthersCutPtDepMiddle", 0.5, "Xi daugh DCA (low<=pt cDCAXiDaugthersCutPtDepUpper{"cDCAXiDaugthersCutPtDepUpper", 0.2, "Xi daugh DCA (pt>=high)"}; - Configurable cCosPACascCutPtDepP0{"cCosPACascCutPtDepP0", 0.2, "Cascade CosPA p0"}; - Configurable cCosPACascCutPtDepP1{"cCosPACascCutPtDepP1", 0.022, "Cascade CosPA p1"}; - Configurable cMaxCascradiuscut{"cMaxCascradiuscut", 200., "Cascade radius max"}; - Configurable cMinCascradiuscut{"cMinCascradiuscut", 1.1, "Cascade radius min"}; - Configurable cMasswindowCasccut{"cMasswindowCasccut", 0.008, "Xi mass window"}; - Configurable cMassXiminus{"cMassXiminus", 1.32171, "Xi mass (GeV/c^2)"}; // PDG + // Selection shared with the Omega(2012) training-table task (plain JSON keys; the mode-B track cuts carry "trk.") + o2::analysis::resonance::EventCuts eventCuts; + XiCuts xiCuts; + K0sCuts k0sCuts; + Xi1530KCuts xi1530KCuts; + o2::analysis::resonance::TrackCuts trackCuts = makeXi1530KTrackCuts(); + PionPidCuts pionPID; + KaonPidCuts kaonPID; + CandidateCuts candidateCuts; // Event Mixing Configurable nEvtMixing{"nEvtMixing", 10, "Number of events to mix"}; ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; ConfigurableAxis cfgMultBins{"cfgMultBins", {VARIABLE_WIDTH, 0.0f, 1.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "Mixing bins - centrality"}; - // Enhanced K0s selections - Configurable cK0sProperLifetimeMax{"cK0sProperLifetimeMax", 20.0, "K0s proper lifetime max (cm/c)"}; - Configurable cK0sArmenterosQtMin{"cK0sArmenterosQtMin", 0.0, "K0s Armenteros qt min"}; - Configurable cK0sArmenterosAlphaCoeff{"cK0sArmenterosAlphaCoeff", 0.2, "K0s Armenteros alpha coefficient"}; - Configurable cK0sDauPosDCAtoPVMin{"cK0sDauPosDCAtoPVMin", 0.05, "K0s positive daughter DCA to PV min"}; - Configurable cK0sDauNegDCAtoPVMin{"cK0sDauNegDCAtoPVMin", 0.05, "K0s negative daughter DCA to PV min"}; - Configurable cK0sRadiusMin{"cK0sRadiusMin", 0.5, "K0s decay radius min"}; - Configurable cK0sRadiusMax{"cK0sRadiusMax", 200.0, "K0s decay radius max"}; - Configurable cK0sCrossMassRejection{"cK0sCrossMassRejection", true, "Enable Lambda mass rejection for K0s"}; - Configurable cK0sCrossMassRejectionWindow{"cK0sCrossMassRejectionWindow", 0.01, "Lambda mass rejection window for K0s (GeV/c^2)"}; - Configurable cK0sDaughterPiTPCNSigmaMax{"cK0sDaughterPiTPCNSigmaMax", 5.0, "Maximum TPC NSigma for K0s daughter pions"}; - Configurable cK0sPosDaughterMinCrossedRows{"cK0sPosDaughterMinCrossedRows", 50, "Minimum TPC crossed rows for K0s positive daughter"}; - Configurable cK0sNegDaughterMinCrossedRows{"cK0sNegDaughterMinCrossedRows", 50, "Minimum TPC crossed rows for K0s negative daughter"}; - - // Pion track selections for 3-body decay - Configurable cPionPtMin{"cPionPtMin", 0.15, "Minimum pion pT"}; - Configurable cPionEtaMax{"cPionEtaMax", 0.8, "Maximum pion |eta|"}; - Configurable cPionDCAxyMax{"cPionDCAxyMax", 0.1, "Maximum pion DCAxy to PV"}; - Configurable cPionDCAzMax{"cPionDCAzMax", 0.2, "Maximum pion DCAz to PV"}; - Configurable cPionTPCNClusMin{"cPionTPCNClusMin", 70, "Minimum TPC clusters for pion"}; - - // Pion PID selections - Configurable cPionTPCNSigmaMax{"cPionTPCNSigmaMax", 3.0, "Maximum TPC NSigma for pion"}; - Configurable cPionTOFNSigmaMax{"cPionTOFNSigmaMax", 3.0, "Maximum TOF NSigma for pion"}; - Configurable cPionUsePtDepPID{"cPionUsePtDepPID", false, "Use pT-dependent PID cuts for pion"}; - Configurable> cPionPIDPtBins{"cPionPIDPtBins", {0.0f, 0.5f, 0.8f, 2.0f, 999.0f}, "pT bin edges for pion PID cuts"}; - Configurable> cPionTPCNSigmaCuts{"cPionTPCNSigmaCuts", {3.0f, 3.0f, 2.0f, 2.0f}, "TPC NSigma cuts per pT bin (pion)"}; - Configurable> cPionTOFNSigmaCuts{"cPionTOFNSigmaCuts", {3.0f, 3.0f, 3.0f, 3.0f}, "TOF NSigma cuts per pT bin (pion)"}; - Configurable> cPionTOFRequired{"cPionTOFRequired", {0, 0, 1, 1}, "Require TOF per pT bin (pion)"}; - - // Xi1530 mass window cut - Configurable cXi1530Mass{"cXi1530Mass", 1.53, "Xi(1530) mass (GeV/c^2)"}; - Configurable cXi1530MassWindow{"cXi1530MassWindow", 0.01, "Xi(1530) mass window (GeV/c^2)"}; - - // PDG masses - double massK0 = MassK0Short; - - // Module-initializer collision tables (matches the resonance-module-initializer producer) + // Module-initializer collision and daughter tables (matches the resonance-module-initializer producer) using ResoCollisions = aod::ResoCollisions_001; using ResoMCCollisions = soa::Join; + using ResoMicroTracks = aod::ResoMicroTracks_001; + using ResoMCCascades = soa::Join; + using ResoMCV0s = soa::Join; + using ResoMCMicroTracks = soa::Join; using BinningTypeVertexContributor = ColumnBinningPolicy; BinningTypeVertexContributor colBinning{{cfgVtxBins, cfgMultBins}}; + Omega2012AnalysisCore core; + void init(InitContext&) { AxisSpec centAxis = {binsCent, "V0M (%)"}; @@ -261,27 +208,6 @@ struct Omega2012Analysis { histos.add("QAbefore/omegaAlphaVsPt", "#alpha_{oa} vs p_{T} before kinematic cuts", kTH2F, {omegaKinPtAxis, openingAngleAxis}); histos.add("QAafter/omegaAlphaVsPt", "#alpha_{oa} vs p_{T} after kinematic cuts", kTH2F, {omegaKinPtAxis, openingAngleAxis}); - // 3-body decay: Xi + pi + K0s - histos.add("omega2012_3body/invmass", "Invariant mass of Omega(2012) → Xi + #pi + K^{0}_{S}", kTH1F, {invMassAxis}); - histos.add("omega2012_3body/massPtCent", "Omega(2012) 3-body mass vs pT vs cent", kTH3F, {invMassAxis, ptAxis, centAxis}); - - // Pion QA histograms for 3-body - histos.add("QAbefore/pionPt", "Pion pT before cuts", kTH1F, {ptAxisQA}); - histos.add("QAbefore/pionEta", "Pion eta before cuts", kTH1F, {{100, -2.0, 2.0, "#eta"}}); - histos.add("QAbefore/pionDCAxy", "Pion DCAxy before cuts", kTH2F, {ptAxisQA, dcaxyAxis}); - histos.add("QAbefore/pionDCAz", "Pion DCAz before cuts", kTH2F, {ptAxisQA, dcazAxis}); - histos.add("QAbefore/pionTPCNcls", "Pion TPC clusters before cuts", kTH1F, {{160, 0, 160, "N_{TPC clusters}"}}); - histos.add("QAbefore/pionTPCNSigma", "Pion TPC NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); - histos.add("QAbefore/pionTOFNSigma", "Pion TOF NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); - - histos.add("QAafter/pionPt", "Pion pT after cuts", kTH1F, {ptAxisQA}); - histos.add("QAafter/pionEta", "Pion eta after cuts", kTH1F, {{100, -2.0, 2.0, "#eta"}}); - histos.add("QAafter/pionDCAxy", "Pion DCAxy after cuts", kTH2F, {ptAxisQA, dcaxyAxis}); - histos.add("QAafter/pionDCAz", "Pion DCAz after cuts", kTH2F, {ptAxisQA, dcazAxis}); - histos.add("QAafter/pionTPCNcls", "Pion TPC clusters after cuts", kTH1F, {{160, 0, 160, "N_{TPC clusters}"}}); - histos.add("QAafter/pionTPCNSigma", "Pion TPC NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); - histos.add("QAafter/pionTOFNSigma", "Pion TOF NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); - // MC truth histograms AxisSpec etaAxis = {100, -2.0, 2.0, "#eta"}; AxisSpec rapidityAxis = {100, -2.0, 2.0, "y"}; @@ -301,590 +227,177 @@ struct Omega2012Analysis { histos.add("MC/hMCRecK0sPt", "MC Reconstructed K0s pT", kTH1F, {ptAxis}); histos.add("MC/hMCTrueXiPt", "MC True Xi pT", kTH1F, {ptAxis}); histos.add("MC/hMCTrueK0sPt", "MC True K0s pT", kTH1F, {ptAxis}); - } - - template - static std::array cascadeDaughterIds(const CascT& xi) - { - auto indices = xi.cascadeIndices(); - return {indices[0], indices[1], indices[2]}; - } - template - static std::array v0DaughterIds(const V0Type& v0) - { - auto indices = v0.indices(); - return {indices[0], indices[1]}; - } - - template - static bool sharesAnyDaughterId(const std::array& first, const std::array& second) - { - for (const auto& firstId : first) { - for (const auto& secondId : second) { - if (firstId == secondId) { - return true; - } + ProcessModes modes; + modes.xiK0s = doprocessData || doprocessMixedEvent || doprocessMC; + modes.xi1530K = doprocessXi1530KMicro || doprocessXi1530KTracks || doprocessXi1530KMCMicro || doprocessXi1530KMixedMicro; + modes.microTracks = doprocessXi1530KMicro || doprocessXi1530KMCMicro || doprocessXi1530KMixedMicro; + modes.mcReco = doprocessMC || doprocessXi1530KMCMicro; + modes.mixing = doprocessMixedEvent || doprocessXi1530KMixedMicro; + core.init(histos, eventCuts, xiCuts, k0sCuts, xi1530KCuts, trackCuts, pionPID, kaonPID, candidateCuts, modes); + + if (modes.xi1530K) { + AxisSpec xiPiMassAxis = {300, 1.4, 1.7, "M_{#Xi#pi} (GeV/#it{c}^{2})"}; + AxisSpec patternAxis = {3, 0.5, 3.5, "charge pattern (1: wrong-sign #pi, 2: wrong-sign K, 3: both)"}; + AxisSpec etaAxisQA = {100, -2.0, 2.0, "#eta"}; + + // Signal charge pattern: Xi- pi+ K- and Xi+ pi- K+ + histos.add("xi1530K/invmass", "Invariant mass of Omega(2012) → #Xi(1530)^{0} K → #Xi #pi K", kTH1F, {invMassAxis}); + histos.add("xi1530K/massPtCent", "Omega(2012) → #Xi(1530)^{0} K mass vs pT vs cent", kTH3F, {invMassAxis, ptAxis, centAxis}); + histos.add("xi1530K/invmass_Mix", "Mixed event invariant mass of Omega(2012) → #Xi(1530)^{0} K", kTH1F, {invMassAxis}); + histos.add("xi1530K/massPtCent_Mix", "Mixed event Omega(2012) → #Xi(1530)^{0} K mass vs pT vs cent", kTH3F, {invMassAxis, ptAxis, centAxis}); + histos.add("xi1530K/massXiPi", "#Xi #pi mass of selected candidates (before the rapidity cut)", kTH1F, {xiPiMassAxis}); + histos.add("xi1530K/massXiPiVsMass", "#Xi #pi mass vs #Xi #pi K mass", kTH2F, {invMassAxis, xiPiMassAxis}); + + // QA of the mode-B inputs (once per object and collision) + histos.add("xi1530K/QAbefore/xiMass", "Xi mass before cuts", kTH1F, {xiMassAxis}); + histos.add("xi1530K/QAbefore/xiPt", "Xi pT before cuts", kTH1F, {ptAxisQA}); + histos.add("xi1530K/QAbefore/xiEta", "Xi eta before cuts", kTH1F, {etaAxisQA}); + histos.add("xi1530K/QAafter/xiMass", "Xi mass after cuts", kTH1F, {xiMassAxis}); + histos.add("xi1530K/QAafter/xiPt", "Xi pT after cuts", kTH1F, {ptAxisQA}); + histos.add("xi1530K/QAafter/xiEta", "Xi eta after cuts", kTH1F, {etaAxisQA}); + histos.add("xi1530K/QAbefore/trackPt", "Track pT before cuts", kTH1F, {ptAxisQA}); + histos.add("xi1530K/QAbefore/trackEta", "Track eta before cuts", kTH1F, {etaAxisQA}); + histos.add("xi1530K/QAbefore/trackDCAxy", "Track DCAxy before cuts", kTH2F, {ptAxisQA, dcaxyAxis}); + histos.add("xi1530K/QAbefore/trackDCAz", "Track DCAz before cuts", kTH2F, {ptAxisQA, dcazAxis}); + histos.add("xi1530K/QAbefore/pionTPCNSigma", "Pion TPC NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAbefore/pionTOFNSigma", "Pion TOF NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAbefore/kaonTPCNSigma", "Kaon TPC NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAbefore/kaonTOFNSigma", "Kaon TOF NSigma before cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAafter/pionPt", "Pion pT after cuts", kTH1F, {ptAxisQA}); + histos.add("xi1530K/QAafter/pionEta", "Pion eta after cuts", kTH1F, {etaAxisQA}); + histos.add("xi1530K/QAafter/pionDCAxy", "Pion DCAxy after cuts", kTH2F, {ptAxisQA, dcaxyAxis}); + histos.add("xi1530K/QAafter/pionDCAz", "Pion DCAz after cuts", kTH2F, {ptAxisQA, dcazAxis}); + histos.add("xi1530K/QAafter/pionTPCNSigma", "Pion TPC NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAafter/pionTOFNSigma", "Pion TOF NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAafter/kaonPt", "Kaon pT after cuts", kTH1F, {ptAxisQA}); + histos.add("xi1530K/QAafter/kaonEta", "Kaon eta after cuts", kTH1F, {etaAxisQA}); + histos.add("xi1530K/QAafter/kaonDCAxy", "Kaon DCAxy after cuts", kTH2F, {ptAxisQA, dcaxyAxis}); + histos.add("xi1530K/QAafter/kaonDCAz", "Kaon DCAz after cuts", kTH2F, {ptAxisQA, dcazAxis}); + histos.add("xi1530K/QAafter/kaonTPCNSigma", "Kaon TPC NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + histos.add("xi1530K/QAafter/kaonTOFNSigma", "Kaon TOF NSigma after cuts", kTH2F, {ptAxisQA, nsigmaAxis}); + + // Charge-pattern controls (never part of the signal) + histos.add("xi1530K_wrongSign/invmassPattern", "Wrong-sign #Xi #pi K mass by charge pattern", kTH2F, {patternAxis, invMassAxis}); + histos.add("xi1530K_wrongSign/massPtPattern", "Wrong-sign #Xi #pi K mass vs pT by charge pattern", kTH3F, {invMassAxis, ptAxis, patternAxis}); + + if (doprocessXi1530KMCMicro) { + histos.add("xi1530K/MC/hMCRecOmega2012Pt", "MC reconstructed Omega(2012) → #Xi(1530)^{0} K pT", kTH1F, {ptAxis}); + histos.add("xi1530K/MC/hMCRecOmega2012PtEta", "MC reconstructed Omega(2012) → #Xi(1530)^{0} K pT vs eta", kTH2F, {ptAxis, etaAxis}); + histos.add("xi1530K/MC/hMCRecMass", "MC reconstructed Omega(2012) → #Xi(1530)^{0} K mass", kTH1F, {invMassAxis}); + histos.add("xi1530K/MC/hMCRecMassXiPi", "MC reconstructed #Xi #pi mass", kTH1F, {xiPiMassAxis}); + histos.add("xi1530K/MC/hMCRecXiPt", "MC reconstructed Xi pT", kTH1F, {ptAxis}); + histos.add("xi1530K/MC/hMCRecPionPt", "MC reconstructed pion pT", kTH1F, {ptAxis}); + histos.add("xi1530K/MC/hMCRecKaonPt", "MC reconstructed kaon pT", kTH1F, {ptAxis}); } } - return false; - } - template - static bool sharesDaughterId(const std::array& daughters, int trackId) - { - for (const auto& daughterId : daughters) { - if (daughterId == trackId) { - return true; - } - } - return false; + LOG(info) << "Size of the histograms in Omega(2012) analysis task"; + histos.print(); } - template - static int trackSourceId(const TrackT& track, const TrackIdsT& trackIds) - { - auto rowIndex = track.globalIndex(); - if (rowIndex >= 0 && rowIndex < trackIds.size()) { - return trackIds.rawIteratorAt(rowIndex).trackId(); - } - return rowIndex; - } - - template - bool kinCuts(const FirstVecT& firstDaughter, const SecondVecT& secondDaughter, const MotherVecT& mother, float& alpha) - { - auto firstP = std::sqrt(firstDaughter.Px() * firstDaughter.Px() + firstDaughter.Py() * firstDaughter.Py() + firstDaughter.Pz() * firstDaughter.Pz()); - auto secondP = std::sqrt(secondDaughter.Px() * secondDaughter.Px() + secondDaughter.Py() * secondDaughter.Py() + secondDaughter.Pz() * secondDaughter.Pz()); - if (firstP < kSmallNumber || secondP < kSmallNumber) { - alpha = 0.f; - return false; - } - - auto cosAlpha = (firstDaughter.Px() * secondDaughter.Px() + firstDaughter.Py() * secondDaughter.Py() + firstDaughter.Pz() * secondDaughter.Pz()) / (firstP * secondP); - if (cosAlpha > 1.) { - cosAlpha = 1.; - } else if (cosAlpha < -1.) { - cosAlpha = -1.; - } - alpha = std::acos(cosAlpha); - - std::vector kinCutsPt = static_cast>(cKinCutsPt); - std::vector kinLowerCutsAlpha = static_cast>(cKinLowerCutsAlpha); - std::vector kinUpperCutsAlpha = static_cast>(cKinUpperCutsAlpha); - - int kinCutsSize = static_cast(kinUpperCutsAlpha.size()); - if (kinCutsSize > static_cast(kinLowerCutsAlpha.size())) { - kinCutsSize = static_cast(kinLowerCutsAlpha.size()); - } - if (kinCutsSize > static_cast(kinCutsPt.size()) - 1) { - kinCutsSize = static_cast(kinCutsPt.size()) - 1; - } - - for (int i = 0; i < kinCutsSize; ++i) { - if ((mother.Pt() > kinCutsPt[i] && mother.Pt() <= kinCutsPt[i + 1]) && (alpha < kinLowerCutsAlpha[i] || alpha > kinUpperCutsAlpha[i])) { - return false; - } - } - - return true; - } - - // Enhanced V0 selection (K0s) with detailed criteria - template - bool v0CutEnhanced(const CollisionType& collision, const V0Type& v0) - { - // Basic kinematic cuts - if (std::abs(v0.eta()) > cMaxV0Etacut) - return false; - if (v0.pt() < cMinPtcut) - return false; - - // Topological cuts - if (v0.v0CosPA() < cK0sMinCosPA) - return false; - if (v0.daughDCA() > cK0sMaxDaughDCA) - return false; - - // Enhanced selections from chk892Flow - // Daughter DCA to PV cuts - if (std::abs(v0.dcapostopv()) < cK0sDauPosDCAtoPVMin) - return false; - if (std::abs(v0.dcanegtopv()) < cK0sDauNegDCAtoPVMin) - return false; - - // Radius cuts - use transRadius instead of v0radius - auto radius = v0.transRadius(); - if (radius < cK0sRadiusMin || radius > cK0sRadiusMax) - return false; - - // DCA to PV - if (std::abs(v0.dcav0topv()) > kMaxDCAV0ToPV) - return false; // max DCA to PV - - // Proper lifetime cut - calculate manually - float dx = v0.decayVtxX() - collision.posX(); - float dy = v0.decayVtxY() - collision.posY(); - float dz = v0.decayVtxZ() - collision.posZ(); - float l = std::sqrt(dx * dx + dy * dy + dz * dz); - float p = std::sqrt(v0.px() * v0.px() + v0.py() * v0.py() + v0.pz() * v0.pz()); - auto properLifetime = (l / (p + kSmallNumber)) * MassK0Short; - if (properLifetime > cK0sProperLifetimeMax) - return false; - - if (v0.qtarm() < cK0sArmenterosQtMin) - return false; - - // Mass window - if (std::abs(v0.mK0Short() - MassK0Short) > cK0sMassWindow) - return false; - - // Competing V0 rejection: remove (Anti)Λ - if (cK0sCrossMassRejection) { - if (std::abs(v0.mLambda() - MassLambda) < cK0sCrossMassRejectionWindow) - return false; - if (std::abs(v0.mAntiLambda() - MassLambda) < cK0sCrossMassRejectionWindow) - return false; - } - - if (std::abs(v0.daughterTPCNSigmaPosPi()) >= cK0sDaughterPiTPCNSigmaMax) - return false; - if (std::abs(v0.daughterTPCNSigmaNegPi()) >= cK0sDaughterPiTPCNSigmaMax) - return false; - - if (v0.nCrossedRowsPos() <= cK0sPosDaughterMinCrossedRows) - return false; - if (v0.nCrossedRowsNeg() <= cK0sNegDaughterMinCrossedRows) - return false; - - if (v0.qtarm() < cK0sArmenterosAlphaCoeff * std::fabs(v0.alpha())) - return false; - - return true; - } - - // Helper function to find pT bin index - int getPtBinIndex(float pt) - { - auto ptBins = static_cast>(cPionPIDPtBins); - for (size_t i = 0; i < ptBins.size() - 1; i++) { - if (pt >= ptBins[i] && pt < ptBins[i + 1]) { - return i; - } - } - return -1; - } - - // Pion PID selection - template - bool pionPidCut(const TrackType& track) - { - float pt = track.pt(); - - if constexpr (IsResoMicrotrack) { - // For ResoMicroTracks - decode PID from flags - float tpcNSigma = o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(track.pidNSigmaPiFlag()); - float tofNSigma = track.hasTOF() ? o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(track.pidNSigmaPiFlag()) : 999.f; - - if (cPionUsePtDepPID) { - int ptBin = getPtBinIndex(pt); - if (ptBin < 0) - return false; - - auto tpcCuts = static_cast>(cPionTPCNSigmaCuts); - auto tofCuts = static_cast>(cPionTOFNSigmaCuts); - auto tofRequired = static_cast>(cPionTOFRequired); - - if (ptBin >= static_cast(tpcCuts.size()) || - ptBin >= static_cast(tofCuts.size()) || - ptBin >= static_cast(tofRequired.size())) { - return false; - } - - if (std::abs(tpcNSigma) >= tpcCuts[ptBin]) - return false; - - if (tofRequired[ptBin] != 0) { - if (!track.hasTOF()) - return false; - if (std::abs(tofNSigma) >= tofCuts[ptBin]) - return false; - } else { - if (track.hasTOF() && std::abs(tofNSigma) >= tofCuts[ptBin]) - return false; - } - - return true; - } else { - bool tpcPass = std::abs(tpcNSigma) < cPionTPCNSigmaMax; - bool tofPass = track.hasTOF() ? std::abs(tofNSigma) < cPionTOFNSigmaMax : true; - return tpcPass && tofPass; - } - } else { - // For ResoTracks - direct access - float tpcNSigma = track.tpcNSigmaPi(); - float tofNSigma = track.hasTOF() ? track.tofNSigmaPi() : 999.f; - - if (cPionUsePtDepPID) { - int ptBin = getPtBinIndex(pt); - if (ptBin < 0) - return false; - - auto tpcCuts = static_cast>(cPionTPCNSigmaCuts); - auto tofCuts = static_cast>(cPionTOFNSigmaCuts); - auto tofRequired = static_cast>(cPionTOFRequired); - - if (ptBin >= static_cast(tpcCuts.size()) || - ptBin >= static_cast(tofCuts.size()) || - ptBin >= static_cast(tofRequired.size())) { - return false; - } - - if (std::abs(tpcNSigma) >= tpcCuts[ptBin]) - return false; - - if (tofRequired[ptBin] != 0) { - if (!track.hasTOF()) - return false; - if (std::abs(tofNSigma) >= tofCuts[ptBin]) - return false; - } else { - if (track.hasTOF() && std::abs(tofNSigma) >= tofCuts[ptBin]) - return false; - } - - return true; - } else { - bool tpcPass = std::abs(tpcNSigma) < cPionTPCNSigmaMax; - bool tofPass = track.hasTOF() ? std::abs(tofNSigma) < cPionTOFNSigmaMax : true; - return tpcPass && tofPass; - } - } - } - - // Pion track selection (for both ResoTracks and ResoMicroTracks) - template - bool pionCut(const TrackType& track) - { - // Basic kinematic cuts - if (track.pt() < cPionPtMin) - return false; - if (std::abs(track.eta()) > cPionEtaMax) - return false; - - // DCA cuts - different access for ResoMicroTracks - if constexpr (IsResoMicrotrack) { - if (o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(track.trackSelectionFlags()) > cPionDCAxyMax) - return false; - if (o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(track.trackSelectionFlags()) > cPionDCAzMax) - return false; - } else { - if (std::abs(track.dcaXY()) > cPionDCAxyMax) - return false; - if (std::abs(track.dcaZ()) > cPionDCAzMax) - return false; - } - - // Track quality cuts - only for ResoTracks - if constexpr (!IsResoMicrotrack) { - if constexpr (requires { track.tpcNClsFound(); }) { - if (track.tpcNClsFound() < cPionTPCNClusMin) - return false; - } - } - - // PID selection - if (!pionPidCut(track)) - return false; - - return true; - } - - // Xi1530 mass window cut - template - bool xi1530MassCut(const XiType& xi, const PionType& pion) - { - // Calculate Xi + pion invariant mass - ROOT::Math::PxPyPzEVector pXi, pPion, pXi1530; - pXi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); - pPion = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(pion.pt(), pion.eta(), pion.phi(), MassPionCharged)); - pXi1530 = pXi + pPion; - - // Check if mass is within Xi(1530) window - float massDiff = std::abs(pXi1530.M() - cXi1530Mass); - return massDiff < cXi1530MassWindow; - } - - // Primary-level cascade kinematics - template - bool cascprimaryTrackCut(const CascT& c) - { - if (std::abs(c.eta()) > cMaxEtaCut) - return false; - if (std::abs(c.pt()) < cMinPtcut) - return false; - if (cDCAxyToPVAsPtForCasc) { - if (std::abs(c.dcaXYCascToPV()) > (cDCAxyToPVByPtCascP0 + cDCAxyToPVByPtCascExp * c.pt())) - return false; - } - if (cDCAzToPVAsPtForCasc) { - if (std::abs(c.dcaZCascToPV()) > (cDCAxyToPVByPtCascP0 + cDCAxyToPVByPtCascExp * std::pow(c.pt(), -1.1f))) - return false; - } - return true; - } - - // Cascade topological selections adapted from xi1530Analysisqa - template - bool casctopCut(const CascT& c) - { - // V0 (Λ) topology inside cascade - if (std::abs(c.daughDCA()) > cDCALambdaDaugtherscut) - return false; - if (std::abs(c.dcav0topv()) < cDCALambdaToPVcut) - return false; - - if (c.sign() < 0) { // Xi- - if (std::abs(c.dcanegtopv()) < cDCAPionToPVcut) - return false; - if (std::abs(c.dcapostopv()) < cDCAProtonToPVcut) - return false; - } else { // Anti-Xi - if (std::abs(c.dcanegtopv()) < cDCAProtonToPVcut) - return false; - if (std::abs(c.dcapostopv()) < cDCAPionToPVcut) - return false; - } - - if (c.v0CosPA() < std::cos(cV0CosPACutPtDepP0 - cV0CosPACutPtDepP1 * c.pt())) - return false; - if (c.transRadius() > cMaxV0radiuscut || c.transRadius() < cMinV0radiuscut) - return false; - if (std::abs(c.mLambda() - MassLambda) > cMasswindowV0cut) - return false; - - // Cascade-level topology - if (std::abs(c.dcabachtopv()) < cDCABachlorToPVcut) - return false; - - if (c.pt() < cDCAXiDaugthersCutPtRangeLower) { - if (c.cascDaughDCA() > cDCAXiDaugthersCutPtDepLower) - return false; - } else if (c.pt() < cDCAXiDaugthersCutPtRangeUpper) { - if (c.cascDaughDCA() > cDCAXiDaugthersCutPtDepMiddle) - return false; - } else { - if (c.cascDaughDCA() > cDCAXiDaugthersCutPtDepUpper) - return false; - } - - if (c.cascCosPA() < std::cos(cCosPACascCutPtDepP0 - cCosPACascCutPtDepP1 * c.pt())) - return false; - if (c.cascTransRadius() > cMaxCascradiuscut || c.cascTransRadius() < cMinCascradiuscut) - return false; - if (std::abs(c.mXi() - cMassXiminus) > cMasswindowCasccut) - return false; - - return true; - } - - template - struct SelectedXiCandidate { - CandidateT candidate; - std::array daughterIds{}; - }; - - template - struct SelectedK0sCandidate { - CandidateT candidate; - std::array daughterIds{}; - }; - - template - auto selectK0sCandidates(const CollisionT& collision, const V0sT& v0s, int& nV0sAfterCuts) - { - using V0Candidate = std::decay_t; - std::vector> selectedK0s; - selectedK0s.reserve(v0s.size()); - - for (const auto& v0 : v0s) { - float dx = v0.decayVtxX() - collision.posX(); - float dy = v0.decayVtxY() - collision.posY(); - float dz = v0.decayVtxZ() - collision.posZ(); - float l = std::sqrt(dx * dx + dy * dy + dz * dz); - float p = std::sqrt(v0.px() * v0.px() + v0.py() * v0.py() + v0.pz() * v0.pz()); - auto properLifetime = (l / (p + kSmallNumber)) * MassK0Short; - - if constexpr (FillQA) { - histos.fill(HIST("QAbefore/k0sMassPt"), v0.pt(), v0.mK0Short()); - histos.fill(HIST("QAbefore/k0sPt"), v0.pt()); - histos.fill(HIST("QAbefore/k0sEta"), v0.eta()); - histos.fill(HIST("QAbefore/k0sCosPA"), v0.pt(), v0.v0CosPA()); - histos.fill(HIST("QAbefore/k0sRadius"), v0.pt(), v0.transRadius()); - histos.fill(HIST("QAbefore/k0sDauDCA"), v0.pt(), v0.daughDCA()); - histos.fill(HIST("QAbefore/k0sDCAtoPV"), v0.pt(), std::abs(v0.dcav0topv())); - histos.fill(HIST("QAbefore/k0sProperLifetime"), v0.pt(), properLifetime); - histos.fill(HIST("QAbefore/k0sArmenteros"), v0.alpha(), v0.qtarm()); - histos.fill(HIST("QAbefore/k0sDauPosDCA"), v0.pt(), std::abs(v0.dcapostopv())); - histos.fill(HIST("QAbefore/k0sDauNegDCA"), v0.pt(), std::abs(v0.dcanegtopv())); - histos.fill(HIST("QAbefore/k0sDauTPCNsigmaPosPi"), v0.pt(), v0.daughterTPCNSigmaPosPi()); - histos.fill(HIST("QAbefore/k0sDauTPCNsigmaNegPi"), v0.pt(), v0.daughterTPCNSigmaNegPi()); - histos.fill(HIST("QAbefore/k0sNCrossedRowsPos"), v0.pt(), v0.nCrossedRowsPos()); - histos.fill(HIST("QAbefore/k0sNCrossedRowsNeg"), v0.pt(), v0.nCrossedRowsNeg()); - } - - if (!v0CutEnhanced(collision, v0)) - continue; - - if constexpr (FillQA) { - nV0sAfterCuts++; - histos.fill(HIST("QAafter/k0sMassPt"), v0.pt(), v0.mK0Short()); - histos.fill(HIST("QAafter/k0sPt"), v0.pt()); - histos.fill(HIST("QAafter/k0sEta"), v0.eta()); - histos.fill(HIST("QAafter/k0sCosPA"), v0.pt(), v0.v0CosPA()); - histos.fill(HIST("QAafter/k0sRadius"), v0.pt(), v0.transRadius()); - histos.fill(HIST("QAafter/k0sDauDCA"), v0.pt(), v0.daughDCA()); - histos.fill(HIST("QAafter/k0sDCAtoPV"), v0.pt(), std::abs(v0.dcav0topv())); - histos.fill(HIST("QAafter/k0sProperLifetime"), v0.pt(), properLifetime); - histos.fill(HIST("QAafter/k0sArmenteros"), v0.alpha(), v0.qtarm()); - histos.fill(HIST("QAafter/k0sDauPosDCA"), v0.pt(), std::abs(v0.dcapostopv())); - histos.fill(HIST("QAafter/k0sDauNegDCA"), v0.pt(), std::abs(v0.dcanegtopv())); - histos.fill(HIST("QAafter/k0sDauTPCNsigmaPosPi"), v0.pt(), v0.daughterTPCNSigmaPosPi()); - histos.fill(HIST("QAafter/k0sDauTPCNsigmaNegPi"), v0.pt(), v0.daughterTPCNSigmaNegPi()); - histos.fill(HIST("QAafter/k0sNCrossedRowsPos"), v0.pt(), v0.nCrossedRowsPos()); - histos.fill(HIST("QAafter/k0sNCrossedRowsNeg"), v0.pt(), v0.nCrossedRowsNeg()); - } - - selectedK0s.push_back({v0, v0DaughterIds(v0)}); - } - - return selectedK0s; - } - - template - auto selectXiCandidates(const CascadesT& cascades, int& nCascAfterCuts) - { - using XiCandidate = std::decay_t; - std::vector> selectedXis; - selectedXis.reserve(cascades.size()); - - for (const auto& xi : cascades) { - if constexpr (FillQA) { - histos.fill(HIST("QAbefore/xiMass"), xi.mXi()); - histos.fill(HIST("QAbefore/xiPt"), xi.pt()); - histos.fill(HIST("QAbefore/xiEta"), xi.eta()); - histos.fill(HIST("QAbefore/xiDCAxy"), xi.pt(), xi.dcaXYCascToPV()); - histos.fill(HIST("QAbefore/xiDCAz"), xi.pt(), xi.dcaZCascToPV()); - histos.fill(HIST("QAbefore/xiV0CosPA"), xi.pt(), xi.v0CosPA()); - histos.fill(HIST("QAbefore/xiCascCosPA"), xi.pt(), xi.cascCosPA()); - histos.fill(HIST("QAbefore/xiV0Radius"), xi.pt(), xi.transRadius()); - histos.fill(HIST("QAbefore/xiCascRadius"), xi.pt(), xi.cascTransRadius()); - histos.fill(HIST("QAbefore/xiV0DauDCA"), xi.pt(), xi.daughDCA()); - histos.fill(HIST("QAbefore/xiCascDauDCA"), xi.pt(), xi.cascDaughDCA()); - } - - if (!cascprimaryTrackCut(xi)) - continue; - if (!casctopCut(xi)) - continue; - - if constexpr (FillQA) { - nCascAfterCuts++; - histos.fill(HIST("QAafter/xiMass"), xi.mXi()); - histos.fill(HIST("QAafter/xiPt"), xi.pt()); - histos.fill(HIST("QAafter/xiEta"), xi.eta()); - histos.fill(HIST("QAafter/xiDCAxy"), xi.pt(), xi.dcaXYCascToPV()); - histos.fill(HIST("QAafter/xiDCAz"), xi.pt(), xi.dcaZCascToPV()); - histos.fill(HIST("QAafter/xiV0CosPA"), xi.pt(), xi.v0CosPA()); - histos.fill(HIST("QAafter/xiCascCosPA"), xi.pt(), xi.cascCosPA()); - histos.fill(HIST("QAafter/xiV0Radius"), xi.pt(), xi.transRadius()); - histos.fill(HIST("QAafter/xiCascRadius"), xi.pt(), xi.cascTransRadius()); - histos.fill(HIST("QAafter/xiV0DauDCA"), xi.pt(), xi.daughDCA()); - histos.fill(HIST("QAafter/xiCascDauDCA"), xi.pt(), xi.cascDaughDCA()); - } - - selectedXis.push_back({xi, cascadeDaughterIds(xi)}); - } - - return selectedXis; - } - - template - void fill(const CollisionT& collision, const CascadesT& cascades, const V0sT& v0s) + // Mode A, same event (data): event QA, Xi and K0S QA, kinematic-cut QA and the Xi K0S mass + template + void fillXiK0s(const CollisionT& collision, const CascadesT& cascades, const V0sT& v0s) { auto cent = collision.cent(); - // Fill event QA histograms (only for same-event) - if constexpr (!IsMix) { - histos.fill(HIST("Event/posZ"), collision.posZ()); - histos.fill(HIST("Event/centrality"), cent); - histos.fill(HIST("Event/posZvsCent"), collision.posZ(), cent); - histos.fill(HIST("Event/nCascades"), cascades.size()); - histos.fill(HIST("Event/nV0s"), v0s.size()); - } + histos.fill(HIST("Event/posZ"), collision.posZ()); + histos.fill(HIST("Event/centrality"), cent); + histos.fill(HIST("Event/posZvsCent"), collision.posZ(), cent); + histos.fill(HIST("Event/nCascades"), cascades.size()); + histos.fill(HIST("Event/nV0s"), v0s.size()); // Count candidates after cuts int nCascAfterCuts = 0; int nV0sAfterCuts = 0; - auto selectedK0s = selectK0sCandidates(collision, v0s, nV0sAfterCuts); - auto selectedXis = selectXiCandidates(cascades, nCascAfterCuts); - - for (const auto& selectedXi : selectedXis) { - const auto& xi = selectedXi.candidate; - - // Build Xi + K0s - for (const auto& selectedK0 : selectedK0s) { - const auto& v0 = selectedK0.candidate; - - if constexpr (!IsMix) { - if (sharesAnyDaughterId(selectedXi.daughterIds, selectedK0.daughterIds)) { - continue; - } - } - - // 4-vectors - ROOT::Math::PxPyPzEVector pXi, pK0s, pRes; - pXi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); - pK0s = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(v0.pt(), v0.eta(), v0.phi(), massK0)); - pRes = pXi + pK0s; - - float alpha = 0.f; - bool kinCutFlag = true; - if (cKinCuts) { - kinCutFlag = kinCuts(pXi, pK0s, pRes, alpha); - if constexpr (!IsMix) { - histos.fill(HIST("QAbefore/omegaAlphaVsPt"), pRes.Pt(), alpha); - } - } - - if (cKinCuts && !kinCutFlag) - continue; - - if constexpr (!IsMix) { - if (cKinCuts) { - histos.fill(HIST("QAafter/omegaAlphaVsPt"), pRes.Pt(), alpha); - } - } - - if (std::abs(pRes.Rapidity()) >= cfgRapidityCut) - continue; - - if constexpr (!IsMix) { - histos.fill(HIST("omega2012/invmass"), pRes.M()); - histos.fill(HIST("omega2012/massPtCent"), pRes.M(), pRes.Pt(), cent); - } else { - histos.fill(HIST("omega2012/invmass_Mix"), pRes.M()); - histos.fill(HIST("omega2012/massPtCent_Mix"), pRes.M(), pRes.Pt(), cent); + auto onK0s = [&](auto const& v0, double properLifetime, bool selected) { + histos.fill(HIST("QAbefore/k0sMassPt"), v0.pt(), v0.mK0Short()); + histos.fill(HIST("QAbefore/k0sPt"), v0.pt()); + histos.fill(HIST("QAbefore/k0sEta"), v0.eta()); + histos.fill(HIST("QAbefore/k0sCosPA"), v0.pt(), v0.v0CosPA()); + histos.fill(HIST("QAbefore/k0sRadius"), v0.pt(), v0.transRadius()); + histos.fill(HIST("QAbefore/k0sDauDCA"), v0.pt(), v0.daughDCA()); + histos.fill(HIST("QAbefore/k0sDCAtoPV"), v0.pt(), std::abs(v0.dcav0topv())); + histos.fill(HIST("QAbefore/k0sProperLifetime"), v0.pt(), properLifetime); + histos.fill(HIST("QAbefore/k0sArmenteros"), v0.alpha(), v0.qtarm()); + histos.fill(HIST("QAbefore/k0sDauPosDCA"), v0.pt(), std::abs(v0.dcapostopv())); + histos.fill(HIST("QAbefore/k0sDauNegDCA"), v0.pt(), std::abs(v0.dcanegtopv())); + histos.fill(HIST("QAbefore/k0sDauTPCNsigmaPosPi"), v0.pt(), v0.daughterTPCNSigmaPosPi()); + histos.fill(HIST("QAbefore/k0sDauTPCNsigmaNegPi"), v0.pt(), v0.daughterTPCNSigmaNegPi()); + histos.fill(HIST("QAbefore/k0sNCrossedRowsPos"), v0.pt(), v0.nCrossedRowsPos()); + histos.fill(HIST("QAbefore/k0sNCrossedRowsNeg"), v0.pt(), v0.nCrossedRowsNeg()); + if (!selected) { + return; + } + nV0sAfterCuts++; + histos.fill(HIST("QAafter/k0sMassPt"), v0.pt(), v0.mK0Short()); + histos.fill(HIST("QAafter/k0sPt"), v0.pt()); + histos.fill(HIST("QAafter/k0sEta"), v0.eta()); + histos.fill(HIST("QAafter/k0sCosPA"), v0.pt(), v0.v0CosPA()); + histos.fill(HIST("QAafter/k0sRadius"), v0.pt(), v0.transRadius()); + histos.fill(HIST("QAafter/k0sDauDCA"), v0.pt(), v0.daughDCA()); + histos.fill(HIST("QAafter/k0sDCAtoPV"), v0.pt(), std::abs(v0.dcav0topv())); + histos.fill(HIST("QAafter/k0sProperLifetime"), v0.pt(), properLifetime); + histos.fill(HIST("QAafter/k0sArmenteros"), v0.alpha(), v0.qtarm()); + histos.fill(HIST("QAafter/k0sDauPosDCA"), v0.pt(), std::abs(v0.dcapostopv())); + histos.fill(HIST("QAafter/k0sDauNegDCA"), v0.pt(), std::abs(v0.dcanegtopv())); + histos.fill(HIST("QAafter/k0sDauTPCNsigmaPosPi"), v0.pt(), v0.daughterTPCNSigmaPosPi()); + histos.fill(HIST("QAafter/k0sDauTPCNsigmaNegPi"), v0.pt(), v0.daughterTPCNSigmaNegPi()); + histos.fill(HIST("QAafter/k0sNCrossedRowsPos"), v0.pt(), v0.nCrossedRowsPos()); + histos.fill(HIST("QAafter/k0sNCrossedRowsNeg"), v0.pt(), v0.nCrossedRowsNeg()); + }; + + auto onXi = [&](auto const& xi, bool selected) { + histos.fill(HIST("QAbefore/xiMass"), xi.mXi()); + histos.fill(HIST("QAbefore/xiPt"), xi.pt()); + histos.fill(HIST("QAbefore/xiEta"), xi.eta()); + histos.fill(HIST("QAbefore/xiDCAxy"), xi.pt(), xi.dcaXYCascToPV()); + histos.fill(HIST("QAbefore/xiDCAz"), xi.pt(), xi.dcaZCascToPV()); + histos.fill(HIST("QAbefore/xiV0CosPA"), xi.pt(), xi.v0CosPA()); + histos.fill(HIST("QAbefore/xiCascCosPA"), xi.pt(), xi.cascCosPA()); + histos.fill(HIST("QAbefore/xiV0Radius"), xi.pt(), xi.transRadius()); + histos.fill(HIST("QAbefore/xiCascRadius"), xi.pt(), xi.cascTransRadius()); + histos.fill(HIST("QAbefore/xiV0DauDCA"), xi.pt(), xi.daughDCA()); + histos.fill(HIST("QAbefore/xiCascDauDCA"), xi.pt(), xi.cascDaughDCA()); + if (!selected) { + return; + } + nCascAfterCuts++; + histos.fill(HIST("QAafter/xiMass"), xi.mXi()); + histos.fill(HIST("QAafter/xiPt"), xi.pt()); + histos.fill(HIST("QAafter/xiEta"), xi.eta()); + histos.fill(HIST("QAafter/xiDCAxy"), xi.pt(), xi.dcaXYCascToPV()); + histos.fill(HIST("QAafter/xiDCAz"), xi.pt(), xi.dcaZCascToPV()); + histos.fill(HIST("QAafter/xiV0CosPA"), xi.pt(), xi.v0CosPA()); + histos.fill(HIST("QAafter/xiCascCosPA"), xi.pt(), xi.cascCosPA()); + histos.fill(HIST("QAafter/xiV0Radius"), xi.pt(), xi.transRadius()); + histos.fill(HIST("QAafter/xiCascRadius"), xi.pt(), xi.cascTransRadius()); + histos.fill(HIST("QAafter/xiV0DauDCA"), xi.pt(), xi.daughDCA()); + histos.fill(HIST("QAafter/xiCascDauDCA"), xi.pt(), xi.cascDaughDCA()); + }; + + const bool kinCutsOn = core.kinCutsEnabled(); + auto onCandidate = [&](auto const& /*xi*/, auto const& /*v0*/, XiK0sCandidateValues const& c) { + if (kinCutsOn) { + histos.fill(HIST("QAbefore/omegaAlphaVsPt"), c.omega.Pt(), c.alpha); + if (!c.passesKinCut) { + return; } + histos.fill(HIST("QAafter/omegaAlphaVsPt"), c.omega.Pt(), c.alpha); } - } + if (!c.inRapidity) { + return; + } + histos.fill(HIST("omega2012/invmass"), c.omega.M()); + histos.fill(HIST("omega2012/massPtCent"), c.omega.M(), c.omega.Pt(), cent); + }; - // Fill event QA for after-cuts counters (only for same-event) - if constexpr (!IsMix) { - histos.fill(HIST("Event/nCascadesAfterCuts"), nCascAfterCuts); - histos.fill(HIST("Event/nV0sAfterCuts"), nV0sAfterCuts); - } + core.forEachXiK0sCandidate(histos, collision, collision, cascades, v0s, false, onXi, onK0s, onCandidate); + + histos.fill(HIST("Event/nCascadesAfterCuts"), nCascAfterCuts); + histos.fill(HIST("Event/nV0sAfterCuts"), nV0sAfterCuts); } void processDummy(aod::ResoCollision const& /*collision*/) @@ -897,10 +410,10 @@ struct Omega2012Analysis { aod::ResoCascades const& resocasc, aod::ResoV0s const& resov0s) { - if (cRecoINELgt0 && !collision.isRecINELgt0()) + if (!core.passesEventCuts(collision)) { return; - - fill(collision, resocasc, resov0s); + } + fillXiK0s(collision, resocasc, resov0s); } PROCESS_SWITCH(Omega2012Analysis, processData, "Process Event for data", false); @@ -908,98 +421,53 @@ struct Omega2012Analysis { aod::ResoCascades const& resocasc, aod::ResoV0s const& resov0s) { - auto cascV0sTuple = std::make_tuple(resocasc, resov0s); Pair pairs{colBinning, nEvtMixing, -1, collisions, cascV0sTuple, &cache}; + const bool kinCutsOn = core.kinCutsEnabled(); for (const auto& [collision1, casc1, collision2, v0s2] : pairs) { - if (cRecoINELgt0 && (!collision1.isRecINELgt0() || !collision2.isRecINELgt0())) + if (!core.passesEventCuts(collision1) || !core.passesEventCuts(collision2)) { continue; - + } auto cent = collision1.cent(); - int unusedXiCount = 0; - int unusedK0sCount = 0; - auto selectedXis = selectXiCandidates(casc1, unusedXiCount); - auto selectedK0s = selectK0sCandidates(collision2, v0s2, unusedK0sCount); - - for (const auto& selectedXi : selectedXis) { - const auto& xi = selectedXi.candidate; - - for (const auto& selectedK0 : selectedK0s) { - const auto& v0 = selectedK0.candidate; - ROOT::Math::PxPyPzEVector pXi, pK0s, pRes; - pXi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); - pK0s = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(v0.pt(), v0.eta(), v0.phi(), massK0)); - pRes = pXi + pK0s; - - float alpha = 0.f; - bool kinCutFlag = true; - if (cKinCuts) { - kinCutFlag = kinCuts(pXi, pK0s, pRes, alpha); - } - - if (cKinCuts && !kinCutFlag) - continue; - - if (std::abs(pRes.Rapidity()) >= cfgRapidityCut) - continue; - - histos.fill(HIST("omega2012/invmass_Mix"), pRes.M()); - histos.fill(HIST("omega2012/massPtCent_Mix"), pRes.M(), pRes.Pt(), cent); + // Xi from collision 1, K0s from collision 2 (selected with the vertex of collision 2) + auto onCandidate = [&](auto const& /*xi*/, auto const& /*v0*/, XiK0sCandidateValues const& c) { + if (kinCutsOn && !c.passesKinCut) { + return; } - } + if (!c.inRapidity) { + return; + } + histos.fill(HIST("omega2012/invmass_Mix"), c.omega.M()); + histos.fill(HIST("omega2012/massPtCent_Mix"), c.omega.M(), c.omega.Pt(), cent); + }; + core.forEachXiK0sCandidate(histos, collision1, collision2, casc1, v0s2, false, nullptr, nullptr, onCandidate); } } PROCESS_SWITCH(Omega2012Analysis, processMixedEvent, "Process Mixed Event", false); // MC reconstructed processing: match reconstructed Xi + K0s pairs to a common Omega(2012) mother void processMC(ResoMCCollisions::iterator const& collision, - soa::Join const& resocasc, - soa::Join const& resov0s) + ResoMCCascades const& resocasc, + ResoMCV0s const& resov0s) { - if (cRecoINELgt0 && !collision.isRecINELgt0()) + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { return; - if (cMCINELgt0 && !collision.isINELgt0()) - return; - if (cMCVtxIn10 && !collision.isVtxIn10()) - return; - - int nCascAfterCuts = 0; - int nV0sAfterCuts = 0; - auto selectedXis = selectXiCandidates(resocasc, nCascAfterCuts); - auto selectedK0s = selectK0sCandidates(collision, resov0s, nV0sAfterCuts); - - for (const auto& selectedXi : selectedXis) { - const auto& xi = selectedXi.candidate; - if (std::abs(xi.pdgCode()) != kXiMinus) - continue; - - for (const auto& selectedK0 : selectedK0s) { - const auto& v0 = selectedK0.candidate; - - if (sharesAnyDaughterId(selectedXi.daughterIds, selectedK0.daughterIds)) - continue; - if (std::abs(v0.pdgCode()) != kK0Short) - continue; - if (xi.motherId() < 0 || xi.motherId() != v0.motherId()) - continue; - if (std::abs(xi.motherPDG()) != kOmega2012Minus || xi.motherPDG() != v0.motherPDG()) - continue; - - ROOT::Math::PxPyPzEVector pXi, pK0s, pRes; - pXi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); - pK0s = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(v0.pt(), v0.eta(), v0.phi(), massK0)); - pRes = pXi + pK0s; - - if (std::abs(pRes.Rapidity()) >= cfgRapidityCut) - continue; - - histos.fill(HIST("MC/hMCRecOmega2012Pt"), pRes.Pt()); - histos.fill(HIST("MC/hMCRecOmega2012PtEta"), pRes.Pt(), pRes.Eta()); - histos.fill(HIST("MC/hMCRecXiPt"), xi.pt()); - histos.fill(HIST("MC/hMCRecK0sPt"), v0.pt()); - } } + // No kinematic cut in the truth-matched spectra (as before the refactoring) + auto onCandidate = [&](auto const& xi, auto const& v0, XiK0sCandidateValues const& c) { + if (classifyXiK0sTruth(xi, v0) != XiK0sTruth::Matched) { + return; + } + if (!c.inRapidity) { + return; + } + histos.fill(HIST("MC/hMCRecOmega2012Pt"), c.omega.Pt()); + histos.fill(HIST("MC/hMCRecOmega2012PtEta"), c.omega.Pt(), c.omega.Eta()); + histos.fill(HIST("MC/hMCRecXiPt"), xi.pt()); + histos.fill(HIST("MC/hMCRecK0sPt"), v0.pt()); + }; + core.forEachXiK0sCandidate(histos, collision, collision, resocasc, resov0s, false, nullptr, nullptr, onCandidate); } PROCESS_SWITCH(Omega2012Analysis, processMC, "Process MC with truth matching", false); @@ -1012,8 +480,9 @@ struct Omega2012Analysis { // Look for Omega(2012) int pdg = mcParticle.pdgCode(); - if (std::abs(pdg) != kOmega2012Minus) + if (std::abs(pdg) != kOmega2012Minus) { continue; + } // Fill generated level histograms auto pt = mcParticle.pt(); @@ -1026,8 +495,9 @@ struct Omega2012Analysis { // Get daughters auto daughters = mcParticle.daughters_as(); - if (daughters.size() != kNumExpectedDaughters) + if (daughters.size() != NumExpectedDaughters) { continue; + } int daughter1PDG = 0, daughter2PDG = 0; ROOT::Math::PxPyPzEVector p1, p2, pMother; @@ -1066,134 +536,155 @@ struct Omega2012Analysis { } PROCESS_SWITCH(Omega2012Analysis, processMCGenerated, "Process MC generated particles", false); - // Fill function for 3-body decay analysis - template - void fillThreeBody(const CollisionT& collision, const CascadesT& cascades, const V0sT& v0s, const TracksT& tracks, const TrackIdsT& trackIds) + // Mode B: Xi(1530)0 K- -> Xi- pi+ K- (and charge conjugate) from one collision, or Xi and tracks from two mixed collisions + template + void fillXi1530K(const CollisionT& collision, float cent, const CascadesT& cascades, const TracksT& tracks, const TrackIdsT& trackIds) { - auto cent = collision.cent(); - int unusedXiCount = 0; - int unusedK0sCount = 0; - auto selectedXis = selectXiCandidates(cascades, unusedXiCount); - auto selectedK0s = selectK0sCandidates(collision, v0s, unusedK0sCount); - - // First loop: xi + pion to check xi1530 mass window - for (const auto& selectedXi : selectedXis) { - const auto& xi = selectedXi.candidate; - - for (const auto& pion : tracks) { - auto pionTrackId = trackSourceId(pion, trackIds); - if (sharesDaughterId(selectedXi.daughterIds, pionTrackId)) { - continue; + // Input QA: same event only (each object once per collision) + auto onXi = [&](auto const& xi, bool selected) { + if constexpr (IsMix) { + return; + } + histos.fill(HIST("xi1530K/QAbefore/xiMass"), xi.mXi()); + histos.fill(HIST("xi1530K/QAbefore/xiPt"), xi.pt()); + histos.fill(HIST("xi1530K/QAbefore/xiEta"), xi.eta()); + if (!selected) { + return; + } + histos.fill(HIST("xi1530K/QAafter/xiMass"), xi.mXi()); + histos.fill(HIST("xi1530K/QAafter/xiPt"), xi.pt()); + histos.fill(HIST("xi1530K/QAafter/xiEta"), xi.eta()); + }; + + auto onTrack = [&](auto const& track, int pionStage, int kaonStage) { + if constexpr (IsMix) { + return; + } + const bool hasTOF = track.hasTOF(); + histos.fill(HIST("xi1530K/QAbefore/trackPt"), track.pt()); + histos.fill(HIST("xi1530K/QAbefore/trackEta"), track.eta()); + histos.fill(HIST("xi1530K/QAbefore/trackDCAxy"), track.pt(), track.dcaXY()); + histos.fill(HIST("xi1530K/QAbefore/trackDCAz"), track.pt(), track.dcaZ()); + histos.fill(HIST("xi1530K/QAbefore/pionTPCNSigma"), track.pt(), track.tpcNSigmaPi()); + histos.fill(HIST("xi1530K/QAbefore/kaonTPCNSigma"), track.pt(), track.tpcNSigmaKa()); + if (hasTOF) { + histos.fill(HIST("xi1530K/QAbefore/pionTOFNSigma"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("xi1530K/QAbefore/kaonTOFNSigma"), track.pt(), track.tofNSigmaKa()); + } + if (pionStage == o2::analysis::resonance::kTrkPID) { + histos.fill(HIST("xi1530K/QAafter/pionPt"), track.pt()); + histos.fill(HIST("xi1530K/QAafter/pionEta"), track.eta()); + histos.fill(HIST("xi1530K/QAafter/pionDCAxy"), track.pt(), track.dcaXY()); + histos.fill(HIST("xi1530K/QAafter/pionDCAz"), track.pt(), track.dcaZ()); + histos.fill(HIST("xi1530K/QAafter/pionTPCNSigma"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("xi1530K/QAafter/pionTOFNSigma"), track.pt(), track.tofNSigmaPi()); } - - // Pion QA before cuts - histos.fill(HIST("QAbefore/pionPt"), pion.pt()); - histos.fill(HIST("QAbefore/pionEta"), pion.eta()); - - if constexpr (IsResoMicrotrack) { - histos.fill(HIST("QAbefore/pionDCAxy"), pion.pt(), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(pion.trackSelectionFlags())); - histos.fill(HIST("QAbefore/pionDCAz"), pion.pt(), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(pion.trackSelectionFlags())); - histos.fill(HIST("QAbefore/pionTPCNSigma"), pion.pt(), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(pion.pidNSigmaPiFlag())); - if (pion.hasTOF()) { - histos.fill(HIST("QAbefore/pionTOFNSigma"), pion.pt(), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(pion.pidNSigmaPiFlag())); - } - } else { - histos.fill(HIST("QAbefore/pionDCAxy"), pion.pt(), pion.dcaXY()); - histos.fill(HIST("QAbefore/pionDCAz"), pion.pt(), pion.dcaZ()); - histos.fill(HIST("QAbefore/pionTPCNSigma"), pion.pt(), pion.tpcNSigmaPi()); - if (pion.hasTOF()) { - histos.fill(HIST("QAbefore/pionTOFNSigma"), pion.pt(), pion.tofNSigmaPi()); - } - if constexpr (requires { pion.tpcNClsFound(); }) { - histos.fill(HIST("QAbefore/pionTPCNcls"), pion.tpcNClsFound()); - } + } + if (kaonStage == o2::analysis::resonance::kTrkPID) { + histos.fill(HIST("xi1530K/QAafter/kaonPt"), track.pt()); + histos.fill(HIST("xi1530K/QAafter/kaonEta"), track.eta()); + histos.fill(HIST("xi1530K/QAafter/kaonDCAxy"), track.pt(), track.dcaXY()); + histos.fill(HIST("xi1530K/QAafter/kaonDCAz"), track.pt(), track.dcaZ()); + histos.fill(HIST("xi1530K/QAafter/kaonTPCNSigma"), track.pt(), track.tpcNSigmaKa()); + if (hasTOF) { + histos.fill(HIST("xi1530K/QAafter/kaonTOFNSigma"), track.pt(), track.tofNSigmaKa()); } + } + }; - if (!pionCut(pion)) - continue; - - // Pion QA after cuts - histos.fill(HIST("QAafter/pionPt"), pion.pt()); - histos.fill(HIST("QAafter/pionEta"), pion.eta()); - - if constexpr (IsResoMicrotrack) { - histos.fill(HIST("QAafter/pionDCAxy"), pion.pt(), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(pion.trackSelectionFlags())); - histos.fill(HIST("QAafter/pionDCAz"), pion.pt(), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(pion.trackSelectionFlags())); - histos.fill(HIST("QAafter/pionTPCNSigma"), pion.pt(), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(pion.pidNSigmaPiFlag())); - if (pion.hasTOF()) { - histos.fill(HIST("QAafter/pionTOFNSigma"), pion.pt(), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(pion.pidNSigmaPiFlag())); - } - } else { - histos.fill(HIST("QAafter/pionDCAxy"), pion.pt(), pion.dcaXY()); - histos.fill(HIST("QAafter/pionDCAz"), pion.pt(), pion.dcaZ()); - histos.fill(HIST("QAafter/pionTPCNSigma"), pion.pt(), pion.tpcNSigmaPi()); - if (pion.hasTOF()) { - histos.fill(HIST("QAafter/pionTOFNSigma"), pion.pt(), pion.tofNSigmaPi()); - } - if constexpr (requires { pion.tpcNClsFound(); }) { - histos.fill(HIST("QAafter/pionTPCNcls"), pion.tpcNClsFound()); + const bool fillWrongSign = core.fillWrongSign(); + auto onCandidate = [&]([[maybe_unused]] auto const& xi, [[maybe_unused]] auto const& pion, [[maybe_unused]] auto const& kaon, Xi1530KCandidateValues const& c) { + if (c.chargePattern != o2::analysis::omega2012ml::kSignalPattern) { + // Charge-pattern controls: same event only, never in the signal histograms + if constexpr (!IsMix) { + if (fillWrongSign && c.inRapidity) { + histos.fill(HIST("xi1530K_wrongSign/invmassPattern"), c.chargePattern, c.omega.M()); + histos.fill(HIST("xi1530K_wrongSign/massPtPattern"), c.omega.M(), c.omega.Pt(), c.chargePattern); } } - - // Check xi1530 mass window cut - if (!xi1530MassCut(xi, pion)) - continue; - - // Second loop: v0 for the selected xi-pion pair - for (const auto& selectedK0 : selectedK0s) { - const auto& v0 = selectedK0.candidate; - if (sharesAnyDaughterId(selectedXi.daughterIds, selectedK0.daughterIds)) { - continue; - } - if (sharesDaughterId(selectedK0.daughterIds, pionTrackId)) { - continue; - } - - // 4-vectors for 3-body decay: Xi + K0s + pion - ROOT::Math::PxPyPzEVector pXi, pK0s, pPion, pRes; - pXi = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(xi.pt(), xi.eta(), xi.phi(), xi.mXi())); - pK0s = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(v0.pt(), v0.eta(), v0.phi(), massK0)); - pPion = ROOT::Math::PxPyPzEVector(ROOT::Math::PtEtaPhiMVector(pion.pt(), pion.eta(), pion.phi(), MassPionCharged)); - - pRes = pXi + pK0s + pPion; - - if (std::abs(pRes.Rapidity()) >= cfgRapidityCut) - continue; - - histos.fill(HIST("omega2012_3body/invmass"), pRes.M()); - histos.fill(HIST("omega2012_3body/massPtCent"), pRes.M(), pRes.Pt(), cent); + return; + } + if constexpr (IsMix) { + if (c.inRapidity) { + histos.fill(HIST("xi1530K/invmass_Mix"), c.omega.M()); + histos.fill(HIST("xi1530K/massPtCent_Mix"), c.omega.M(), c.omega.Pt(), cent); } + return; } - } + histos.fill(HIST("xi1530K/massXiPi"), c.massXiPi); + if (!c.inRapidity) { + return; + } + histos.fill(HIST("xi1530K/invmass"), c.omega.M()); + histos.fill(HIST("xi1530K/massPtCent"), c.omega.M(), c.omega.Pt(), cent); + histos.fill(HIST("xi1530K/massXiPiVsMass"), c.omega.M(), c.massXiPi); + if constexpr (IsMC) { + if (classifyXi1530KTruth(xi, pion, kaon) != Xi1530KTruth::Matched) { + return; + } + histos.fill(HIST("xi1530K/MC/hMCRecOmega2012Pt"), c.omega.Pt()); + histos.fill(HIST("xi1530K/MC/hMCRecOmega2012PtEta"), c.omega.Pt(), c.omega.Eta()); + histos.fill(HIST("xi1530K/MC/hMCRecMass"), c.omega.M()); + histos.fill(HIST("xi1530K/MC/hMCRecMassXiPi"), c.massXiPi); + histos.fill(HIST("xi1530K/MC/hMCRecXiPt"), xi.pt()); + histos.fill(HIST("xi1530K/MC/hMCRecPionPt"), pion.pt()); + histos.fill(HIST("xi1530K/MC/hMCRecKaonPt"), kaon.pt()); + } + }; + + core.forEachXi1530KCandidate(histos, collision, cascades, tracks, trackIds, false, onXi, onTrack, onCandidate); } - // 3-body decay analysis: Xi + pi + K0s with ResoTracks - void processThreeBodyWithTracks(ResoCollisions::iterator const& collision, - aod::ResoCascades const& resocasc, - aod::ResoV0s const& resov0s, - aod::ResoTracks const& resotracks, - aod::ResoTrackTracks const& resotrackids) + void processXi1530KMicro(ResoCollisions::iterator const& collision, + aod::ResoCascades const& resocasc, + ResoMicroTracks const& resomicrotracks) { - if (cRecoINELgt0 && !collision.isRecINELgt0()) + if (!core.passesEventCuts(collision)) { return; + } + fillXi1530K(collision, collision.cent(), resocasc, resomicrotracks, nullptr); + } + PROCESS_SWITCH(Omega2012Analysis, processXi1530KMicro, "Process Xi(1530)0 K mode with ResoMicroTracks_001", false); - fillThreeBody(collision, resocasc, resov0s, resotracks, resotrackids); + void processXi1530KTracks(ResoCollisions::iterator const& collision, + aod::ResoCascades const& resocasc, + aod::ResoTracks const& resotracks, + aod::ResoTrackTracks const& resotrackids) + { + if (!core.passesEventCuts(collision)) { + return; + } + fillXi1530K(collision, collision.cent(), resocasc, resotracks, resotrackids); } - PROCESS_SWITCH(Omega2012Analysis, processThreeBodyWithTracks, "Process 3-body decay with ResoTracks", false); - - // 3-body decay analysis: Xi + pi + K0s with ResoMicroTracks - void processThreeBodyWithMicroTracks(ResoCollisions::iterator const& collision, - aod::ResoCascades const& resocasc, - aod::ResoV0s const& resov0s, - aod::ResoMicroTracks const& resomicrotracks, - aod::ResoMicroTrackTracks const& resomicrotrackids) + PROCESS_SWITCH(Omega2012Analysis, processXi1530KTracks, "Process Xi(1530)0 K mode with ResoTracks", false); + + void processXi1530KMCMicro(ResoMCCollisions::iterator const& collision, + ResoMCCascades const& resocasc, + ResoMCMicroTracks const& resomicrotracks) { - if (cRecoINELgt0 && !collision.isRecINELgt0()) + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { return; + } + fillXi1530K(collision, collision.cent(), resocasc, resomicrotracks, nullptr); + } + PROCESS_SWITCH(Omega2012Analysis, processXi1530KMCMicro, "Process Xi(1530)0 K mode with truth matching (ResoMicroTracks_001)", false); - fillThreeBody(collision, resocasc, resov0s, resomicrotracks, resomicrotrackids); + // Mixed events: Xi from collision 1, pion and kaon from collision 2 + void processXi1530KMixedMicro(ResoCollisions const& collisions, + aod::ResoCascades const& resocasc, + ResoMicroTracks const& resomicrotracks) + { + auto cascTracksTuple = std::make_tuple(resocasc, resomicrotracks); + Pair pairs{colBinning, nEvtMixing, -1, collisions, cascTracksTuple, &cache}; + for (const auto& [collision1, casc1, collision2, tracks2] : pairs) { + if (!core.passesEventCuts(collision1) || !core.passesEventCuts(collision2)) { + continue; + } + fillXi1530K(collision1, collision1.cent(), casc1, tracks2, nullptr); + } } - PROCESS_SWITCH(Omega2012Analysis, processThreeBodyWithMicroTracks, "Process 3-body decay with ResoMicroTracks", false); + PROCESS_SWITCH(Omega2012Analysis, processXi1530KMixedMicro, "Process mixed events of the Xi(1530)0 K mode (ResoMicroTracks_001)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Tasks/Resonances/omega2012TrainingTable.cxx b/PWGLF/Tasks/Resonances/omega2012TrainingTable.cxx new file mode 100644 index 00000000000..18502ccedb1 --- /dev/null +++ b/PWGLF/Tasks/Resonances/omega2012TrainingTable.cxx @@ -0,0 +1,404 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file omega2012TrainingTable.cxx +/// \brief Derived Omega(2012) training tables, separately for the XiK0s and Xi1530K decay modes +/// \author Bong-Hwi Lim +/// +/// Mode A (XiK0s): processXiK0s, processXiK0sMC. Mode B (Xi1530K, Xi- pi+ K- with micro tracks): +/// processXi1530KMicro, processXi1530KMCMicro. Generated parents: processMCTrue. +/// Each candidate is written only to the tables of its own mode (PWGLF/DataModel/LFOmega2012MlTables.h); +/// the selection is the one of omega2012Analysis.cxx (PWGLF/Core/Omega2012AnalysisCore.h, same JSON keys). + +#include "PWGLF/Core/Omega2012AnalysisCore.h" +#include "PWGLF/Core/Omega2012MlFeatures.h" +#include "PWGLF/Core/ResoAnalysisSelectionCore.h" +#include "PWGLF/DataModel/LFOmega2012MlTables.h" +#include "PWGLF/DataModel/LFResonanceTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::analysis::omega2012; + +namespace ml = o2::analysis::omega2012ml; + +static_assert(std::extent_v == ml::XiK0sNMasterFeatures, + "OmXiK0sFeatures column size must match the XiK0s feature contract"); +static_assert(std::extent_v == ml::Xi1530KNMasterFeatures, + "OmXi1530KFeatures column size must match the Xi1530K feature contract"); + +struct Omega2012TrainingTable { + using ResoCollisions = aod::ResoCollisions_001; + using ResoMCCols = soa::Join; + using ResoMicroTracks = aod::ResoMicroTracks_001; + using ResoMCMicroTracks = soa::Join; + using ResoMCCascades = soa::Join; + using ResoMCV0s = soa::Join; + using ResoMCParents = aod::ResoMCParents_001; + + // FNV-1a over the bit patterns of the master features, for parity logs + static constexpr uint64_t FnvOffsetBasis = 14695981039346656037ULL; + static constexpr uint64_t FnvPrime = 1099511628211ULL; + static constexpr unsigned int BitsPerByte = 8; + static constexpr unsigned int BitsPerFloat = 32; + static constexpr uint32_t ByteMask = 0xffU; + static constexpr int NBuildStatus = 6; + + Produces mlEvents; + Produces mlCascades; + Produces mlV0s; + Produces mlTracks; + Produces mlXiK0sCandidates; + Produces mlXiK0sInputs; + Produces mlXiK0sTruth; + Produces mlXi1530KCandidates; + Produces mlXi1530KInputs; + Produces mlXi1530KTruth; + Produces mlGenAudit; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Selection shared with the Omega(2012) histogram task (same JSON keys) + o2::analysis::resonance::EventCuts eventCuts; + XiCuts xiCuts; + K0sCuts k0sCuts; + Xi1530KCuts xi1530KCuts; + o2::analysis::resonance::TrackCuts trackCuts = makeXi1530KTrackCuts(); + PionPidCuts pionPID; + KaonPidCuts kaonPID; + CandidateCuts candidateCuts; + + Configurable omega2012MlExportStage{"omega2012MlExportStage", "loose", "Candidate export stage: loose (pass bits >= 1) or selected (all pass bits of the mode)"}; + Configurable omega2012MlLooseAudit{"omega2012MlLooseAudit", true, "Record the loose cut flow and mass/pT/activity spectrum per mode"}; + Configurable omega2012MlParityRows{"omega2012MlParityRows", 0, "Log feature hashes of the first N exported candidates of each mode"}; + Configurable omega2012MlExportWrongSign{"omega2012MlExportWrongSign", false, "Also export the wrong-sign charge-pattern controls of the Xi1530K mode"}; + + Omega2012AnalysisCore core; + int64_t mlEventRow = -1; + std::unordered_map mlCascadeRows; + std::unordered_map mlV0Rows; + std::unordered_map mlTrackRows; + int parityLoggedXiK0s = 0; + int parityLoggedXi1530K = 0; + + void init(InitContext&) + { + const bool xiK0s = doprocessXiK0s || doprocessXiK0sMC; + const bool xi1530K = doprocessXi1530KMicro || doprocessXi1530KMCMicro; + if ((doprocessXiK0s && doprocessXiK0sMC) || (doprocessXi1530KMicro && doprocessXi1530KMCMicro)) { + LOG(fatal) << "Enable at most one of processXiK0s/processXiK0sMC and one of processXi1530KMicro/processXi1530KMCMicro"; + } + if (!xiK0s && !xi1530K) { + LOG(fatal) << "Enable at least one candidate process function"; + } + if (doprocessMCTrue && !doprocessXiK0sMC && !doprocessXi1530KMCMicro) { + LOG(fatal) << "processMCTrue requires processXiK0sMC or processXi1530KMCMicro"; + } + if (omega2012MlParityRows < 0) { + LOG(fatal) << "omega2012MlParityRows must not be negative"; + } + if (omega2012MlExportStage.value != "loose" && omega2012MlExportStage.value != "selected") { + LOG(fatal) << "omega2012MlExportStage must be loose or selected"; + } + + ProcessModes modes; + modes.xiK0s = xiK0s; + modes.xi1530K = xi1530K; + modes.microTracks = xi1530K; + modes.mcReco = doprocessXiK0sMC || doprocessXi1530KMCMicro; + modes.mcGen = doprocessMCTrue; + LooseStageOptions looseOptions; + looseOptions.audit = omega2012MlLooseAudit; + looseOptions.exportSelected = omega2012MlExportStage.value == "selected"; + core.init(histos, eventCuts, xiCuts, k0sCuts, xi1530KCuts, trackCuts, pionPID, kaonPID, candidateCuts, modes, looseOptions); + + // Candidates that violate the canonical or feature contract are skipped, never written. + const std::array statusLabels{"Ok", "InvalidChargePattern", "ReusedDaughter", "InvalidMomentum", "InvalidKinematics", "InvalidContract"}; + if (xiK0s) { + auto skipped = histos.add("ML/xiK0s/exportSkipped", "XiK0s mode skipped candidates;build status;candidates", HistType::kTH1D, {{NBuildStatus, -0.5, NBuildStatus - 0.5}}); + for (std::size_t i = 0; i < statusLabels.size(); ++i) { + skipped->GetXaxis()->SetBinLabel(i + 1, statusLabels[i]); + } + } + if (xi1530K) { + auto skipped = histos.add("ML/xi1530K/exportSkipped", "Xi1530K mode skipped candidates;build status;candidates", HistType::kTH1D, {{NBuildStatus, -0.5, NBuildStatus - 0.5}}); + for (std::size_t i = 0; i < statusLabels.size(); ++i) { + skipped->GetXaxis()->SetBinLabel(i + 1, statusLabels[i]); + } + histos.add("ML/xi1530K/chargePattern", "Xi1530K mode candidates handed to the export;charge pattern (0 signal);candidates", HistType::kTH1D, {{4, -0.5, 3.5}}); + } + + LOG(info) << "Size of the histograms in Omega(2012) training table task"; + histos.print(); + } + + template + static uint64_t featureHash(std::array const& features) + { + uint64_t hash = FnvOffsetBasis; + for (const auto& value : features) { + const auto bits = std::bit_cast(value); + for (unsigned int shift = 0; shift < BitsPerFloat; shift += BitsPerByte) { + hash = (hash ^ ((bits >> shift) & ByteMask)) * FnvPrime; + } + } + return hash; + } + + template + void writeEvent(Collision const& collision, DecayMode mode) + { + mlCascadeRows.clear(); + mlV0Rows.clear(); + mlTrackRows.clear(); + // ResoCollisions_001 carries no run number or BC; the reduced collision row identifies the event within its DF. + mlEvents(static_cast(mode), static_cast(collision.globalIndex()), + collision.posZ(), collision.bMagField(), collision.cent(), collision.multiplicity(), collision.isRecINELgt0()); + mlEventRow = mlEvents.lastIndex(); + } + + template + int64_t writeCascade(Cascade const& xi) + { + const auto id = static_cast(xi.globalIndex()); + if (auto it = mlCascadeRows.find(id); it != mlCascadeRows.end()) { + return it->second; + } + const auto indices = xi.cascadeIndices(); + const std::array daughterIds{indices[0], indices[1], indices[2]}; + mlCascades(mlEventRow, id, xi.px(), xi.py(), xi.pz(), static_cast(xi.sign()), xi.mXi(), xi.mLambda(), + xi.v0CosPA(), xi.cascCosPA(), xi.daughDCA(), xi.cascDaughDCA(), + xi.dcapostopv(), xi.dcanegtopv(), xi.dcabachtopv(), xi.dcav0topv(), xi.dcaXYCascToPV(), xi.dcaZCascToPV(), + xi.transRadius(), xi.cascTransRadius(), xi.decayVtxX(), xi.decayVtxY(), xi.decayVtxZ(), + xi.daughterTPCNSigmaPosPi10(), xi.daughterTPCNSigmaPosKa10(), xi.daughterTPCNSigmaPosPr10(), + xi.daughterTPCNSigmaNegPi10(), xi.daughterTPCNSigmaNegKa10(), xi.daughterTPCNSigmaNegPr10(), + xi.daughterTPCNSigmaBachPi10(), xi.daughterTPCNSigmaBachKa10(), xi.daughterTPCNSigmaBachPr10(), + xi.daughterTOFNSigmaPosPi10(), xi.daughterTOFNSigmaPosKa10(), xi.daughterTOFNSigmaPosPr10(), + xi.daughterTOFNSigmaNegPi10(), xi.daughterTOFNSigmaNegKa10(), xi.daughterTOFNSigmaNegPr10(), + xi.daughterTOFNSigmaBachPi10(), xi.daughterTOFNSigmaBachKa10(), xi.daughterTOFNSigmaBachPr10(), + xi.nCrossedRowsPos(), xi.nCrossedRowsNeg(), xi.nCrossedRowsBach(), daughterIds.data()); + const auto row = mlCascades.lastIndex(); + mlCascadeRows.emplace(id, row); + return row; + } + + template + int64_t writeV0(V0 const& v0) + { + const auto id = static_cast(v0.globalIndex()); + if (auto it = mlV0Rows.find(id); it != mlV0Rows.end()) { + return it->second; + } + const auto indices = v0.indices(); + const std::array daughterIds{indices[0], indices[1]}; + mlV0s(mlEventRow, id, v0.px(), v0.py(), v0.pz(), v0.mK0Short(), v0.mLambda(), v0.mAntiLambda(), + v0.v0CosPA(), v0.daughDCA(), v0.dcapostopv(), v0.dcanegtopv(), v0.dcav0topv(), v0.transRadius(), + v0.decayVtxX(), v0.decayVtxY(), v0.decayVtxZ(), v0.alpha(), v0.qtarm(), + v0.daughterTPCNSigmaPosPi10(), v0.daughterTPCNSigmaPosKa10(), v0.daughterTPCNSigmaPosPr10(), + v0.daughterTPCNSigmaNegPi10(), v0.daughterTPCNSigmaNegKa10(), v0.daughterTPCNSigmaNegPr10(), + v0.daughterTOFNSigmaPosPi10(), v0.daughterTOFNSigmaPosKa10(), v0.daughterTOFNSigmaPosPr10(), + v0.daughterTOFNSigmaNegPi10(), v0.daughterTOFNSigmaNegKa10(), v0.daughterTOFNSigmaNegPr10(), + v0.nCrossedRowsPos(), v0.nCrossedRowsNeg(), daughterIds.data()); + const auto row = mlV0s.lastIndex(); + mlV0Rows.emplace(id, row); + return row; + } + + template + int64_t writeTrack(Track const& track) + { + const auto id = static_cast(track.globalIndex()); + if (auto it = mlTrackRows.find(id); it != mlTrackRows.end()) { + return it->second; + } + mlTracks(mlEventRow, static_cast(track.trackId()), track.px(), track.py(), track.pz(), + track.pidNSigmaPiFlag(), track.pidNSigmaKaFlag(), track.pidNSigmaPrFlag(), + track.trackSelectionFlags(), track.trackFlags(), track.tpcNClsCrossedRows(), track.itsClusterMap()); + const auto row = mlTracks.lastIndex(); + mlTrackRows.emplace(id, row); + return row; + } + + template + void writeXiK0sCandidate(Collision const& collision, Cascade const& xi, V0 const& v0, XiK0sCandidateValues const& values, uint16_t passBits) + { + const auto canonical = ml::canonicalizeXiK0s(ml::makeCascadeSnapshot(collision, xi), ml::makeV0Snapshot(collision, v0)); + if (canonical.status != ml::BuildStatus::Ok) { + histos.fill(HIST("ML/xiK0s/exportSkipped"), static_cast(canonical.status)); + return; + } + const auto pack = ml::buildXiK0sFeatures(canonical.candidate); + if (pack.status != ml::BuildStatus::Ok) { + histos.fill(HIST("ML/xiK0s/exportSkipped"), static_cast(pack.status)); + return; + } + if (parityLoggedXiK0s < omega2012MlParityRows) { + LOGP(info, "OMEGA2012MLPARITY mode=XiK0s collision={} cascade={} v0={} passBits={} featureHash={}", + collision.globalIndex(), xi.globalIndex(), v0.globalIndex(), passBits, featureHash(pack.master)); + ++parityLoggedXiK0s; + } + const auto cascadeRow = writeCascade(xi); + const auto v0Row = writeV0(v0); + mlXiK0sCandidates(mlEventRow, cascadeRow, v0Row, + static_cast(values.omega.M()), static_cast(values.omega.Pt()), + static_cast(values.omega.Rapidity()), static_cast(values.omega.Eta()), + static_cast(values.omega.Phi()), values.alpha, static_cast(xi.sign()), passBits); + const auto row = mlXiK0sCandidates.lastIndex(); + mlXiK0sInputs(row, pack.master.data(), static_cast(pack.status)); + if constexpr (IsMC) { + const bool matched = classifyXiK0sTruth(xi, v0) == XiK0sTruth::Matched; + mlXiK0sTruth(row, matched ? uint8_t{1} : uint8_t{2}, matched ? xi.motherPDG() : 0, + matched ? static_cast(xi.motherId()) : int64_t{-1}, v0.motherPDG(), xi.motherPDG()); + } else { + mlXiK0sTruth(row, uint8_t{0}, 0, int64_t{-1}, 0, 0); + } + } + + template + void writeXi1530KCandidate(Collision const& collision, Cascade const& xi, Track const& pion, Track const& kaon, + Xi1530KCandidateValues const& values, uint16_t passBits) + { + histos.fill(HIST("ML/xi1530K/chargePattern"), values.chargePattern); + if (values.chargePattern != ml::kSignalPattern && !omega2012MlExportWrongSign) { + return; // charge-pattern controls are exported only on request + } + const auto canonical = ml::canonicalizeXi1530K(ml::makeCascadeSnapshot(collision, xi), ml::makeTrackSnapshot(pion), ml::makeTrackSnapshot(kaon)); + if (canonical.status != ml::BuildStatus::Ok) { + histos.fill(HIST("ML/xi1530K/exportSkipped"), static_cast(canonical.status)); + return; + } + const auto pack = ml::buildXi1530KFeatures(canonical.candidate); + if (pack.status != ml::BuildStatus::Ok) { + histos.fill(HIST("ML/xi1530K/exportSkipped"), static_cast(pack.status)); + return; + } + if (parityLoggedXi1530K < omega2012MlParityRows) { + LOGP(info, "OMEGA2012MLPARITY mode=Xi1530K collision={} cascade={} tracks={},{} chargePattern={} passBits={} featureHash={}", + collision.globalIndex(), xi.globalIndex(), canonical.candidate.pion.sourceTrackId, canonical.candidate.kaon.sourceTrackId, + values.chargePattern, passBits, featureHash(pack.master)); + ++parityLoggedXi1530K; + } + const auto cascadeRow = writeCascade(xi); + const auto pionRow = writeTrack(pion); + const auto kaonRow = writeTrack(kaon); + mlXi1530KCandidates(mlEventRow, cascadeRow, pionRow, kaonRow, + static_cast(values.omega.M()), values.massXiPi, values.massXiK, values.massPiK, + static_cast(values.omega.Pt()), static_cast(values.omega.Rapidity()), + static_cast(values.omega.Eta()), static_cast(values.omega.Phi()), + values.chargePattern, passBits); + const auto row = mlXi1530KCandidates.lastIndex(); + mlXi1530KInputs(row, pack.master.data(), static_cast(pack.status)); + if constexpr (IsMC) { + const bool matched = classifyXi1530KTruth(xi, pion, kaon) == Xi1530KTruth::Matched; + mlXi1530KTruth(row, matched ? uint8_t{1} : uint8_t{2}, matched ? kaon.motherPDG() : 0, + matched ? static_cast(kaon.motherId()) : int64_t{-1}, static_cast(xi.motherId())); + } else { + mlXi1530KTruth(row, uint8_t{0}, 0, int64_t{-1}, int64_t{-1}); + } + } + + template + void exportXiK0s(Collision const& collision, Cascades const& cascades, V0s const& v0s) + { + writeEvent(collision, DecayMode::XiK0s); + // Selection only: the analysis histograms belong to the Omega(2012) analysis task + core.forEachXiK0sCandidate(histos, collision, collision, cascades, v0s, true, nullptr, nullptr, nullptr, + [&](auto const& coll, auto const& xi, auto const& v0, XiK0sCandidateValues const& values, uint16_t passBits) { + writeXiK0sCandidate(coll, xi, v0, values, passBits); + }); + } + + template + void exportXi1530K(Collision const& collision, Cascades const& cascades, Tracks const& tracks) + { + writeEvent(collision, DecayMode::Xi1530K); + core.forEachXi1530KCandidate(histos, collision, cascades, tracks, nullptr, true, nullptr, nullptr, nullptr, + [&](auto const& coll, auto const& xi, auto const& pion, auto const& kaon, + Xi1530KCandidateValues const& values, uint16_t passBits) { + writeXi1530KCandidate(coll, xi, pion, kaon, values, passBits); + }); + } + + void processXiK0s(ResoCollisions::iterator const& collision, aod::ResoCascades const& cascades, aod::ResoV0s const& v0s) + { + if (!core.passesEventCuts(collision)) { + return; + } + exportXiK0s(collision, cascades, v0s); + } + PROCESS_SWITCH(Omega2012TrainingTable, processXiK0s, "Write XiK0s-mode candidates from data", true); + + void processXiK0sMC(ResoMCCols::iterator const& collision, ResoMCCascades const& cascades, ResoMCV0s const& v0s) + { + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + exportXiK0s(collision, cascades, v0s); + } + PROCESS_SWITCH(Omega2012TrainingTable, processXiK0sMC, "Write XiK0s-mode candidates with truth from reconstructed MC", false); + + void processXi1530KMicro(ResoCollisions::iterator const& collision, aod::ResoCascades const& cascades, ResoMicroTracks const& tracks) + { + if (!core.passesEventCuts(collision)) { + return; + } + exportXi1530K(collision, cascades, tracks); + } + PROCESS_SWITCH(Omega2012TrainingTable, processXi1530KMicro, "Write Xi1530K-mode candidates from data micro v001 tables", false); + + void processXi1530KMCMicro(ResoMCCols::iterator const& collision, ResoMCCascades const& cascades, ResoMCMicroTracks const& tracks) + { + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + exportXi1530K(collision, cascades, tracks); + } + PROCESS_SWITCH(Omega2012TrainingTable, processXi1530KMCMicro, "Write Xi1530K-mode candidates with truth from reconstructed MC micro v001 tables", false); + + void processMCTrue(ResoMCCols::iterator const& collision, ResoMCParents const& resoParents) + { + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + core.forEachGeneratedOmega2012(histos, resoParents, [&](auto const& part, GeneratedChannel channel) { + mlGenAudit(collision.globalIndex(), static_cast(part.originalMcParticleId()), + part.pdgCode(), part.daughterPDG1(), part.daughterPDG2(), static_cast(channel), part.pt(), part.y()); + }); + } + PROCESS_SWITCH(Omega2012TrainingTable, processMCTrue, "Write generated Omega(2012) parents of selected reconstructed MC events", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}