From aeb0e7fc46df0d27b76a1e8c945aec5e8adbff5a Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Wed, 22 Jul 2026 18:55:23 +0200 Subject: [PATCH 1/5] updated tasks to run on derived data for phi and kstar --- PWGLF/DataModel/LFCKSSpinalignmentTables.h | 113 ++- .../Resonances/cksspinalignment.cxx | 689 +++++++------ PWGLF/Tasks/Resonances/cksspinalignder.cxx | 939 +++++------------- PWGLF/Tasks/Resonances/phiflowder.cxx | 299 +++++- 4 files changed, 976 insertions(+), 1064 deletions(-) diff --git a/PWGLF/DataModel/LFCKSSpinalignmentTables.h b/PWGLF/DataModel/LFCKSSpinalignmentTables.h index 7988f76da31..0571c108b4e 100644 --- a/PWGLF/DataModel/LFCKSSpinalignmentTables.h +++ b/PWGLF/DataModel/LFCKSSpinalignmentTables.h @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file LFCKSSpinalignmentTables.h -/// \brief DataModel for Charged KStar spin alignment +/// \brief DataModel for reduced K0s and pion tables for charged KStar analysis /// /// \author Prottay Das @@ -20,7 +20,6 @@ #include #include -#include #include namespace o2::aod @@ -29,77 +28,89 @@ namespace kshortpionevent { DECLARE_SOA_COLUMN(Cent, cent, float); DECLARE_SOA_COLUMN(Posz, posz, float); -// DECLARE_SOA_COLUMN(CollIndex, collIndex, float); DECLARE_SOA_COLUMN(PsiFT0C, psiFT0C, float); DECLARE_SOA_COLUMN(PsiFT0A, psiFT0A, float); DECLARE_SOA_COLUMN(PsiTPC, psiTPC, float); } // namespace kshortpionevent + DECLARE_SOA_TABLE(KShortpionEvents, "AOD", "KSHORTPIONEVENT", o2::soa::Index<>, kshortpionevent::Cent, kshortpionevent::Posz, - // kshortpionevent::CollIndex, kshortpionevent::PsiFT0C, kshortpionevent::PsiFT0A, - kshortpionevent::PsiTPC) + kshortpionevent::PsiTPC); + using KShortpionEvent = KShortpionEvents::iterator; -namespace kshortpionpair +// ------------------------------------------------------------ +// K0s candidate table +// ------------------------------------------------------------ +namespace kshorttrack { DECLARE_SOA_INDEX_COLUMN(KShortpionEvent, kshortpionevent); -DECLARE_SOA_COLUMN(V0Cospa, v0Cospa, float); //! V0 Cospa -DECLARE_SOA_COLUMN(V0Radius, v0Radius, float); //! V0 Radius -DECLARE_SOA_COLUMN(DcaPositive, dcaPositive, float); //! DCA Positive -DECLARE_SOA_COLUMN(DcaNegative, dcaNegative, float); //! DCA Negative -DECLARE_SOA_COLUMN(DcaBetweenDaughter, dcaBetweenDaughter, float); //! DCA between daughters -DECLARE_SOA_COLUMN(V0Lifetime, v0Lifetime, float); //! KShort lifetime -DECLARE_SOA_COLUMN(KShortPx, kShortPx, float); //! KShort Px -DECLARE_SOA_COLUMN(KShortPy, kShortPy, float); //! KShort Py -DECLARE_SOA_COLUMN(KShortPz, kShortPz, float); //! KShort Pz -DECLARE_SOA_COLUMN(KShortMass, kShortMass, float); //! KShort Mass -DECLARE_SOA_COLUMN(PionBachPx, pionBachPx, float); //! Bachelor Pion Px -DECLARE_SOA_COLUMN(PionBachPy, pionBachPy, float); //! Bachelor Pion Py -DECLARE_SOA_COLUMN(PionBachPz, pionBachPz, float); //! Bachelor Pion Pz -/*DECLARE_SOA_COLUMN(PionBachTPC, pionBachTPC, float); //! Bachelor Pion nsigmatpc -DECLARE_SOA_COLUMN(PionBachTOFHit, pionBachTOFHit, int); //! Bachelor Pion tof hit availability -DECLARE_SOA_COLUMN(PionBachTOF, pionBachTOF, float); //! Bachelor Pion nsigmatof*/ -DECLARE_SOA_COLUMN(PionBachIndex, pionBachIndex, int); //! Bachelor Pion index -DECLARE_SOA_COLUMN(PionIndex1, pionIndex1, int); //! Daughter Pion index1 -DECLARE_SOA_COLUMN(PionIndex2, pionIndex2, int); //! Daughter Pion index2 -DECLARE_SOA_COLUMN(PionPidMask, pionPidMask, uint8_t); //! bitmask for PID selections -} // namespace kshortpionpair + +DECLARE_SOA_COLUMN(V0Cospa, v0Cospa, float); +DECLARE_SOA_COLUMN(V0Radius, v0Radius, float); +DECLARE_SOA_COLUMN(DcaPositive, dcaPositive, float); +DECLARE_SOA_COLUMN(DcaNegative, dcaNegative, float); +DECLARE_SOA_COLUMN(DcaBetweenDaughter, dcaBetweenDaughter, float); +// DECLARE_SOA_COLUMN(V0Lifetime, v0Lifetime, float); + +DECLARE_SOA_COLUMN(KShortPx, kShortPx, float); +DECLARE_SOA_COLUMN(KShortPy, kShortPy, float); +DECLARE_SOA_COLUMN(KShortPz, kShortPz, float); +DECLARE_SOA_COLUMN(KShortMass, kShortMass, float); + +DECLARE_SOA_COLUMN(PionIndex1, pionIndex1, int64_t); +DECLARE_SOA_COLUMN(PionIndex2, pionIndex2, int64_t); +} // namespace kshorttrack + DECLARE_SOA_TABLE(KShortTracks, "AOD", "KSHORTTRACK", o2::soa::Index<>, - kshortpionpair::KShortpionEventId, - kshortpionpair::V0Cospa, - kshortpionpair::V0Radius, - kshortpionpair::DcaPositive, - kshortpionpair::DcaNegative, - kshortpionpair::DcaBetweenDaughter, - kshortpionpair::V0Lifetime, - // kshortpionpair::Armenteros, - kshortpionpair::KShortPx, - kshortpionpair::KShortPy, - kshortpionpair::KShortPz, - kshortpionpair::KShortMass, - kshortpionpair::PionIndex1, - kshortpionpair::PionIndex2); + kshorttrack::KShortpionEventId, + kshorttrack::V0Cospa, + kshorttrack::V0Radius, + kshorttrack::DcaPositive, + kshorttrack::DcaNegative, + kshorttrack::DcaBetweenDaughter, + // kshorttrack::V0Lifetime, + kshorttrack::KShortPx, + kshorttrack::KShortPy, + kshorttrack::KShortPz, + kshorttrack::KShortMass, + kshorttrack::PionIndex1, + kshorttrack::PionIndex2); using KShortTrack = KShortTracks::iterator; +// ------------------------------------------------------------ +// Bachelor pion table +// ------------------------------------------------------------ +namespace piontrack +{ +DECLARE_SOA_INDEX_COLUMN(KShortpionEvent, kshortpionevent); + +DECLARE_SOA_COLUMN(PionBachPx, pionBachPx, float); +DECLARE_SOA_COLUMN(PionBachPy, pionBachPy, float); +DECLARE_SOA_COLUMN(PionBachPz, pionBachPz, float); +DECLARE_SOA_COLUMN(Charge, charge, int8_t); +DECLARE_SOA_COLUMN(PionBachIndex, pionBachIndex, int64_t); +DECLARE_SOA_COLUMN(PionPidMask, pionPidMask, uint8_t); +} // namespace piontrack + DECLARE_SOA_TABLE(PionTracks, "AOD", "PIONTRACK", o2::soa::Index<>, - kshortpionpair::KShortpionEventId, - kshortpionpair::PionBachPx, - kshortpionpair::PionBachPy, - kshortpionpair::PionBachPz, - // kshortpionpair::PionBachSign, - // kshortpionpair::PionBachTPC, - // kshortpionpair::PionBachTOFHit, - // kshortpionpair::PionBachTOF, - kshortpionpair::PionBachIndex, - kshortpionpair::PionPidMask); + piontrack::KShortpionEventId, + piontrack::PionBachPx, + piontrack::PionBachPy, + piontrack::PionBachPz, + piontrack::Charge, + piontrack::PionBachIndex, + piontrack::PionPidMask); using PionTrack = PionTracks::iterator; + } // namespace o2::aod + #endif // PWGLF_DATAMODEL_LFCKSSPINALIGNMENTTABLES_H_ diff --git a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx index 014002da82b..a92d1a94329 100644 --- a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx +++ b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file cksspinalignment.cxx -/// \brief Table producer for Charged KStar spin alignment +/// \brief Reduced table producer for K0s and charged pions for later charged K* reconstruction /// /// \author prottay.das@cern.ch @@ -28,7 +28,6 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" -#include #include #include #include @@ -41,13 +40,11 @@ #include #include -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) -#include +#include -#include +#include #include #include -#include #include using namespace o2; @@ -61,8 +58,6 @@ struct cksspinalignment { Produces kshortTrack; Produces pionTrack; - Service ccdb; - struct : ConfigurableGroup { Configurable requireRCTFlagChecker{"requireRCTFlagChecker", true, "Check event quality in run condition table"}; Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; @@ -74,215 +69,268 @@ struct cksspinalignment { Configurable useGoodITSLayersAll{"useGoodITSLayersAll", true, "Apply kIsGoodITSLayersAll selection bit"}; Configurable cfgCutOccupancy{"cfgCutOccupancy", 2000, "Occupancy cut"}; - // events Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; - Configurable cfgCutCentralityMax{"cfgCutCentralityMax", 80.0f, "Accepted maximum Centrality"}; - Configurable cfgCutCentralityMin{"cfgCutCentralityMin", 0.0f, "Accepted minimum Centrality"}; + Configurable cfgCutCentralityMax{"cfgCutCentralityMax", 80.0f, "Accepted maximum centrality"}; + Configurable cfgCutCentralityMin{"cfgCutCentralityMin", 0.0f, "Accepted minimum centrality"}; - // Configs for track - Configurable cfgCutPt{"cfgCutPt", 0.2, "Pt cut on daughter track"}; - Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; + Configurable cfgCutPt{"cfgCutPt", 0.2f, "Prefilter minimum pT for bachelor pion"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Prefilter maximum eta for bachelor pion"}; - // Configs for pion struct : ConfigurableGroup { - Configurable itsPIDSelection{"itsPIDSelection", true, "PID ITS"}; - Configurable lowITSPIDNsigma{"lowITSPIDNsigma", -3.0, "lower cut on PID nsigma for ITS"}; - Configurable highITSPIDNsigma{"highITSPIDNsigma", 3.0, "higher cut on PID nsigma for ITS"}; - Configurable itsclusterPiMeson{"itsclusterPiMeson", 5, "Minimum number of ITS cluster for pi meson track"}; - Configurable tpcCrossedRowsPiMeson{"tpcCrossedRowsPiMeson", 80, "Minimum number of TPC Crossed Rows for pi meson track"}; - Configurable cutDCAxyPiMeson{"cutDCAxyPiMeson", 0.1, "Maximum DCAxy for pi meson track"}; - Configurable cutDCAzPiMeson{"cutDCAzPiMeson", 0.1, "Maximum DCAz for pi meson track"}; - Configurable cutEtaPiMeson{"cutEtaPiMeson", 0.8, "Maximum eta for pi meson track"}; - Configurable cutPTPiMeson{"cutPTPiMeson", 0.8, "Maximum pt for pi meson track"}; - Configurable usePID{"usePID", true, "Flag for using PID selection for pi meson track"}; - Configurable nsigmaCutTPCPiMeson{"nsigmaCutTPCPiMeson", 3.0, "Maximum nsigma cut TPC for pi meson track"}; - Configurable nsigmaCutTOFPiMeson{"nsigmaCutTOFPiMeson", 3.0, "Maximum nsigma cut TOF for pi meson track"}; - Configurable cutTOFBetaPiMeson{"cutTOFBetaPiMeson", 3.0, "Maximum beta cut for pi meson track"}; + Configurable itsPIDSelection{"itsPIDSelection", true, "Apply ITS PID for bachelor pion"}; + Configurable lowITSPIDNsigma{"lowITSPIDNsigma", -3.0f, "Lower ITS n-sigma pion cut"}; + Configurable highITSPIDNsigma{"highITSPIDNsigma", 3.0f, "Upper ITS n-sigma pion cut"}; + Configurable itsclusterPiMeson{"itsclusterPiMeson", 5, "Minimum ITS clusters for bachelor pion"}; + Configurable tpcCrossedRowsPiMeson{"tpcCrossedRowsPiMeson", 80, "Minimum TPC crossed rows for bachelor pion"}; + Configurable cutDCAxyPiMeson{"cutDCAxyPiMeson", 0.1f, "Maximum DCAxy for bachelor pion"}; + Configurable cutDCAzPiMeson{"cutDCAzPiMeson", 0.1f, "Maximum DCAz for bachelor pion"}; + Configurable cutEtaPiMeson{"cutEtaPiMeson", 0.8f, "Maximum eta for bachelor pion"}; + Configurable cutPTPiMeson{"cutPTPiMeson", 0.2f, "Minimum pT for bachelor pion"}; + Configurable usePID{"usePID", true, "Apply pion PID"}; + Configurable nsigmaCutTPCPiMeson{"nsigmaCutTPCPiMeson", 3.0f, "Maximum TPC n-sigma pion cut"}; + Configurable nsigmaCutTOFPiMeson{"nsigmaCutTOFPiMeson", 3.0f, "Maximum TOF n-sigma pion cut"}; + Configurable cutTOFBetaPiMeson{"cutTOFBetaPiMeson", 0.5f, "Minimum TOF beta for bachelor pion"}; + Configurable pSwitchPID{"pSwitchPID", 0.5f, "pT switch for pT-dependent pion PID"}; } grpPion; enum PionPidBits : uint8_t { - kPID1 = 1u << 0, // selectionPID - kPID2 = 1u << 1, // selectionPID2 - kPID3 = 1u << 2, // selectionPID3 - kPID4 = 1u << 3 // selectionPID4 + kPID1 = 1u << 0, + kPID2 = 1u << 1, + kPID3 = 1u << 2, + kPID4 = 1u << 3 }; - // Configs for V0 - Configurable confV0PtMin{"confV0PtMin", 0.f, "Minimum transverse momentum of V0"}; - Configurable confV0PtMax{"confV0PtMax", 1000.f, "Maximum transverse momentum of V0"}; - Configurable confV0Rap{"confV0Rap", 0.8f, "Rapidity range of V0"}; - Configurable confV0DCADaughMax{"confV0DCADaughMax", 1.0f, "Maximum DCA between the V0 daughters"}; - Configurable confV0CPAMin{"confV0CPAMin", 0.9998f, "Minimum CPA of V0"}; - Configurable confV0TranRadV0Min{"confV0TranRadV0Min", 1.5f, "Minimum transverse radius"}; - Configurable confV0TranRadV0Max{"confV0TranRadV0Max", 100.f, "Maximum transverse radius"}; - Configurable cMaxV0DCA{"cMaxV0DCA", 1.2, "Maximum V0 DCA to PV"}; - Configurable cMinV0DCAPi{"cMinV0DCAPi", 0.05, "Minimum V0 daughters DCA to PV for Pi"}; - Configurable cMaxV0LifeTime{"cMaxV0LifeTime", 50, "Maximum V0 life time"}; - Configurable qtArmenterosMin{"qtArmenterosMin", 0.2, "Minimum armenteros cut for K0s"}; - // config for V0 daughters - Configurable confDaughEta{"confDaughEta", 0.8f, "V0 Daugh sel: max eta"}; - Configurable cfgDaughPiPt{"cfgDaughPiPt", 0.2, "minimum daughter pion pt"}; - Configurable confDaughTPCnclsMin{"confDaughTPCnclsMin", 50.f, "V0 Daugh sel: Min. nCls TPC"}; - Configurable confDaughTPCncrwsMin{"confDaughTPCncrwsMin", 70.f, "V0 Daugh sel: Min. nCrws TPC"}; - Configurable confDaughPIDCuts{"confDaughPIDCuts", 3, "PID selections for Kshortpion daughters"}; - - // configurable for chargedkstar - Configurable cfgKstarMassMin{"cfgKstarMassMin", 0.75f, "K* mass min"}; - Configurable cfgKstarMassMax{"cfgKstarMassMax", 1.05f, "K* mass max"}; - - Configurable iMNbins{"iMNbins", 50, "Number of bins in invariant mass"}; - Configurable lbinIM{"lbinIM", 1.09, "lower bin value in IM histograms"}; - Configurable hbinIM{"hbinIM", 1.14, "higher bin value in IM histograms"}; + Configurable confV0PtMin{"confV0PtMin", 0.0f, "Minimum K0s pT"}; + Configurable confV0PtMax{"confV0PtMax", 1000.0f, "Maximum K0s pT"}; + Configurable confV0Rap{"confV0Rap", 0.8f, "Maximum K0s rapidity"}; + Configurable confV0DCADaughMax{"confV0DCADaughMax", 1.0f, "Maximum DCA between V0 daughters"}; + Configurable confV0CPAMin{"confV0CPAMin", 0.9998, "Minimum K0s CPA"}; + Configurable confV0TranRadV0Min{"confV0TranRadV0Min", 1.5f, "Minimum K0s transverse radius"}; + Configurable confV0TranRadV0Max{"confV0TranRadV0Max", 100.0f, "Maximum K0s transverse radius"}; + Configurable cMaxV0DCA{"cMaxV0DCA", 1.2, "Maximum K0s DCA to PV"}; + Configurable cMinV0DCAPi{"cMinV0DCAPi", 0.05f, "Minimum V0 daughter DCA to PV"}; + Configurable cMaxV0LifeTime{"cMaxV0LifeTime", 50.0f, "Maximum K0s lifetime"}; + Configurable qtArmenterosMin{"qtArmenterosMin", 0.2f, "Minimum Armenteros qT/|alpha| for K0s"}; + + Configurable confDaughEta{"confDaughEta", 0.8f, "Maximum V0 daughter eta"}; + Configurable cfgDaughPiPt{"cfgDaughPiPt", 0.2f, "Minimum V0 daughter pion pT"}; + Configurable confDaughTPCnclsMin{"confDaughTPCnclsMin", 50.0f, "Minimum V0 daughter TPC clusters"}; + Configurable confDaughTPCncrwsMin{"confDaughTPCncrwsMin", 70.0f, "Minimum V0 daughter TPC crossed rows"}; + Configurable confDaughPIDCuts{"confDaughPIDCuts", 3.0f, "TPC pion PID cut for V0 daughters"}; + + Configurable cfgK0sMassMin{"cfgK0sMassMin", 0.45f, "Minimum K0s invariant mass"}; + Configurable cfgK0sMassMax{"cfgK0sMassMax", 0.55f, "Maximum K0s invariant mass"}; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; RCTFlagsChecker rctChecker; + void init(o2::framework::InitContext&) { rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, rctCut.cfgEvtRCTFlagCheckerZDCCheck, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad); - AxisSpec thnAxisInvMass{iMNbins, lbinIM, hbinIM, "#it{M} (GeV/#it{c}^{2})"}; - histos.add("hCent", "hCent", kTH1F, {{8, 0, 80.0}}); - histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{5, 0, 5.0}}); - histos.add("hTrkSelInfo", "hTrkSelInfo", kTH1F, {{5, 0, 5.0}}); - histos.add("hKShortMass", "hKShortMass", kTH1F, {thnAxisInvMass}); - histos.add("hV0Info", "hV0Info", kTH1F, {{5, 0, 5.0}}); + histos.add("hCent", "hCent", kTH1F, {{16, 0.0f, 80.0f}}); + histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{6, 0.0f, 6.0f}}); + histos.add("hTrkSelInfo", "hTrkSelInfo", kTH1F, {{5, 0.0f, 5.0f}}); + histos.add("hV0Info", "hV0Info", kTH1F, {{5, 0.0f, 5.0f}}); + histos.add("hKShortMass", "hKShortMass;M_{#pi^{+}#pi^{-}} (GeV/#it{c}^{2});Counts", kTH1F, {{200, 0.4f, 0.6f}}); + histos.add("hNStoredK0s", "hNStoredK0s;N_{K^{0}_{S}};Events", kTH1F, {{100, 0.0f, 100.0f}}); + histos.add("hNStoredPions", "hNStoredPions;N_{#pi};Events", kTH1F, {{500, 0.0f, 500.0f}}); } template bool selectionTrack(const T& candidate) { - if (candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() >= grpPion.itsclusterPiMeson && candidate.tpcNClsCrossedRows() > grpPion.tpcCrossedRowsPiMeson && std::abs(candidate.dcaXY()) <= grpPion.cutDCAxyPiMeson && std::abs(candidate.dcaZ()) <= grpPion.cutDCAzPiMeson && std::abs(candidate.eta()) <= grpPion.cutEtaPiMeson && candidate.pt() >= grpPion.cutPTPiMeson) { - return true; - } - return false; + return candidate.isGlobalTrack() && + candidate.isPVContributor() && + candidate.itsNCls() >= grpPion.itsclusterPiMeson && + candidate.tpcNClsCrossedRows() >= grpPion.tpcCrossedRowsPiMeson && + std::abs(candidate.dcaXY()) <= grpPion.cutDCAxyPiMeson && + std::abs(candidate.dcaZ()) <= grpPion.cutDCAzPiMeson && + std::abs(candidate.eta()) <= grpPion.cutEtaPiMeson && + candidate.pt() >= grpPion.cutPTPiMeson; } template bool selectionPID(const T& candidate) { - if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < grpPion.nsigmaCutTPCPiMeson) { - return true; + const float nTPC = candidate.tpcNSigmaPi(); + + if (!candidate.hasTOF()) { + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson; } - if (candidate.hasTOF() && candidate.beta() > grpPion.cutTOFBetaPiMeson && std::abs(candidate.tpcNSigmaPi()) < grpPion.nsigmaCutTPCPiMeson && std::abs(candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson) { - return true; + + if (candidate.beta() <= grpPion.cutTOFBetaPiMeson) { + return false; } - return false; + + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson && + std::abs(candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson; } template bool selectionPID2(const T& candidate) { - if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < grpPion.nsigmaCutTPCPiMeson) { - return true; + const float nTPC = candidate.tpcNSigmaPi(); + + if (!candidate.hasTOF()) { + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson; } - if (candidate.hasTOF() && candidate.beta() > grpPion.cutTOFBetaPiMeson && std::sqrt(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson) { - return true; + + if (candidate.beta() <= grpPion.cutTOFBetaPiMeson) { + return false; } - return false; + + const float nTOF = candidate.tofNSigmaPi(); + const float nCombined = std::sqrt(nTPC * nTPC + nTOF * nTOF); + + return nCombined < grpPion.nsigmaCutTOFPiMeson; } template bool selectionPID3(const T& candidate) { - auto px = candidate.px(); - auto py = candidate.py(); - // auto pz = candidate.pz(); - auto pt = std::sqrt(px * px + py * py); - float lowmom = 0.5; - if (pt < lowmom) { - if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < grpPion.nsigmaCutTPCPiMeson) { - return true; - } else if (candidate.hasTOF() && std::sqrt(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson) { - return true; - } - } else if (candidate.hasTOF() && std::sqrt(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson) { - return true; + const float pt = candidate.pt(); + const float nTPC = candidate.tpcNSigmaPi(); + + if (pt < grpPion.pSwitchPID) { + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson; + } + + if (!candidate.hasTOF()) { + return false; + } + + if (candidate.beta() <= grpPion.cutTOFBetaPiMeson) { + return false; } - return false; + + const float nTOF = candidate.tofNSigmaPi(); + const float nCombined = std::sqrt(nTPC * nTPC + nTOF * nTOF); + + return nCombined < grpPion.nsigmaCutTOFPiMeson; } template bool selectionPID4(const T& candidate) { - // Use total momentum p (as you said). If you really want pT, replace with sqrt(px^2+py^2). - const float px = candidate.px(); - const float py = candidate.py(); - // const float pz = candidate.pz(); - const float pt = std::sqrt(px * px + py * py); - - constexpr float pSwitch = 0.5f; // GeV/c - + const float pt = candidate.pt(); const float nTPC = candidate.tpcNSigmaPi(); - // Low momentum: TPC-only, TOF not required and not used - if (pt < pSwitch) { - return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson; // e.g. 3 + if (pt < grpPion.pSwitchPID) { + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson; } - // High momentum: TOF hit mandatory + separate 3σ cuts if (!candidate.hasTOF()) { return false; } + if (candidate.beta() <= grpPion.cutTOFBetaPiMeson) { + return false; + } + const float nTOF = candidate.tofNSigmaPi(); - return (std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson) && - (std::abs(nTOF) < grpPion.nsigmaCutTOFPiMeson); + + return std::abs(nTPC) < grpPion.nsigmaCutTPCPiMeson && + std::abs(nTOF) < grpPion.nsigmaCutTOFPiMeson; } template uint8_t pionPidMask(const T& trk) { - uint8_t m = 0; - if (selectionPID(trk)) - m |= kPID1; - if (selectionPID2(trk)) - m |= kPID2; - if (selectionPID3(trk)) - m |= kPID3; - if (selectionPID4(trk)) - m |= kPID4; - return m; + uint8_t mask = 0; + + if (selectionPID(trk)) { + mask |= kPID1; + } + + if (selectionPID2(trk)) { + mask |= kPID2; + } + + if (selectionPID3(trk)) { + mask |= kPID3; + } + + if (selectionPID4(trk)) { + mask |= kPID4; + } + + return mask; } + struct StoredK0s { + float px; + float py; + float pz; + float mass; + float cospa; + float radius; + float dcaPositive; + float dcaNegative; + float dcaBetweenDaughters; + float lifetime; + int64_t positiveIndex; + int64_t negativeIndex; + }; + + struct StoredPion { + float px; + float py; + float pz; + int8_t charge; + int64_t index; + uint8_t pidMask; + }; + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter centralityFilter = (aod::cent::centFT0C < cfgCutCentralityMax && aod::cent::centFT0C >= cfgCutCentralityMin); + Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && aod::track::pt > cfgCutPt); + + using EventCandidates = soa::Filtered>; + using AllTrackCandidates = soa::Filtered>; + using ResoV0s = aod::V0Datas; + template bool selectionV0(Collision const& collision, V0 const& candidate) { - if (std::abs(candidate.dcav0topv()) > cMaxV0DCA) { - return false; - } const float pT = candidate.pt(); const float tranRad = candidate.v0radius(); const float dcaDaughv0 = std::abs(candidate.dcaV0daughters()); const float cpav0 = candidate.v0cosPA(); - float ctauKShort = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * (o2::constants::physics::MassK0); + const float ctauKShort = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0; - if (pT < confV0PtMin) { + if (std::abs(candidate.dcav0topv()) > cMaxV0DCA) { return false; } - if (pT > confV0PtMax) { + + if (pT < confV0PtMin || pT > confV0PtMax) { return false; } + if (dcaDaughv0 > confV0DCADaughMax) { return false; } + if (cpav0 < confV0CPAMin) { return false; } - if (tranRad < confV0TranRadV0Min) { - return false; - } - if (tranRad > confV0TranRadV0Max) { + + if (tranRad < confV0TranRadV0Min || tranRad > confV0TranRadV0Max) { return false; } + if (std::abs(ctauKShort) > cMaxV0LifeTime) { return false; } + if ((candidate.qtarm() / std::abs(candidate.alpha())) < qtArmenterosMin) { return false; } - if (std::abs(candidate.yK0Short()) > confV0Rap) { // use full rapidity 0.8 for K0s + + if (std::abs(candidate.yK0Short()) > confV0Rap) { return false; } + return true; } @@ -292,241 +340,236 @@ struct cksspinalignment { auto postrack = candidate.template posTrack_as(); auto negtrack = candidate.template negTrack_as(); - const auto ncrfc = 0.8; + constexpr float minCrossedRowsOverFindable = 0.8f; if (postrack.sign() < 0 || negtrack.sign() > 0) { return false; } - if (postrack.tpcNClsCrossedRows() < confDaughTPCncrwsMin || negtrack.tpcNClsCrossedRows() < confDaughTPCncrwsMin) { + + if (postrack.tpcNClsCrossedRows() < confDaughTPCncrwsMin || + negtrack.tpcNClsCrossedRows() < confDaughTPCncrwsMin) { return false; } - if (postrack.tpcNClsFound() < confDaughTPCnclsMin || negtrack.tpcNClsFound() < confDaughTPCnclsMin) { + + if (postrack.tpcNClsFound() < confDaughTPCnclsMin || + negtrack.tpcNClsFound() < confDaughTPCnclsMin) { return false; } - if (postrack.tpcCrossedRowsOverFindableCls() < ncrfc || negtrack.tpcCrossedRowsOverFindableCls() < ncrfc) { + + if (postrack.tpcCrossedRowsOverFindableCls() < minCrossedRowsOverFindable || + negtrack.tpcCrossedRowsOverFindableCls() < minCrossedRowsOverFindable) { return false; } - if (std::abs(postrack.tpcNSigmaPi()) > confDaughPIDCuts || std::abs(negtrack.tpcNSigmaPi()) > confDaughPIDCuts) { + + if (std::abs(postrack.tpcNSigmaPi()) > confDaughPIDCuts || + std::abs(negtrack.tpcNSigmaPi()) > confDaughPIDCuts) { return false; } - if (candidate.positivept() < cfgDaughPiPt || candidate.negativept() < cfgDaughPiPt) { + + if (candidate.positivept() < cfgDaughPiPt || + candidate.negativept() < cfgDaughPiPt) { return false; } - if (std::abs(candidate.positiveeta()) > confDaughEta || std::abs(candidate.negativeeta()) > confDaughEta) { + + if (std::abs(candidate.positiveeta()) > confDaughEta || + std::abs(candidate.negativeeta()) > confDaughEta) { return false; } - if (std::abs(candidate.dcapostopv()) < cMinV0DCAPi || std::abs(candidate.dcanegtopv()) < cMinV0DCAPi) { + + if (std::abs(candidate.dcapostopv()) < cMinV0DCAPi || + std::abs(candidate.dcanegtopv()) < cMinV0DCAPi) { return false; } return true; } - std::tuple getK0sTags(const auto& v0, const auto& collision) + template + bool selectionK0s(Collision const& collision, V0 const& v0) { - // auto postrack = v0.template posTrack_as(); - // auto negtrack = v0.template negTrack_as(); - - int kshortTag = 0; - - if (isSelectedV0Daughter(v0) && v0.mK0Short() > lbinIM && v0.mK0Short() < hbinIM) { - kshortTag = 1; + if (!isSelectedV0Daughter(v0)) { + return false; } - if (!kshortTag) { - return {0, false}; // No valid tags + if (!selectionV0(collision, v0)) { + return false; } - if (!selectionV0(collision, v0)) { - return {0, false}; // Fails selection + if (v0.mK0Short() < cfgK0sMassMin || v0.mK0Short() > cfgK0sMassMax) { + return false; } - return {kshortTag, true}; // Valid candidate + return true; } - ROOT::Math::PxPyPzMVector kshort, pion, pionbach, antiPion; - ROOT::Math::PxPyPzMVector kshortDummy, pionDummy; - - Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter centralityFilter = (nabs(aod::cent::centFT0C) < cfgCutCentralityMax && nabs(aod::cent::centFT0C) > cfgCutCentralityMin); - Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && (aod::track::pt) > cfgCutPt); - - using EventCandidates = soa::Filtered>; - using AllTrackCandidates = soa::Filtered>; - using ResoV0s = aod::V0Datas; - void processData(EventCandidates::iterator const& collision, AllTrackCandidates const& tracks, ResoV0s const& V0s) + void processData(EventCandidates::iterator const& collision, + AllTrackCandidates const& tracks, + ResoV0s const& v0s) { o2::aod::ITSResponse itsResponse; - std::vector kshortMother, pionBachelor; - std::vector v0Cospa = {}; - std::vector v0Radius = {}; - std::vector dcaPositive = {}; - std::vector dcaNegative = {}; - std::vector positiveIndex = {}; - std::vector negativeIndex = {}; - std::vector dcaBetweenDaughter = {}; - std::vector v0Lifetime = {}; - // std::vector armenteros = {}; - std::vector pionBachelorIndex = {}; - // std::vector pionBachelorSign = {}; - // std::vector pionBachelorTPC = {}; - // std::vector pionBachelorTOF = {}; - // std::vector pionBachelorTOFHit = {}; - std::vector pionBachelorPidMask = {}; - - int numbV0 = 0; - auto centrality = collision.centFT0C(); - auto vz = collision.posZ(); - int occupancy = collision.trackOccupancyInTimeRange(); - auto psiFT0C = collision.psiFT0C(); - auto psiFT0A = collision.psiFT0A(); - auto psiTPC = collision.psiTPC(); - // auto psiTPCR = collision.psiTPCR(); - // auto psiTPCL = collision.psiTPCL(); + + const float centrality = collision.centFT0C(); + const float vz = collision.posZ(); + const int occupancy = collision.trackOccupancyInTimeRange(); + + const float psiFT0C = collision.psiFT0C(); + const float psiFT0A = collision.psiFT0A(); + const float psiTPC = collision.psiTPC(); + histos.fill(HIST("hEvtSelInfo"), 0.5); - // if ((!rctCut.requireRCTFlagChecker || rctChecker(collision)) && collision.selection_bit(aod::evsel::kNoSameBunchPileup) && collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && collision.selection_bit(aod::evsel::kNoITSROFrameBorder) && collision.sel8() && collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll) && occupancy < cfgCutOccupancy) { - if ((!rctCut.requireRCTFlagChecker || rctChecker(collision)) && collision.selection_bit(aod::evsel::kNoSameBunchPileup) && collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && collision.selection_bit(aod::evsel::kNoITSROFrameBorder) && (!useNoCollInTimeRangeStandard || collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) && collision.sel8() && (!useGoodITSLayersAll || collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) && occupancy < cfgCutOccupancy) { - histos.fill(HIST("hEvtSelInfo"), 1.5); - if (collision.triggereventep()) { - histos.fill(HIST("hEvtSelInfo"), 2.5); - histos.fill(HIST("hCent"), centrality); - - for (const auto& track1 : tracks) { - histos.fill(HIST("hTrkSelInfo"), 0.5); - if (!selectionTrack(track1)) { - continue; - } - histos.fill(HIST("hTrkSelInfo"), 1.5); - - if (grpPion.itsPIDSelection && !(itsResponse.nSigmaITS(track1) > grpPion.lowITSPIDNsigma && itsResponse.nSigmaITS(track1) < grpPion.highITSPIDNsigma)) { - continue; - } - histos.fill(HIST("hTrkSelInfo"), 2.5); - /* - if (grpPion.usePID && !selectionPID2(track1)) { - continue; - }*/ - const uint8_t mask = pionPidMask(track1); - if (grpPion.usePID && mask == 0) { - continue; - } - histos.fill(HIST("hTrkSelInfo"), 3.5); - - auto track1ID = track1.globalIndex(); - auto track1sign = track1.sign(); - if (track1sign == 0) - continue; - /*auto track1nsigTPC = track1.tpcNSigmaPi(); - auto track1nsigTOF = -999.9; - auto track1TOFHit = -1; - if (track1.hasTOF()) { - track1TOFHit = 1; - track1nsigTOF = track1.tofNSigmaPi(); - histos.fill(HIST("hTrkSelInfo"), 4.5); - }*/ - pionbach = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), o2::constants::physics::MassPionCharged); - pionBachelor.push_back(pionbach); - pionBachelorIndex.push_back(track1ID); - // pionBachelorSign.push_back(track1sign); - // pionBachelorTPC.push_back(track1nsigTPC); - // pionBachelorTOF.push_back(track1nsigTOF); - // pionBachelorTOFHit.push_back(track1TOFHit); - pionBachelorPidMask.push_back(mask); - } - for (const auto& v0 : V0s) { - histos.fill(HIST("hV0Info"), 0.5); - auto [kshortTag, isValid] = getK0sTags(v0, collision); - if (kshortTag && isValid) { - histos.fill(HIST("hV0Info"), 1.5); - auto postrack1 = v0.template posTrack_as(); - auto negtrack1 = v0.template negTrack_as(); - positiveIndex.push_back(postrack1.globalIndex()); - negativeIndex.push_back(negtrack1.globalIndex()); - v0Cospa.push_back(v0.v0cosPA()); - v0Radius.push_back(v0.v0radius()); - dcaPositive.push_back(std::abs(v0.dcapostopv())); - dcaNegative.push_back(std::abs(v0.dcanegtopv())); - dcaBetweenDaughter.push_back(std::abs(v0.dcaV0daughters())); - v0Lifetime.push_back(v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * (o2::constants::physics::MassK0)); - // armenteros.push_back((v0.qtarm() / std::abs(v0.alpha()))); - - pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), o2::constants::physics::MassPionCharged); - antiPion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), o2::constants::physics::MassPionCharged); - kshort = pion + antiPion; - // chargedkstar = kshort + pionbach; - kshortMother.push_back(kshort); - // chargedkstarMother.push_back(chargedkstar); - histos.fill(HIST("hKShortMass"), kshort.M()); - } - numbV0 = numbV0 + 1; - } - if (numbV0 > 1 && v0Cospa.size() > 1 && !kshortMother.empty() && !pionBachelor.empty()) { - - std::vector useK0(kshortMother.size(), 0); - std::vector usePi(pionBachelor.size(), 0); - bool anyPair = false; - for (size_t ik0 = 0; ik0 < kshortMother.size(); ++ik0) { - const auto& k0 = kshortMother[ik0]; - const int posId = positiveIndex[ik0]; - const int negId = negativeIndex[ik0]; - - for (size_t ipi = 0; ipi < pionBachelor.size(); ++ipi) { - const int piId = pionBachelorIndex[ipi]; - - // avoid self-combination: bachelor pion can't be a V0 daughter - if (piId == posId || piId == negId) { - continue; - } - - const auto kstar = k0 + pionBachelor[ipi]; - const float m = kstar.M(); - - if (m < cfgKstarMassMin || m > cfgKstarMassMax) { - continue; - } - - useK0[ik0] = 1; - usePi[ipi] = 1; - anyPair = true; - } - } - - // only write event + tables if at least one K* candidate exists in window - if (anyPair) { - histos.fill(HIST("hEvtSelInfo"), 3.5); - // kshortpionEvent(centrality, vz, collision.index(), psiFT0C, psiFT0A, psiTPC); - kshortpionEvent(centrality, vz, psiFT0C, psiFT0A, psiTPC); - auto indexEvent = kshortpionEvent.lastIndex(); - // write only used K0s - for (size_t ik0 = 0; ik0 < kshortMother.size(); ++ik0) { - if (!useK0[ik0]) { - continue; - } - kshortDummy = kshortMother[ik0]; - kshortTrack(indexEvent, v0Cospa[ik0], v0Radius[ik0], dcaPositive[ik0], dcaNegative[ik0], dcaBetweenDaughter[ik0], v0Lifetime[ik0], kshortDummy.Px(), kshortDummy.Py(), kshortDummy.Pz(), kshortDummy.M(), positiveIndex[ik0], negativeIndex[ik0]); - } - // write only used pions - for (size_t ipi = 0; ipi < pionBachelor.size(); ++ipi) { - if (!usePi[ipi]) { - continue; - } - pionDummy = pionBachelor[ipi]; - /* - pionTrack(indexEvent, - pionDummy.Px(), pionDummy.Py(), pionDummy.Pz(), - pionBachelorTPC[ipi], - pionBachelorTOFHit[ipi], - pionBachelorTOF[ipi], - pionBachelorIndex[ipi]);*/ - pionTrack(indexEvent, pionDummy.Px(), pionDummy.Py(), pionDummy.Pz(), pionBachelorIndex[ipi], pionBachelorPidMask[ipi]); - } - } - } + + if ((rctCut.requireRCTFlagChecker && !rctChecker(collision)) || + !collision.selection_bit(aod::evsel::kNoSameBunchPileup) || + !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || + !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || + !collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || + (useNoCollInTimeRangeStandard && + !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) || + !collision.sel8() || + (useGoodITSLayersAll && + !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) || + occupancy >= cfgCutOccupancy) { + return; + } + + histos.fill(HIST("hEvtSelInfo"), 1.5); + + if (!collision.triggereventep()) { + return; + } + + histos.fill(HIST("hEvtSelInfo"), 2.5); + histos.fill(HIST("hCent"), centrality); + + std::vector selectedK0s; + std::vector selectedPions; + + for (const auto& track : tracks) { + histos.fill(HIST("hTrkSelInfo"), 0.5); + + if (!selectionTrack(track)) { + continue; + } + + histos.fill(HIST("hTrkSelInfo"), 1.5); + + const float nSigmaITS = + itsResponse.nSigmaITS(track); + + if (grpPion.itsPIDSelection && + (nSigmaITS <= grpPion.lowITSPIDNsigma || + nSigmaITS >= grpPion.highITSPIDNsigma)) { + continue; + } + + histos.fill(HIST("hTrkSelInfo"), 2.5); + + const uint8_t mask = pionPidMask(track); + + if (grpPion.usePID && mask == 0) { + continue; + } + + const int8_t charge = static_cast(track.sign()); + + if (charge == 0) { + continue; + } + + histos.fill(HIST("hTrkSelInfo"), 3.5); + + selectedPions.push_back({track.px(), + track.py(), + track.pz(), + charge, + track.globalIndex(), + mask}); + } + + for (const auto& v0 : v0s) { + histos.fill(HIST("hV0Info"), 0.5); + + if (!selectionK0s(collision, v0)) { + continue; } + + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + ROOT::Math::PxPyPzMVector pionPos(v0.pxpos(), v0.pypos(), v0.pzpos(), o2::constants::physics::MassPionCharged); + ROOT::Math::PxPyPzMVector pionNeg(v0.pxneg(), v0.pyneg(), v0.pzneg(), o2::constants::physics::MassPionCharged); + ROOT::Math::PxPyPzMVector k0s = pionPos + pionNeg; + + // const float lifetime = + // v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * + // o2::constants::physics::MassK0; + + selectedK0s.push_back({static_cast(k0s.Px()), + static_cast(k0s.Py()), + static_cast(k0s.Pz()), + static_cast(k0s.M()), + static_cast(v0.v0cosPA()), + static_cast(v0.v0radius()), + static_cast(std::abs(v0.dcapostopv())), + static_cast(std::abs(v0.dcanegtopv())), + static_cast(std::abs(v0.dcaV0daughters())), + // lifetime, + postrack.globalIndex(), + negtrack.globalIndex()}); + + histos.fill(HIST("hV0Info"), 1.5); + histos.fill(HIST("hKShortMass"), k0s.M()); + } + + histos.fill(HIST("hNStoredK0s"), selectedK0s.size()); + histos.fill(HIST("hNStoredPions"), selectedPions.size()); + + if (selectedK0s.empty() && selectedPions.empty()) { + return; + } + + histos.fill(HIST("hEvtSelInfo"), 3.5); + + kshortpionEvent(centrality, + vz, + psiFT0C, + psiFT0A, + psiTPC); + + const int64_t indexEvent = kshortpionEvent.lastIndex(); + + for (const auto& k0s : selectedK0s) { + kshortTrack(indexEvent, + k0s.cospa, + k0s.radius, + k0s.dcaPositive, + k0s.dcaNegative, + k0s.dcaBetweenDaughters, + // k0s.lifetime, + k0s.px, + k0s.py, + k0s.pz, + k0s.mass, + k0s.positiveIndex, + k0s.negativeIndex); + } + + for (const auto& pion : selectedPions) { + pionTrack(indexEvent, + pion.px, + pion.py, + pion.pz, + pion.charge, + pion.index, + pion.pidMask); } + + histos.fill(HIST("hEvtSelInfo"), 4.5); } + PROCESS_SWITCH(cksspinalignment, processData, "Process data", true); }; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; diff --git a/PWGLF/Tasks/Resonances/cksspinalignder.cxx b/PWGLF/Tasks/Resonances/cksspinalignder.cxx index 830d4764f71..54d2958bbd9 100644 --- a/PWGLF/Tasks/Resonances/cksspinalignder.cxx +++ b/PWGLF/Tasks/Resonances/cksspinalignder.cxx @@ -10,16 +10,12 @@ // or submit itself to any jurisdiction. /// \file cksspinalignder.cxx -/// \brief Analysis task for ChargedKStar spin alignment +/// \brief Derived task for charged KStar spin alignment using reduced K0s and pion tables /// /// \author prottay.das@cern.ch #include "PWGLF/DataModel/LFCKSSpinalignmentTables.h" -#include "Common/Core/RecoDecay.h" - -#include -#include #include #include #include @@ -35,22 +31,13 @@ #include #include -#include -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) -#include -#include +#include +#include #include -#include -#include // for std::fabs -#include +#include #include -#include #include -#include -#include -#include -#include using namespace o2; using namespace o2::framework; @@ -59,108 +46,65 @@ using namespace o2::soa; struct cksspinalignder { - struct : ConfigurableGroup { - Configurable cfgURL{"cfgURL", "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"}; - Configurable ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "Latest acceptable timestamp of creation for the object"}; - } cfgCcdbParam; - - // event sel///////// - Configurable centMin{"centMin", 0, "Minimum Centrality"}; - Configurable centMax{"centMax", 80, "Maximum Centrality"}; - // V0 selection //////////// - Configurable cosPA{"cosPA", 0.995, "Cosine Pointing Angle"}; - Configurable radiusMin{"radiusMin", 1.2, "Minimum V0 radius"}; - Configurable radiusMax{"radiusMax", 100, "Maximum V0 radius"}; - Configurable dcaPion{"dcaPion", 0.1, "DCA Pion"}; - Configurable dcaDaughters{"dcaDaughters", 1.0, "DCA between daughters"}; - Configurable lifetimeMax{"lifetimeMax", 20, "Maximum V0 lifetime"}; - Configurable ptMin{"ptMin", 0.5, "V0 Pt minimum"}; - Configurable ptMax{"ptMax", 10.0, "V0 Pt maximum"}; - Configurable v0eta{"v0eta", 0.8, "Eta cut on lambda"}; - // pion sel///////// - Configurable usePID{"usePID", false, "Enable PID selection using stored bitmask"}; - // Choose which PID bit to require (1..4). If pidExclusive=true -> require ONLY that PID. - Configurable pidChoice{"pidChoice", 2, "PID choice: 1,2,3,4"}; - Configurable pidExclusive{"pidExclusive", false, "If true require ONLY the chosen PID (mask == bit). If false require (mask & bit)!=0"}; - - // Event Mixing - Configurable nEvtMixing{"nEvtMixing", 5, "Number of events to mix"}; - ConfigurableAxis cfgVtxBins{"cfgVtxBins", {10, -10, 10}, "Mixing bins - z-vertex"}; - ConfigurableAxis cfgMultBins{"cfgMultBins", {8, 0.0, 80}, "Mixing bins - centrality"}; - Configurable ptMix{"ptMix", 0.2f, "ME: pT bin width"}; - Configurable etaMix{"etaMix", 0.2f, "ME: eta bin width"}; - Configurable phiMix{"phiMix", 0.3f, "ME: phi bin width"}; - - // THnsparse bining - ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {50, 1.09, 1.14}, "#it{M} (GeV/#it{c}^{2})"}; - ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.0, 10.0}, "#it{p}_{T}"}; - ConfigurableAxis configThnAxisSA{"configThnAxisSA", {20, -1.0, 1.0}, "cos#it{#theta *}"}; - ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {8, 0.0, 80.0}, "Centrality"}; - - // Enable access to the CCDB for the offset and correction constants and save them in dedicated variables. - Service ccdb; - o2::ccdb::CcdbApi ccdbApi; + // Event selection + Configurable centMin{"centMin", 0.0f, "Minimum centrality"}; + Configurable centMax{"centMax", 80.0f, "Maximum centrality"}; + + // K0s selection from stored reduced table + Configurable applyK0sSelection{"applyK0sSelection", true, "Apply additional K0s selection in derived task"}; + Configurable cosPA{"cosPA", 0.995f, "Minimum V0 cosine pointing angle"}; + Configurable radiusMin{"radiusMin", 1.2f, "Minimum V0 radius"}; + Configurable radiusMax{"radiusMax", 100.0f, "Maximum V0 radius"}; + Configurable dcaPion{"dcaPion", 0.1f, "Minimum V0 daughter DCA to PV"}; + Configurable dcaDaughters{"dcaDaughters", 1.0f, "Maximum DCA between V0 daughters"}; + Configurable k0sPtMin{"k0sPtMin", 0.0f, "Minimum K0s pT"}; + Configurable k0sPtMax{"k0sPtMax", 10.0f, "Maximum K0s pT"}; + Configurable k0sEtaMax{"k0sEtaMax", 0.8f, "Maximum K0s eta"}; + Configurable k0sMassMin{"k0sMassMin", 0.45f, "Minimum K0s mass"}; + Configurable k0sMassMax{"k0sMassMax", 0.55f, "Maximum K0s mass"}; + + // Bachelor pion PID + Configurable usePID{"usePID", true, "Apply bachelor pion PID using stored bitmask"}; + Configurable pidChoice{"pidChoice", 4, "PID choice: 1, 2, 3, or 4"}; + Configurable pidExclusive{"pidExclusive", false, "If true require only selected PID bit. If false require selected bit among possible bits"}; + + // Charge selection + Configurable chargeSelection{"chargeSelection", 0, "Bachelor pion charge selection: 0 = both, +1 = K*+, -1 = K*-"}; + + // Event-plane choice + Configurable epChoice{"epChoice", 0, "Event plane: 0 = FT0C, 1 = FT0A, 2 = TPC"}; + + // Event mixing + Configurable nEvtMixing{"nEvtMixing", 5, "Number of mixed events"}; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {10, -10.0f, 10.0f}, "Mixing bins in z vertex"}; + ConfigurableAxis cfgCentBins{"cfgCentBins", {8, 0.0f, 80.0f}, "Mixing bins in centrality"}; + + // Sparse axes + ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {150, 0.75f, 1.05f}, "#it{M}_{K^{0}_{S}#pi} (GeV/#it{c}^{2})"}; + ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.0f, 10.0f}, "#it{p}_{T}^{K^{0}_{S}#pi} (GeV/#it{c})"}; + ConfigurableAxis configThnAxisSA{"configThnAxisSA", {20, -1.0f, 1.0f}, "cos#it{#theta}^{*}"}; + ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {8, 0.0f, 80.0f}, "Centrality (%)"}; + ConfigurableAxis configThnAxisCharge{"configThnAxisCharge", {2, -1.5f, 1.5f}, "Bachelor pion charge"}; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(o2::framework::InitContext&) { - - histos.add("hCentrality", "Centrality distribution", kTH1F, {{80, 0, 80}}); - histos.add("hKShortMass", "hKShortMass", kTH1F, {{100, 0.45, 0.55}}); - - histos.add("hSparsesame", "hSparsesame", kTHnSparseF, {configThnAxisInvMass, configThnAxisPt, configThnAxisSA, configThnAxisCentrality}); - histos.add("hSparsemix", "hSparsemix", kTHnSparseF, {configThnAxisInvMass, configThnAxisPt, configThnAxisSA, configThnAxisCentrality}); - - ccdb->setURL(cfgCcdbParam.cfgURL); - ccdbApi.init("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + histos.add("hCentrality", "Centrality distribution;Centrality (%);Events", kTH1F, {{80, 0.0f, 80.0f}}); + histos.add("hKShortMass", "K^{0}_{S} mass;M_{#pi^{+}#pi^{-}} (GeV/#it{c}^{2});Counts", kTH1F, {{200, 0.4f, 0.6f}}); + + histos.add("hSparsesame", + "Same-event K^{0}_{S}#pi;M_{K^{0}_{S}#pi};p_{T};cos#theta^{*};centrality;charge", + kTHnSparseF, + {configThnAxisInvMass, configThnAxisPt, configThnAxisSA, configThnAxisCentrality, configThnAxisCharge}); + + histos.add("hSparsemix", + "Mixed-event K^{0}_{S}#pi;M_{K^{0}_{S}#pi};p_{T};cos#theta^{*};centrality;charge", + kTHnSparseF, + {configThnAxisInvMass, configThnAxisPt, configThnAxisSA, configThnAxisCentrality, configThnAxisCharge}); } - template - bool selectionV0(T const& candidate) - { - auto kshortPt = std::sqrt(candidate.kShortPx() * candidate.kShortPx() + candidate.kShortPy() * candidate.kShortPy()); - auto px = candidate.kShortPx(); - auto py = candidate.kShortPy(); - auto pz = candidate.kShortPz(); - auto p = std::sqrt(px * px + py * py + pz * pz); - auto kshortEta = 0.5 * std::log((p + pz) / (p - pz)); - - if (std::abs(kshortEta) > v0eta) { - return false; - } - if (candidate.v0Cospa() < cosPA) { - return false; - } - if (candidate.v0Radius() > radiusMax) { - return false; - } - if (candidate.v0Radius() < radiusMin) { - return false; - } - if (candidate.dcaBetweenDaughter() > dcaDaughters) { - return false; - } - if (std::abs(candidate.dcaPositive()) < dcaPion && std::abs(candidate.dcaNegative()) < dcaPion) { - return false; - } - if (candidate.v0Lifetime() > lifetimeMax) { - return false; - } - if (kshortPt < ptMin) { - return false; - } - if (kshortPt > ptMax) { - return false; - } - return true; - } - - // ----- PID via stored bitmask (no selectionPID function needed) ----- - // Assumes your PionTracks table has a uint8_t column named pionPidMask(). - // Bit mapping: PID1=1, PID2=2, PID3=4, PID4=8. + // PID bit mapping from producer static constexpr uint8_t kPID1 = 1u << 0; static constexpr uint8_t kPID2 = 1u << 1; static constexpr uint8_t kPID3 = 1u << 2; @@ -168,8 +112,7 @@ struct cksspinalignder { uint8_t requiredPidBit() const { - const int choice = pidChoice.value; - switch (choice) { + switch (pidChoice.value) { case 1: return kPID1; case 2: @@ -179,8 +122,8 @@ struct cksspinalignder { case 4: return kPID4; default: - LOGF(warn, "pidChoice=%d is invalid. Using PID2.", choice); - return kPID2; + LOGF(warn, "pidChoice=%d is invalid. Using PID4.", pidChoice.value); + return kPID4; } } @@ -190,640 +133,300 @@ struct cksspinalignder { if (!usePID.value) { return true; } - const uint8_t bit = requiredPidBit(); + const uint8_t mask = pion.pionPidMask(); + const uint8_t bit = requiredPidBit(); if (pidExclusive.value) { - return mask == bit; // ONLY chosen PID + return mask == bit; } - return (mask & bit) != 0; // chosen PID (maybe others too) - } - // --------------------------------------------------------------- - std::tuple computePtEtaPhi(float px, float py, float pz) - { - float pt = std::sqrt(px * px + py * py); - float p = std::sqrt(px * px + py * py + pz * pz); - float eta = (p != std::abs(pz)) ? 0.5 * std::log((p + pz) / (p - pz)) : 0.0f; // avoid division by zero - float phi = RecoDecay::constrainAngle(std::atan2(py, px)); - return {pt, eta, phi}; + return (mask & bit) != 0; } - ROOT::Math::PtEtaPhiMVector kshort, pion, chkstar; - ROOT::Math::PtEtaPhiMVector kshortmix, pionmix, chkstarmix; - ROOT::Math::PxPyPzMVector fourVecDauCM, fourVecDauCMmix; - ROOT::Math::XYZVector threeVecDauCM, eventplaneVecNorm, threeVecDauCMmix, eventplaneVecNormmix; - - Filter centralityFilter = (nabs(aod::kshortpionevent::cent) < centMax && nabs(aod::kshortpionevent::cent) > centMin); - - using EventCandidates = soa::Filtered; - - void processSameData(EventCandidates::iterator const& collision, aod::KShortTracks const& V0s, aod::PionTracks const& piontracks) + template + bool passChargeSelection(const TPion& pion) const { - auto centrality = collision.cent(); - histos.fill(HIST("hCentrality"), centrality); - - auto psiFT0C = collision.psiFT0C(); + const int charge = static_cast(pion.charge()); - for (const auto& v0 : V0s) { - if (!selectionV0(v0)) { - continue; - } - auto [kshortPt, kshortEta, kshortPhi] = computePtEtaPhi(v0.kShortPx(), v0.kShortPy(), v0.kShortPz()); - kshort = ROOT::Math::PtEtaPhiMVector(kshortPt, kshortEta, kshortPhi, v0.kShortMass()); - histos.fill(HIST("hKShortMass"), kshort.M()); - - for (const auto& piontrack : piontracks) { - - // PID selection via stored mask - if (!passSelectedPID(piontrack)) { - continue; - } - auto [pionPt, pionEta, pionPhi] = computePtEtaPhi(piontrack.pionBachPx(), piontrack.pionBachPy(), piontrack.pionBachPz()); - pion = ROOT::Math::PtEtaPhiMVector(pionPt, pionEta, pionPhi, o2::constants::physics::MassPionCharged); - - if (piontrack.pionBachIndex() == v0.pionIndex1() || piontrack.pionBachIndex() == v0.pionIndex2()) - continue; // checking if bachelor pion is khort daughter or not -> skip further processing if such is the case + if (charge == 0) { + return false; + } - chkstar = kshort + pion; + if (chargeSelection.value == 0) { + return true; + } - ROOT::Math::Boost boost{chkstar.BoostToCM()}; - fourVecDauCM = boost(kshort); - threeVecDauCM = fourVecDauCM.Vect(); - eventplaneVecNorm = ROOT::Math::XYZVector(std::sin(2.0 * psiFT0C), -std::cos(2.0 * psiFT0C), 0); - auto cosThetaStar = eventplaneVecNorm.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(eventplaneVecNorm.Mag2()); + return charge == chargeSelection.value; + } - histos.fill(HIST("hSparsesame"), chkstar.M(), chkstar.Pt(), cosThetaStar, centrality); - } + template + bool selectionV0(const TK0s& candidate) const + { + if (!applyK0sSelection.value) { + return true; } - } - PROCESS_SWITCH(cksspinalignder, processSameData, "Process same data", true); - // Processing Event Mixing - SliceCache cache; - using BinningType = ColumnBinningPolicy; - BinningType colBinning{{cfgVtxBins, cfgMultBins}, true}; - Preslice tracksPerCollisionV0 = aod::kshortpionpair::kshortpioneventId; - Preslice tracksPerCollisionBach = aod::kshortpionpair::kshortpioneventId; + ROOT::Math::PxPyPzMVector k0s(candidate.kShortPx(), + candidate.kShortPy(), + candidate.kShortPz(), + candidate.kShortMass()); - void processMixedData(EventCandidates const& collisions, aod::KShortTracks const& V0s, aod::PionTracks const& piontracks) - { - for (const auto& [collision1, collision2] : selfCombinations(colBinning, nEvtMixing, -1, collisions, collisions)) { - if (collision1.index() == collision2.index()) { - continue; - } - auto centrality = collision1.cent(); - auto psiFT0Cmix = collision1.psiFT0C(); + const float k0sPt = k0s.Pt(); + const float k0sEta = k0s.Eta(); - auto groupV0 = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); - auto groupPion = piontracks.sliceBy(tracksPerCollisionBach, collision2.index()); - for (const auto& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(groupV0, groupPion))) { - if (!selectionV0(t1)) - continue; - auto [kshortPtmix, kshortEtamix, kshortPhimix] = computePtEtaPhi(t1.kShortPx(), t1.kShortPy(), t1.kShortPz()); - kshortmix = ROOT::Math::PtEtaPhiMVector(kshortPtmix, kshortEtamix, kshortPhimix, t1.kShortMass()); + if (std::abs(k0sEta) > k0sEtaMax.value) { + return false; + } - // PID selection via stored mask (on the bachelor pion from the mixed event) - if (!passSelectedPID(t2)) { - continue; - } - auto [pionPtmix, pionEtamix, pionPhimix] = computePtEtaPhi(t2.pionBachPx(), t2.pionBachPy(), t2.pionBachPz()); - pionmix = ROOT::Math::PtEtaPhiMVector(pionPtmix, pionEtamix, pionPhimix, o2::constants::physics::MassPionCharged); + if (candidate.kShortMass() < k0sMassMin.value || + candidate.kShortMass() > k0sMassMax.value) { + return false; + } - chkstarmix = kshortmix + pionmix; + if (candidate.v0Cospa() < cosPA.value) { + return false; + } - ROOT::Math::Boost boost{chkstarmix.BoostToCM()}; - fourVecDauCMmix = boost(kshortmix); - threeVecDauCMmix = fourVecDauCMmix.Vect(); - eventplaneVecNormmix = ROOT::Math::XYZVector(std::sin(2.0 * psiFT0Cmix), -std::cos(2.0 * psiFT0Cmix), 0); - auto cosThetaStarmix = eventplaneVecNormmix.Dot(threeVecDauCMmix) / std::sqrt(threeVecDauCMmix.Mag2()) / std::sqrt(eventplaneVecNormmix.Mag2()); + if (candidate.v0Radius() < radiusMin.value || + candidate.v0Radius() > radiusMax.value) { + return false; + } - histos.fill(HIST("hSparsemix"), chkstarmix.M(), chkstarmix.Pt(), cosThetaStarmix, centrality); - } + if (candidate.dcaBetweenDaughter() > dcaDaughters.value) { + return false; } - } - PROCESS_SWITCH(cksspinalignder, processMixedData, "Process mixed data", true); - - // ===================================================================================== - // MEV2-style mixing for Charged K* : KEEP K0s (from event E) FIXED, MIX bachelor pion - // pion is taken from other events E' but matched in (event class + pion pt/eta/phi) - // ===================================================================================== - - // One status bin (unused dimension, kept for compatibility with linearKey) - static constexpr int N_STATUS = 1; - - // ---------- Binner for bachelor-pion kinematics ---------- - struct MixBinnerPi { - float ptMin, ptMax, ptStep; - float etaMin, etaMax, etaStep; - float phiMin, phiMax, phiStep; - - // Dummy "mass" axis (we keep 1 bin to keep the same 6D infrastructure) - static constexpr float mMin = 0.f; - static constexpr float mMax = 1.f; - static constexpr int nM_ = 1; - static constexpr float mStep = (mMax - mMin) / nM_; - - int nPt_{1}, nEta_{1}, nPhi_{1}; - - MixBinnerPi(float ptMin_, float ptMax_, float ptStep_, - float etaAbsMax_, float etaStep_, - float phiStep_) - : ptMin(ptMin_), ptMax(ptMax_), ptStep(ptStep_), etaMin(-etaAbsMax_), etaMax(+etaAbsMax_), etaStep(etaStep_), phiMin(0.f), phiMax(static_cast(2.0 * TMath::Pi())), phiStep(phiStep_) - { - ptStep = (ptStep > 0.f ? ptStep : 0.2f); - etaStep = (etaStep > 0.f ? etaStep : 0.2f); - phiStep = (phiStep > 0.f ? phiStep : 0.2f); - - nPt_ = std::max(1, static_cast(std::floor((ptMax - ptMin) / ptStep + 0.5f))); - nEta_ = std::max(1, static_cast(std::floor((etaMax - etaMin) / etaStep + 0.5f))); - nPhi_ = std::max(1, static_cast(std::ceil((phiMax - phiMin) / phiStep))); + + if (candidate.dcaPositive() < dcaPion.value || + candidate.dcaNegative() < dcaPion.value) { + return false; } - inline int nPt() const { return nPt_; } - inline int nEta() const { return nEta_; } - inline int nPhi() const { return nPhi_; } - inline int nM() const { return nM_; } - - inline int binFromValue(float v, float vmin, float step, int nBins) const - { - if (!std::isfinite(v)) - return -1; - const float x = (v - vmin) / step; - int b = static_cast(std::floor(x + 1e-6f)); - if (b < 0) - return -1; - if (b >= nBins) - b = nBins - 1; - return b; + if (k0sPt < k0sPtMin.value || + k0sPt > k0sPtMax.value) { + return false; } - inline int ptBin(float pt) const { return binFromValue(pt, ptMin, ptStep, nPt_); } - inline int etaBin(float eta) const { return binFromValue(eta, etaMin, etaStep, nEta_); } - inline int phiBin(float phi) const { return binFromValue(phi, phiMin, phiStep, nPhi_); } - inline int massBin(float m) const { return binFromValue(m, mMin, mStep, nM_); } // always 0 if finite - }; - - // Buffer candidate: one pion stored in a bin - struct BufferCandPi { - int64_t collisionIdx; // to enforce different event - int64_t rowIndex; // row id in PionTracks (for iteratorAt) - uint16_t ptBin, etaBin, phiBin, mBin; - }; - - struct MatchRef { - int64_t collisionIdx; - int64_t rowIndex; - }; - - // Key: (colBin, stat, pt, eta, phi, m) - static inline size_t linearKey(int colBin, int statBin, - int ptBin, int etaBin, int phiBin, int mBin, - int nStat, int nPt, int nEta, int nPhi, int nM) - { - return ((((((static_cast(colBin) * nStat + statBin) * nPt + ptBin) * nEta + etaBin) * nPhi + phiBin) * nM + mBin)); + return true; } - // ===================================================================================== - // Mixed-event method: - // PASS 1: build buffer of pions by (event-class bin) + (pi pt/eta/phi bins) - // PASS 2: loop same-event (K0, piSame) pairs; replace piSame by piX from buffer - // keep K0 from current event fixed; use psiFT0C of current event. - // ===================================================================================== - void processMixedDataMEV2_Pion(EventCandidates const& collisions, - aod::KShortTracks const& V0s, - aod::PionTracks const& piontracks) + template + bool k0sBelongsToEvent(const TK0s& k0s, const TCollision& collision) const { - static thread_local std::mt19937 rng(0xBADC0FFE); - - // Build pion binner from your configurables - // Use your own pt range config for bachelor pions if you have it; otherwise reuse ptMin/ptMax - MixBinnerPi mb{ - ptMin, ptMax, /*ptStep*/ ptMix.value, - /*|eta|max*/ v0eta.value, /*etaStep*/ etaMix.value, - /*phiStep*/ phiMix.value}; - - const int nCol = colBinning.getAllBinsCount(); // event-class bins (vz, cent) - const int nStat = N_STATUS; // 1 - const int nPt = mb.nPt(); - const int nEta = mb.nEta(); - const int nPhi = mb.nPhi(); - const int nM = mb.nM(); - - const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; - std::vector> buffer(nKeys); - - // ---------------- PASS 1: fill 6D buffer with bachelor pions ---------------- - for (auto const& col : collisions) { - const int colBin = colBinning.getBin(std::make_tuple(col.posz(), col.cent())); - if (colBin < 0) - continue; - - auto slicePi = piontracks.sliceBy(tracksPerCollisionBach, col.index()); - - for (auto const& pi : slicePi) { + return k0s.kshortpioneventId() == collision.globalIndex(); + } - // apply configurable PID selection stored in bitmask (your passSelectedPID) - if (!passSelectedPID(pi)) - continue; + template + bool pionBelongsToEvent(const TPion& pion, const TCollision& collision) const + { + return pion.kshortpioneventId() == collision.globalIndex(); + } - auto [pt, eta, phi] = computePtEtaPhi(pi.pionBachPx(), pi.pionBachPy(), pi.pionBachPz()); - phi = RecoDecay::constrainAngle(phi, 0.0F); + template + float getEventPlane(const TCollision& collision) const + { + switch (epChoice.value) { + case 0: + return collision.psiFT0C(); + case 1: + return collision.psiFT0A(); + case 2: + return collision.psiTPC(); + default: + LOGF(warn, "epChoice=%d is invalid. Using FT0C.", epChoice.value); + return collision.psiFT0C(); + } + } - const int ptB = mb.ptBin(pt); - const int etaB = mb.etaBin(eta); - const int phiB = mb.phiBin(phi); - const int mB = 0; // dummy axis has 1 bin - if (ptB < 0 || etaB < 0 || phiB < 0) - continue; + bool computeCosThetaStar(const ROOT::Math::PxPyPzMVector& k0s, + const ROOT::Math::PxPyPzMVector& kstar, + float psi, + float& cosThetaStar) const + { + ROOT::Math::Boost boost{kstar.BoostToCM()}; + ROOT::Math::PxPyPzMVector k0sCM = boost(k0s); + ROOT::Math::XYZVector k0sCMVec = k0sCM.Vect(); - const int stat = 0; - const size_t key = linearKey(colBin, stat, ptB, etaB, phiB, mB, nStat, nPt, nEta, nPhi, nM); + ROOT::Math::XYZVector eventPlaneVector(std::sin(2.0f * psi), + -std::cos(2.0f * psi), + 0.0f); - const int64_t row = static_cast(pi.globalIndex()); + const double denom = + std::sqrt(k0sCMVec.Mag2()) * + std::sqrt(eventPlaneVector.Mag2()); - buffer[key].push_back(BufferCandPi{ - .collisionIdx = static_cast(col.index()), - .rowIndex = row, - .ptBin = static_cast(ptB), - .etaBin = static_cast(etaB), - .phiBin = static_cast(phiB), - .mBin = static_cast(mB)}); - } + if (denom <= 0.0 || !std::isfinite(denom)) { + return false; } - // ---------------- PASS 2: same-event (K0,piSame), replace ONLY pion by piX ---------------- - for (auto const& col1 : collisions) { - const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); - if (colBin < 0) - continue; + cosThetaStar = + static_cast(eventPlaneVector.Dot(k0sCMVec) / denom); - auto poolK0 = V0s.sliceBy(tracksPerCollisionV0, col1.index()); - auto poolPi = piontracks.sliceBy(tracksPerCollisionBach, col1.index()); - if (poolK0.size() == 0 || poolPi.size() == 0) - continue; + cosThetaStar = std::clamp(cosThetaStar, -1.0f, 1.0f); - const float centrality = col1.cent(); - const float psiFT0Cmix = col1.psiFT0C(); + return true; + } - // Loop over K0s from the *current* event (fixed object) - for (auto const& k0 : poolK0) { - if (!selectionV0(k0)) - continue; + template + bool fillKstarPair(const TK0s& k0sCand, + const TPion& pionCand, + const TCollision& collision, + bool isMixed) + { + if (!selectionV0(k0sCand)) { + return false; + } - // Build K0 4-vector once - auto [kpt, keta, kphi] = computePtEtaPhi(k0.kShortPx(), k0.kShortPy(), k0.kShortPz()); - kphi = RecoDecay::constrainAngle(kphi, 0.0F); - kshortmix = ROOT::Math::PtEtaPhiMVector(kpt, keta, kphi, k0.kShortMass()); + if (!passSelectedPID(pionCand)) { + return false; + } - // Loop over same-event pions to define the "same-event pair" base - for (auto const& piSame : poolPi) { - if (!passSelectedPID(piSame)) - continue; + if (!passChargeSelection(pionCand)) { + return false; + } - // avoid self-combination: bachelor pion can't be a K0 daughter - if (piSame.pionBachIndex() == k0.pionIndex1() || piSame.pionBachIndex() == k0.pionIndex2()) { - continue; - } + // Same-event self-combination rejection: + // bachelor pion must not be one of the K0s daughters. + if (!isMixed) { + if (pionCand.pionBachIndex() == k0sCand.pionIndex1() || + pionCand.pionBachIndex() == k0sCand.pionIndex2()) { + return false; + } + } - // Compute the kinematic bin of THIS pion (this defines where to pick mixed pions) - auto [ppt, peta, pphi] = computePtEtaPhi(piSame.pionBachPx(), piSame.pionBachPy(), piSame.pionBachPz()); - pphi = RecoDecay::constrainAngle(pphi, 0.0F); + ROOT::Math::PxPyPzMVector k0s(k0sCand.kShortPx(), + k0sCand.kShortPy(), + k0sCand.kShortPz(), + k0sCand.kShortMass()); - const int ptB = mb.ptBin(ppt); - const int etaB = mb.etaBin(peta); - const int phiB = mb.phiBin(pphi); - const int mB = 0; - if (ptB < 0 || etaB < 0 || phiB < 0) - continue; + ROOT::Math::PxPyPzMVector pion(pionCand.pionBachPx(), + pionCand.pionBachPy(), + pionCand.pionBachPz(), + o2::constants::physics::MassPionCharged); - const int stat = 0; - const size_t key = linearKey(colBin, stat, ptB, etaB, phiB, mB, nStat, nPt, nEta, nPhi, nM); - auto const& binVec = buffer[key]; - if (binVec.empty()) - continue; + const auto kstar = k0s + pion; - // collect partners from same bin but different collision - std::vector matches; - matches.reserve(binVec.size()); - const int64_t curColIdx = static_cast(col1.index()); + const float mass = kstar.M(); + const float pt = kstar.Pt(); + const float centrality = collision.cent(); + const float charge = static_cast(pionCand.charge()); + const float psi = getEventPlane(collision); - for (auto const& bc : binVec) { - if (bc.collisionIdx == curColIdx) - continue; - matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); - } - if (matches.empty()) - continue; + float cosThetaStar = 0.0f; - // random unique sampling (cap = nEvtMixing) - const int cap = nEvtMixing.value; - if (cap > 0 && cap < static_cast(matches.size())) { - std::uniform_int_distribution dist(0, static_cast(matches.size()) - 1); - std::unordered_set chosen; - chosen.reserve(static_cast(cap) * 2); - while (static_cast(chosen.size()) < cap) - chosen.insert(dist(rng)); - - std::vector subset; - subset.reserve(cap); - for (const int& idx : chosen) - subset.push_back(matches[idx]); - matches.swap(subset); - } else { - std::shuffle(matches.begin(), matches.end(), rng); - } + if (!computeCosThetaStar(k0s, kstar, psi, cosThetaStar)) { + return false; + } - // const float wBase = 1.0f / static_cast(matches.size()); - const float wBase = 1.0; + if (isMixed) { + histos.fill(HIST("hSparsemix"), mass, pt, cosThetaStar, centrality, charge); + } else { + histos.fill(HIST("hSparsesame"), mass, pt, cosThetaStar, centrality, charge); + } - // Replace pion by mixed pion piX, keep k0 fixed - for (auto const& m : matches) { - auto piX = piontracks.iteratorAt(m.rowIndex); // must match PASS 1 rowIndex convention + return true; + } - // safety: PID selection again - if (!passSelectedPID(piX)) - continue; + Filter centralityFilter = + (aod::kshortpionevent::cent >= centMin) && + (aod::kshortpionevent::cent < centMax); - // OPTIONAL keep: also avoid piX being a K0 daughter (same k0 from current event) - if (piX.pionBachIndex() == k0.pionIndex1() || piX.pionBachIndex() == k0.pionIndex2()) { - continue; - } + using EventCandidates = soa::Filtered; - auto [xpt, xeta, xphi] = computePtEtaPhi(piX.pionBachPx(), piX.pionBachPy(), piX.pionBachPz()); - xphi = RecoDecay::constrainAngle(xphi, 0.0F); - pionmix = ROOT::Math::PtEtaPhiMVector(xpt, xeta, xphi, o2::constants::physics::MassPionCharged); + void processSameData(EventCandidates::iterator const& collision, + aod::KShortTracks const& k0sTracks, + aod::PionTracks const& pionTracks) + { + const float centrality = collision.cent(); - chkstarmix = kshortmix + pionmix; + histos.fill(HIST("hCentrality"), centrality); - ROOT::Math::Boost boost{chkstarmix.BoostToCM()}; - fourVecDauCMmix = boost(kshortmix); - threeVecDauCMmix = fourVecDauCMmix.Vect(); + for (const auto& k0sCand : k0sTracks) { + if (!k0sBelongsToEvent(k0sCand, collision)) { + continue; + } - eventplaneVecNormmix = ROOT::Math::XYZVector(std::sin(2.f * psiFT0Cmix), - -std::cos(2.f * psiFT0Cmix), 0.f); + if (!selectionV0(k0sCand)) { + continue; + } - const float cosThetaStarmix = - eventplaneVecNormmix.Dot(threeVecDauCMmix) / - std::sqrt(threeVecDauCMmix.Mag2()) / - std::sqrt(eventplaneVecNormmix.Mag2()); + histos.fill(HIST("hKShortMass"), k0sCand.kShortMass()); - histos.fill(HIST("hSparsemix"), - chkstarmix.M(), chkstarmix.Pt(), cosThetaStarmix, centrality, wBase); - } + for (const auto& pionCand : pionTracks) { + if (!pionBelongsToEvent(pionCand, collision)) { + continue; } + + fillKstarPair(k0sCand, pionCand, collision, false); } } } - PROCESS_SWITCH(cksspinalignder, processMixedDataMEV2_Pion, "Process mixed data MEV4 (keep K0 fixed, replace pion with mixed pion)", true); - - /* -// ---- Configuration knobs (already in your task or add similarly) ---- -// ptMix: step for K0 pT binning, etaMix step for eta, phiMix step for phi -// configThnAxisInvMass etc are your histogram axes (not used for binning here) - -// We can keep one status bin for now -static constexpr int N_STATUS = 1; - -// Binner for K0 kinematics -struct MixBinnerK0 { - float ptMin, ptMax, ptStep; - float etaMin, etaMax, etaStep; - float phiMin, phiMax, phiStep; - - // K0 mass binning (you can widen if you want) - static constexpr float mMin = 0.45f; - static constexpr float mMax = 0.55f; // exclusive-ish - static constexpr int nM_ = 1; - static constexpr float mStep = (mMax - mMin) / nM_; - - int nPt_, nEta_, nPhi_; + PROCESS_SWITCH(cksspinalignder, processSameData, "Process same-event charged KStar", true); - MixBinnerK0(float ptMin_, float ptMax_, float ptStep_, - float etaAbsMax_, float etaStep_, - float phiStep_) - : ptMin(ptMin_), ptMax(ptMax_), ptStep(ptStep_), - etaMin(-etaAbsMax_), etaMax(+etaAbsMax_), etaStep(etaStep_), - phiMin(0.f), phiMax(static_cast(2.0 * TMath::Pi())), phiStep(phiStep_) - { - ptStep = (ptStep > 0.f ? ptStep : 0.1f); - etaStep = (etaStep > 0.f ? etaStep : 0.1f); - phiStep = (phiStep > 0.f ? phiStep : 0.1f); - - nPt_ = std::max(1, static_cast(std::floor((ptMax - ptMin) / ptStep + 0.5f))); - nEta_ = std::max(1, static_cast(std::floor((etaMax - etaMin) / etaStep + 0.5f))); - nPhi_ = std::max(1, static_cast(std::ceil((phiMax - phiMin) / phiStep))); - } - - inline int nPt() const { return nPt_; } - inline int nEta() const { return nEta_; } - inline int nPhi() const { return nPhi_; } - inline int nM() const { return nM_; } - - inline int binFromValue(float v, float vmin, float step, int nBins) const - { - if (!std::isfinite(v)) return -1; - const float x = (v - vmin) / step; - int b = static_cast(std::floor(x + 1e-6f)); - if (b < 0) return -1; - if (b >= nBins) b = nBins - 1; - return b; - } - - inline int ptBin(float pt) const { return binFromValue(pt, ptMin, ptStep, nPt_); } - inline int etaBin(float eta) const { return binFromValue(eta, etaMin, etaStep, nEta_); } - inline int phiBin(float phi) const { return binFromValue(phi, phiMin, phiStep, nPhi_); } - inline int massBin(float m) const { return binFromValue(m, mMin, mStep, nM_); } -}; - - struct BufferCandK0 { - int64_t collisionIdx; // to enforce different event - int64_t rowIndex; // row id in KShortTracks (for iteratorAt) - uint16_t ptBin, etaBin, phiBin, mBin; - }; - - struct MatchRef { - int64_t collisionIdx; - int64_t rowIndex; - }; - - // Key: (colBin, stat, pt, eta, phi, m) - static inline size_t linearKey(int colBin, int statBin, - int ptBin, int etaBin, int phiBin, int mBin, - int nStat, int nPt, int nEta, int nPhi, int nM) - { - return ((((((static_cast(colBin) * nStat + statBin) * nPt + ptBin) * nEta + etaBin) * nPhi + phiBin) * nM + mBin)); - } + SliceCache cache; + using BinningType = ColumnBinningPolicy; - void processMixedDataMEV4(EventCandidates const& collisions, - aod::KShortTracks const& V0s, - aod::PionTracks const& piontracks) + template + void fillMixedPairs(const TK0sGroup& k0sGroup, + const TPionGroup& pionGroup, + const TCollision& referenceCollision) { - static thread_local std::mt19937 rng(0xBADC0FFE); - - // Build K0 binner from your configurables - MixBinnerK0 mb{ - ptMin.value, ptMax.value, ptMix.value, - v0eta.value, etaMix.value, - phiMix.value - }; - - const int nCol = colBinning.getAllBinsCount(); - const int nStat = N_STATUS; - const int nPt = mb.nPt(); - const int nEta = mb.nEta(); - const int nPhi = mb.nPhi(); - const int nM = mb.nM(); - - const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; - std::vector> buffer(nKeys); - - // ---------------- PASS 1: fill 6D buffer with K0s ---------------- - for (auto const& col : collisions) { - const int colBin = colBinning.getBin(std::make_tuple(col.posz(), col.cent())); - if (colBin < 0) continue; - - auto sliceK0 = V0s.sliceBy(tracksPerCollisionV0, col.index()); - - for (auto const& k0 : sliceK0) { - if (!selectionV0(k0)) continue; - - auto [kpt, keta, kphi] = computePtEtaPhi(k0.kShortPx(), k0.kShortPy(), k0.kShortPz()); - kphi = RecoDecay::constrainAngle(kphi, 0.0F); - - const int ptB = mb.ptBin(kpt); - const int etaB = mb.etaBin(keta); - const int phiB = mb.phiBin(kphi); - const int mB = mb.massBin(k0.kShortMass()); - if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) continue; - - const int stat = 0; - const size_t key = linearKey(colBin, stat, ptB, etaB, phiB, mB, nStat, nPt, nEta, nPhi, nM); - - // IMPORTANT: rowIndex must match iteratorAt accessor below. - const int64_t row = static_cast(k0.globalIndex()); // if this fails -> use k0.index() consistently - - buffer[key].push_back(BufferCandK0{ - .collisionIdx = static_cast(col.index()), - .rowIndex = row, - .ptBin = static_cast(ptB), - .etaBin = static_cast(etaB), - .phiBin = static_cast(phiB), - .mBin = static_cast(mB) - }); + for (const auto& k0sCand : k0sGroup) { + for (const auto& pionCand : pionGroup) { + fillKstarPair(k0sCand, pionCand, referenceCollision, true); } } - - // ---------------- PASS 2: build same-event (K0,π) pairs, replace only K0 by K0X ---------------- - for (auto const& col1 : collisions) { - const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); - if (colBin < 0) continue; - - auto poolK0 = V0s.sliceBy(tracksPerCollisionV0, col1.index()); - auto poolPi = piontracks.sliceBy(tracksPerCollisionBach, col1.index()); - - if (poolK0.size() == 0 || poolPi.size() == 0) continue; // no .empty() in ASoA slice - - const float centrality = col1.cent(); - const float psiFT0Cmix = col1.psiFT0C(); - - // same-event pair loop - for (auto const& k0 : poolK0) { - if (!selectionV0(k0)) continue; - - auto [kpt, keta, kphi] = computePtEtaPhi(k0.kShortPx(), k0.kShortPy(), k0.kShortPz()); - kphi = RecoDecay::constrainAngle(kphi, 0.0F); - - const int ptB = mb.ptBin(kpt); - const int etaB = mb.etaBin(keta); - const int phiB = mb.phiBin(kphi); - const int mB = mb.massBin(k0.kShortMass()); - if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) continue; - - const int stat = 0; - const size_t key = linearKey(colBin, stat, ptB, etaB, phiB, mB, nStat, nPt, nEta, nPhi, nM); - auto const& binVec = buffer[key]; - if (binVec.empty()) continue; - - // pre-collect eligible K0 partners from other events (same bin, different collision) - std::vector matches; - matches.reserve(binVec.size()); - const int64_t curColIdx = static_cast(col1.index()); - - for (auto const& bc : binVec) { - if (bc.collisionIdx == curColIdx) continue; - matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); - } - if (matches.empty()) continue; - - // choose random unique subset of K0 partners - const int cap = nEvtMixing.value; - if (cap > 0 && cap < static_cast(matches.size())) { - std::uniform_int_distribution dist(0, static_cast(matches.size()) - 1); - std::unordered_set chosen; - chosen.reserve(static_cast(cap) * 2); - while (static_cast(chosen.size()) < cap) chosen.insert(dist(rng)); - - std::vector subset; - subset.reserve(cap); - for (const int& idx : chosen) subset.push_back(matches[idx]); - matches.swap(subset); - } else { - std::shuffle(matches.begin(), matches.end(), rng); } - const float wBase = 1.0f / static_cast(matches.size()); - - // Now keep π from SAME event fixed (this is your requirement) - for (auto const& piSame : poolPi) { - if (!passSelectedPID(piSame)) continue; - - // avoid self-combination if you stored K0 daughter indices - if (piSame.pionBachIndex() == k0.pionIndex1() || piSame.pionBachIndex() == k0.pionIndex2()) { - continue; - } - - auto [ppt, peta, pphi] = computePtEtaPhi(piSame.pionBachPx(), piSame.pionBachPy(), piSame.pionBachPz()); - pphi = RecoDecay::constrainAngle(pphi, 0.0F); - pionmix = ROOT::Math::PtEtaPhiMVector(ppt, peta, pphi, o2::constants::physics::MassPionCharged); + void processMixedData(EventCandidates const& collisions, + aod::KShortTracks const& k0sTracks, + aod::PionTracks const& pionTracks) + { + BinningType colBinning{{cfgVtxBins, cfgCentBins}, true}; - // replace K0 by mixed K0X from other event, but keep piSame - for (auto const& m : matches) { - auto k0X = V0s.iteratorAt(m.rowIndex); // if globalIndex fails -> use index consistently in PASS1 + for (const auto& [collision1, collision2] : + selfCombinations(colBinning, nEvtMixing.value, -1, collisions, collisions)) { - if (!selectionV0(k0X)) continue; + if (collision1.globalIndex() == collision2.globalIndex()) { + continue; + } - auto [xpt, xeta, xphi] = computePtEtaPhi(k0X.kShortPx(), k0X.kShortPy(), k0X.kShortPz()); - xphi = RecoDecay::constrainAngle(xphi, 0.0F); - kshortmix = ROOT::Math::PtEtaPhiMVector(xpt, xeta, xphi, k0X.kShortMass()); + // Direction 1: K0s from event 1 + pion from event 2. + for (const auto& k0sCand : k0sTracks) { + if (!k0sBelongsToEvent(k0sCand, collision1)) { + continue; + } - chkstarmix = kshortmix + pionmix; + for (const auto& pionCand : pionTracks) { + if (!pionBelongsToEvent(pionCand, collision2)) { + continue; + } - ROOT::Math::Boost boost{chkstarmix.BoostToCM()}; - fourVecDauCMmix = boost(kshortmix); - threeVecDauCMmix = fourVecDauCMmix.Vect(); + fillKstarPair(k0sCand, pionCand, collision1, true); + } + } - eventplaneVecNormmix = ROOT::Math::XYZVector(std::sin(2.f * psiFT0Cmix), - -std::cos(2.f * psiFT0Cmix), 0.f); + // Direction 2: K0s from event 2 + pion from event 1. + for (const auto& k0sCand : k0sTracks) { + if (!k0sBelongsToEvent(k0sCand, collision2)) { + continue; + } - const float cosThetaStarmix = - eventplaneVecNormmix.Dot(threeVecDauCMmix) / - std::sqrt(threeVecDauCMmix.Mag2()) / - std::sqrt(eventplaneVecNormmix.Mag2()); + for (const auto& pionCand : pionTracks) { + if (!pionBelongsToEvent(pionCand, collision1)) { + continue; + } - histos.fill(HIST("hSparsemix"), - chkstarmix.M(), chkstarmix.Pt(), cosThetaStarmix, centrality, wBase); - } - } + fillKstarPair(k0sCand, pionCand, collision2, true); + } } } } - PROCESS_SWITCH(cksspinalignder, processMixedDataMEV4, - "Process mixed data MEV4 (replace K0, keep same-event pion)", true); - */ + PROCESS_SWITCH(cksspinalignder, processMixedData, "Process mixed-event charged KStar", true); }; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; diff --git a/PWGLF/Tasks/Resonances/phiflowder.cxx b/PWGLF/Tasks/Resonances/phiflowder.cxx index 1293ceb8b3b..0435db8bf70 100644 --- a/PWGLF/Tasks/Resonances/phiflowder.cxx +++ b/PWGLF/Tasks/Resonances/phiflowder.cxx @@ -17,24 +17,18 @@ #include "PWGLF/DataModel/LFPhiFlowTables.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 using namespace o2; @@ -59,9 +53,18 @@ struct phiflowder { ConfigurableAxis axisInvMass{"axisInvMass", {120, 0.9f, 1.2f}, "#it{M}_{K^{+}K^{-}} (GeV/#it{c}^{2})"}; ConfigurableAxis axisPhiPt{"axisPhiPt", {100, 0.f, 10.f}, "#it{p}_{T}^{K^{+}K^{-}} (GeV/#it{c})"}; ConfigurableAxis axisCent{"axisCent", {80, 0.f, 80.f}, "Centrality (%)"}; + ConfigurableAxis axisV1{"axisV1", {1000, -1.0f, 1.0f}, "v1"}; + ConfigurableAxis axisEta{"axisEta", {8, -0.8f, 0.8f}, "Eta"}; ConfigurableAxis axisNKaons{"axisNKaons", {300, 0.f, 300.f}, "Number of stored kaons per event"}; ConfigurableAxis axisNPairs{"axisNPairs", {500, 0.f, 5000.f}, "Number of K^{+}K^{-} pairs per event"}; + Configurable nEvtMixing{"nEvtMixing", 5, "Number of events to mix"}; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {10, -10, 10}, "Mixing bins - z-vertex"}; + ConfigurableAxis cfgCentBins{"cfgCentBins", {8, 0.0, 80}, "Mixing bins - centrality"}; + /*Configurable ptMix{"ptMix", 0.2f, "ME: pT bin width"}; + Configurable etaMix{"etaMix", 0.2f, "ME: eta bin width"}; + Configurable phiMix{"phiMix", 0.3f, "ME: phi bin width"};*/ + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; static constexpr uint8_t kPID1 = 1u << 0; @@ -74,12 +77,13 @@ struct phiflowder { void init(InitContext const&) { histos.add("hCentrality", "Centrality distribution", kTH1F, {axisCent}); - histos.add("hPhiMassSame", "Same-event K^{+}K^{-} invariant mass", kTH1F, {axisInvMass}); - histos.add("hPhiPtSame", "Same-event K^{+}K^{-} pT", kTH1F, {axisPhiPt}); - histos.add("hNplusPerEvent", "Number of stored K^{+} per event", kTH1F, {axisNKaons}); - histos.add("hNminusPerEvent", "Number of stored K^{-} per event", kTH1F, {axisNKaons}); - histos.add("hNphiSamePerEvent", "Number of same-event K^{+}K^{-} pairs per event", kTH1F, {axisNPairs}); - histos.add("hSparseSame", "Same-event K^{+}K^{-};M;pT;centrality", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent}); + histos.add("hSparseSame", "Same-event K^{+}K^{-};M;pT;centrality", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}); + histos.add("hSparseMix", "Mixed-event K^{+}K^{-};M;pT;centrality", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}); + histos.add("hpQxytpvscent", "hpQxytpvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); + histos.add("hpQxpvscent", "hpQxpvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); + histos.add("hpQxtvscent", "hpQxtvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); + histos.add("hpQypvscent", "hpQypvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); + histos.add("hpQytvscent", "hpQytvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); } uint8_t getRequiredPidBit() const @@ -135,6 +139,74 @@ struct phiflowder { return angle >= deepAngleCut.value; } + template + int fillMixedPairs(const TPosKaons& posGroup, + const TNegKaons& negGroup, + float centrality, float qxZDCA, float qxZDCC, float qyZDCA, float qyZDCC) + { + int nMixedPairs = 0; + + for (const auto& kPlus : posGroup) { + if (!passStoredPID(kPlus)) { + continue; + } + + ROOT::Math::PxPyPzMVector plusVec( + kPlus.px(), + kPlus.py(), + kPlus.pz(), + o2::constants::physics::MassKPlus); + + for (const auto& kMinus : negGroup) { + if (!passStoredPID(kMinus)) { + continue; + } + + if (!passDeepAngle(kPlus.px(), + kPlus.py(), + kPlus.pz(), + kMinus.px(), + kMinus.py(), + kMinus.pz())) { + continue; + } + + ROOT::Math::PxPyPzMVector minusVec( + kMinus.px(), + kMinus.py(), + kMinus.pz(), + o2::constants::physics::MassKPlus); + + const auto phiCandidate = plusVec + minusVec; + + const float phiMass = phiCandidate.M(); + const float phiPt = phiCandidate.Pt(); + float etaCand = phiCandidate.Eta(); + double phiCand = phiCandidate.Phi(); + if (phiCand < 0.0) { + phiCand += 2.0 * TMath::Pi(); + } + double const cosNPhi = std::cos(phiCand); + double const sinNPhi = std::sin(phiCand); + + auto ux = cosNPhi; // real part of candidate q vector + auto uy = sinNPhi; // imaginary part of candidate q vector + auto oddv1 = ux * (qxZDCA - qxZDCC) + uy * (qyZDCA - qyZDCC); + // auto evenv1 = ux * (qxZDCA + qxZDCC) + uy * (qyZDCA + qyZDCC); + + if (phiMass < massMin.value || phiMass > massMax.value) { + continue; + } + + histos.fill(HIST("hSparseMix"), phiMass, phiPt, centrality, etaCand, oddv1); + + nMixedPairs++; + } + } + + return nMixedPairs; + } + Filter centralityFilter = (aod::kaonevent::cent >= centMin) && (aod::kaonevent::cent < centMax); using EventCandidates = soa::Filtered; @@ -145,12 +217,27 @@ struct phiflowder { void processSameData(EventCandidates::iterator const& collision, aod::KaonTracks const& /*kaontracks*/) { const float centrality = collision.cent(); + auto qxZDCA = collision.qxA(); + auto qyZDCA = collision.qyA(); + auto qxZDCC = collision.qxC(); + auto qyZDCC = collision.qyC(); + + auto qxtQxp = qxZDCC * qxZDCA; + auto qytQyp = qyZDCC * qyZDCA; + auto qxytp = qxtQxp + qytQyp; + auto qxpQyt = qxZDCA * qyZDCC; + auto qxtQyp = qxZDCC * qyZDCA; histos.fill(HIST("hCentrality"), centrality); + histos.fill(HIST("hpQxytpvscent"), centrality, qxytp); + histos.fill(HIST("hpQxpvscent"), centrality, qxZDCA); + histos.fill(HIST("hpQxtvscent"), centrality, qxZDCC); + histos.fill(HIST("hpQypvscent"), centrality, qyZDCA); + histos.fill(HIST("hpQytvscent"), centrality, qyZDCC); auto posThisColl = posKaons->sliceByCached(aod::kaonpair::kaoneventId, collision.globalIndex(), cache); auto negThisColl = negKaons->sliceByCached(aod::kaonpair::kaoneventId, collision.globalIndex(), cache); - + /* int nPlus = 0; int nMinus = 0; int nPhiSame = 0; @@ -169,7 +256,7 @@ struct phiflowder { histos.fill(HIST("hNplusPerEvent"), nPlus); histos.fill(HIST("hNminusPerEvent"), nMinus); - + */ for (const auto& kPlus : posThisColl) { if (!passStoredPID(kPlus)) { continue; @@ -192,23 +279,191 @@ struct phiflowder { const float phiMass = phiCandidate.M(); const float phiPt = phiCandidate.Pt(); + const float etaCand = phiCandidate.Eta(); + double phiCand = phiCandidate.Phi(); + if (phiCand < 0.0) { + phiCand += 2.0 * TMath::Pi(); + } + double const cosNPhi = std::cos(phiCand); + double const sinNPhi = std::sin(phiCand); + + auto ux = cosNPhi; // real part of candidate q vector + auto uy = sinNPhi; // imaginary part of candidate q vector + auto oddv1 = ux * (qxZDCA - qxZDCC) + uy * (qyZDCA - qyZDCC); + // auto evenv1 = ux * (qxZDCA + qxZDCC) + uy * (qyZDCA + qyZDCC); if (phiMass < massMin || phiMass > massMax) { continue; } - histos.fill(HIST("hPhiMassSame"), phiMass); - histos.fill(HIST("hPhiPtSame"), phiPt); - histos.fill(HIST("hSparseSame"), phiMass, phiPt, centrality); + histos.fill(HIST("hSparseSame"), phiMass, phiPt, centrality, etaCand, oddv1); + // nPhiSame++; + } + } + } + + PROCESS_SWITCH(phiflowder, processSameData, "Process same-event K+K- pairs", true); + + void processLikeSign(EventCandidates::iterator const& collision, aod::KaonTracks const& /*kaontracks*/) + { + const float centrality = collision.cent(); + + auto posThisColl = posKaons->sliceByCached(aod::kaonpair::kaoneventId, collision.globalIndex(), cache); + auto negThisColl = negKaons->sliceByCached(aod::kaonpair::kaoneventId, collision.globalIndex(), cache); + /* + int nPlus = 0; + int nMinus = 0; + int nPhiLike = 0; - nPhiSame++; + for (const auto& kPlus : posThisColl) { + if (passStoredPID(kPlus)) { + nPlus++; } } - histos.fill(HIST("hNphiSamePerEvent"), nPhiSame); + for (const auto& kMinus : negThisColl) { + if (passStoredPID(kMinus)) { + nMinus++; + } + } + + histos.fill(HIST("hNplusLikePerEvent"), nPlus); + histos.fill(HIST("hNminusLikePerEvent"), nMinus); + */ + + // K+K+ pairs + for (const auto& kPlus1 : posThisColl) { + if (!passStoredPID(kPlus1)) { + continue; + } + + ROOT::Math::PxPyPzMVector plusVec1(kPlus1.px(), kPlus1.py(), kPlus1.pz(), o2::constants::physics::MassKPlus); + + for (const auto& kPlus2 : posThisColl) { + if (!passStoredPID(kPlus2)) { + continue; + } + + // Reject self-pairs and count every unique K+K+ pair only once. + if (kPlus2.kaonIndex() <= kPlus1.kaonIndex()) { + continue; + } + + if (!passDeepAngle(kPlus1.px(), kPlus1.py(), kPlus1.pz(), kPlus2.px(), kPlus2.py(), kPlus2.pz())) { + continue; + } + + ROOT::Math::PxPyPzMVector plusVec2(kPlus2.px(), kPlus2.py(), kPlus2.pz(), o2::constants::physics::MassKPlus); + + const auto likeCandidate = plusVec1 + plusVec2; + const float likeMass = likeCandidate.M(); + const float likePt = likeCandidate.Pt(); + + if (likeMass < massMin.value || likeMass > massMax.value) { + continue; + } + histos.fill(HIST("hSparseLike"), likeMass, likePt, centrality); + // nPhiLike++; + } + } + + // K-K- pairs + for (const auto& kMinus1 : negThisColl) { + if (!passStoredPID(kMinus1)) { + continue; + } + + ROOT::Math::PxPyPzMVector minusVec1(kMinus1.px(), kMinus1.py(), kMinus1.pz(), o2::constants::physics::MassKPlus); + + for (const auto& kMinus2 : negThisColl) { + if (!passStoredPID(kMinus2)) { + continue; + } + + // Reject self-pairs and count every unique K-K- pair only once. + if (kMinus2.kaonIndex() <= kMinus1.kaonIndex()) { + continue; + } + + if (!passDeepAngle(kMinus1.px(), kMinus1.py(), kMinus1.pz(), kMinus2.px(), kMinus2.py(), kMinus2.pz())) { + continue; + } + + ROOT::Math::PxPyPzMVector minusVec2(kMinus2.px(), kMinus2.py(), kMinus2.pz(), o2::constants::physics::MassKPlus); + + const auto likeCandidate = minusVec1 + minusVec2; + const float likeMass = likeCandidate.M(); + const float likePt = likeCandidate.Pt(); + + if (likeMass < massMin.value || likeMass > massMax.value) { + continue; + } + + histos.fill(HIST("hSparseLike"), likeMass, likePt, centrality); + // nPhiLike++; + } + } } - PROCESS_SWITCH(phiflowder, processSameData, "Process same-event K+K- pairs", true); + PROCESS_SWITCH(phiflowder, processLikeSign, "Process like-sign K+K+ and K-K- pairs", true); + + // Processing Event Mixing + using BinningType = ColumnBinningPolicy; + BinningType colBinning{{cfgVtxBins, cfgCentBins}, true}; + + void processMixedData(EventCandidates const& collisions, + aod::KaonTracks const& /*kaontracks*/) + { + BinningType colBinning{{cfgVtxBins, cfgCentBins}, true}; + for (const auto& [collision1, collision2] : selfCombinations(colBinning, nEvtMixing.value, -1, collisions, collisions)) { + + if (collision1.globalIndex() == collision2.globalIndex()) { + continue; + } + + const float centrality1 = collision1.cent(); + const float centrality2 = collision2.cent(); + + auto qxZDCA = collision1.qxA(); + auto qyZDCA = collision1.qyA(); + auto qxZDCC = collision1.qxC(); + auto qyZDCC = collision1.qyC(); + + // K+ from event 1 and K- from event 2 + auto posGroup1 = posKaons->sliceByCached( + aod::kaonpair::kaoneventId, + collision1.globalIndex(), + cache); + + auto negGroup2 = negKaons->sliceByCached( + aod::kaonpair::kaoneventId, + collision2.globalIndex(), + cache); + + // K+ from event 2 and K- from event 1 + auto posGroup2 = posKaons->sliceByCached( + aod::kaonpair::kaoneventId, + collision2.globalIndex(), + cache); + + auto negGroup1 = negKaons->sliceByCached( + aod::kaonpair::kaoneventId, + collision1.globalIndex(), + cache); + + int nMixedPairs = 0; + + nMixedPairs += fillMixedPairs(posGroup1, + negGroup2, + centrality1, qxZDCA, qxZDCC, qyZDCA, qyZDCC); + + nMixedPairs += fillMixedPairs(posGroup2, + negGroup1, + centrality2, qxZDCA, qxZDCC, qyZDCA, qyZDCC); + } + } + + PROCESS_SWITCH(phiflowder, processMixedData, "Process mixed-event K+K- pairs", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 588e274c46cbe222864531f56dfa9cbf4c1cc698 Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Thu, 23 Jul 2026 10:16:37 +0200 Subject: [PATCH 2/5] corrected the errors --- PWGLF/TableProducer/Resonances/cksspinalignment.cxx | 2 +- PWGLF/Tasks/Resonances/phiflowder.cxx | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx index a92d1a94329..3796e09ff0f 100644 --- a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx +++ b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx @@ -268,7 +268,7 @@ struct cksspinalignment { float dcaPositive; float dcaNegative; float dcaBetweenDaughters; - float lifetime; + // float lifetime; int64_t positiveIndex; int64_t negativeIndex; }; diff --git a/PWGLF/Tasks/Resonances/phiflowder.cxx b/PWGLF/Tasks/Resonances/phiflowder.cxx index 0435db8bf70..58e85c9d28f 100644 --- a/PWGLF/Tasks/Resonances/phiflowder.cxx +++ b/PWGLF/Tasks/Resonances/phiflowder.cxx @@ -225,8 +225,8 @@ struct phiflowder { auto qxtQxp = qxZDCC * qxZDCA; auto qytQyp = qyZDCC * qyZDCA; auto qxytp = qxtQxp + qytQyp; - auto qxpQyt = qxZDCA * qyZDCC; - auto qxtQyp = qxZDCC * qyZDCA; + // auto qxpQyt = qxZDCA * qyZDCC; + // auto qxtQyp = qxZDCC * qyZDCA; histos.fill(HIST("hCentrality"), centrality); histos.fill(HIST("hpQxytpvscent"), centrality, qxytp); From d358c363ee3123f597d81688a3d9373e362c7c5f Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Thu, 23 Jul 2026 10:30:53 +0200 Subject: [PATCH 3/5] updated --- .../Resonances/cksspinalignment.cxx | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx index 3796e09ff0f..8166cda6b06 100644 --- a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx +++ b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx @@ -470,7 +470,7 @@ struct cksspinalignment { continue; } - const int8_t charge = static_cast(track.sign()); + const auto charge = static_cast(track.sign()); if (charge == 0) { continue; @@ -496,6 +496,9 @@ struct cksspinalignment { auto postrack = v0.template posTrack_as(); auto negtrack = v0.template negTrack_as(); + const auto posId = static_cast(postrack.globalIndex()); + const auto negId = static_cast(negtrack.globalIndex()); + ROOT::Math::PxPyPzMVector pionPos(v0.pxpos(), v0.pypos(), v0.pzpos(), o2::constants::physics::MassPionCharged); ROOT::Math::PxPyPzMVector pionNeg(v0.pxneg(), v0.pyneg(), v0.pzneg(), o2::constants::physics::MassPionCharged); ROOT::Math::PxPyPzMVector k0s = pionPos + pionNeg; @@ -504,18 +507,18 @@ struct cksspinalignment { // v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * // o2::constants::physics::MassK0; - selectedK0s.push_back({static_cast(k0s.Px()), - static_cast(k0s.Py()), - static_cast(k0s.Pz()), - static_cast(k0s.M()), - static_cast(v0.v0cosPA()), - static_cast(v0.v0radius()), - static_cast(std::abs(v0.dcapostopv())), - static_cast(std::abs(v0.dcanegtopv())), - static_cast(std::abs(v0.dcaV0daughters())), + selectedK0s.push_back({k0s.Px(), + k0s.Py(), + k0s.Pz(), + k0s.M(), + v0.v0cosPA(), + v0.v0radius(), + std::abs(v0.dcapostopv()), + std::abs(v0.dcanegtopv()), + std::abs(v0.dcaV0daughters()), // lifetime, - postrack.globalIndex(), - negtrack.globalIndex()}); + posId, + negId}); histos.fill(HIST("hV0Info"), 1.5); histos.fill(HIST("hKShortMass"), k0s.M()); From 46050fa8f8e6a8b653e4a25dd10dc7d2f519fb94 Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Thu, 23 Jul 2026 11:25:58 +0200 Subject: [PATCH 4/5] corrected unused var --- PWGLF/Tasks/Resonances/phiflowder.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/phiflowder.cxx b/PWGLF/Tasks/Resonances/phiflowder.cxx index 58e85c9d28f..017ed45cc93 100644 --- a/PWGLF/Tasks/Resonances/phiflowder.cxx +++ b/PWGLF/Tasks/Resonances/phiflowder.cxx @@ -55,7 +55,6 @@ struct phiflowder { ConfigurableAxis axisCent{"axisCent", {80, 0.f, 80.f}, "Centrality (%)"}; ConfigurableAxis axisV1{"axisV1", {1000, -1.0f, 1.0f}, "v1"}; ConfigurableAxis axisEta{"axisEta", {8, -0.8f, 0.8f}, "Eta"}; - ConfigurableAxis axisNKaons{"axisNKaons", {300, 0.f, 300.f}, "Number of stored kaons per event"}; ConfigurableAxis axisNPairs{"axisNPairs", {500, 0.f, 5000.f}, "Number of K^{+}K^{-} pairs per event"}; Configurable nEvtMixing{"nEvtMixing", 5, "Number of events to mix"}; @@ -84,6 +83,7 @@ struct phiflowder { histos.add("hpQxtvscent", "hpQxtvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); histos.add("hpQypvscent", "hpQypvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); histos.add("hpQytvscent", "hpQytvscent", HistType::kTHnSparseF, {axisCent, axisV1}, true); + histos.add("hMixpairs", "hMixpairs", HistType::kTHnSparseF, {axisNPairs}, true); } uint8_t getRequiredPidBit() const @@ -461,6 +461,8 @@ struct phiflowder { negGroup1, centrality2, qxZDCA, qxZDCC, qyZDCA, qyZDCC); } + + histos.fill(HIST("hMixpairs"), nMixedPairs); } PROCESS_SWITCH(phiflowder, processMixedData, "Process mixed-event K+K- pairs", true); From c73660d2a86e1e5f7f781ead8fb5f1b1e293558e Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Thu, 23 Jul 2026 13:16:23 +0200 Subject: [PATCH 5/5] updated 2 --- .../Resonances/cksspinalignment.cxx | 18 +++++++++--------- PWGLF/Tasks/Resonances/phiflowder.cxx | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx index 8166cda6b06..296f7a93582 100644 --- a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx +++ b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx @@ -507,15 +507,15 @@ struct cksspinalignment { // v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * // o2::constants::physics::MassK0; - selectedK0s.push_back({k0s.Px(), - k0s.Py(), - k0s.Pz(), - k0s.M(), - v0.v0cosPA(), - v0.v0radius(), - std::abs(v0.dcapostopv()), - std::abs(v0.dcanegtopv()), - std::abs(v0.dcaV0daughters()), + selectedK0s.push_back({static_cast(k0s.Px()), + static_cast(k0s.Py()), + static_cast(k0s.Pz()), + static_cast(k0s.M()), + static_cast(v0.v0cosPA()), + static_cast(v0.v0radius()), + static_cast(std::abs(v0.dcapostopv())), + static_cast(std::abs(v0.dcanegtopv())), + static_cast(std::abs(v0.dcaV0daughters())), // lifetime, posId, negId}); diff --git a/PWGLF/Tasks/Resonances/phiflowder.cxx b/PWGLF/Tasks/Resonances/phiflowder.cxx index 017ed45cc93..4b6b044ce7f 100644 --- a/PWGLF/Tasks/Resonances/phiflowder.cxx +++ b/PWGLF/Tasks/Resonances/phiflowder.cxx @@ -460,9 +460,9 @@ struct phiflowder { nMixedPairs += fillMixedPairs(posGroup2, negGroup1, centrality2, qxZDCA, qxZDCC, qyZDCA, qyZDCC); - } - histos.fill(HIST("hMixpairs"), nMixedPairs); + histos.fill(HIST("hMixpairs"), nMixedPairs); + } } PROCESS_SWITCH(phiflowder, processMixedData, "Process mixed-event K+K- pairs", true);