diff --git a/PWGLF/DataModel/LFCKSSpinalignmentTables.h b/PWGLF/DataModel/LFCKSSpinalignmentTables.h index 6593d22fe6b..74abf2f0a82 100644 --- a/PWGLF/DataModel/LFCKSSpinalignmentTables.h +++ b/PWGLF/DataModel/LFCKSSpinalignmentTables.h @@ -27,6 +27,7 @@ #include "Framework/AnalysisDataModel.h" #include +#include namespace o2::aod { @@ -34,7 +35,7 @@ namespace kshortpionevent { DECLARE_SOA_COLUMN(Cent, cent, float); DECLARE_SOA_COLUMN(Posz, posz, float); -DECLARE_SOA_COLUMN(CollIndex, collIndex, 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); @@ -43,7 +44,7 @@ DECLARE_SOA_TABLE(KShortpionEvents, "AOD", "KSHORTPIONEVENT", o2::soa::Index<>, kshortpionevent::Cent, kshortpionevent::Posz, - kshortpionevent::CollIndex, + // kshortpionevent::CollIndex, kshortpionevent::PsiFT0C, kshortpionevent::PsiFT0A, kshortpionevent::PsiTPC) @@ -65,12 +66,13 @@ DECLARE_SOA_COLUMN(KShortMass, kShortMass, float); //! KShort Ma 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(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(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_TABLE(KShortTracks, "AOD", "KSHORTTRACK", o2::soa::Index<>, @@ -98,10 +100,11 @@ DECLARE_SOA_TABLE(PionTracks, "AOD", "PIONTRACK", kshortpionpair::PionBachPy, kshortpionpair::PionBachPz, // kshortpionpair::PionBachSign, - kshortpionpair::PionBachTPC, - kshortpionpair::PionBachTOFHit, - kshortpionpair::PionBachTOF, - kshortpionpair::PionBachIndex); + // kshortpionpair::PionBachTPC, + // kshortpionpair::PionBachTOFHit, + // kshortpionpair::PionBachTOF, + kshortpionpair::PionBachIndex, + kshortpionpair::PionPidMask); using PionTrack = PionTracks::iterator; } // namespace o2::aod diff --git a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx index b76c9382798..f897dcb1dc8 100644 --- a/PWGLF/TableProducer/Resonances/cksspinalignment.cxx +++ b/PWGLF/TableProducer/Resonances/cksspinalignment.cxx @@ -73,6 +73,8 @@ struct cksspinalignment { Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", true, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; } rctCut; + Configurable useNoCollInTimeRangeStandard{"useNoCollInTimeRangeStandard", false, "Apply kNoCollInTimeRangeStandard selection bit"}; + Configurable useGoodITSLayersAll{"useGoodITSLayersAll", true, "Apply kIsGoodITSLayersAll selection bit"}; Configurable cfgCutOccupancy{"cfgCutOccupancy", 2000, "Occupancy cut"}; // events @@ -101,6 +103,13 @@ struct cksspinalignment { Configurable cutTOFBetaPiMeson{"cutTOFBetaPiMeson", 3.0, "Maximum beta cut for pi meson track"}; } grpPion; + enum PionPidBits : uint8_t { + kPID1 = 1u << 0, // selectionPID + kPID2 = 1u << 1, // selectionPID2 + kPID3 = 1u << 2, // selectionPID3 + kPID4 = 1u << 3 // selectionPID4 + }; + // Configs for V0 Configurable confV0PtMin{"confV0PtMin", 0.f, "Minimum transverse momentum of V0"}; Configurable confV0PtMax{"confV0PtMax", 1000.f, "Maximum transverse momentum of V0"}; @@ -120,6 +129,10 @@ struct cksspinalignment { 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"}; @@ -130,6 +143,8 @@ struct cksspinalignment { { 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}); @@ -157,6 +172,81 @@ struct cksspinalignment { return false; } + template + bool selectionPID2(const T& candidate) + { + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < grpPion.nsigmaCutTPCPiMeson) { + return true; + } + if (candidate.hasTOF() && candidate.beta() > grpPion.cutTOFBetaPiMeson && std::sqrt(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) < grpPion.nsigmaCutTOFPiMeson) { + return true; + } + return false; + } + + 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; + } + return false; + } + + 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 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 + } + + // High momentum: TOF hit mandatory + separate 3σ cuts + if (!candidate.hasTOF()) { + return false; + } + + const float nTOF = candidate.tofNSigmaPi(); + 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; + } + template bool selectionV0(Collision const& collision, V0 const& candidate) { @@ -280,11 +370,12 @@ struct cksspinalignment { std::vector dcaBetweenDaughter = {}; std::vector v0Lifetime = {}; // std::vector armenteros = {}; - std::vector pionBachelorIndex = {}; + std::vector pionBachelorIndex = {}; // std::vector pionBachelorSign = {}; - std::vector pionBachelorTPC = {}; - std::vector pionBachelorTOF = {}; - std::vector pionBachelorTOFHit = {}; + // std::vector pionBachelorTPC = {}; + // std::vector pionBachelorTOF = {}; + // std::vector pionBachelorTOFHit = {}; + std::vector pionBachelorPidMask = {}; int numbV0 = 0; auto centrality = collision.centFT0C(); @@ -296,10 +387,12 @@ struct cksspinalignment { // auto psiTPCR = collision.psiTPCR(); // auto psiTPCL = collision.psiTPCL(); 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.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) && 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) && 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); @@ -312,8 +405,12 @@ struct cksspinalignment { continue; } histos.fill(HIST("hTrkSelInfo"), 2.5); - - if (grpPion.usePID && !selectionPID(track1)) { + /* + if (grpPion.usePID && !selectionPID2(track1)) { + continue; + }*/ + const uint8_t mask = pionPidMask(track1); + if (grpPion.usePID && mask == 0) { continue; } histos.fill(HIST("hTrkSelInfo"), 3.5); @@ -322,21 +419,22 @@ struct cksspinalignment { auto track1sign = track1.sign(); if (track1sign == 0) continue; - auto track1nsigTPC = track1.tpcNSigmaPi(); + /*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); + // 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); @@ -365,23 +463,66 @@ struct cksspinalignment { } numbV0 = numbV0 + 1; } - if (numbV0 > 1 && v0Cospa.size() > 1) { - histos.fill(HIST("hEvtSelInfo"), 3.5); - kshortpionEvent(centrality, vz, collision.index(), psiFT0C, psiFT0A, psiTPC); - auto indexEvent = kshortpionEvent.lastIndex(); - //// Fill track table for Charged KStar////////////////// - for (auto icks = kshortMother.begin(); icks != kshortMother.end(); ++icks) { - auto iter = std::distance(kshortMother.begin(), icks); - kshortDummy = kshortMother.at(iter); - // chargedkstarDummy = chargedkstarMother.at(iter); - - kshortTrack(indexEvent, v0Cospa.at(iter), v0Radius.at(iter), dcaPositive.at(iter), dcaNegative.at(iter), dcaBetweenDaughter.at(iter), v0Lifetime.at(iter), kshortDummy.Px(), kshortDummy.Py(), kshortDummy.Pz(), kshortDummy.M(), positiveIndex.at(iter), negativeIndex.at(iter)); + 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; + } } - for (auto ipi = pionBachelor.begin(); ipi != pionBachelor.end(); ++ipi) { - auto iterpi = std::distance(pionBachelor.begin(), ipi); - pionDummy = pionBachelor.at(iterpi); - pionTrack(indexEvent, pionDummy.Px(), pionDummy.Py(), pionDummy.Pz(), pionBachelorTPC.at(iterpi), pionBachelorTOFHit.at(iterpi), pionBachelorTOF.at(iterpi), pionBachelorIndex.at(iterpi)); + // 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]); + } } } } diff --git a/PWGLF/Tasks/Resonances/cksspinalignder.cxx b/PWGLF/Tasks/Resonances/cksspinalignder.cxx index ae29ab6a6fc..9546b96d9d0 100644 --- a/PWGLF/Tasks/Resonances/cksspinalignder.cxx +++ b/PWGLF/Tasks/Resonances/cksspinalignder.cxx @@ -38,12 +38,15 @@ #include // for std::fabs #include // #include +#include #include +#include #include // <<< CHANGED: for dedup sets #include #include #include #include // <<< CHANGED: for seenMap +#include #include #include @@ -77,14 +80,18 @@ struct cksspinalignder { Configurable ptMax{"ptMax", 10.0, "V0 Pt maximum"}; Configurable v0eta{"v0eta", 0.8, "Eta cut on lambda"}; // pion sel///////// - Configurable nsigmaCutTPC{"nsigmaCutTPC", 3.0, "N sigma TPC cut for bachelor pions"}; - Configurable nsigmaCutTOF{"nsigmaCutTOF", 3.0, "N sigma TOF cut for bachelor pions"}; - Configurable usePID{"usePID", false, "Flag for using PID selection"}; + 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})"}; @@ -153,25 +160,47 @@ struct cksspinalignder { return true; } - template - bool selectionPID(const T& candidate) + // ----- 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. + static constexpr uint8_t kPID1 = 1u << 0; + static constexpr uint8_t kPID2 = 1u << 1; + static constexpr uint8_t kPID3 = 1u << 2; + static constexpr uint8_t kPID4 = 1u << 3; + + uint8_t requiredPidBit() const { - auto px = candidate.pionBachPx(); - auto py = candidate.pionBachPy(); - auto pz = candidate.pionBachPz(); - auto p = std::sqrt(px * px + py * py + pz * pz); - float lowmom = 0.5; - if (p < lowmom) { - if (!candidate.pionBachTOFHit() && std::abs(candidate.pionBachTPC()) < nsigmaCutTPC) { - return true; - } else if (candidate.pionBachTOFHit() && std::sqrt(candidate.pionBachTPC() * candidate.pionBachTPC() + candidate.pionBachTOF() * candidate.pionBachTOF()) < nsigmaCutTOF) { - return true; - } - } else if (candidate.pionBachTOFHit() && std::sqrt(candidate.pionBachTPC() * candidate.pionBachTPC() + candidate.pionBachTOF() * candidate.pionBachTOF()) < nsigmaCutTOF) { + const int choice = pidChoice.value; + switch (choice) { + case 1: + return kPID1; + case 2: + return kPID2; + case 3: + return kPID3; + case 4: + return kPID4; + default: + LOGF(warn, "pidChoice=%d is invalid. Using PID2.", choice); + return kPID2; + } + } + + template + bool passSelectedPID(const TPion& pion) const + { + if (!usePID.value) { return true; } - return false; + const uint8_t bit = requiredPidBit(); + const uint8_t mask = pion.pionPidMask(); + + if (pidExclusive.value) { + return mask == bit; // ONLY chosen PID + } + return (mask & bit) != 0; // chosen PID (maybe others too) } + // --------------------------------------------------------------- std::tuple computePtEtaPhi(float px, float py, float pz) { @@ -207,15 +236,17 @@ struct cksspinalignder { 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 (usePID && !selectionPID(piontrack)) - continue; // checking PID - chkstar = kshort + pion; ROOT::Math::Boost boost{chkstar.BoostToCM()}; @@ -254,9 +285,11 @@ struct cksspinalignder { auto [kshortPtmix, kshortEtamix, kshortPhimix] = computePtEtaPhi(t1.kShortPx(), t1.kShortPy(), t1.kShortPz()); kshortmix = ROOT::Math::PtEtaPhiMVector(kshortPtmix, kshortEtamix, kshortPhimix, t1.kShortMass()); + // 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()); - if (usePID && !selectionPID(t2)) - continue; // checking PID pionmix = ROOT::Math::PtEtaPhiMVector(pionPtmix, pionEtamix, pionPhimix, o2::constants::physics::MassPionCharged); chkstarmix = kshortmix + pionmix; @@ -272,6 +305,526 @@ struct cksspinalignder { } } 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))); + } + + 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_); } // 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)); + } + + // ===================================================================================== + // 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) + { + 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) { + + // apply configurable PID selection stored in bitmask (your passSelectedPID) + if (!passSelectedPID(pi)) + continue; + + auto [pt, eta, phi] = computePtEtaPhi(pi.pionBachPx(), pi.pionBachPy(), pi.pionBachPz()); + phi = RecoDecay::constrainAngle(phi, 0.0F); + + 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; + + const int stat = 0; + const size_t key = linearKey(colBin, stat, ptB, etaB, phiB, mB, nStat, nPt, nEta, nPhi, nM); + + const int64_t row = static_cast(pi.globalIndex()); + + 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)}); + } + } + + // ---------------- 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; + + auto poolK0 = V0s.sliceBy(tracksPerCollisionV0, col1.index()); + auto poolPi = piontracks.sliceBy(tracksPerCollisionBach, col1.index()); + if (poolK0.size() == 0 || poolPi.size() == 0) + continue; + + const float centrality = col1.cent(); + const float psiFT0Cmix = col1.psiFT0C(); + + // Loop over K0s from the *current* event (fixed object) + for (auto const& k0 : poolK0) { + if (!selectionV0(k0)) + continue; + + // 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()); + + // Loop over same-event pions to define the "same-event pair" base + for (auto const& piSame : poolPi) { + if (!passSelectedPID(piSame)) + continue; + + // avoid self-combination: bachelor pion can't be a K0 daughter + if (piSame.pionBachIndex() == k0.pionIndex1() || piSame.pionBachIndex() == k0.pionIndex2()) { + continue; + } + + // 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); + + 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; + + 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; + + // collect partners from same bin but 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; + + // 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); + } + + // const float wBase = 1.0f / static_cast(matches.size()); + const float wBase = 1.0; + + // 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 + + // safety: PID selection again + if (!passSelectedPID(piX)) + continue; + + // OPTIONAL keep: also avoid piX being a K0 daughter (same k0 from current event) + if (piX.pionBachIndex() == k0.pionIndex1() || piX.pionBachIndex() == k0.pionIndex2()) { + continue; + } + + 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); + + chkstarmix = kshortmix + pionmix; + + ROOT::Math::Boost boost{chkstarmix.BoostToCM()}; + fourVecDauCMmix = boost(kshortmix); + threeVecDauCMmix = fourVecDauCMmix.Vect(); + + eventplaneVecNormmix = ROOT::Math::XYZVector(std::sin(2.f * psiFT0Cmix), + -std::cos(2.f * psiFT0Cmix), 0.f); + + const float cosThetaStarmix = + eventplaneVecNormmix.Dot(threeVecDauCMmix) / + std::sqrt(threeVecDauCMmix.Mag2()) / + std::sqrt(eventplaneVecNormmix.Mag2()); + + histos.fill(HIST("hSparsemix"), + chkstarmix.M(), chkstarmix.Pt(), cosThetaStarmix, centrality, wBase); + } + } + } + } + } + + 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_; + + 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)); + } + + void processMixedDataMEV4(EventCandidates const& collisions, + aod::KShortTracks const& V0s, + aod::PionTracks const& piontracks) + { + 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) + }); + } + } + + // ---------------- 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); + + // 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 + + if (!selectionV0(k0X)) 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()); + + chkstarmix = kshortmix + pionmix; + + ROOT::Math::Boost boost{chkstarmix.BoostToCM()}; + fourVecDauCMmix = boost(kshortmix); + threeVecDauCMmix = fourVecDauCMmix.Vect(); + + eventplaneVecNormmix = ROOT::Math::XYZVector(std::sin(2.f * psiFT0Cmix), + -std::cos(2.f * psiFT0Cmix), 0.f); + + const float cosThetaStarmix = + eventplaneVecNormmix.Dot(threeVecDauCMmix) / + std::sqrt(threeVecDauCMmix.Mag2()) / + std::sqrt(eventplaneVecNormmix.Mag2()); + + histos.fill(HIST("hSparsemix"), + chkstarmix.M(), chkstarmix.Pt(), cosThetaStarmix, centrality, wBase); + } + } + } + } + } + + PROCESS_SWITCH(cksspinalignder, processMixedDataMEV4, + "Process mixed data MEV4 (replace K0, keep same-event pion)", true); + */ }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {