diff --git a/PWGLF/Core/K1AnalysisMicroCore.h b/PWGLF/Core/K1AnalysisMicroCore.h new file mode 100644 index 00000000000..4e2cd0e720c --- /dev/null +++ b/PWGLF/Core/K1AnalysisMicroCore.h @@ -0,0 +1,741 @@ +// 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 K1AnalysisMicroCore.h +/// \brief Shared K1(1270) selection, truth classification and candidate enumeration of the K1 resonance tasks +/// \author Su-Jeong Ji , Bong-Hwi Lim +/// +/// The core owns the selection and its cut-flow instrumentation (CutFlow/*, ML/*). The tasks own their +/// output histograms and fill them from the pair and candidate hooks of forEachCandidate(). + +#ifndef PWGLF_CORE_K1ANALYSISMICROCORE_H_ +#define PWGLF_CORE_K1ANALYSISMICROCORE_H_ + +#include "PWGLF/Core/K1MlFeatures.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::k1micro +{ + +using o2::analysis::resonance::EventCuts; +using o2::analysis::resonance::isCutEnabled; +using o2::analysis::resonance::isInRange; +using o2::analysis::resonance::isInWindow; +using o2::analysis::resonance::PIDCutConfig; +using o2::analysis::resonance::ResoAnalysisSelectionCore; +using o2::analysis::resonance::TrackCuts; +using o2::analysis::resonance::TrackStage; + +enum class K1TruthChannel { + None = 0, + RhoK = 1, + KStarPi = 2 +}; + +// The value is the species index of the PID configuration in ResoAnalysisSelectionCore. +enum class Species : int { + Pion = 0, + Kaon = 1 +}; + +// Cumulative selection bits of an unlike-sign candidate handed to the export hook. +enum CandidatePassBit : uint16_t { + kPassLoose = 1, // valid canonical candidate inside the K1 rapidity window + kPassQuality = 2, // track quality of all three tracks + kPassPID = 4, // TOF requirement and PID of all three tracks + kPassPair = 8, // pion-pair pT and secondary mass window + kPassCandidate = 16 // candidate cuts +}; +inline constexpr uint16_t PassBitsSelected = kPassLoose | kPassQuality | kPassPID | kPassPair | kPassCandidate; + +inline constexpr double MassRho770 = 0.77526; // PDG 2024, not available in o2::constants::physics +inline constexpr int NCandidateStages = 12; +inline constexpr int NTruthChannels = 3; + +// Configurable groups without prefix: the JSON keys are the plain configurable names. + +/// Pion PID selection +struct PionPidCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 3.0, "TPC nSigma cut for Pion (-999: off)"}; // TPC + o2::framework::Configurable cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion (-999: off)"}; // TOF + o2::framework::Configurable nsigmaCutCombinedPion{"nsigmaCutCombinedPion", -999, "Combined nSigma cut for Pion"}; // Combined + o2::framework::Configurable cUseOnlyTOFTrackPi{"cUseOnlyTOFTrackPi", false, "Use only TOF track for PID selection"}; // Use only TOF track for Pion PID selection + 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 selection +struct KaonPidCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cMaxTPCnSigmaKaon{"cMaxTPCnSigmaKaon", 3.0, "TPC nSigma cut for Kaon (-999: off)"}; // TPC + o2::framework::Configurable cMaxTOFnSigmaKaon{"cMaxTOFnSigmaKaon", 3.0, "TOF nSigma cut for Kaon (-999: off)"}; // TOF + o2::framework::Configurable nsigmaCutCombinedKaon{"nsigmaCutCombinedKaon", -999, "Combined nSigma cut for Kaon"}; // Combined + o2::framework::Configurable cUseOnlyTOFTrackKa{"cUseOnlyTOFTrackKa", false, "Use only TOF track for PID selection"}; // Use only TOF track for Kaon PID selection + 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)"}; +}; + +/// Secondary selection (-999 switches a cut off; the values it needs are then not computed) +struct SecondaryCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cMinSecondaryPtCut{"cMinSecondaryPtCut", 0.5, "Min pT cut for secondary selection"}; + o2::framework::Configurable cfgModeK892orRho{"cfgModeK892orRho", false, "Secondary scenario for K892 (true) or Rho (false)"}; + o2::framework::Configurable cSecondaryMasswindow{"cSecondaryMasswindow", -999, "Secondary inv mass selection window"}; + o2::framework::Configurable cMinAnotherSecondaryMassCut{"cMinAnotherSecondaryMassCut", -999, "Min inv. mass selection of another secondary scenario"}; + o2::framework::Configurable cMaxAnotherSecondaryMassCut{"cMaxAnotherSecondaryMassCut", -999, "MAx inv. mass selection of another secondary scenario"}; + o2::framework::Configurable cMinPiKaMassCut{"cMinPiKaMassCut", -999, "bPion-Kaon pair inv mass selection minimum"}; + o2::framework::Configurable cMaxPiKaMassCut{"cMaxPiKaMassCut", -999, "bPion-Kaon pair inv mass selection maximum"}; + o2::framework::Configurable cMinAngle{"cMinAngle", -999, "Minimum angle between the secondary resonance and the bachelor"}; + o2::framework::Configurable cMaxAngle{"cMaxAngle", -999, "Maximum angle between the secondary resonance and the bachelor"}; + o2::framework::Configurable cMinPairAsym{"cMinPairAsym", -999, "Minimum pair asymmetry"}; + o2::framework::Configurable cMaxPairAsym{"cMaxPairAsym", -999, "Maximum pair asymmetry"}; +}; + +/// Common TOF switch and K1 selection +struct CandidateCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cByPassTOF{"cByPassTOF", false, "Bypass the TOF nSigma selection"}; + o2::framework::Configurable cK1MaxRap{"cK1MaxRap", 0.5, "K1 maximum rapidity"}; + o2::framework::Configurable cK1MinRap{"cK1MinRap", -0.5, "K1 minimum rapidity"}; +}; + +/// Process functions enabled in the task; they decide the configuration checks and the registered histograms. +struct ProcessModes { + bool microTracks = false; // any process function reading micro tracks (quantised DCA and nSigma) + bool mcReco = false; // reconstructed MC with full or micro tracks + bool mcRecoMicro = false; // reconstructed MC with micro tracks + bool mcGen = false; // generated K1 parents in selected reconstructed events +}; + +/// Loose-stage traversal of unlike-sign micro candidates for the export hook. +/// With the defaults and without an export hook, the candidate loop 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 +}; + +/// Selected (pion, pion, kaon) triplet handed to the candidate hook. +/// mass13, mass23, angle and pairAsym are computed only inside the K1 rapidity window, and there only +/// when a candidate cut needs them or forEachCandidate() was asked for them; otherwise they are 0. +struct K1CandidateValues { + ROOT::Math::PxPyPzMVector k1; // pion1 + pion2 + kaon + ROOT::Math::PxPyPzMVector secondary; // pion1 + pion2 + double mass13 = 0.; // pion1 + kaon (K*0 candidate in the unlike-sign case) + double mass23 = 0.; // pion2 + kaon + double angle = 0.; // opening angle between the secondary resonance and the bachelor + double pairAsym = 0.; // energy asymmetry between the secondary resonance and the bachelor + bool isUnlikeSign = false; // opposite-sign pion pair + bool inRapidity = false; // K1 rapidity window + bool passesCandidateCuts = false; // secondary mass, pi-K mass, angle and asymmetry cuts (inside the rapidity window) +}; + +// Truth classification from the immediate mothers and the sibling IDs of the reconstructed daughters. +template +bool hasSibling(const Track& directDaughter, int resonanceId) +{ + if (resonanceId < 0) { + return false; + } + const auto siblings = directDaughter.siblingIds(); + return siblings[0] == resonanceId || siblings[1] == resonanceId; +} + +template +bool matchesKStarPi(const Track& resonancePion, const Track& directPion, const Kaon& kaon) +{ + const int charge = kaon.pdgCode() > 0 ? 1 : -1; + if (resonancePion.motherId() != kaon.motherId() || + resonancePion.motherId() == directPion.motherId()) { + return false; + } + if (resonancePion.motherPDG() != charge * o2::constants::physics::Pdg::kK0Star892 || kaon.motherPDG() != charge * o2::constants::physics::Pdg::kK0Star892) { + return false; + } + if (resonancePion.pdgCode() != -charge * kPiPlus || directPion.pdgCode() != charge * kPiPlus || + directPion.motherPDG() != charge * o2::constants::physics::Pdg::kK1_1270Plus) { + return false; + } + return hasSibling(directPion, kaon.motherId()); +} + +template +K1TruthChannel classifyK1Truth(const Track& pion1, const Track& pion2, const Kaon& kaon) +{ + if (std::abs(pion1.pdgCode()) != kPiPlus || std::abs(pion2.pdgCode()) != kPiPlus || + std::abs(kaon.pdgCode()) != kKPlus) { + return K1TruthChannel::None; + } + if (pion1.motherId() < 0 || pion2.motherId() < 0 || kaon.motherId() < 0) { + return K1TruthChannel::None; + } + const int charge = kaon.pdgCode() > 0 ? 1 : -1; + const bool rhoPions = pion1.motherId() == pion2.motherId() && + pion1.motherPDG() == kRho770_0 && pion2.motherPDG() == kRho770_0 && + pion1.pdgCode() == -pion2.pdgCode(); + if (rhoPions && kaon.motherPDG() == charge * o2::constants::physics::Pdg::kK1_1270Plus && + kaon.motherId() != pion1.motherId() && hasSibling(kaon, pion1.motherId())) { + return K1TruthChannel::RhoK; + } + if (matchesKStarPi(pion1, pion2, kaon) || matchesKStarPi(pion2, pion1, kaon)) { + return K1TruthChannel::KStarPi; + } + return K1TruthChannel::None; +} + +// Immediate decay channel of a generated K1 from the PDG codes of its two daughters. +inline K1TruthChannel classifyGeneratedK1(int charge, int daughter1, int daughter2) +{ + if ((daughter1 == kRho770_0 && daughter2 == charge * kKPlus) || + (daughter2 == kRho770_0 && daughter1 == charge * kKPlus)) { + return K1TruthChannel::RhoK; + } + if ((daughter1 == charge * o2::constants::physics::Pdg::kK0Star892 && daughter2 == charge * kPiPlus) || + (daughter2 == charge * o2::constants::physics::Pdg::kK0Star892 && daughter1 == charge * kPiPlus)) { + return K1TruthChannel::KStarPi; + } + return K1TruthChannel::None; +} + +/// K1 selection and candidate enumeration shared by the K1 tasks. +/// The task owns the configurable groups and the histogram registry and passes them in init(). +class K1AnalysisMicroCore +{ + public: + void init(o2::framework::HistogramRegistry& histos, + EventCuts const& eventCuts, TrackCuts const& trackCuts, + PionPidCuts const& pionPidCuts, KaonPidCuts const& kaonPidCuts, + SecondaryCuts const& secondaryCuts, CandidateCuts const& candidateCuts, + ProcessModes const& modes, LooseStageOptions const& looseOptions = {}) + { + mSecondaryCuts = secondaryCuts; + 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.cMaxTPCnSigmaPion.value; + pion.maxTOFnSigma = pionPidCuts.cMaxTOFnSigmaPion.value; + pion.combinedNSigma = pionPidCuts.nsigmaCutCombinedPion.value; + pion.onlyTOFTracks = pionPidCuts.cUseOnlyTOFTrackPi.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 = "cMaxTPCnSigmaPion"; + pion.maxTOFName = "cMaxTOFnSigmaPion"; + pion.tpcCutsName = "cPionTPCNSigmaCuts"; + pion.tofCutsName = "cPionTOFNSigmaCuts"; + auto& kaon = pid[static_cast(Species::Kaon)]; + kaon.species = "Kaon"; + kaon.maxTPCnSigma = kaonPidCuts.cMaxTPCnSigmaKaon.value; + kaon.maxTOFnSigma = kaonPidCuts.cMaxTOFnSigmaKaon.value; + kaon.combinedNSigma = kaonPidCuts.nsigmaCutCombinedKaon.value; + kaon.onlyTOFTracks = kaonPidCuts.cUseOnlyTOFTrackKa.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 = "cMaxTPCnSigmaKaon"; + kaon.maxTOFName = "cMaxTOFnSigmaKaon"; + kaon.tpcCutsName = "cKaonTPCNSigmaCuts"; + kaon.tofCutsName = "cKaonTOFNSigmaCuts"; + mSelection.init(eventCuts, trackCuts, std::move(pid), mCandidateCuts.cByPassTOF, modes.microTracks); + + mSecondaryWindowOn = isCutEnabled(mSecondaryCuts.cSecondaryMasswindow); + mAnotherMassCutOn = isCutEnabled(mSecondaryCuts.cMinAnotherSecondaryMassCut) || isCutEnabled(mSecondaryCuts.cMaxAnotherSecondaryMassCut); + mPiKaMassCutOn = isCutEnabled(mSecondaryCuts.cMinPiKaMassCut) || isCutEnabled(mSecondaryCuts.cMaxPiKaMassCut); + mAngleCutOn = isCutEnabled(mSecondaryCuts.cMinAngle) || isCutEnabled(mSecondaryCuts.cMaxAngle); + mPairAsymCutOn = isCutEnabled(mSecondaryCuts.cMinPairAsym) || isCutEnabled(mSecondaryCuts.cMaxPairAsym); + + registerHistograms(histos, modes); + } + + template + bool passesEventCuts(const CollisionType& collision) + { + return mSelection.passesEventCuts(collision); + } + + template + bool passesMCEventCuts(const CollisionType& collision) + { + return mSelection.passesMCEventCuts(collision); + } + + // Full selection stage of a track (quality, TOF requirement, PID) + template + int trackSelectionStage(const TrackType& track) + { + const int qualityStage = mSelection.trackQualityStage(track); + 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 constexpr (S == Species::Pion) { + if (hasTOF) { + tofNSigma = track.tofNSigmaPi(); + } + } else { + if (hasTOF) { + tofNSigma = track.tofNSigmaKa(); + } + } + if (!mSelection.passesPID(SpeciesIndex, track.pt(), hasTOF, tpcNSigma, tofNSigma)) { + return TrackStage::kTrkTOFRequired; + } + return TrackStage::kTrkPID; + } + + // Unordered (pion, pion, kaon) candidate enumeration of one collision (or one mixed pair of collisions). + // dTracks1: bachelor kaons, dTracks2: pions. The core fills the cut-flow instrumentation; the hooks + // (nullptr to skip) receive + // - onPair(trk1, trk2, secondary, passesPairPt): every pion pair passing the pion selection, + // trk1 being the pion with the lower index, + // - onCandidate(kaon, pion1, pion2, K1CandidateValues): every triplet passing the pion, pair and kaon + // selection, with the canonical pion roles (pion1: opposite sign to the kaon in the unlike-sign case), + // - onExport(collision, kaon, same-sign pion, opposite-sign pion, truth channel, pass bits): the unlike-sign + // micro same-event candidates at the loose or selected stage configured by LooseStageOptions. + // computeValues requests mass13, mass23, angle and pairAsym for every candidate in the rapidity window. + template + void forEachCandidate(o2::framework::HistogramRegistry& histos, const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2, + bool computeValues, PairHook onPair = nullptr, CandidateHook onCandidate = nullptr, ExportHook onExport = nullptr) + { + if (dTracks1.size() == 0 || dTracks2.size() == 0) { + return; + } + constexpr bool HasPairHook = !std::is_same_v; + constexpr bool HasCandidateHook = !std::is_same_v; + constexpr bool HasExportHook = !std::is_same_v; + // Sets are local to this reconstructed collision: IDs cannot leak across DFs. + // Source-file/DF deduplication across split collisions belongs in the audit. + std::array, NTruthChannels> matchedMothers; + + // Selection cache: every track is selected once, not once per pair x bachelor. + // dTracks1: bachelor kaons, dTracks2: pions (different collisions in mixed events). + constexpr bool FillCutFlow = IsResoMicrotrack && !IsMix; + const int64_t firstKaonIndex = dTracks1.begin().index(); + const int64_t firstPionIndex = dTracks2.begin().index(); + const auto kaonSelected = buildSelectionCache(histos, dTracks1, firstKaonIndex); + const auto pionSelected = buildSelectionCache(histos, dTracks2, firstPionIndex); + + // Only micro same-event ML work needs traversal before conventional cuts. + // The canonical-candidate validity is also required before a selected-stage export. + bool visitLoose = false; + bool checkValidity = false; + if constexpr (IsResoMicrotrack && !IsMix) { + visitLoose = mLooseOptions.audit || (!mLooseOptions.exportSelected && HasExportHook); + checkValidity = visitLoose || (mLooseOptions.exportSelected && HasExportHook); + } + std::vector kaonQuality(dTracks1.size(), 0), pionQuality(dTracks2.size(), 0); + if (visitLoose) { + for (const auto& track : dTracks1) { + kaonQuality[getCacheIndex(track, firstKaonIndex, kaonQuality.size())] = mSelection.trackQualityStage(track) == TrackStage::kTrkClusters; + } + for (const auto& track : dTracks2) { + pionQuality[getCacheIndex(track, firstPionIndex, pionQuality.size())] = mSelection.trackQualityStage(track) == TrackStage::kTrkClusters; + } + } + + // Values needed only by switched-on cuts or by the task are computed only then + const bool isK892Mode = mSecondaryCuts.cfgModeK892orRho; + const bool needAngle = computeValues || mAngleCutOn; + const bool needPairAsym = computeValues || mPairAsymCutOn; + // K892 mode: the K* candidate is (trk1, K), rho mode: the rho is (trk1, trk2) + const bool needMass13 = computeValues || (isK892Mode ? mSecondaryWindowOn : mAnotherMassCutOn) || (isK892Mode && (needAngle || needPairAsym)); + const bool needMass23 = computeValues || mPiKaMassCutOn; + const bool rhoWindowOn = mSecondaryWindowOn && !isK892Mode; + + ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResonanceSecondary, lDecayDaughter_bach, lResonanceK1, lPair13, lPair23; + // Unordered pion pairs: each (pion, pion, kaon) triplet is visited once. + // Here trk1 is the pion with the lower index; the roles are assigned once the bachelor is known. + for (const auto& [trk1, trk2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(dTracks2, dTracks2))) { + // trk1: pion, trk2: pion, bTrack: kaon + const bool pionsSelected = pionSelected[getCacheIndex(trk1, firstPionIndex, pionSelected.size())] && pionSelected[getCacheIndex(trk2, firstPionIndex, pionSelected.size())]; + bool pairPt = false; + bool rhoWindow = true; + if (pionsSelected || visitLoose) { + // Resonance reconstruction + lDecayDaughter1.SetCoordinates(trk1.px(), trk1.py(), trk1.pz(), o2::constants::physics::MassPionCharged); + lDecayDaughter2.SetCoordinates(trk2.px(), trk2.py(), trk2.pz(), o2::constants::physics::MassPionCharged); + lResonanceSecondary = lDecayDaughter1 + lDecayDaughter2; + pairPt = !(lResonanceSecondary.Pt() < mSecondaryCuts.cMinSecondaryPtCut); + rhoWindow = !rhoWindowOn || isInWindow(lResonanceSecondary.M(), MassRho770, mSecondaryCuts.cSecondaryMasswindow); + } + if constexpr (FillCutFlow) { + // Early stages count potential triplets: each pair carries N bachelor trials. + // This preserves the pair-first reconstruction and avoids a new cubic data loop. + // Distinct pion IDs are guaranteed by the strictly upper index policy (stage 1 is always passed). + const int lastStage = !pionsSelected ? 1 : !pairPt ? 3 + : !rhoWindow ? 4 + : 5; + for (int stage = 0; stage <= lastStage; ++stage) { + histos.fill(HIST("CutFlow/candidates"), stage, 0, static_cast(dTracks1.size())); + } + if constexpr (IsMC) { + // Match before rejecting quality/pT so both channels have an upstream numerator. + if (std::abs(trk1.pdgCode()) == kPiPlus && trk1.pdgCode() == -trk2.pdgCode()) { + for (const auto& bachelor : dTracks1) { + const auto channel = classifyK1Truth(trk1, trk2, bachelor); + if (channel != K1TruthChannel::None) { + for (int stage = 0; stage <= lastStage; ++stage) { + histos.fill(HIST("CutFlow/candidates"), stage, static_cast(channel)); + } + } + } + } + } + } + if (!pionsSelected && !visitLoose) { + continue; + } + if constexpr (HasPairHook) { + if (pionsSelected) { + onPair(trk1, trk2, lResonanceSecondary, pairPt); + } + } + if (!pairPt && !visitLoose) { + continue; + } + // Secondary mass window (rho mode): the bachelor loop is skipped for rejected pairs + if (!rhoWindow && !visitLoose) { + continue; + } + + for (const auto& bTrack : dTracks1) { + if (bTrack.index() == trk1.index() || bTrack.index() == trk2.index()) { + continue; + } + K1TruthChannel flowChannel = K1TruthChannel::None; + if constexpr (IsMC && IsResoMicrotrack && !IsMix) { + flowChannel = classifyK1Truth(trk1, trk2, bTrack); + } + auto countCandidate = [&](int stage) { + if constexpr (FillCutFlow) { + histos.fill(HIST("CutFlow/candidates"), stage, 0); + if constexpr (IsMC) { + if (flowChannel != K1TruthChannel::None) { + histos.fill(HIST("CutFlow/candidates"), stage, static_cast(flowChannel)); + } + } + } + }; + const bool pairSelected = pionsSelected && pairPt && rhoWindow; + const bool bachelorSelected = kaonSelected[getCacheIndex(bTrack, firstKaonIndex, kaonSelected.size())]; + if (pairSelected) { + countCandidate(6); + } + if ((!pairSelected || !bachelorSelected) && !visitLoose) { + continue; + } + const bool tripletSelected = pairSelected && bachelorSelected; + if (tripletSelected) { + countCandidate(7); + } + + // Canonical assignment of the pion roles, once the bachelor is known. + // Unlike-sign pair: the pion with the sign opposite to the kaon is pion 1 (K*0 partner), the other is pion 2. + // Like-sign pair (the rule is ambiguous): the pion with the lower index is pion 1. + const bool isUnlikeSign = trk1.sign() * trk2.sign() < 0; + const bool swapPions = isUnlikeSign && trk1.sign() == bTrack.sign(); + const auto& pion1 = swapPions ? trk2 : trk1; + const auto& pion2 = swapPions ? trk1 : trk2; + const auto& lPion1 = swapPions ? lDecayDaughter2 : lDecayDaughter1; + const auto& lPion2 = swapPions ? lDecayDaughter1 : lDecayDaughter2; + + // K1 reconstruction + lDecayDaughter_bach.SetCoordinates(bTrack.px(), bTrack.py(), bTrack.pz(), o2::constants::physics::MassKaonCharged); + lResonanceK1 = lResonanceSecondary + lDecayDaughter_bach; + K1CandidateValues values; + values.k1 = lResonanceK1; + values.secondary = lResonanceSecondary; + values.isUnlikeSign = isUnlikeSign; + + auto countMl = [&](int stage) { + if (mLooseOptions.audit && isUnlikeSign) { + histos.fill(HIST("ML/looseCutflow"), stage, 0); + if (flowChannel != K1TruthChannel::None) { + const int stratum = 2 * static_cast(flowChannel) - (bTrack.sign() > 0 ? 1 : 0); + histos.fill(HIST("ML/looseCutflow"), stage, stratum); + } + } + }; + bool validLoose = false; + if constexpr (IsResoMicrotrack && !IsMix) { + if (checkValidity && isUnlikeSign) { + const auto canonical = o2::analysis::k1ml::canonicalizeUS(o2::analysis::k1ml::makeTrackSnapshot(bTrack), o2::analysis::k1ml::makeTrackSnapshot(pion2), o2::analysis::k1ml::makeTrackSnapshot(pion1)); + validLoose = canonical.status == o2::analysis::k1ml::BuildStatus::Ok && + std::isfinite(lResonanceK1.M()) && std::isfinite(lResonanceK1.Pt()) && + std::isfinite(lResonanceK1.Rapidity()) && std::isfinite(lResonanceK1.Eta()) && + std::isfinite(lResonanceK1.Phi()) && + o2::analysis::k1ml::buildMasterFeatures(canonical.candidate).status == o2::analysis::k1ml::BuildStatus::Ok; + if (validLoose) { + countMl(0); + } + } + } + + // Stage L common acceptance uses the existing inclusive rapidity window. + values.inRapidity = lResonanceK1.Rapidity() >= mCandidateCuts.cK1MinRap && lResonanceK1.Rapidity() <= mCandidateCuts.cK1MaxRap; + if (!values.inRapidity) { + if constexpr (HasCandidateHook) { + if (tripletSelected) { + onCandidate(bTrack, pion1, pion2, values); + } + } + continue; + } + if (tripletSelected) { + countCandidate(8); + } + + if (needMass13) { + lPair13 = lPion1 + lDecayDaughter_bach; + values.mass13 = lPair13.M(); + } + if (needMass23) { + lPair23 = lPion2 + lDecayDaughter_bach; + values.mass23 = lPair23.M(); + } + // Rho mode: secondary = (trk1, trk2) against the bachelor. K892 mode: secondary = (trk1, K) against trk2. + if (needAngle) { + values.angle = isK892Mode ? ROOT::Math::VectorUtil::Angle(lPair13, lPion2) : ROOT::Math::VectorUtil::Angle(lResonanceSecondary, lDecayDaughter_bach); + } + if (needPairAsym) { + values.pairAsym = isK892Mode ? (lPair13.E() - lPion2.E()) / (lPair13.E() + lPion2.E()) + : (lResonanceSecondary.E() - lDecayDaughter_bach.E()) / (lResonanceSecondary.E() + lDecayDaughter_bach.E()); + } + + // Candidate cuts (each one is evaluated only if switched on) + values.passesCandidateCuts = + (!isK892Mode || !mSecondaryWindowOn || (isInWindow(values.mass13, o2::constants::physics::MassK0Star892, mSecondaryCuts.cSecondaryMasswindow) && pion1.sign() != bTrack.sign())) && + (!mAnotherMassCutOn || isInRange(isK892Mode ? lResonanceSecondary.M() : values.mass13, mSecondaryCuts.cMinAnotherSecondaryMassCut, mSecondaryCuts.cMaxAnotherSecondaryMassCut)) && + (!mPiKaMassCutOn || isInRange(values.mass23, mSecondaryCuts.cMinPiKaMassCut, mSecondaryCuts.cMaxPiKaMassCut)) && + (!mAngleCutOn || isInRange(values.angle, mSecondaryCuts.cMinAngle, mSecondaryCuts.cMaxAngle)) && + (!mPairAsymCutOn || isInRange(values.pairAsym, mSecondaryCuts.cMinPairAsym, mSecondaryCuts.cMaxPairAsym)); + auto exportCandidate = [&](uint16_t passBits) { + if constexpr (IsResoMicrotrack && !IsMix && HasExportHook) { + onExport(collision, bTrack, pion2, pion1, flowChannel, passBits); + } else { + static_cast(passBits); + } + }; + if (visitLoose && validLoose) { + countMl(1); + if (mLooseOptions.audit) { + histos.fill(HIST("ML/looseMassPtActivity"), lResonanceK1.M(), lResonanceK1.Pt(), collision.cent()); + } + const bool qualityPass = kaonQuality[getCacheIndex(bTrack, firstKaonIndex, kaonQuality.size())] && + pionQuality[getCacheIndex(trk1, firstPionIndex, pionQuality.size())] && + pionQuality[getCacheIndex(trk2, firstPionIndex, pionQuality.size())]; + uint16_t passBits = kPassLoose; + if (qualityPass) { + passBits |= kPassQuality; + countMl(2); + if (pionsSelected && bachelorSelected) { + passBits |= kPassPID; + countMl(3); + if (pairPt && rhoWindow) { + passBits |= kPassPair; + countMl(4); + if (values.passesCandidateCuts) { + passBits |= kPassCandidate; + countMl(5); + } + } + } + } + if (!mLooseOptions.exportSelected) { + exportCandidate(passBits); + } + } + // Stage C retains the frozen conventional selections. + if (!tripletSelected) { + continue; + } + if (values.passesCandidateCuts) { + countCandidate(9); + countCandidate(isUnlikeSign ? 10 : 11); + if (isUnlikeSign && mLooseOptions.exportSelected && validLoose) { + exportCandidate(PassBitsSelected); + } + if constexpr (IsMC && IsResoMicrotrack && !IsMix) { + if (flowChannel != K1TruthChannel::None) { + const int mother = flowChannel == K1TruthChannel::RhoK ? bTrack.motherId() : std::abs(pion1.motherPDG()) == o2::constants::physics::Pdg::kK1_1270Plus ? pion1.motherId() + : pion2.motherId(); + if (matchedMothers[static_cast(flowChannel)].insert(mother).second) { + histos.fill(HIST("CutFlow/uniqueMothersPerCollision"), static_cast(flowChannel)); + } + } + } + } + if constexpr (HasCandidateHook) { + onCandidate(bTrack, pion1, pion2, values); + } + } // bTrack + } + } // forEachCandidate + + // Generated K1 parents of a selected reconstructed MC collision. The optional callback receives + // (parent, immediate channel) for the parents inside the K1 rapidity window. + // Parents belong to selected reconstructed events; split reco collisions + // repeat parent sets. This is not an unconditional generated denominator. + template + void forEachGeneratedK1(o2::framework::HistogramRegistry& histos, const ParentsType& resoParents, Callback callback = nullptr) + { + for (const auto& part : resoParents) { + if (std::abs(part.pdgCode()) != o2::constants::physics::Pdg::kK1_1270Plus) { + continue; + } + const int charge = part.pdgCode() > 0 ? 1 : -1; + const K1TruthChannel channel = classifyGeneratedK1(charge, part.daughterPDG1(), part.daughterPDG2()); + histos.fill(HIST("CutFlow/generated"), 0, static_cast(channel)); + if (part.y() < mCandidateCuts.cK1MinRap || part.y() > mCandidateCuts.cK1MaxRap) { + continue; + } + histos.fill(HIST("CutFlow/generated"), 1, static_cast(channel)); + if constexpr (!std::is_same_v) { + callback(part, channel); + } + } + } + + private: + // Selection cache of one track slice. The row number of a grouped slice is global, hence the offset. + 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); + } + + template + std::vector buildSelectionCache(o2::framework::HistogramRegistry& histos, const TracksType& tracks, int64_t firstIndex) + { + std::vector selected(tracks.size(), 0); + for (const auto& track : tracks) { + const int stage = trackSelectionStage(track); + selected[getCacheIndex(track, firstIndex, selected.size())] = (stage == TrackStage::kTrkPID) ? 1 : 0; + if constexpr (FillCutFlow) { + for (int i = 0; i <= stage; ++i) { + histos.fill(HIST("CutFlow/tracks"), i, static_cast(S)); + } + } + } + return selected; + } + + // Cut-flow instrumentation of the selection; the output histograms belong to the tasks. + void registerHistograms(o2::framework::HistogramRegistry& histos, ProcessModes const& modes) + { + using o2::framework::HistType; + if (mLooseOptions.audit) { + auto flow = histos.add("ML/looseCutflow", "US triplets;stage;signal stratum", HistType::kTH2D, {{6, -0.5, 5.5}, {5, -0.5, 4.5}}); + const std::array labels{"structural US", "loose acceptance", "track quality", "TOF + PID", "pair requirements", "selected US"}; + const std::array strata{"all US", "rhoK+", "rhoK-", "KstarPi+", "KstarPi-"}; + for (std::size_t i = 0; i < labels.size(); ++i) { + flow->GetXaxis()->SetBinLabel(i + 1, labels[i]); + } + for (std::size_t i = 0; i < strata.size(); ++i) { + flow->GetYaxis()->SetBinLabel(i + 1, strata[i]); + } + histos.add("ML/looseMassPtActivity", "Loose US;mass (GeV/c^{2});pT (GeV/c);FT0M percentile", HistType::kTH3D, + {{300, 0.7, 3.7}, {{0., 0.5, 1., 2., 3., 5., 8., 15., 30., 100.}, "pT"}, {{0., 10., 30., 50., 70., 100., 110.}, "FT0M percentile"}}); + } + + // Micro-only instrumentation: category 0 includes all combinations, not just unmatched. + constexpr int NTrackStages = TrackStage::kTrkNStages; + auto trackFlow = histos.add("CutFlow/tracks", "Micro 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/candidates", "Unordered micro triplets;stage;category", HistType::kTH2D, {{NCandidateStages, -0.5, NCandidateStages - 0.5}, {3, -0.5, 2.5}}); + const std::array candidateLabels{"input unordered triplets", "distinct pion IDs", "pion selection (quality+PID)", "pion pair constructed", "pair pT", "secondary mass window (rho mode)", "three distinct IDs", "kaon selection (quality+PID)", "K1 rapidity", "candidate cuts", "final US", "final LS"}; + 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, "rho K"); + candidateFlow->GetYaxis()->SetBinLabel(3, "K* pi"); + if (modes.mcRecoMicro) { + auto mothers = histos.add("CutFlow/uniqueMothersPerCollision", "Final unique K1 IDs summed over reconstructed collisions (not globally deduplicated)", HistType::kTH1D, {{2, 0.5, 2.5}}); + mothers->GetXaxis()->SetBinLabel(1, "rho K"); + mothers->GetXaxis()->SetBinLabel(2, "K* pi"); + } + if (modes.mcGen) { + auto generated = histos.add("CutFlow/generated", "K1 parent rows conditional on selected reconstructed events;stage;immediate channel", HistType::kTH2D, {{2, -0.5, 1.5}, {3, -0.5, 2.5}}); + generated->GetXaxis()->SetBinLabel(1, "all K1 parent rows"); + generated->GetXaxis()->SetBinLabel(2, "K1 rapidity window"); + generated->GetYaxis()->SetBinLabel(1, "other / unresolved"); + generated->GetYaxis()->SetBinLabel(2, "rho K"); + generated->GetYaxis()->SetBinLabel(3, "K* pi"); + } + } + + ResoAnalysisSelectionCore mSelection; + SecondaryCuts mSecondaryCuts; + CandidateCuts mCandidateCuts; + LooseStageOptions mLooseOptions; + + // Derived once in init(): which candidate cuts are switched on. + bool mSecondaryWindowOn = false; + bool mAnotherMassCutOn = false; + bool mPiKaMassCutOn = false; + bool mAngleCutOn = false; + bool mPairAsymCutOn = false; +}; + +} // namespace o2::analysis::k1micro + +#endif // PWGLF_CORE_K1ANALYSISMICROCORE_H_ diff --git a/PWGLF/Core/K1MlFeatures.h b/PWGLF/Core/K1MlFeatures.h new file mode 100644 index 00000000000..262fb0feff7 --- /dev/null +++ b/PWGLF/Core/K1MlFeatures.h @@ -0,0 +1,873 @@ +// 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 K1MlFeatures.h +/// \brief Canonical micro001 K1 candidate and feature contract helpers +/// \author Bong-Hwi Lim +/// + +#ifndef PWGLF_CORE_K1MLFEATURES_H_ +#define PWGLF_CORE_K1MLFEATURES_H_ + +#include "PWGLF/DataModel/LFResonanceTables.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::k1ml +{ +inline constexpr std::size_t NItsLayers = 7; +inline constexpr std::size_t NCandidateTracks = 3; +inline constexpr std::size_t NMasterFeatures = 125; + +enum class Role : uint8_t { Kaon, + PionSame, + PionOpp }; +enum class Profile : uint8_t { DetectorV1, + RelationalV1, + SubstructureV1 }; +enum class BuildStatus : uint8_t { + Ok, + InvalidChargePattern, + ReusedTrack, + InvalidMomentum, + InvalidKinematics, + InvalidContract +}; + +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 CandidateSnapshot { + std::array tracks{}; // K, same-charge pion, opposite-charge pion +}; + +struct CanonicalizationResult { + CandidateSnapshot candidate{}; + BuildStatus status = BuildStatus::InvalidChargePattern; + explicit operator bool() const { return status == BuildStatus::Ok; } +}; + +struct KinematicAudit { + float massKPiPi = 0.f; + float massPiPi = 0.f; + float massKaonPiSame = 0.f; + float massKaonPiOpp = 0.f; + float candidatePt = 0.f; + float candidateEta = 0.f; + float candidatePhi = 0.f; + float scalarSumPt = 0.f; +}; + +struct FeaturePack { + std::array master{}; + BuildStatus status = BuildStatus::InvalidContract; + KinematicAudit kinematics{}; +}; + +// Adapt the real ResoMicroTracks_001 public columns and dynamic accessors. +// Keeping the raw bytes alongside decoded accessors makes the encoding auditable. +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(TrackSnapshot const& track) +{ + return std::isfinite(track.px) && std::isfinite(track.py) && std::isfinite(track.pz) && + std::hypot(track.px, track.py) > 0.f; +} + +inline CanonicalizationResult canonicalizeUS(TrackSnapshot const& kaon, + TrackSnapshot const& pionA, + TrackSnapshot const& pionB) +{ + CanonicalizationResult result; + if (kaon.sourceTrackId == pionA.sourceTrackId || kaon.sourceTrackId == pionB.sourceTrackId || + pionA.sourceTrackId == pionB.sourceTrackId) { + result.status = BuildStatus::ReusedTrack; + return result; + } + if ((kaon.charge != 1 && kaon.charge != -1) || (pionA.charge != 1 && pionA.charge != -1) || + (pionB.charge != 1 && pionB.charge != -1) || pionA.charge == pionB.charge || + (pionA.charge != kaon.charge && pionB.charge != kaon.charge)) { + result.status = BuildStatus::InvalidChargePattern; + return result; + } + if (!hasFiniteMomentum(kaon) || !hasFiniteMomentum(pionA) || !hasFiniteMomentum(pionB)) { + result.status = BuildStatus::InvalidMomentum; + return result; + } + const auto& same = pionA.charge == kaon.charge ? pionA : pionB; + const auto& opposite = pionA.charge == kaon.charge ? pionB : pionA; + result.candidate.tracks = {kaon, same, opposite}; + result.status = BuildStatus::Ok; + return result; +} + +namespace detail +{ +struct EncodedValue { + float value; + float valid; + float overflow; +}; + +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) ? -3.5f : 3.5f, .valid = 1.f, .overflow = 1.f}; + } + return {.value = decoded, .valid = 1.f, .overflow = 0.f}; +} + +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}; +} + +struct Kinematics { + double pt; + double eta; + double phi; + double energy; +}; + +inline bool getKinematics(TrackSnapshot const& track, double mass, Kinematics& out) +{ + if (!hasFiniteMomentum(track)) { + return false; + } + const double px = track.px, py = track.py, pz = track.pz; + out.pt = std::hypot(px, py); + const double p = std::hypot(out.pt, pz); + out.eta = std::asinh(pz / out.pt); + out.phi = std::atan2(py, px); + out.energy = std::sqrt(p * p + mass * mass); + return std::isfinite(out.pt) && std::isfinite(out.eta) && std::isfinite(out.phi) && std::isfinite(out.energy); +} + +inline float invariantMass(Kinematics const& a, Kinematics const& b) +{ + const auto pxa = a.pt * std::cos(a.phi), pya = a.pt * std::sin(a.phi); + const auto pxb = b.pt * std::cos(b.phi), pyb = b.pt * std::sin(b.phi); + const double e = a.energy + b.energy; + const double px = pxa + pxb, py = pya + pyb; + // pz is recovered from pt*sinh(eta), matching the stored three-momentum. + const double pz = a.pt * std::sinh(a.eta) + b.pt * std::sinh(b.eta); + const double m2 = e * e - px * px - py * py - pz * pz; + return static_cast(std::sqrt(std::max(0.0, m2))); +} + +inline float invariantMass(std::array const& tracks, + std::array const& indices) +{ + double e = 0., px = 0., py = 0., pz = 0.; + for (const auto& i : indices) { + const auto& track = tracks[i]; + e += track.energy; + px += track.pt * std::cos(track.phi); + py += track.pt * std::sin(track.phi); + pz += track.pt * std::sinh(track.eta); + } + return static_cast(std::sqrt(std::max(0.0, e * e - px * px - py * py - pz * pz))); +} + +inline void append(std::array& out, std::size_t& index, EncodedValue value) +{ + out[index++] = value.value; + out[index++] = value.valid; + out[index++] = value.overflow; +} + +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 and projection indices are generated from feature_contract_v1.json. +inline constexpr std::array MasterFeatureNames{ + "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", + "kaon.pt_fraction", + "pion_same.tpc_nsigma_pi", + "pion_same.tpc_nsigma_pi_valid", + "pion_same.tpc_nsigma_pi_overflow", + "pion_same.tpc_nsigma_ka", + "pion_same.tpc_nsigma_ka_valid", + "pion_same.tpc_nsigma_ka_overflow", + "pion_same.tpc_nsigma_pr", + "pion_same.tpc_nsigma_pr_valid", + "pion_same.tpc_nsigma_pr_overflow", + "pion_same.tof_nsigma_pi", + "pion_same.tof_nsigma_pi_valid", + "pion_same.tof_nsigma_pi_overflow", + "pion_same.tof_nsigma_ka", + "pion_same.tof_nsigma_ka_valid", + "pion_same.tof_nsigma_ka_overflow", + "pion_same.tof_nsigma_pr", + "pion_same.tof_nsigma_pr_valid", + "pion_same.tof_nsigma_pr_overflow", + "pion_same.abs_dca_xy", + "pion_same.abs_dca_xy_valid", + "pion_same.abs_dca_xy_overflow", + "pion_same.abs_dca_z", + "pion_same.abs_dca_z_valid", + "pion_same.abs_dca_z_overflow", + "pion_same.passed_ptdep_dca_xy", + "pion_same.passed_ptdep_dca_z", + "pion_same.has_tof", + "pion_same.tpc_crossed_rows", + "pion_same.its_hit_l0", + "pion_same.its_hit_l1", + "pion_same.its_hit_l2", + "pion_same.its_hit_l3", + "pion_same.its_hit_l4", + "pion_same.its_hit_l5", + "pion_same.its_hit_l6", + "pion_same.is_pv_contributor", + "pion_same.pt_fraction", + "pion_opp.tpc_nsigma_pi", + "pion_opp.tpc_nsigma_pi_valid", + "pion_opp.tpc_nsigma_pi_overflow", + "pion_opp.tpc_nsigma_ka", + "pion_opp.tpc_nsigma_ka_valid", + "pion_opp.tpc_nsigma_ka_overflow", + "pion_opp.tpc_nsigma_pr", + "pion_opp.tpc_nsigma_pr_valid", + "pion_opp.tpc_nsigma_pr_overflow", + "pion_opp.tof_nsigma_pi", + "pion_opp.tof_nsigma_pi_valid", + "pion_opp.tof_nsigma_pi_overflow", + "pion_opp.tof_nsigma_ka", + "pion_opp.tof_nsigma_ka_valid", + "pion_opp.tof_nsigma_ka_overflow", + "pion_opp.tof_nsigma_pr", + "pion_opp.tof_nsigma_pr_valid", + "pion_opp.tof_nsigma_pr_overflow", + "pion_opp.abs_dca_xy", + "pion_opp.abs_dca_xy_valid", + "pion_opp.abs_dca_xy_overflow", + "pion_opp.abs_dca_z", + "pion_opp.abs_dca_z_valid", + "pion_opp.abs_dca_z_overflow", + "pion_opp.passed_ptdep_dca_xy", + "pion_opp.passed_ptdep_dca_z", + "pion_opp.has_tof", + "pion_opp.tpc_crossed_rows", + "pion_opp.its_hit_l0", + "pion_opp.its_hit_l1", + "pion_opp.its_hit_l2", + "pion_opp.its_hit_l3", + "pion_opp.its_hit_l4", + "pion_opp.its_hit_l5", + "pion_opp.its_hit_l6", + "pion_opp.is_pv_contributor", + "pion_opp.pt_fraction", + "kaon__pion_same.delta_eta", + "kaon__pion_same.sin_delta_phi", + "kaon__pion_same.cos_delta_phi", + "kaon__pion_same.z_pt", + "kaon__pion_opp.delta_eta", + "kaon__pion_opp.sin_delta_phi", + "kaon__pion_opp.cos_delta_phi", + "kaon__pion_opp.z_pt", + "pion_same__pion_opp.delta_eta", + "pion_same__pion_opp.sin_delta_phi", + "pion_same__pion_opp.cos_delta_phi", + "pion_same__pion_opp.z_pt", + "mass_pi_pi", + "mass_kaon_pion_opp"}; +inline constexpr std::array DetectorV1Projection{ + 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, + 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, + 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}; +inline constexpr std::array RelationalV1Projection{ + 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}; +inline constexpr std::array SubstructureV1Projection{ + 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}; +static_assert(MasterFeatureNames.size() == NMasterFeatures, "master feature name count must match NMasterFeatures"); +static_assert(detail::isStrictlyIncreasingBelow(DetectorV1Projection, NMasterFeatures), "DetectorV1 projection indices must be strictly increasing and below NMasterFeatures"); +static_assert(detail::isStrictlyIncreasingBelow(RelationalV1Projection, NMasterFeatures), "RelationalV1 projection indices must be strictly increasing and below NMasterFeatures"); +static_assert(detail::isStrictlyIncreasingBelow(SubstructureV1Projection, NMasterFeatures), "SubstructureV1 projection indices must be strictly increasing and below NMasterFeatures"); +inline constexpr std::string_view FeatureContractSha256 = "39f38ece001581d8ebf57392fad045759a63ba49f57c411b9b566d7c1cc58b8a"; + +inline FeaturePack buildMasterFeatures(CandidateSnapshot const& candidate) +{ + FeaturePack pack; + for (auto const& track : candidate.tracks) { + if ((track.charge != 1 && track.charge != -1) || !hasFiniteMomentum(track)) { + pack.status = BuildStatus::InvalidMomentum; + return pack; + } + } + if (candidate.tracks[0].sourceTrackId == candidate.tracks[1].sourceTrackId || + candidate.tracks[0].sourceTrackId == candidate.tracks[2].sourceTrackId || + candidate.tracks[1].sourceTrackId == candidate.tracks[2].sourceTrackId || + candidate.tracks[0].charge != candidate.tracks[1].charge || + candidate.tracks[0].charge == candidate.tracks[2].charge) { + pack.status = BuildStatus::InvalidChargePattern; + return pack; + } + + std::array kin{}; + for (std::size_t i = 0; i < NCandidateTracks; ++i) { + if (!detail::getKinematics(candidate.tracks[i], i == 0 ? o2::constants::physics::MassKaonCharged : o2::constants::physics::MassPionCharged, kin[i])) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + } + const double sumPt = kin[0].pt + kin[1].pt + kin[2].pt; + if (!std::isfinite(sumPt) || sumPt <= 0.) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + const double totalPx = static_cast(candidate.tracks[0].px) + candidate.tracks[1].px + candidate.tracks[2].px; + const double totalPy = static_cast(candidate.tracks[0].py) + candidate.tracks[1].py + candidate.tracks[2].py; + const double totalPz = static_cast(candidate.tracks[0].pz) + candidate.tracks[1].pz + candidate.tracks[2].pz; + const double candidatePt = std::hypot(totalPx, totalPy); + if (!std::isfinite(candidatePt)) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + pack.kinematics.scalarSumPt = static_cast(sumPt); + pack.kinematics.candidatePt = static_cast(candidatePt); + pack.kinematics.candidateEta = candidatePt > 0. ? static_cast(std::asinh(totalPz / candidatePt)) : 0.f; + pack.kinematics.candidatePhi = static_cast(std::atan2(totalPy, totalPx)); + pack.kinematics.massKPiPi = detail::invariantMass(kin, std::array{0, 1, 2}); + pack.kinematics.massPiPi = detail::invariantMass(kin[1], kin[2]); + pack.kinematics.massKaonPiSame = detail::invariantMass(kin[0], kin[1]); + pack.kinematics.massKaonPiOpp = detail::invariantMass(kin[0], kin[2]); + + std::size_t index = 0; + for (std::size_t i = 0; i < NCandidateTracks; ++i) { + const auto& track = candidate.tracks[i]; + for (const float& decoded : track.tpcNSigma) { + detail::append(pack.master, index, detail::encodePID(decoded)); + } + for (const float& decoded : track.tofNSigma) { + detail::append(pack.master, index, track.hasTOF ? detail::encodePID(decoded) : detail::EncodedValue{.value = 0.f, .valid = 0.f, .overflow = 0.f}); + } + detail::append(pack.master, index, detail::encodeDCA(track.dcaXY)); + detail::append(pack.master, index, detail::encodeDCA(track.dcaZ)); + pack.master[index++] = track.passedPtDependentDCAxy ? 1.f : 0.f; + pack.master[index++] = track.passedPtDependentDCAz ? 1.f : 0.f; + pack.master[index++] = track.hasTOF ? 1.f : 0.f; + pack.master[index++] = static_cast(track.tpcNClsCrossedRows); + for (const bool& hit : track.itsHit) { + pack.master[index++] = hit ? 1.f : 0.f; + } + pack.master[index++] = track.isPVContributor ? 1.f : 0.f; + pack.master[index++] = static_cast(kin[i].pt / sumPt); + } + + constexpr std::array, 3> PairIndices{{{{0, 1}}, {{0, 2}}, {{1, 2}}}}; + for (auto const& pair : PairIndices) { + const auto i = pair[0], j = pair[1]; + const double dphi = kin[i].phi - kin[j].phi; + const double zDenominator = kin[i].pt + kin[j].pt; + if (!std::isfinite(dphi) || zDenominator <= 0.) { + pack.status = BuildStatus::InvalidKinematics; + return pack; + } + pack.master[index++] = static_cast(kin[i].eta - kin[j].eta); + pack.master[index++] = static_cast(std::sin(dphi)); + pack.master[index++] = static_cast(std::cos(dphi)); + pack.master[index++] = static_cast(std::min(kin[i].pt, kin[j].pt) / zDenominator); + } + pack.master[index++] = pack.kinematics.massPiPi; + pack.master[index++] = pack.kinematics.massKaonPiOpp; + const bool allFinite = std::all_of(pack.master.begin(), pack.master.end(), [](float x) { return std::isfinite(x); }); + if (index != pack.master.size() || !allFinite) { + pack.status = BuildStatus::InvalidContract; + return pack; + } + pack.status = BuildStatus::Ok; + return pack; +} + +inline std::vector projectFeatures(FeaturePack const& pack, Profile profile) +{ + if (pack.status != BuildStatus::Ok) { + return {}; + } + std::vector projected; + switch (profile) { + case Profile::DetectorV1: + projected.reserve(DetectorV1Projection.size()); + for (const auto& i : DetectorV1Projection) { + projected.push_back(pack.master[i]); + } + break; + case Profile::RelationalV1: + projected.reserve(RelationalV1Projection.size()); + for (const auto& i : RelationalV1Projection) { + projected.push_back(pack.master[i]); + } + break; + case Profile::SubstructureV1: + projected.reserve(SubstructureV1Projection.size()); + for (const auto& i : SubstructureV1Projection) { + projected.push_back(pack.master[i]); + } + break; + default: + return {}; + } + return projected; +} +} // namespace o2::analysis::k1ml + +#endif // PWGLF_CORE_K1MLFEATURES_H_ diff --git a/PWGLF/Core/ResoAnalysisSelectionCore.h b/PWGLF/Core/ResoAnalysisSelectionCore.h new file mode 100644 index 00000000000..225f65124b7 --- /dev/null +++ b/PWGLF/Core/ResoAnalysisSelectionCore.h @@ -0,0 +1,462 @@ +// 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 ResoAnalysisSelectionCore.h +/// \brief Event, track-quality and PID selection of resonance daughters from the reduced v001 resonance tables +/// \author Bong-Hwi Lim +/// +/// The same selection code serves full ResoTracks (exact values) and ResoMicroTracks (quantised DCA and nSigma, +/// see LFResonanceTables.h); only the comparisons are quantisation aware. A cut set to DisabledCut is off. + +#ifndef PWGLF_CORE_RESOANALYSISSELECTIONCORE_H_ +#define PWGLF_CORE_RESOANALYSISSELECTIONCORE_H_ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::resonance +{ + +inline constexpr float DisabledCut = -999.f; // an optional cut with this value is off and not evaluated +inline constexpr double DCAGridStep = 0.025; // v001 micro DCA encoding, lower-inclusive bins up to DCAGridMax +inline constexpr double DCAGridMax = 0.15; +inline constexpr double PIDGridStart = 2.0; // v001 micro nSigma encoding: 0.25 bins in [2.0, 3.5] +inline constexpr double PIDGridStep = 0.25; +inline constexpr double PIDGridMax = 3.5; +inline constexpr double GridTolerance = 1e-4; +inline constexpr std::size_t MinPtBinEdges = 2; // a pT dependent PID table needs at least one bin +inline constexpr float ProducerDCAPtP0 = 0.004f; // resonanceModuleInitializer cfgTightDCAOffset default +inline constexpr float ProducerDCAPtCoeff = 0.013f; // resonanceModuleInitializer cfgTightDCAPtCoefficient default +inline constexpr float ProducerDCAPtPower = 1.f; // resonanceModuleInitializer cfgTightDCAPtPower default +inline constexpr float ConfigTolerance = 1e-6f; + +// Last stage passed by a track; a cut-flow histogram can be filled directly from this value. +enum TrackStage : int { + kTrkInput = 0, + kTrkPt, + kTrkEta, + kTrkDCAxy, + kTrkDCAz, + kTrkFlags, + kTrkClusters, + kTrkTOFRequired, + kTrkPID, + kTrkNStages +}; + +// Resolved PID cut of one species at a given pT. +struct PIDCut { + double tpcMax = 0.; + double tofMax = 0.; + double combined = 0.; + bool tofRequired = false; +}; + +// A cut is on unless it carries the disabled value (tolerant to the float parsing of the JSON value). +inline bool isCutEnabled(float value) +{ + return value > DisabledCut + 1.f; +} + +// v001 micro values are lower-inclusive bin edges: a maximum cut on the grid keeps bins below it. +inline bool passesBinnedMax(double decoded, double cut) +{ + return decoded < cut - o2::constants::math::Epsilon; +} + +// Minimum cut on the grid keeps the bin starting at the cut. +inline bool passesBinnedMin(double decoded, double cut) +{ + return decoded >= cut - o2::constants::math::Epsilon; +} + +template +bool passesMax(double value, double cut) +{ + if constexpr (IsResoMicrotrack) { + return passesBinnedMax(value, cut); + } else { + return value < cut; + } +} + +inline bool isInRange(double value, double minimum, double maximum) +{ + if (isCutEnabled(minimum) && value < minimum) { + return false; + } + if (isCutEnabled(maximum) && value > maximum) { + return false; + } + return true; +} + +inline bool isInWindow(double value, double center, double width) +{ + return std::abs(value - center) < width; +} + +// Preserve pT-bin membership [low, high). +inline int getPtBinIndex(float pt, const std::vector& ptBins) +{ + for (std::size_t i = 1; i < ptBins.size(); ++i) { + if (pt >= ptBins[i - 1] && pt < ptBins[i]) { + return static_cast(i - 1); + } + } + return -1; +} + +// Configurable groups without prefix: the JSON keys are the plain configurable names. + +/// Event selection +struct EventCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cRecoINELgt0{"cRecoINELgt0", false, "Apply reconstructed INEL>0 selection"}; + o2::framework::Configurable cMCINELgt0{"cMCINELgt0", false, "Require generator INEL>0 in MC processes"}; + o2::framework::Configurable cMCVtxIn10{"cMCVtxIn10", false, "Require generator |vz| < 10 cm in MC processes"}; +}; + +/// Track selections (common for all daughter species, -999 switches an optional cut off) +struct TrackCuts : o2::framework::ConfigurableGroup { + o2::framework::Configurable cMinPtcut{"cMinPtcut", 0.15, "Track minium pt cut"}; + o2::framework::Configurable cMaxEtacut{"cMaxEtacut", -999.f, "Track maximum |eta| cut (-999: off)"}; + // DCAr to PV + o2::framework::Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.1, "Track DCAr cut to PV Maximum"}; + // DCAz to PV + o2::framework::Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 0.1, "Track DCAz cut to PV Maximum"}; + o2::framework::Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; + o2::framework::Configurable cfgUsePtDepDCA{"cfgUsePtDepDCA", false, "Use pT dependent DCA cut instead of the fixed maximum"}; + o2::framework::Configurable cDCAToPVByPtP0{"cDCAToPVByPtP0", 0.004f, "pT dependent DCA cut = P0 + coefficient / pT^power (cm)"}; + o2::framework::Configurable cDCAToPVByPtCoeff{"cDCAToPVByPtCoeff", 0.013f, "Coefficient in the pT dependent DCA cut"}; + o2::framework::Configurable cDCAToPVByPtPower{"cDCAToPVByPtPower", 1.f, "Power in the pT dependent DCA cut"}; + o2::framework::Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + o2::framework::Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + o2::framework::Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + o2::framework::Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor + o2::framework::Configurable cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; + o2::framework::Configurable cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"}; + o2::framework::Configurable cfgTPCcluster{"cfgTPCcluster", 0, "Number of TPC cluster (found clusters, ResoTracks only)"}; + o2::framework::Configurable cfgTPCCrossedRowsMin{"cfgTPCCrossedRowsMin", 0, "Minimum number of TPC crossed rows"}; + o2::framework::Configurable cfgITSNClsMin{"cfgITSNClsMin", 0, "Minimum number of ITS clusters (ResoMicroTracks only)"}; + o2::framework::Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; +}; + +/// PID selection of one daughter species, filled by the task from its own (species-named) configurables. +/// The configurable names are used only in the configuration messages. +struct PIDCutConfig { + std::string species; // e.g. "Pion", used in messages + double maxTPCnSigma = DisabledCut; + double maxTOFnSigma = DisabledCut; + double combinedNSigma = DisabledCut; // combined TPC-TOF cut, on when > 0 + bool onlyTOFTracks = false; // require a TOF signal + bool usePtDependent = false; // use the pT binned cuts below instead of the fixed maxima + std::vector ptBins; // bin edges; the other vectors have one entry per bin + std::vector tpcNSigmaCuts; + std::vector tofNSigmaCuts; + std::vector tofRequired; + std::string maxTPCName; // configurable names for the messages + std::string maxTOFName; + std::string tpcCutsName; + std::string tofCutsName; +}; + +/// Event, track-quality, TOF-requirement and PID selection of resonance daughters. +/// The species index is the position of its PIDCutConfig in init(). +class ResoAnalysisSelectionCore +{ + public: + // byPassTOF skips the TOF nSigma cut and the pT binned TOF requirement; microTracks enables the grid checks. + void init(EventCuts const& eventCuts, TrackCuts const& trackCuts, std::vector pidCuts, bool byPassTOF, bool microTracks) + { + mEventCuts = eventCuts; + mTrackCuts = trackCuts; + mPID = std::move(pidCuts); + mByPassTOF = byPassTOF; + checkConfiguration(microTracks); + } + + template + bool passesEventCuts(const CollisionType& collision) + { + return !(mEventCuts.cRecoINELgt0 && !collision.isRecINELgt0()); + } + + template + bool passesMCEventCuts(const CollisionType& collision) + { + if (mEventCuts.cMCINELgt0 && !collision.isINELgt0()) { + return false; + } + if (mEventCuts.cMCVtxIn10 && !collision.isVtxIn10()) { + return false; + } + return true; + } + + // Resolve the PID cut of one species at a given pT; false if the pT is outside all pT-dependent bins. + bool getPIDCut(int species, float pt, PIDCut& cut) + { + const auto& config = mPID[species]; + cut.tpcMax = config.maxTPCnSigma; + cut.tofMax = config.maxTOFnSigma; + cut.combined = config.combinedNSigma; + cut.tofRequired = false; + if (config.usePtDependent) { + const int ptBin = getPtBinIndex(pt, config.ptBins); + if (ptBin < 0) { + return false; + } + const auto bin = static_cast(ptBin); + cut.tpcMax = config.tpcNSigmaCuts[bin]; + cut.tofMax = config.tofNSigmaCuts[bin]; + cut.tofRequired = config.tofRequired[bin] != 0; + } + return true; + } + + // Track quality selection shared by all species. Returns the last stage that was passed. + // Full tracks store exact values; micro tracks store quantised DCA (see LFResonanceTables.h). + template + int trackQualityStage(const TrackType& track) + { + const double pt = track.pt(); + const double dcaXY = track.dcaXY(); + const double dcaZ = track.dcaZ(); + // Invalid micro DCA codes decode to NaN + if (!std::isfinite(pt) || !std::isfinite(track.eta()) || !std::isfinite(dcaXY) || !std::isfinite(dcaZ)) { + return kTrkInput; + } + if (std::abs(pt) < mTrackCuts.cMinPtcut) { + return kTrkInput; + } + if (isCutEnabled(mTrackCuts.cMaxEtacut) && !(std::abs(track.eta()) < mTrackCuts.cMaxEtacut)) { + return kTrkPt; + } + + if (mTrackCuts.cfgUsePtDepDCA) { + if constexpr (IsResoMicrotrack) { + if (!track.passedPtDependentDCAxy()) { + return kTrkEta; + } + if (!track.passedPtDependentDCAz()) { + return kTrkDCAxy; + } + } else { + const double dcaPtCut = mTrackCuts.cDCAToPVByPtP0 + mTrackCuts.cDCAToPVByPtCoeff * std::pow(pt, -static_cast(mTrackCuts.cDCAToPVByPtPower)); + if (!(std::abs(dcaXY) < dcaPtCut)) { + return kTrkEta; + } + if (!(std::abs(dcaZ) < dcaPtCut)) { + return kTrkDCAxy; + } + } + } else { + if (isCutEnabled(mTrackCuts.cMaxDCArToPVcut)) { + if constexpr (IsResoMicrotrack) { + if (!passesBinnedMax(dcaXY, mTrackCuts.cMaxDCArToPVcut)) { + return kTrkEta; + } + } else { + if (!(std::abs(dcaXY) <= mTrackCuts.cMaxDCArToPVcut)) { + return kTrkEta; + } + } + } + if (isCutEnabled(mTrackCuts.cMaxDCAzToPVcut)) { + if constexpr (IsResoMicrotrack) { + if (!passesBinnedMax(dcaZ, mTrackCuts.cMaxDCAzToPVcut)) { + return kTrkDCAxy; + } + } else { + if (!(std::abs(dcaZ) <= mTrackCuts.cMaxDCAzToPVcut)) { + return kTrkDCAxy; + } + } + } + } + if (isCutEnabled(mTrackCuts.cMinDCAzToPVcut)) { + if constexpr (IsResoMicrotrack) { + if (!passesBinnedMin(dcaZ, mTrackCuts.cMinDCAzToPVcut)) { + return kTrkDCAxy; + } + } else { + if (!(std::abs(dcaZ) >= mTrackCuts.cMinDCAzToPVcut)) { + return kTrkDCAxy; + } + } + } + + // Track flags + if ((mTrackCuts.cfgPrimaryTrack && !track.isPrimaryTrack()) || + (mTrackCuts.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) || + (mTrackCuts.cfgGlobalTrack && !track.isGlobalTrack()) || + (mTrackCuts.cfgPVContributor && !track.isPVContributor()) || + (mTrackCuts.cfgUseITSRefit && !track.passedITSRefit()) || + (mTrackCuts.cfgUseTPCRefit && !track.passedTPCRefit())) { + return kTrkDCAz; + } + + // Clusters: found clusters exist only in ResoTracks, ITS clusters only in ResoMicroTracks + if constexpr (!IsResoMicrotrack) { + if constexpr (requires { track.tpcNClsFound(); }) { + if (track.tpcNClsFound() < mTrackCuts.cfgTPCcluster) { + return kTrkFlags; + } + } + } + if constexpr (requires { track.tpcNClsCrossedRows(); }) { + if (track.tpcNClsCrossedRows() < mTrackCuts.cfgTPCCrossedRowsMin) { + return kTrkFlags; + } + } + if constexpr (IsResoMicrotrack) { + if constexpr (requires { track.itsNCls(); }) { + if (track.itsNCls() < mTrackCuts.cfgITSNClsMin) { + return kTrkFlags; + } + } + } + return kTrkClusters; + } + + // TOF signal requirement of the track (global, per species, or per pT bin) + template + bool passesTOFRequired(int species, const TrackType& track) + { + bool required = mTrackCuts.cfgHasTOF || mPID[species].onlyTOFTracks; + PIDCut cut; + // A pT outside all bins is rejected by passesPID + if (!mByPassTOF && getPIDCut(species, track.pt(), cut) && cut.tofRequired) { + required = true; + } + return !required || track.hasTOF(); + } + + // PID selection from the nSigma values of the species; tofNSigma is used only with hasTOF. + template + bool passesPID(int species, float pt, bool hasTOF, double tpcNSigma, double tofNSigma) + { + PIDCut cut; + if (!getPIDCut(species, pt, cut)) { + return false; + } + if (isCutEnabled(cut.tpcMax) && !passesMax(std::abs(tpcNSigma), cut.tpcMax)) { + return false; + } + // Missing TOF is handled by passesTOFRequired; here the TPC alone decides + if (mByPassTOF || !hasTOF) { + return true; + } + bool tofPassed = !isCutEnabled(cut.tofMax) || passesMax(std::abs(tofNSigma), cut.tofMax); + if (!tofPassed && cut.combined > 0 && tpcNSigma * tpcNSigma + tofNSigma * tofNSigma < cut.combined * cut.combined) { + tofPassed = true; + } + return tofPassed; + } + + private: + void checkConfiguration(bool microTracks) + { + // Consistency of the pT dependent PID configuration + bool anyOnlyTOF = false; + for (const auto& config : mPID) { + if (config.usePtDependent) { + const auto& bins = config.ptBins; + if (bins.size() < MinPtBinEdges || config.tpcNSigmaCuts.size() != bins.size() - 1 || + config.tofNSigmaCuts.size() != bins.size() - 1 || config.tofRequired.size() != bins.size() - 1) { + LOG(fatal) << config.species << " pT dependent PID vectors must have (number of pT bin edges - 1) entries"; + } + } + anyOnlyTOF = anyOnlyTOF || config.onlyTOFTracks; + } + if (mByPassTOF && anyOnlyTOF) { + LOG(warning) << "cByPassTOF skips the TOF nSigma selection, but cUseOnlyTOFTrack* still requires a TOF signal"; + } + + // Micro tracks store quantised DCA and nSigma: a cut off the grid would silently act as a different cut. + if (!microTracks) { + return; + } + auto checkDCAGrid = [](const char* name, double cut) { + const double nearest = std::min(std::max(std::round(cut / DCAGridStep) * DCAGridStep, 0.), DCAGridMax); + if (std::abs(cut - nearest) > GridTolerance) { + LOG(fatal) << name << " = " << cut << " is not on the quantised DCA grid (multiples of " << DCAGridStep << " up to " << DCAGridMax << "); nearest value: " << nearest; + } + }; + auto checkPIDGrid = [](const std::string& name, double cut) { + const double nearest = std::min(std::max(PIDGridStart + std::round((cut - PIDGridStart) / PIDGridStep) * PIDGridStep, PIDGridStart), PIDGridMax); + if (std::abs(cut - nearest) > GridTolerance) { + LOG(fatal) << name << " = " << cut << " is not on the quantised nSigma grid ({2.0, 2.25, ..., 3.5}); nearest value: " << nearest; + } + }; + if (mTrackCuts.cfgUsePtDepDCA) { + LOG(info) << "Micro tracks use the producer pT dependent DCA flags (0.004 + 0.013 / pT); cDCAToPVByPt* are ignored"; + if (std::abs(mTrackCuts.cDCAToPVByPtP0 - ProducerDCAPtP0) > ConfigTolerance || std::abs(mTrackCuts.cDCAToPVByPtCoeff - ProducerDCAPtCoeff) > ConfigTolerance || std::abs(mTrackCuts.cDCAToPVByPtPower - ProducerDCAPtPower) > ConfigTolerance) { + LOG(warning) << "cDCAToPVByPt* differ from the producer defaults, but micro tracks always use the producer formula"; + } + } else { + if (isCutEnabled(mTrackCuts.cMaxDCArToPVcut)) { + checkDCAGrid("cMaxDCArToPVcut", mTrackCuts.cMaxDCArToPVcut); + } + if (isCutEnabled(mTrackCuts.cMaxDCAzToPVcut)) { + checkDCAGrid("cMaxDCAzToPVcut", mTrackCuts.cMaxDCAzToPVcut); + } + } + if (isCutEnabled(mTrackCuts.cMinDCAzToPVcut)) { + checkDCAGrid("cMinDCAzToPVcut", mTrackCuts.cMinDCAzToPVcut); + } + bool anyCombined = false; + for (const auto& config : mPID) { + if (isCutEnabled(config.maxTPCnSigma) && !config.usePtDependent) { + checkPIDGrid(config.maxTPCName, config.maxTPCnSigma); + } + if (isCutEnabled(config.maxTOFnSigma) && !config.usePtDependent) { + checkPIDGrid(config.maxTOFName, config.maxTOFnSigma); + } + anyCombined = anyCombined || config.combinedNSigma > 0; + } + for (const auto& config : mPID) { + if (!config.usePtDependent) { + continue; + } + for (const auto& cut : config.tpcNSigmaCuts) { + if (isCutEnabled(cut)) { + checkPIDGrid(config.tpcCutsName, cut); + } + } + for (const auto& cut : config.tofNSigmaCuts) { + if (isCutEnabled(cut)) { + checkPIDGrid(config.tofCutsName, cut); + } + } + } + if (anyCombined) { + LOG(warning) << "nsigmaCutCombined* on micro tracks uses quantised nSigma values (approximate)"; + } + } + + EventCuts mEventCuts; + TrackCuts mTrackCuts; + std::vector mPID; + bool mByPassTOF = false; +}; + +} // namespace o2::analysis::resonance + +#endif // PWGLF_CORE_RESOANALYSISSELECTIONCORE_H_ diff --git a/PWGLF/DataModel/LFK1MlTables.h b/PWGLF/DataModel/LFK1MlTables.h new file mode 100644 index 00000000000..892d2adefe0 --- /dev/null +++ b/PWGLF/DataModel/LFK1MlTables.h @@ -0,0 +1,125 @@ +// 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 LFK1MlTables.h +/// \brief Derived K1 microtrack training and audit tables +/// \author Bong-Hwi Lim +/// +#ifndef PWGLF_DATAMODEL_LFK1MLTABLES_H_ +#define PWGLF_DATAMODEL_LFK1MLTABLES_H_ + +#include +#include + +#include + +namespace o2::aod +{ +namespace k1ml +{ +// Persisted relations target only these derived tables, so AO2D merge tools +// can relocate them. Source AO2D row numbers remain scalar audit values. +DECLARE_SOA_COLUMN(K1RecoCollisionId, k1RecoCollisionId, int64_t); //! row index of the reduced ResoCollisions_001 collision in the input DF +DECLARE_SOA_COLUMN(K1PosZ, k1PosZ, float); +DECLARE_SOA_COLUMN(K1BField, k1BField, float); +DECLARE_SOA_COLUMN(K1Centrality, k1Centrality, float); +DECLARE_SOA_COLUMN(K1Multiplicity, k1Multiplicity, float); +DECLARE_SOA_COLUMN(K1RecINELgt0, k1RecINELgt0, bool); +DECLARE_SOA_COLUMN(K1SourceTrackId, k1SourceTrackId, int64_t); //! trackId of the reduced micro track +DECLARE_SOA_COLUMN(K1Px, k1Px, float); +DECLARE_SOA_COLUMN(K1Py, k1Py, float); +DECLARE_SOA_COLUMN(K1Pz, k1Pz, float); +DECLARE_SOA_COLUMN(K1PidPi, k1PidPi, uint8_t); +DECLARE_SOA_COLUMN(K1PidKa, k1PidKa, uint8_t); +DECLARE_SOA_COLUMN(K1PidPr, k1PidPr, uint8_t); +DECLARE_SOA_COLUMN(K1SelectionFlags, k1SelectionFlags, uint8_t); +DECLARE_SOA_COLUMN(K1TrackFlags, k1TrackFlags, uint8_t); +DECLARE_SOA_COLUMN(K1CrossedRows, k1CrossedRows, uint8_t); +DECLARE_SOA_COLUMN(K1ItsClusterMap, k1ItsClusterMap, uint8_t); +DECLARE_SOA_COLUMN(K1Mass, k1Mass, float); +DECLARE_SOA_COLUMN(K1MassPiPi, k1MassPiPi, float); +DECLARE_SOA_COLUMN(K1MassKaPiSame, k1MassKaPiSame, float); +DECLARE_SOA_COLUMN(K1MassKaPiOpp, k1MassKaPiOpp, float); +DECLARE_SOA_COLUMN(K1ScalarSumPt, k1ScalarSumPt, float); +DECLARE_SOA_COLUMN(K1PiPiPt, k1PiPiPt, float); +DECLARE_SOA_COLUMN(K1Pt, k1Pt, float); +DECLARE_SOA_COLUMN(K1Y, k1Y, float); +DECLARE_SOA_COLUMN(K1Eta, k1Eta, float); +DECLARE_SOA_COLUMN(K1Phi, k1Phi, float); +DECLARE_SOA_COLUMN(K1Charge, k1Charge, int8_t); +// Cumulative selection bits of the unlike-sign candidate: +// 1 = valid canonical candidate inside the K1 rapidity window (loose stage), +// 2 = track quality of all three tracks, +// 4 = TOF requirement and PID of all three tracks, +// 8 = pion-pair pT and secondary mass window, +// 16 = candidate cuts. +// Candidates written at the "selected" export stage carry 31. +DECLARE_SOA_COLUMN(K1BaselinePassBits, k1BaselinePassBits, uint16_t); +DECLARE_SOA_COLUMN(K1MasterFeatures, k1MasterFeatures, float[125]); //! 125 master features of the K1 ML feature contract (PWGLF/Core/K1MlFeatures.h) +DECLARE_SOA_COLUMN(K1FeatureStatus, k1FeatureStatus, uint8_t); //! o2::analysis::k1ml::BuildStatus; only Ok (0) rows are written +DECLARE_SOA_COLUMN(K1TruthStatus, k1TruthStatus, uint8_t); //! 0 data, 1 matched, 2 unmatched +DECLARE_SOA_COLUMN(K1TruthChannel, k1TruthChannel, uint8_t); //! 0 none, 1 rho K, 2 K* pi +DECLARE_SOA_COLUMN(K1MotherPdg, k1MotherPdg, int32_t); +DECLARE_SOA_COLUMN(K1MotherId, k1MotherId, int64_t); +DECLARE_SOA_COLUMN(K1GeneratedPt, k1GeneratedPt, float); //! NaN in K1MlTruth: not filled for reconstructed candidates +DECLARE_SOA_COLUMN(K1GeneratedY, k1GeneratedY, float); //! NaN in K1MlTruth: not filled for reconstructed candidates +DECLARE_SOA_COLUMN(K1OriginalMcParticleId, k1OriginalMcParticleId, int64_t); +DECLARE_SOA_COLUMN(K1DaughterPdg1, k1DaughterPdg1, int32_t); +DECLARE_SOA_COLUMN(K1DaughterPdg2, k1DaughterPdg2, int32_t); +DECLARE_SOA_COLUMN(K1GenSelectedRecoEvent, k1GenSelectedRecoEvent, bool); //! true: parents are taken from selected reconstructed events +} // namespace k1ml + +DECLARE_SOA_TABLE(K1MlEvents, "AOD", "K1MLEVENT", + o2::soa::Index<>, k1ml::K1RecoCollisionId, k1ml::K1PosZ, k1ml::K1BField, + k1ml::K1Centrality, k1ml::K1Multiplicity, k1ml::K1RecINELgt0); +namespace k1ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(K1MlEvent, k1MlEvent, int, K1MlEvents, ""); +} // namespace k1ml +DECLARE_SOA_TABLE(K1MlTracks, "AOD", "K1MLTRACK", + o2::soa::Index<>, k1ml::K1MlEventId, k1ml::K1SourceTrackId, + k1ml::K1Px, k1ml::K1Py, k1ml::K1Pz, + k1ml::K1PidPi, k1ml::K1PidKa, k1ml::K1PidPr, + k1ml::K1SelectionFlags, k1ml::K1TrackFlags, + k1ml::K1CrossedRows, k1ml::K1ItsClusterMap); +namespace k1ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(K1MlKaonTrack, k1MlKaonTrack, int, K1MlTracks, "_Kaon"); +DECLARE_SOA_INDEX_COLUMN_FULL(K1MlSamePionTrack, k1MlSamePionTrack, int, K1MlTracks, "_Same"); +DECLARE_SOA_INDEX_COLUMN_FULL(K1MlOppPionTrack, k1MlOppPionTrack, int, K1MlTracks, "_Opp"); +} // namespace k1ml +DECLARE_SOA_TABLE(K1MlCandidates, "AOD", "K1MLCANDIDATE", + o2::soa::Index<>, k1ml::K1MlEventId, + k1ml::K1MlKaonTrackId, k1ml::K1MlSamePionTrackId, k1ml::K1MlOppPionTrackId, + k1ml::K1Mass, k1ml::K1MassPiPi, k1ml::K1MassKaPiSame, k1ml::K1MassKaPiOpp, + k1ml::K1ScalarSumPt, k1ml::K1PiPiPt, + k1ml::K1Pt, k1ml::K1Y, k1ml::K1Eta, k1ml::K1Phi, k1ml::K1Charge, + k1ml::K1BaselinePassBits); +namespace k1ml +{ +DECLARE_SOA_INDEX_COLUMN_FULL(K1MlCandidate, k1MlCandidate, int, K1MlCandidates, ""); +} // namespace k1ml +DECLARE_SOA_TABLE(K1MlInputs, "AOD", "K1MLINPUT", + o2::soa::Index<>, k1ml::K1MlCandidateId, + k1ml::K1MasterFeatures, k1ml::K1FeatureStatus); +DECLARE_SOA_TABLE(K1MlTruth, "AOD", "K1MLTRUTH", + o2::soa::Index<>, k1ml::K1MlCandidateId, + k1ml::K1TruthStatus, k1ml::K1TruthChannel, + k1ml::K1MotherPdg, k1ml::K1MotherId, + k1ml::K1GeneratedPt, k1ml::K1GeneratedY); +DECLARE_SOA_TABLE(K1MlGenAudit, "AOD", "K1MLGENAUDIT", + o2::soa::Index<>, k1ml::K1RecoCollisionId, k1ml::K1OriginalMcParticleId, + k1ml::K1MotherPdg, k1ml::K1DaughterPdg1, k1ml::K1DaughterPdg2, + k1ml::K1TruthChannel, k1ml::K1GeneratedPt, k1ml::K1GeneratedY, + k1ml::K1GenSelectedRecoEvent); +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFK1MLTABLES_H_ diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index f8a6669861e..b1f09db8d9a 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -79,6 +79,11 @@ o2physics_add_dpl_workflow(k1analysismicro PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(k1-training-table + SOURCES k1TrainingTable.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(phianalysisrun3 SOURCES phianalysisrun3.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Resonances/k1AnalysisMicro.cxx b/PWGLF/Tasks/Resonances/k1AnalysisMicro.cxx index 0cbab0d427c..9d6877ca7e5 100644 --- a/PWGLF/Tasks/Resonances/k1AnalysisMicro.cxx +++ b/PWGLF/Tasks/Resonances/k1AnalysisMicro.cxx @@ -14,12 +14,11 @@ /// \author Su-Jeong Ji , Bong-Hwi Lim /// +#include "PWGLF/Core/K1AnalysisMicroCore.h" +#include "PWGLF/Core/ResoAnalysisSelectionCore.h" #include "PWGLF/DataModel/LFResonanceTables.h" -#include -#include #include -#include #include #include #include @@ -28,129 +27,134 @@ #include #include #include +#include #include #include #include -#include -#include // FIXME +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#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::constants::math; -; +using namespace o2::analysis::resonance; +using namespace o2::analysis::k1micro; + +/// Histogram binning, QA and debug output +struct HistogramOptions : ConfigurableGroup { + Configurable cNbinsDiv{"cNbinsDiv", 1, "Integer to divide the number of bins"}; + Configurable additionalQAplots{"additionalQAplots", true, "Additional QA plots"}; + Configurable cfgTruthDebug{"cfgTruthDebug", 0, "Maximum logged matched candidates per truth channel"}; +}; + +enum BinAnti : unsigned int { + kNormal = 0, + kAnti, + kNAEnd +}; + +enum BinType : unsigned int { + kK1P = 0, + kK1N, + kK1P_Mix, + kK1N_Mix, + kK1P_GenINEL10, + kK1N_GenINEL10, + kK1P_GenINELgt10, + kK1N_GenINELgt10, + kK1P_GenTrig10, + kK1N_GenTrig10, + kK1P_GenEvtSel, + kK1N_GenEvtSel, + kK1P_Rec, + kK1N_Rec, + kTYEnd +}; + +enum class QAFolder { + Before, // QA/*: before the candidate cuts + After, // QAcut/*: after the candidate cuts + MC // QAMC/*: matched K1 truth candidates +}; struct K1AnalysisMicro { - enum BinAnti : unsigned int { - kNormal = 0, - kAnti, - kNAEnd - }; - enum BinType : unsigned int { - kK1P = 0, - kK1N, - kK1P_Mix, - kK1N_Mix, - kK1P_GenINEL10, - kK1N_GenINEL10, - kK1P_GenINELgt10, - kK1N_GenINELgt10, - kK1P_GenTrig10, - kK1N_GenTrig10, - kK1P_GenEvtSel, - kK1N_GenEvtSel, - kK1P_Rec, - kK1N_Rec, - kTYEnd - }; + // Module-initializer v001 tables; full tracks keep their unversioned schema as a fallback. + using ResoCollisions = aod::ResoCollisions_001; + using ResoMCCols = soa::Join; + using ResoTracks = aod::ResoTracks; // no v001 exists; K1 does not need ResoTrackTracks (trackId unused) + using ResoMicroTracks = aod::ResoMicroTracks_001; + using ResoMCTracks = soa::Join; + using ResoMCMicroTracks = soa::Join; + using ResoMCParents = aod::ResoMCParents_001; + SliceCache cache; - Preslice perRCol = aod::resodaughter::resoCollisionId; - Preslice perCollision = aod::track::collisionId; + // Registered only to enable the slice cache that SameKindPair (event mixing) needs, as in Xi1820Analysis + Preslice perResoCollisionTrack = aod::resodaughter::resoCollisionId; + Preslice perResoCollisionMicroTrack = aod::resodaughter::resoCollisionId; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - using ResoMCCols = soa::Join; - //// Configurables - Configurable cNbinsDiv{"cNbinsDiv", 1, "Integer to divide the number of bins"}; + // Selection shared with the K1 training-table task (plain JSON keys, no group prefix) + EventCuts eventCuts; + TrackCuts trackCuts; + PionPidCuts pionPID; + KaonPidCuts kaonPID; + SecondaryCuts secondaryCuts; + CandidateCuts candidateCuts; + HistogramOptions histogramOptions; + /// Event Mixing Configurable nEvtMixing{"nEvtMixing", 5, "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, 20.0f, 40.0f, 60.0f, 80.0f, 100.0f, 200.0f, 99999.f}, "Mixing bins - multiplicity"}; - /// Pre-selection cuts - Configurable cMinPtcut{"cMinPtcut", 0.15, "Track minium pt cut"}; - - /// DCA Selections - // DCAr to PV - Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.1, "Track DCAr cut to PV Maximum"}; - // DCAz to PV - Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 0.1, "Track DCAz cut to PV Maximum"}; - Configurable cMinDCAzToPVcut{"cMinDCAzToPVcut", 0.0, "Track DCAz cut to PV Minimum"}; - - /// PID Selections - Configurable cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 3.0, "TPC nSigma cut for Pion"}; // TPC - Configurable cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion"}; // TOF - Configurable nsigmaCutCombinedPion{"nsigmaCutCombinedPion", -999, "Combined nSigma cut for Pion"}; // Combined - Configurable cTOFVeto{"cTOFVeto", true, "TOF Veto, if false, TOF is nessessary for PID selection"}; // TOF Veto - Configurable cUseOnlyTOFTrackPi{"cUseOnlyTOFTrackPi", false, "Use only TOF track for PID selection"}; // Use only TOF track for Pion PID selection - // Kaon - Configurable cMaxTPCnSigmaKaon{"cMaxTPCnSigmaKaon", 3.0, "TPC nSigma cut for Kaon"}; // TPC - Configurable cMaxTOFnSigmaKaon{"cMaxTOFnSigmaKaon", 3.0, "TOF nSigma cut for Kaon"}; // TOF - Configurable nsigmaCutCombinedKaon{"nsigmaCutCombinedKaon", -999, "Combined nSigma cut for Kaon"}; // Combined - Configurable cUseOnlyTOFTrackKa{"cUseOnlyTOFTrackKa", false, "Use only TOF track for PID selection"}; // Use only TOF track for Kaon PID selection - // Track selections - Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz - Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) - Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz - Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor - Configurable additionalQAplots{"additionalQAplots", true, "Additional QA plots"}; - Configurable additionalEvsel{"additionalEvsel", true, "Additional event selcection"}; - Configurable cfgTPCcluster{"cfgTPCcluster", 0, "Number of TPC cluster"}; - Configurable cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; - Configurable cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"}; - Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; - - // Secondary selection - Configurable cMinSecondaryPtCut{"cMinSecondaryPtCut", 0.5, "Min pT cut for secondary selection"}; - /* - Configurable cfgModeK892orRho{"cfgModeK892orRho", false, "Secondary scenario for K892 (true) or Rho (false)"}; - Configurable cSecondaryMasswindow{"cSecondaryMasswindow", 0.1, "Secondary inv mass selection window"}; - Configurable cMinAnotherSecondaryMassCut{"cMinAnotherSecondaryMassCut", 0, "Min inv. mass selection of another secondary scenario"}; - Configurable cMaxAnotherSecondaryMassCut{"cMaxAnotherSecondaryMassCut", 999, "MAx inv. mass selection of another secondary scenario"}; - Configurable cMinPiKaMassCut{"cMinPiKaMassCut", 0, "bPion-Kaon pair inv mass selection minimum"}; - Configurable cMaxPiKaMassCut{"cMaxPiKaMassCut", 999, "bPion-Kaon pair inv mass selection maximum"}; - Configurable cMinAngle{"cMinAngle", 0, "Minimum angle between K(892)0 and bachelor pion"}; - Configurable cMaxAngle{"cMaxAngle", 4, "Maximum angle between K(892)0 and bachelor pion"}; - Configurable cMinPairAsym{"cMinPairAsym", -1, "Minimum pair asymmetry"}; - Configurable cMaxPairAsym{"cMaxPairAsym", 1, "Maximum pair asymmetry"}; -*/ - - // K1 selection - Configurable cK1MaxRap{"cK1MaxRap", 0.5, "K1 maximum rapidity"}; - Configurable cK1MinRap{"cK1MinRap", -0.5, "K1 minimum rapidity"}; - - void init(o2::framework::InitContext&) + + K1AnalysisMicroCore core; + std::array truthDebugCounts{}; + + void init(InitContext&) { + const int sameEventModes = static_cast(doprocessResoTracks) + static_cast(doprocessResoMicroTracks) + + static_cast(doprocessMC) + static_cast(doprocessMCMicro); + const int mixedEventModes = static_cast(doprocessME) + static_cast(doprocessMEMicro); + if (sameEventModes > 1 || mixedEventModes > 1) { + LOG(fatal) << "Enable at most one same-event mode and one mixing mode"; + } + + ProcessModes modes; + modes.microTracks = doprocessResoMicroTracks || doprocessMCMicro || doprocessMEMicro; + modes.mcReco = doprocessMC || doprocessMCMicro; + modes.mcRecoMicro = doprocessMCMicro; + modes.mcGen = doprocessMCTrue; + core.init(histos, eventCuts, trackCuts, pionPID, kaonPID, secondaryCuts, candidateCuts, modes); + registerHistograms(modes); + + // Print output histograms statistics + LOG(info) << "Size of the histograms in K1 Analysis Task"; + histos.print(); + } + + void registerHistograms(ProcessModes const& modes) + { + const int nBinsDiv = histogramOptions.cNbinsDiv; std::vector centBinning = {0., 1., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 80., 90., 100., 200.}; AxisSpec centAxis = {centBinning, "T0M (%)"}; AxisSpec ptAxis = {150, 0, 15, "#it{p}_{T} (GeV/#it{c})"}; AxisSpec dcaxyAxis = {300, 0, 3, "DCA_{#it{xy}} (cm)"}; AxisSpec dcazAxis = {500, 0, 5, "DCA_{#it{z}} (cm)"}; - AxisSpec invMassAxisK892 = {1400 / cNbinsDiv, 0.6, 2.0, "Invariant Mass (GeV/#it{c}^2)"}; // K(892)0 - AxisSpec invMassAxisRho = {2000 / cNbinsDiv, 0.0, 2.0, "Invariant Mass (GeV/#it{c}^2)"}; // rho - AxisSpec invMassAxisReso = {1600 / cNbinsDiv, 0.9f, 2.5f, "Invariant Mass (GeV/#it{c}^2)"}; // K1 - AxisSpec invMassAxisScan = {250, 0, 2.5, "Invariant Mass (GeV/#it{c}^2)"}; // For selection + AxisSpec invMassAxisK892 = {1400 / nBinsDiv, 0.6, 2.0, "Invariant Mass (GeV/#it{c}^2)"}; // K(892)0 + AxisSpec invMassAxisRho = {2000 / nBinsDiv, 0.0, 2.0, "Invariant Mass (GeV/#it{c}^2)"}; // rho + AxisSpec invMassAxisReso = {1600 / nBinsDiv, 0.9f, 2.5f, "Invariant Mass (GeV/#it{c}^2)"}; // K1 AxisSpec pidQAAxis = {130, -6.5, 6.5}; - AxisSpec dataTypeAxis = {9, 0, 9, "Histogram types"}; - AxisSpec mcTypeAxis = {4, 0, 4, "Histogram types"}; // THnSparse AxisSpec axisAnti = {BinAnti::kNAEnd, 0, BinAnti::kNAEnd, "Type of bin: Normal or Anti"}; AxisSpec axisType = {BinType::kTYEnd, 0, BinType::kTYEnd, "Type of bin with charge and mix"}; - AxisSpec mcLabelAxis = {5, -0.5, 4.5, "MC Label"}; // DCA QA // Primary pion @@ -223,7 +227,16 @@ struct K1AnalysisMicro { histos.add("k1invmass_Mix", "Invariant mass of K1(1270) (ME)", HistType::kTH1F, {invMassAxisReso}); // MC - if (doprocessMC) { + if (modes.mcReco) { + AxisSpec channelAxis = {3, -0.5, 2.5, "0: non-K1, 1: rho K, 2: K* pi"}; + histos.add("MCReco/collisions", "Selected reconstructed MC collisions", HistType::kTH1D, {{1, 0, 1}}); + histos.add("MCReco/microTracks", "Input micro tracks in selected MC collisions", HistType::kTH1D, {{1, 0, 1}}); + histos.add("MCReco/channel", "All selected pi-pi-K combinations by truth channel", HistType::kTH1D, {channelAxis}); + histos.add("MCReco/mass", "Reconstructed mass by truth channel", HistType::kTH2D, {channelAxis, invMassAxisReso}); + histos.add("MCReco/pt", "Reconstructed pT by truth channel", HistType::kTH2D, {channelAxis, ptAxis}); + histos.add("MCReco/piPiMass", "pi-pi mass by truth channel", HistType::kTH2D, {channelAxis, invMassAxisRho}); + histos.add("MCReco/pi1KMass", "First pion-kaon mass by truth channel", HistType::kTH2D, {channelAxis, invMassAxisK892}); + histos.add("MCReco/pi2KMass", "Second pion-kaon mass by truth channel", HistType::kTH2D, {channelAxis, invMassAxisK892}); histos.add("k1invmass_MC", "Invariant mass of K1(1270)", HistType::kTH1F, {invMassAxisReso}); histos.add("k1invmass_MC_noK1", "Invariant mass of K1(1270)", HistType::kTH1F, {invMassAxisReso}); @@ -254,494 +267,324 @@ struct K1AnalysisMicro { histos.add("QAMC/hInvmassSecon_PiKa", "Invariant mass of secondary resonance vs pion-kaon", HistType::kTH2F, {invMassAxisRho, invMassAxisK892}); histos.add("QAMC/hInvmassSecon", "Invariant mass of secondary resonance", HistType::kTH1F, {invMassAxisRho}); histos.add("QAMC/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1F, {ptAxis}); - } // doprocessMC - // Print output histograms statistics - LOG(info) << "Size of the histograms in K1 Analysis Task"; - histos.print(); - } // init - - // PDG code - int kPDGRho770 = 113; - int kK1Plus = 10323; - - template - bool trackCut(const TrackType& track) - { - if constexpr (!IsResoMicrotrack) { - // basic track cuts - if (std::abs(track.pt()) < cMinPtcut) - return false; - if (std::abs(track.dcaXY()) > cMaxDCArToPVcut) - return false; - if (std::abs(track.dcaZ()) > cMaxDCAzToPVcut) - return false; - if (track.tpcNClsFound() < cfgTPCcluster) - return false; - if (cfgHasTOF && !track.hasTOF()) - return false; - if (cfgUseITSRefit && !track.passedITSRefit()) - return false; - if (cfgUseTPCRefit && !track.passedTPCRefit()) - return false; - if (cfgPVContributor && !track.isPVContributor()) - return false; - if (cfgPrimaryTrack && !track.isPrimaryTrack()) - return false; - if (cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) - return false; - if (cfgGlobalTrack && !track.isGlobalTrack()) - return false; - } else { - if (std::abs(track.pt()) < cMinPtcut) - return false; - if (o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(track.trackSelectionFlags()) > cMaxDCArToPVcut - Epsilon) - return false; - if (o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(track.trackSelectionFlags()) > cMaxDCAzToPVcut - Epsilon) - return false; - if (cfgPrimaryTrack && !track.isPrimaryTrack()) - return false; - if (cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) - return false; - if (cfgPVContributor && !track.isPVContributor()) - return false; + } // mcReco + if (modes.mcGen) { + AxisSpec channelAxis = {3, -0.5, 2.5, "0: other/unresolved, 1: rho K, 2: K* pi"}; + histos.add("MCGen/chargeChannel", "K1 parents in selected reconstructed events, inside the K1 rapidity window", HistType::kTH2D, {{2, -1.5, 1.5, "K1 charge"}, channelAxis}); + histos.add("MCGen/ptChannel", "Generated K1 pT by immediate decay channel", HistType::kTH2D, {channelAxis, ptAxis}); } - return true; } - // Pion PID selection tools - template - bool selectionPIDpion(const T& candidate) + // Track QA of a pion; isPrimary selects the trkppion (first) or trkspion (second) histograms + template + void fillPionQA(const TrackType& track, bool isPrimary) { - if constexpr (!IsResoMicrotrack) { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - if (std::abs(candidate.tpcNSigmaPi()) < cMaxTPCnSigmaPion) { - tpcPIDPassed = true; - } else { - return false; - } - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaPi()) < cMaxTOFnSigmaPion) { - tofPIDPassed = true; + const bool hasTOF = track.hasTOF(); + if (isPrimary) { + if constexpr (Folder == QAFolder::Before) { + histos.fill(HIST("QA/trkppionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QA/trkppionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QA/trkppionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); } - if ((nsigmaCutCombinedPion > 0) && (candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi() < nsigmaCutCombinedPion * nsigmaCutCombinedPion)) { - tofPIDPassed = true; + histos.fill(HIST("QA/trkppionpT"), track.pt()); + histos.fill(HIST("QA/trkppionDCAxy"), track.dcaXY()); + histos.fill(HIST("QA/trkppionDCAz"), track.dcaZ()); + } else if constexpr (Folder == QAFolder::After) { + histos.fill(HIST("QAcut/trkppionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QAcut/trkppionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QAcut/trkppionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); } + histos.fill(HIST("QAcut/trkppionpT"), track.pt()); + histos.fill(HIST("QAcut/trkppionDCAxy"), track.dcaXY()); + histos.fill(HIST("QAcut/trkppionDCAz"), track.dcaZ()); } else { - if (!cTOFVeto) { - return false; + histos.fill(HIST("QAMC/trkppionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QAMC/trkppionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QAMC/trkppionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); } - tofPIDPassed = true; - } - if (tpcPIDPassed && tofPIDPassed) { - return true; + histos.fill(HIST("QAMC/trkppionpT"), track.pt()); + histos.fill(HIST("QAMC/trkppionDCAxy"), track.dcaXY()); + histos.fill(HIST("QAMC/trkppionDCAz"), track.dcaZ()); } } else { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - tpcPIDPassed = std::abs(o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(candidate.pidNSigmaPiFlag())) < cMaxTPCnSigmaPion + Epsilon; - tofPIDPassed = candidate.hasTOF() ? std::abs(o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(candidate.pidNSigmaPiFlag())) < cMaxTOFnSigmaPion + Epsilon : true; - if (tpcPIDPassed && tofPIDPassed) { - return true; + if constexpr (Folder == QAFolder::Before) { + histos.fill(HIST("QA/trkspionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QA/trkspionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QA/trkspionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + } + histos.fill(HIST("QA/trkspionpT"), track.pt()); + histos.fill(HIST("QA/trkspionDCAxy"), track.dcaXY()); + histos.fill(HIST("QA/trkspionDCAz"), track.dcaZ()); + } else if constexpr (Folder == QAFolder::After) { + histos.fill(HIST("QAcut/trkspionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QAcut/trkspionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QAcut/trkspionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + } + histos.fill(HIST("QAcut/trkspionpT"), track.pt()); + histos.fill(HIST("QAcut/trkspionDCAxy"), track.dcaXY()); + histos.fill(HIST("QAcut/trkspionDCAz"), track.dcaZ()); + } else { + histos.fill(HIST("QAMC/trkspionTPCPID"), track.pt(), track.tpcNSigmaPi()); + if (hasTOF) { + histos.fill(HIST("QAMC/trkspionTOFPID"), track.pt(), track.tofNSigmaPi()); + histos.fill(HIST("QAMC/trkspionTPCTOFPID"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + } + histos.fill(HIST("QAMC/trkspionpT"), track.pt()); + histos.fill(HIST("QAMC/trkspionDCAxy"), track.dcaXY()); + histos.fill(HIST("QAMC/trkspionDCAz"), track.dcaZ()); } } - return false; } - // Kaon PID selection tools - template - bool selectionPIDkaon(const T& candidate) + // Track QA of the bachelor kaon + template + void fillKaonQA(const TrackType& track) { - if constexpr (!IsResoMicrotrack) { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - if (std::abs(candidate.tpcNSigmaKa()) < cMaxTPCnSigmaKaon) { - tpcPIDPassed = true; - } else { - return false; + const bool hasTOF = track.hasTOF(); + if constexpr (Folder == QAFolder::Before) { + histos.fill(HIST("QA/trkkaonTPCPID"), track.pt(), track.tpcNSigmaKa()); + if (hasTOF) { + histos.fill(HIST("QA/trkkaonTOFPID"), track.pt(), track.tofNSigmaKa()); + histos.fill(HIST("QA/trkkaonTPCTOFPID"), track.tpcNSigmaKa(), track.tofNSigmaKa()); } - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaKa()) < cMaxTOFnSigmaKaon) { - tofPIDPassed = true; - } - if ((nsigmaCutCombinedKaon > 0) && (candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa() < nsigmaCutCombinedKaon * nsigmaCutCombinedKaon)) { - tofPIDPassed = true; - } - } else { - if (!cTOFVeto) { - return false; - } - tofPIDPassed = true; - } - if (tpcPIDPassed && tofPIDPassed) { - return true; + histos.fill(HIST("QA/trkkaonpT"), track.pt()); + histos.fill(HIST("QA/trkkaonDCAxy"), track.dcaXY()); + histos.fill(HIST("QA/trkkaonDCAz"), track.dcaZ()); + } else if constexpr (Folder == QAFolder::After) { + histos.fill(HIST("QAcut/trkkaonTPCPID"), track.pt(), track.tpcNSigmaKa()); + if (hasTOF) { + histos.fill(HIST("QAcut/trkkaonTOFPID"), track.pt(), track.tofNSigmaKa()); + histos.fill(HIST("QAcut/trkkaonTPCTOFPID"), track.tpcNSigmaKa(), track.tofNSigmaKa()); } + histos.fill(HIST("QAcut/trkkaonpT"), track.pt()); + histos.fill(HIST("QAcut/trkkaonDCAxy"), track.dcaXY()); + histos.fill(HIST("QAcut/trkkaonDCAz"), track.dcaZ()); } else { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - tpcPIDPassed = std::abs(o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(candidate.pidNSigmaKaFlag())) < cMaxTPCnSigmaKaon + Epsilon; - tofPIDPassed = candidate.hasTOF() ? std::abs(o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(candidate.pidNSigmaKaFlag())) < cMaxTOFnSigmaKaon + Epsilon : true; - if (tpcPIDPassed && tofPIDPassed) { - return true; + histos.fill(HIST("QAMC/trkkaonTPCPID"), track.pt(), track.tpcNSigmaKa()); + if (hasTOF) { + histos.fill(HIST("QAMC/trkkaonTOFPID"), track.pt(), track.tofNSigmaKa()); + histos.fill(HIST("QAMC/trkkaonTPCTOFPID"), track.tpcNSigmaKa(), track.tofNSigmaKa()); } + histos.fill(HIST("QAMC/trkkaonpT"), track.pt()); + histos.fill(HIST("QAMC/trkkaonDCAxy"), track.dcaXY()); + histos.fill(HIST("QAMC/trkkaonDCAz"), track.dcaZ()); } - return false; } - template - bool isTrueK1(const T& trk1, const T& trk2, const T2& bTrack) - { - if (std::abs(trk1.pdgCode()) != kPiPlus || std::abs(trk2.pdgCode()) != kPiPlus) - return false; - if (std::abs(bTrack.pdgCode()) != kKPlus) - return false; - auto mother1 = trk1.motherId(); - auto mother2 = trk2.motherId(); - if (mother1 != mother2) - return false; - if (((std::abs(trk1.motherPDG()) && std::abs(trk2.motherPDG()) != kPDGRho770) && (std::abs(bTrack.motherPDG()) != kK1Plus)) || (std::abs(trk1.motherPDG()) && std::abs(bTrack.motherPDG()) != kK0Star892 && (std::abs(trk2.motherPDG()) != kK1Plus)) || (std::abs(trk2.motherPDG()) && std::abs(bTrack.motherPDG()) != kK0Star892 && (std::abs(trk1.motherPDG()) != kK1Plus))) - return false; - auto siblings = bTrack.siblingIds(); - if (siblings[0] != mother1 && siblings[1] != mother2) - return false; - return true; - } // isTrueK1 - - template - bool isTrueK892(const T& trk1, const T& trk2) - { - if (std::abs(trk1.pdgCode()) != kPiPlus || std::abs(trk2.pdgCode()) != kKPlus) - return false; - auto mother1 = trk1.motherId(); - auto mother2 = trk2.motherId(); - if (mother1 != mother2) - return false; - if (std::abs(trk1.motherPDG()) != kK0Star892) - return false; - return true; - } - - template - bool isTrueRho(const T& trk1, const T& trk2) - { - if (std::abs(trk1.pdgCode()) != kPiPlus || std::abs(trk2.pdgCode()) != kPiPlus) - return false; - auto mother1 = trk1.motherId(); - auto mother2 = trk2.motherId(); - if (mother1 != mother2) - return false; - if (std::abs(trk1.motherPDG()) != kPDGRho770) - return false; - return true; - } + // Histograms of the selected pion pairs and (pion, pion, kaon) candidates of one collision + // (or one mixed pair of collisions). The selection itself is the shared K1 core. + // dTracks1: bachelor kaons, dTracks2: pions. template void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2) { - auto multiplicity = collision.cent(); - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonanceSecondary, lDecayDaughter_bach, lResonanceK1; - for (const auto& [trk1, trk2] : combinations(CombinationsFullIndexPolicy(dTracks2, dTracks2))) { - // Full index policy is needed to consider all possible combinations - if (trk1.index() == trk2.index()) - continue; // We need to run (0,1), (1,0) pairs too. But the same id pairs are not needed. - // trk1: pion, trk2: pion, bTrack: kaon - if (!trackCut(trk1) || !trackCut(trk2)) - continue; + const bool fillQA = !IsMix && histogramOptions.additionalQAplots; + const auto multiplicity = collision.cent(); + + // Pion pair passing the pion selection; trk1 is the pion with the lower index + auto onPair = [&](auto const& trk1, auto const& trk2, ROOT::Math::PxPyPzMVector const& secondary, bool passesPairPt) { + if (fillQA) { + fillPionQA(trk1, true); + fillPionQA(trk2, false); + } + if (!passesPairPt) { + return; + } + if (fillQA) { + histos.fill(HIST("QA/hInvmassSecon"), secondary.M()); + } + if constexpr (IsMC) { + histos.fill(HIST("QAMC/hpT_Secondary"), secondary.Pt()); + } + }; - auto trk1pt = trk1.pt(); - auto trk2pt = trk2.pt(); - auto isTrk1hasTOF = trk1.hasTOF(); - auto isTrk2hasTOF = trk2.hasTOF(); - - if constexpr (!IsResoMicrotrack) { - auto trk1NSigmaPiTPC = trk1.tpcNSigmaPi(); - auto trk1NSigmaPiTOF = (isTrk1hasTOF) ? trk1.tofNSigmaPi() : -999.; - auto trk2NSigmaPiTPC = trk2.tpcNSigmaPi(); - auto trk2NSigmaPiTOF = (isTrk2hasTOF) ? trk2.tofNSigmaPi() : -999.; - - if (cUseOnlyTOFTrackPi && !isTrk1hasTOF) - continue; - if (!selectionPIDpion(trk1) || !selectionPIDpion(trk2)) - continue; - - if constexpr (!IsMix) { - - histos.fill(HIST("QA/trkppionTPCPID"), trk1pt, trk1NSigmaPiTPC); - if (isTrk1hasTOF) { - histos.fill(HIST("QA/trkppionTOFPID"), trk1pt, trk1NSigmaPiTOF); - histos.fill(HIST("QA/trkppionTPCTOFPID"), trk1NSigmaPiTPC, trk1NSigmaPiTOF); - } - histos.fill(HIST("QA/trkppionpT"), trk1pt); - histos.fill(HIST("QA/trkppionDCAxy"), trk1.dcaXY()); - histos.fill(HIST("QA/trkppionDCAz"), trk1.dcaZ()); - - histos.fill(HIST("QA/trkspionTPCPID"), trk2pt, trk2NSigmaPiTPC); - if (isTrk2hasTOF) { - histos.fill(HIST("QA/trkspionTOFPID"), trk2pt, trk2NSigmaPiTOF); - histos.fill(HIST("QA/trkspionTPCTOFPID"), trk2NSigmaPiTPC, trk2NSigmaPiTOF); - } - histos.fill(HIST("QA/trkspionpT"), trk2pt); - histos.fill(HIST("QA/trkspionDCAxy"), trk2.dcaXY()); - histos.fill(HIST("QA/trkspionDCAz"), trk2.dcaZ()); - } - } else { - histos.fill(HIST("QA/trkppionTPCPID"), trk1pt, o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk1.pidNSigmaPiFlag())); - if (isTrk1hasTOF) { - histos.fill(HIST("QA/trkppionTOFPID"), trk1pt, o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk1.pidNSigmaPiFlag())); - histos.fill(HIST("QA/trkppionTPCTOFPID"), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk1.pidNSigmaPiFlag()), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk1.pidNSigmaPiFlag())); - } - histos.fill(HIST("QA/trkppionpT"), trk1pt); - histos.fill(HIST("QA/trkppionDCAxy"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(trk1.trackSelectionFlags())); - histos.fill(HIST("QA/trkppionDCAz"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(trk1.trackSelectionFlags())); - - histos.fill(HIST("QA/trkspionTPCPID"), trk2pt, o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk2.pidNSigmaPiFlag())); - if (isTrk2hasTOF) { - histos.fill(HIST("QA/trkspionTOFPID"), trk2pt, o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk2.pidNSigmaPiFlag())); - histos.fill(HIST("QA/trkspionTPCTOFPID"), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk2.pidNSigmaPiFlag()), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk2.pidNSigmaPiFlag())); - } - histos.fill(HIST("QA/trkspionpT"), trk2pt); - histos.fill(HIST("QA/trkspionDCAxy"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(trk2.trackSelectionFlags())); - histos.fill(HIST("QA/trkspionDCAz"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(trk2.trackSelectionFlags())); + // Candidate passing the pion, pair and kaon selection; pion1 is the K*0 partner in the unlike-sign case + auto onCandidate = [&](auto const& kaon, auto const& pion1, auto const& pion2, K1CandidateValues const& c) { + if (fillQA) { + fillKaonQA(kaon); + } + if (!c.inRapidity) { + return; } - // Resonance reconstruction - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), MassPionCharged); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), MassPionCharged); - lResonanceSecondary = lDecayDaughter1 + lDecayDaughter2; + // QA histogram before the candidate cuts + if (fillQA) { + histos.fill(HIST("QA/K1OA"), c.angle); + histos.fill(HIST("QA/K1PairAsym"), c.pairAsym); + histos.fill(HIST("QA/hInvmassK892_Rho"), c.mass13, c.secondary.M()); + histos.fill(HIST("QA/hInvmassSecon_PiKa"), c.secondary.M(), c.mass23); + histos.fill(HIST("QA/hpT_Secondary"), c.secondary.Pt()); + } + if (!c.passesCandidateCuts) { + return; + } - if (lResonanceSecondary.Pt() < cMinSecondaryPtCut) - continue; + // QA histograms after the candidate cuts + if (fillQA) { + fillPionQA(pion1, true); + fillPionQA(pion2, false); + fillKaonQA(kaon); + histos.fill(HIST("QAcut/K1OA"), c.angle); + histos.fill(HIST("QAcut/K1PairAsym"), c.pairAsym); + histos.fill(HIST("QAcut/hInvmassK892_Rho"), c.mass13, c.secondary.M()); + histos.fill(HIST("QAcut/hInvmassSecon_PiKa"), c.secondary.M(), c.mass23); + histos.fill(HIST("QAcut/hInvmassSecon"), c.secondary.M()); + histos.fill(HIST("QAcut/hpT_Secondary"), c.secondary.Pt()); + } - if constexpr (!IsMix) { - histos.fill(HIST("QA/hInvmassSecon"), lResonanceSecondary.M()); + const unsigned int typeNormal = BinAnti::kNormal; + if constexpr (IsMix) { + const unsigned int typeK1 = kaon.sign() > 0 ? BinType::kK1P_Mix : BinType::kK1N_Mix; + histos.fill(HIST("hInvmass_K1_Mix"), typeNormal, typeK1, multiplicity, c.k1.Pt(), c.k1.M()); + histos.fill(HIST("k1invmass_Mix"), c.k1.M()); + return; } - if constexpr (IsMC) { - /* - if (isTrueK892(trk1, trk2)) - histos.fill(HIST("QAMC/hpT_Secondary"), lResonanceSecondary.Pt()); - } else { - if (isTrueRho(trk1, trk2)) - histos.fill(HIST("QAMC/hpT_Secondary"), lResonanceSecondary.Pt()); - } - */ - histos.fill(HIST("QAMC/hpT_Secondary"), lResonanceSecondary.Pt()); + unsigned int typeK1 = kaon.sign() > 0 ? BinType::kK1P : BinType::kK1N; + if (c.isUnlikeSign) { + histos.fill(HIST("k1invmass"), c.k1.M()); + histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, c.k1.Pt(), c.k1.M()); + } else { + histos.fill(HIST("k1invmass_LS"), c.k1.M()); + histos.fill(HIST("hInvmass_K1_LS"), typeNormal, typeK1, multiplicity, c.k1.Pt(), c.k1.M()); } - // Mass Window cut is removed - - for (const auto& bTrack : dTracks1) { - if (bTrack.index() == trk1.index() || bTrack.index() == trk2.index()) - continue; - if (!trackCut(bTrack)) - continue; - if (!selectionPIDkaon(bTrack)) - continue; - - // K1 reconstruction - lDecayDaughter_bach.SetXYZM(bTrack.px(), bTrack.py(), bTrack.pz(), MassKaonCharged); - lResonanceK1 = lResonanceSecondary + lDecayDaughter_bach; - - // Cuts - if (lResonanceK1.Rapidity() > cK1MaxRap || lResonanceK1.Rapidity() < cK1MinRap) - continue; - - auto lK1Angle = lResonanceSecondary.Angle(lDecayDaughter_bach.Vect()); - auto lPairAsym = (lResonanceSecondary.E() - lDecayDaughter_bach.E()) / (lResonanceSecondary.E() + lDecayDaughter_bach.E()); - - TLorentzVector temp13 = lDecayDaughter1 + lDecayDaughter_bach; - TLorentzVector temp23 = lDecayDaughter2 + lDecayDaughter_bach; - - // QA histogram - if constexpr (!IsMix) { - histos.fill(HIST("QA/K1OA"), lK1Angle); - histos.fill(HIST("QA/K1PairAsym"), lPairAsym); - histos.fill(HIST("QA/hInvmassK892_Rho"), temp13.M(), lResonanceSecondary.M()); - histos.fill(HIST("QA/hInvmassSecon_PiKa"), lResonanceSecondary.M(), temp23.M()); - histos.fill(HIST("QA/hpT_Secondary"), lResonanceSecondary.Pt()); + + if constexpr (IsMC) { + const auto channel = classifyK1Truth(pion1, pion2, kaon); + const int channelBin = static_cast(channel); + histos.fill(HIST("MCReco/channel"), channelBin); + histos.fill(HIST("MCReco/mass"), channelBin, c.k1.M()); + histos.fill(HIST("MCReco/pt"), channelBin, c.k1.Pt()); + histos.fill(HIST("MCReco/piPiMass"), channelBin, c.secondary.M()); + histos.fill(HIST("MCReco/pi1KMass"), channelBin, c.mass13); + histos.fill(HIST("MCReco/pi2KMass"), channelBin, c.mass23); + if (channel == K1TruthChannel::None) { + histos.fill(HIST("k1invmass_MC_noK1"), c.k1.M()); + return; } - // Selection cuts are removed - // QA histograms after the cuts are removed as no cuts are applied - - if constexpr (!IsMix) { - unsigned int typeK1 = bTrack.sign() > 0 ? BinType::kK1P : BinType::kK1N; - unsigned int typeNormal = BinAnti::kNormal; - if (trk1.sign() * trk2.sign() < 0) { - histos.fill(HIST("k1invmass"), lResonanceK1.M()); - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, lResonanceK1.Pt(), lResonanceK1.M()); - } else { - histos.fill(HIST("k1invmass_LS"), lResonanceK1.M()); - histos.fill(HIST("hInvmass_K1_LS"), typeNormal, typeK1, multiplicity, lResonanceK1.Pt(), lResonanceK1.M()); - } - - if constexpr (IsMC) { - if (isTrueK1(trk1, trk2, bTrack)) { - typeK1 = bTrack.sign() > 0 ? BinType::kK1P_Rec : BinType::kK1N_Rec; - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, lResonanceK1.Pt(), lResonanceK1.M()); - histos.fill(HIST("k1invmass_MC"), lResonanceK1.M()); - histos.fill(HIST("QAMC/K1OA"), lK1Angle); - histos.fill(HIST("QAMC/K1PairAsym"), lPairAsym); - histos.fill(HIST("QAMC/hInvmassK892_Rho"), temp13.M(), lResonanceSecondary.M()); - histos.fill(HIST("QAMC/hInvmassSecon_PiKa"), lResonanceSecondary.M(), temp23.M()); - histos.fill(HIST("QAMC/hInvmassSecon"), lResonanceSecondary.M()); - histos.fill(HIST("QAMC/hpT_Seocondary"), lResonanceSecondary.Pt()); - - if constexpr (!IsResoMicrotrack) { - - auto trk1NSigmaPiTPC = trk1.tpcNSigmaPi(); - auto trk1NSigmaPiTOF = (isTrk1hasTOF) ? trk1.tofNSigmaPi() : -999.; - auto trk2NSigmaPiTPC = trk2.tpcNSigmaPi(); - auto trk2NSigmaPiTOF = (isTrk2hasTOF) ? trk2.tofNSigmaPi() : -999.; - - // PID QA primary pion - histos.fill(HIST("QAMC/trkppionTPCPID"), trk1pt, trk1NSigmaPiTPC); - if (isTrk1hasTOF) { - histos.fill(HIST("QAMC/trkppionTOFPID"), trk1pt, trk1NSigmaPiTOF); - histos.fill(HIST("QAMC/trkppionTPCTOFPID"), trk1NSigmaPiTPC, trk1NSigmaPiTOF); - } - histos.fill(HIST("QAMC/trkppionpT"), trk1pt); - histos.fill(HIST("QAMC/trkppionDCAxy"), trk1.dcaXY()); - histos.fill(HIST("QAMC/trkppionDCAz"), trk1.dcaZ()); - - // PID QA secondary pion - histos.fill(HIST("QAMC/trkspionTPCPID"), trk2pt, trk2NSigmaPiTPC); - if (isTrk2hasTOF) { - histos.fill(HIST("QAMC/trkspionTOFPID"), trk2pt, trk2NSigmaPiTOF); - histos.fill(HIST("QAMC/trkspionTPCTOFPID"), trk2NSigmaPiTPC, trk2NSigmaPiTOF); - } - histos.fill(HIST("QAMC/trkspionpT"), trk2pt); - histos.fill(HIST("QAMC/trkspionDCAxy"), trk2.dcaXY()); - histos.fill(HIST("QAMC/trkspionDCAz"), trk2.dcaZ()); - - } else { - - histos.fill(HIST("QAMC/trkppionTPCPID"), trk1pt, o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk1.pidNSigmaSelectionFlags())); - if (isTrk1hasTOF) { - histos.fill(HIST("QAMC/trkppionTOFPID"), trk1pt, o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk1.pidNSigmaSelectionFlags())); - histos.fill(HIST("QAMC/trkppionTPCTOFPID"), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk1.pidNSigmaSelectionFlags()), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk1.pidNSigmaSelectionFlags())); - } - histos.fill(HIST("QAMC/trkppionpT"), trk1pt); - histos.fill(HIST("QAMC/trkppionDCAxy"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(trk1.trackSelectionFlags())); - histos.fill(HIST("QAMC/trkppionDCAz"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(trk1.trackSelectionFlags())); - - // PID QA secondary pion - histos.fill(HIST("QAMC/trkspionTPCPID"), trk2pt, o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk2.pidNSigmaSelectionFlags())); - if (isTrk2hasTOF) { - histos.fill(HIST("QAMC/trkspionTOFPID"), trk2pt, o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk2.pidNSigmaSelectionFlags())); - histos.fill(HIST("QAMC/trkspionTPCTOFPID"), o2::aod::resomicrodaughter::PidNSigma::getTPCnSigma(trk2.pidNSigmaSelectionFlags()), o2::aod::resomicrodaughter::PidNSigma::getTOFnSigma(trk2.pidNSigmaSelectionFlags())); - } - histos.fill(HIST("QAMC/trkspionpT"), trk2pt); - histos.fill(HIST("QAMC/trkspionDCAxy"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAxy(trk2.trackSelectionFlags())); - histos.fill(HIST("QAMC/trkspionDCAz"), o2::aod::resomicrodaughter::ResoMicroTrackSelFlag::decodeDCAz(trk2.trackSelectionFlags())); - } - } else { - histos.fill(HIST("k1invmass_MC_noK1"), lResonanceK1.M()); - } - } // IsMC - } else { - unsigned int typeK1 = bTrack.sign() > 0 ? BinType::kK1P_Mix : BinType::kK1N_Mix; - unsigned int typeNormal = BinAnti::kNormal; - histos.fill(HIST("hInvmass_K1_Mix"), typeNormal, typeK1, multiplicity, lResonanceK1.Pt(), lResonanceK1.M()); - histos.fill(HIST("k1invmass_Mix"), lResonanceK1.M()); + if (truthDebugCounts[channelBin] < histogramOptions.cfgTruthDebug) { + ++truthDebugCounts[channelBin]; + LOGP(info, "K1Truth channel={} collision={} tracks=({},{},{}) pdg=({},{},{}) mothers=({},{},{}) motherPDG=({},{},{}) siblings=(({},{}),({},{}),({},{}))", + channelBin, collision.globalIndex(), + pion1.globalIndex(), pion2.globalIndex(), kaon.globalIndex(), + pion1.pdgCode(), pion2.pdgCode(), kaon.pdgCode(), pion1.motherId(), pion2.motherId(), kaon.motherId(), + pion1.motherPDG(), pion2.motherPDG(), kaon.motherPDG(), + pion1.siblingIds()[0], pion1.siblingIds()[1], pion2.siblingIds()[0], pion2.siblingIds()[1], kaon.siblingIds()[0], kaon.siblingIds()[1]); } - } // bTrack - } - } // fillHistograms + typeK1 = kaon.sign() > 0 ? BinType::kK1P_Rec : BinType::kK1N_Rec; + histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, c.k1.Pt(), c.k1.M()); + histos.fill(HIST("k1invmass_MC"), c.k1.M()); + histos.fill(HIST("QAMC/K1OA"), c.angle); + histos.fill(HIST("QAMC/K1PairAsym"), c.pairAsym); + histos.fill(HIST("QAMC/hInvmassK892_Rho"), c.mass13, c.secondary.M()); + histos.fill(HIST("QAMC/hInvmassSecon_PiKa"), c.secondary.M(), c.mass23); + histos.fill(HIST("QAMC/hInvmassSecon"), c.secondary.M()); + histos.fill(HIST("QAMC/hpT_Secondary"), c.secondary.Pt()); + + // PID QA primary and secondary pion + fillPionQA(pion1, true); + fillPionQA(pion2, false); + fillKaonQA(kaon); + } + }; + + core.forEachCandidate(histos, collision, dTracks1, dTracks2, IsMC || fillQA, onPair, onCandidate); + } - void processResoTracks(aod::ResoCollision const& collision, - aod::ResoTracks const& resotracks) + void processResoTracks(ResoCollisions::iterator const& collision, + ResoTracks const& resotracks) { + if (!core.passesEventCuts(collision)) { + return; + } fillHistograms(collision, resotracks, resotracks); } PROCESS_SWITCH(K1AnalysisMicro, processResoTracks, "Process ResoTracks", false); - void processResoMicroTracks(aod::ResoCollision const& collision, - aod::ResoMicroTracks const& resomicrotracks) + void processResoMicroTracks(ResoCollisions::iterator const& collision, + ResoMicroTracks const& resomicrotracks) { + if (!core.passesEventCuts(collision)) { + return; + } fillHistograms(collision, resomicrotracks, resomicrotracks); } PROCESS_SWITCH(K1AnalysisMicro, processResoMicroTracks, "Process ResoMicroTracks", true); - void processMC(aod::ResoCollision const& collision, - soa::Join const& resotracks) + void processMC(ResoMCCols::iterator const& collision, + ResoMCTracks const& resotracks) { + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + histos.fill(HIST("MCReco/collisions"), 0.5); fillHistograms(collision, resotracks, resotracks); } PROCESS_SWITCH(K1AnalysisMicro, processMC, "Process Event for MC", false); - void processMCTrue(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents) + void processMCMicro(ResoMCCols::iterator const& collision, ResoMCMicroTracks const& tracks) { - auto multiplicity = collision.cent(); - for (const auto& part : resoParents) { - if (std::abs(part.pdgCode()) != kK1Plus) - continue; - if (std::abs(part.y()) > 0.5) { - continue; - } - bool pass1 = false; - bool pass2 = false; - bool pass3 = false; - bool pass4 = false; - if (std::abs(part.daughterPDG1()) == 313 || std::abs(part.daughterPDG2()) == 313) { // At least one decay into K892 - pass2 = true; - } - if (std::abs(part.daughterPDG1()) == kPiPlus || std::abs(part.daughterPDG2()) == kPiPlus) { // At lest one decay into pion - pass1 = true; - } - if (std::abs(part.daughterPDG1()) == kPDGRho770 || std::abs(part.daughterPDG2()) == kPDGRho770) { - pass4 = true; - } - if (std::abs(part.daughterPDG1()) == kKPlus || std::abs(part.daughterPDG2()) == kKPlus) { - pass3 = true; - } - if (!pass1 || !pass2 || !pass3 || !pass4) // If we have both decay products - continue; - auto typeNormal = part.pdgCode() > 0 ? BinAnti::kNormal : BinAnti::kAnti; - if (collision.isVtxIn10()) // INEL>10 - { - auto typeK1 = part.pdgCode() > 0 ? BinType::kK1P_GenINEL10 : BinType::kK1N_GenINEL10; - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, part.pt(), 1); - } - if (collision.isVtxIn10() && collision.isInSel8()) // INEL>10, vtx10 - { - auto typeK1 = part.pdgCode() > 0 ? BinType::kK1P_GenINELgt10 : BinType::kK1N_GenINELgt10; - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, part.pt(), 1); - } - if (collision.isVtxIn10() && collision.isTriggerTVX()) // vtx10, TriggerTVX - { - auto typeK1 = part.pdgCode() > 0 ? BinType::kK1P_GenTrig10 : BinType::kK1N_GenTrig10; - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, part.pt(), 1); - } - if (collision.isInAfterAllCuts()) // after all event selection - { - auto typeK1 = part.pdgCode() > 0 ? BinType::kK1P_GenEvtSel : BinType::kK1N_GenEvtSel; - histos.fill(HIST("hInvmass_K1"), typeNormal, typeK1, multiplicity, part.pt(), 1); - } + // The modular producer already selected these reconstructed collisions. + // Apply precisely the same reconstruction loop as the data baseline. + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; } + histos.fill(HIST("MCReco/collisions"), 0.5); + histos.fill(HIST("MCReco/microTracks"), 0.5, tracks.size()); + fillHistograms(collision, tracks, tracks); } - PROCESS_SWITCH(K1AnalysisMicro, processMCTrue, "Process Event for MC", false); + PROCESS_SWITCH(K1AnalysisMicro, processMCMicro, "Process reconstructed MC with micro v001 tables", false); + + void processMCTrue(ResoMCCols::iterator const& collision, ResoMCParents const& resoParents) + { + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + // Keep other/unresolved immediate decays too; never require both pairs. + core.forEachGeneratedK1(histos, resoParents, [&](auto const& part, K1TruthChannel channel) { + const int charge = part.pdgCode() > 0 ? 1 : -1; + histos.fill(HIST("MCGen/chargeChannel"), charge, static_cast(channel)); + histos.fill(HIST("MCGen/ptChannel"), static_cast(channel), part.pt()); + }); + } + PROCESS_SWITCH(K1AnalysisMicro, processMCTrue, "Process generated K1 in selected events with v001 parents", false); // Processing Event Mixing using BinningTypeVtxZT0M = ColumnBinningPolicy; - void processME(o2::aod::ResoCollisions const& collisions, aod::ResoTracks const& resotracks) + void processME(ResoCollisions const& collisions, ResoTracks const& resotracks) { auto tracksTuple = std::make_tuple(resotracks); BinningTypeVtxZT0M colBinning{{cfgVtxBins, cfgMultBins}, true}; - SameKindPair pairs{colBinning, nEvtMixing, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + SameKindPair pairs{colBinning, nEvtMixing, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { + if (!core.passesEventCuts(collision1) || !core.passesEventCuts(collision2)) { + continue; + } fillHistograms(collision1, tracks1, tracks2); } }; PROCESS_SWITCH(K1AnalysisMicro, processME, "Process EventMixing light without partition", false); // Processing Event Mixing -- Micro - // using BinningTypeVtxZT0M = ColumnBinningPolicy; - void processMEMicro(o2::aod::ResoCollisions const& collisions, aod::ResoMicroTracks const& resomicrotracks) + void processMEMicro(ResoCollisions const& collisions, ResoMicroTracks const& resomicrotracks) { auto tracksTuple = std::make_tuple(resomicrotracks); BinningTypeVtxZT0M colBinning{{cfgVtxBins, cfgMultBins}, true}; - SameKindPair pairs{colBinning, nEvtMixing, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + SameKindPair pairs{colBinning, nEvtMixing, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { + if (!core.passesEventCuts(collision1) || !core.passesEventCuts(collision2)) { + continue; + } fillHistograms(collision1, tracks1, tracks2); } }; - PROCESS_SWITCH(K1AnalysisMicro, processMEMicro, "Process EventMixing light without partition", true); -}; // struct + PROCESS_SWITCH(K1AnalysisMicro, processMEMicro, "Process EventMixing light without partition", false); +}; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { diff --git a/PWGLF/Tasks/Resonances/k1TrainingTable.cxx b/PWGLF/Tasks/Resonances/k1TrainingTable.cxx new file mode 100644 index 00000000000..8b070a97be1 --- /dev/null +++ b/PWGLF/Tasks/Resonances/k1TrainingTable.cxx @@ -0,0 +1,281 @@ +// 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 k1TrainingTable.cxx +/// \brief Derived K1 training table workflow +/// \author Bong-Hwi Lim +/// + +#include "PWGLF/Core/K1AnalysisMicroCore.h" +#include "PWGLF/Core/K1MlFeatures.h" +#include "PWGLF/Core/ResoAnalysisSelectionCore.h" +#include "PWGLF/DataModel/LFK1MlTables.h" +#include "PWGLF/DataModel/LFResonanceTables.h" + +#include +#include +#include +#include +#include +#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 + +using namespace o2; +using namespace o2::framework; +using namespace o2::constants::physics; +using namespace o2::analysis::resonance; +using namespace o2::analysis::k1micro; + +static_assert(std::extent_v == o2::analysis::k1ml::NMasterFeatures, + "K1MasterFeatures column size must match the K1 ML feature contract"); + +/// Writes the unlike-sign K1 micro candidates of the shared K1 selection, with the +/// canonical (kaon, same-sign pion, opposite-sign pion) tracks and the master features. +struct K1TrainingTable { + using ResoCollisions = aod::ResoCollisions_001; + using ResoMCCols = soa::Join; + using ResoMicroTracks = aod::ResoMicroTracks_001; + using ResoMCMicroTracks = 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; + + Produces k1MlEvents; + Produces k1MlTracks; + Produces k1MlCandidates; + Produces k1MlInputs; + Produces k1MlTruth; + Produces k1MlGenAudit; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Selection shared with the K1 histogram task (plain JSON keys, no group prefix) + EventCuts eventCuts; + TrackCuts trackCuts; + PionPidCuts pionPID; + KaonPidCuts kaonPID; + SecondaryCuts secondaryCuts; + CandidateCuts candidateCuts; + + Configurable k1MlExportStage{"k1MlExportStage", "loose", "Candidate export stage: loose (pass bits >= 1) or selected (pass bits = 31)"}; + Configurable k1MlLooseAudit{"k1MlLooseAudit", true, "Record loose US cutflow and mass/pT/activity spectrum"}; + Configurable k1MlParityRows{"k1MlParityRows", 0, "Log feature hashes of the first N exported candidates"}; + + K1AnalysisMicroCore core; + int64_t k1MlEventRow = -1; + std::unordered_map k1MlTrackRows; + int k1MlParityLogged = 0; + + void init(InitContext&) + { + if (static_cast(doprocessResoMicroTracks) + static_cast(doprocessMCMicro) != 1 || + (doprocessMCTrue && !doprocessMCMicro)) { + LOG(fatal) << "K1 training table requires exactly one of processResoMicroTracks and processMCMicro; processMCTrue requires processMCMicro"; + } + if (k1MlParityRows < 0) { + LOG(fatal) << "k1MlParityRows must not be negative"; + } + if (k1MlExportStage.value != "loose" && k1MlExportStage.value != "selected") { + LOG(fatal) << "k1MlExportStage must be loose or selected"; + } + + ProcessModes modes; + modes.microTracks = true; + modes.mcReco = doprocessMCMicro; + modes.mcRecoMicro = doprocessMCMicro; + modes.mcGen = doprocessMCTrue; + LooseStageOptions looseOptions; + looseOptions.audit = k1MlLooseAudit; + looseOptions.exportSelected = k1MlExportStage.value == "selected"; + core.init(histos, eventCuts, trackCuts, pionPID, kaonPID, secondaryCuts, candidateCuts, modes, looseOptions); + if (doprocessMCMicro) { + histos.add("MCReco/collisions", "Selected reconstructed MC collisions", HistType::kTH1D, {{1, 0, 1}}); + histos.add("MCReco/microTracks", "Input micro tracks in selected MC collisions", HistType::kTH1D, {{1, 0, 1}}); + } + + // Candidates that violate the canonical or feature contract are skipped, never written. + auto skipped = histos.add("ML/exportSkipped", "Skipped candidates;K1 ML build status;candidates", HistType::kTH1D, {{6, -0.5, 5.5}}); + const std::array statusLabels{"Ok", "InvalidChargePattern", "ReusedTrack", "InvalidMomentum", "InvalidKinematics", "InvalidContract"}; + for (std::size_t i = 0; i < statusLabels.size(); ++i) { + skipped->GetXaxis()->SetBinLabel(i + 1, statusLabels[i]); + } + + LOG(info) << "Size of the histograms in K1 training table task"; + histos.print(); + } + + template + int64_t writeK1MlTrack(Track const& track) + { + const auto id = static_cast(track.globalIndex()); + if (auto it = k1MlTrackRows.find(id); it != k1MlTrackRows.end()) { + return it->second; + } + k1MlTracks(k1MlEventRow, 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 int64_t row = k1MlTracks.lastIndex(); + k1MlTrackRows.emplace(id, row); + return row; + } + + template + void writeK1MlEvent(Collision const& collision) + { + k1MlTrackRows.clear(); + // ResoCollisions_001 carries no run number or BC; the reduced collision row identifies the event within its DF. + k1MlEvents(collision.globalIndex(), + collision.posZ(), collision.bMagField(), collision.cent(), collision.multiplicity(), collision.isRecINELgt0()); + k1MlEventRow = k1MlEvents.lastIndex(); + } + + template + void logParity(Collision const& collision, o2::analysis::k1ml::CandidateSnapshot const& candidate, o2::analysis::k1ml::FeaturePack const& pack) + { + if (k1MlParityLogged >= k1MlParityRows) { + return; + } + uint64_t hash = FnvOffsetBasis; + for (const auto& value : pack.master) { + const auto bits = std::bit_cast(value); + for (unsigned int shift = 0; shift < BitsPerFloat; shift += BitsPerByte) { + hash = (hash ^ ((bits >> shift) & ByteMask)) * FnvPrime; + } + } + LOGP(info, "K1MLPARITY collision={} tracks={},{},{} featureHash={}", + collision.globalIndex(), candidate.tracks[0].sourceTrackId, candidate.tracks[1].sourceTrackId, + candidate.tracks[2].sourceTrackId, hash); + ++k1MlParityLogged; + } + + template + void writeK1MlCandidate(Collision const& collision, Kaon const& kaon, Pion const& samePion, Pion const& oppPion, + K1TruthChannel channel, uint16_t passBits) + { + using o2::analysis::k1ml::BuildStatus; + const auto canonical = o2::analysis::k1ml::canonicalizeUS(o2::analysis::k1ml::makeTrackSnapshot(kaon), o2::analysis::k1ml::makeTrackSnapshot(samePion), o2::analysis::k1ml::makeTrackSnapshot(oppPion)); + if (canonical.status != BuildStatus::Ok) { + histos.fill(HIST("ML/exportSkipped"), static_cast(canonical.status)); + return; + } + const auto pack = o2::analysis::k1ml::buildMasterFeatures(canonical.candidate); + if (pack.status != BuildStatus::Ok) { + histos.fill(HIST("ML/exportSkipped"), static_cast(pack.status)); + return; + } + logParity(collision, canonical.candidate, pack); + const auto kaonRow = writeK1MlTrack(kaon); + const auto sameRow = writeK1MlTrack(samePion); + const auto oppRow = writeK1MlTrack(oppPion); + ROOT::Math::PxPyPzMVector k{kaon.px(), kaon.py(), kaon.pz(), MassKaonCharged}; + ROOT::Math::PxPyPzMVector s{samePion.px(), samePion.py(), samePion.pz(), MassPionCharged}; + ROOT::Math::PxPyPzMVector o{oppPion.px(), oppPion.py(), oppPion.pz(), MassPionCharged}; + const auto mother = k + s + o; + k1MlCandidates(k1MlEventRow, kaonRow, sameRow, oppRow, + static_cast(mother.M()), static_cast((s + o).M()), + static_cast((k + s).M()), static_cast((k + o).M()), + pack.kinematics.scalarSumPt, static_cast((s + o).Pt()), + static_cast(mother.Pt()), static_cast(mother.Rapidity()), + static_cast(mother.Eta()), static_cast(mother.Phi()), + static_cast(kaon.sign()), passBits); + const int64_t row = k1MlCandidates.lastIndex(); + k1MlInputs(row, pack.master.data(), static_cast(pack.status)); + if constexpr (IsMC) { + const bool matched = channel != K1TruthChannel::None; + int motherId = -1; + if (matched) { + motherId = channel == K1TruthChannel::RhoK ? kaon.motherId() : std::abs(samePion.motherPDG()) == Pdg::kK1_1270Plus ? samePion.motherId() + : oppPion.motherId(); + } + // Immediate-mother information alone does not prove a physical UID. + k1MlTruth(row, matched ? uint8_t{1} : uint8_t{2}, static_cast(channel), + matched ? kaon.sign() * Pdg::kK1_1270Plus : 0, static_cast(motherId), + std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); + } else { + k1MlTruth(row, uint8_t{0}, uint8_t{0}, 0, int64_t{-1}, + std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); + } + } + + void processResoMicroTracks(ResoCollisions::iterator const& collision, ResoMicroTracks const& tracks) + { + if (!core.passesEventCuts(collision)) { + return; + } + writeK1MlEvent(collision); + // Selection only: the K1 analysis histograms belong to the K1 analysis task + core.forEachCandidate(histos, collision, tracks, tracks, false, nullptr, nullptr, + [this](auto const& coll, auto const& kaon, auto const& samePion, auto const& oppPion, + K1TruthChannel channel, uint16_t passBits) { + writeK1MlCandidate(coll, kaon, samePion, oppPion, channel, passBits); + }); + } + PROCESS_SWITCH(K1TrainingTable, processResoMicroTracks, "Write K1 candidates from data micro v001 tables", true); + + void processMCMicro(ResoMCCols::iterator const& collision, ResoMCMicroTracks const& tracks) + { + // The modular producer already selected these reconstructed collisions. + if (!core.passesEventCuts(collision) || !core.passesMCEventCuts(collision)) { + return; + } + histos.fill(HIST("MCReco/collisions"), 0.5); + histos.fill(HIST("MCReco/microTracks"), 0.5, tracks.size()); + writeK1MlEvent(collision); + core.forEachCandidate(histos, collision, tracks, tracks, false, nullptr, nullptr, + [this](auto const& coll, auto const& kaon, auto const& samePion, auto const& oppPion, + K1TruthChannel channel, uint16_t passBits) { + writeK1MlCandidate(coll, kaon, samePion, oppPion, channel, passBits); + }); + } + PROCESS_SWITCH(K1TrainingTable, processMCMicro, "Write K1 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.forEachGeneratedK1(histos, resoParents, [&](auto const& part, K1TruthChannel channel) { + k1MlGenAudit(collision.globalIndex(), static_cast(part.originalMcParticleId()), + part.pdgCode(), part.daughterPDG1(), part.daughterPDG2(), static_cast(channel), + part.pt(), part.y(), true); + }); + } + PROCESS_SWITCH(K1TrainingTable, processMCTrue, "Write generated K1 parents of selected reconstructed MC events", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}