diff --git a/PWGLF/TableProducer/Resonances/doublephitable.cxx b/PWGLF/TableProducer/Resonances/doublephitable.cxx index 57744ddf693..b2819b76cef 100644 --- a/PWGLF/TableProducer/Resonances/doublephitable.cxx +++ b/PWGLF/TableProducer/Resonances/doublephitable.cxx @@ -720,6 +720,12 @@ struct doublephitable { const auto phi2Fit = k3Fit + k4Fit; const auto pairFit = phi1Fit + phi2Fit; + if (pairFit.Pt() < 6.0) { + continue; + } + if (pairFit.M() < 2.5 || pairFit.M() > 3.2) { + continue; + } // No pair-pT or pair-mass cut here. PhiPhiPairPayload pair; diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index 043d798ab62..a14647b525a 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -67,13 +67,14 @@ struct doublephimeson { Configurable additionalEvsel{"additionalEvsel", true, "Additional event selection"}; Configurable isDeep{"isDeep", true, "Store deep angle"}; Configurable cutMinNsigmaTPC{"cutMinNsigmaTPC", -2.5, "nsigma cut TPC"}; - Configurable cutNsigmaTPC{"cutNsigmaTPC", 3.0, "nsigma cut TPC"}; - Configurable cutNsigmaTOF{"cutNsigmaTOF", 3.0, "nsigma cut TOF"}; + Configurable cutNsigmaTPC{"cutNsigmaTPC", 2.5, "nsigma cut TPC"}; + Configurable cutNsigmaTOF{"cutNsigmaTOF", 2.5, "nsigma cut TOF"}; Configurable momTOFCut{"momTOFCut", 1.8, "minimum pT cut for madnatory TOF"}; Configurable maxKaonPt{"maxKaonPt", 100.0, "maximum kaon pt cut"}; Configurable cfgCrossPhiLow{"cfgCrossPhiLow", 1.01, "Lower edge of phi mass window for cross-pairing (ghost) veto"}; Configurable cfgCrossPhiHigh{"cfgCrossPhiHigh", 1.03, "Upper edge of phi mass window for cross-pairing (ghost) veto"}; Configurable useParametrized{"useParametrized", false, "Use pT dependent mass peak and width"}; + Configurable useCrossPairRejection{"useCrossPairRejection", true, "Use cross pair phi signal compatibilaty"}; // ------------------------------------------------------------ // pT-dependent phi mass peak and width from single-phi BW fits // @@ -337,17 +338,17 @@ struct doublephimeson { } if (PIDStrategy == 1003) { - if (ptcand < 0.5 && TOFHit != 1 && nsigmaTPC > -2.0 && nsigmaTPC < 2.0) { + if (ptcand < 0.5 && TOFHit != 1 && std::abs(nsigmaTPC) < cutNsigmaTPC) { return true; } - if (ptcand < 0.5 && TOFHit == 1 && std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < 2.5) { + if (ptcand < 0.5 && TOFHit == 1 && std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < cutNsigmaTPC) { return true; } if (ptcand >= 0.5) { - if (TOFHit != 1 && nsigmaTPC > -2.0 && nsigmaTPC < 2.0) { + if (TOFHit != 1 && nsigmaTPC > -2.0 && nsigmaTPC < cutNsigmaTPC) { return true; } - if (TOFHit == 1 && std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < 2.5) { + if (TOFHit == 1 && std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < cutNsigmaTOF) { return true; } } @@ -1695,7 +1696,6 @@ struct doublephimeson { if (id2 <= id1) { continue; } - const double kplus2pt = std::hypot(t2.phid1Px(), t2.phid1Py()); const double kminus2pt = std::hypot(t2.phid2Px(), t2.phid2Py()); @@ -1736,25 +1736,20 @@ struct doublephimeson { t1.phid1Index() == t2.phid2Index() || t1.phid2Index() == t2.phid1Index() || t1.phid2Index() == t2.phid2Index()) { - // LOGF(info,"track share",t1.phid1Index(),t1.phid2Index(),t2.phid1Index(),t2.phid2Index()); + LOGF(info, "track share %d %d %d %d", t1.phid1Index(), t1.phid2Index(), t2.phid1Index(), t2.phid2Index()); continue; } - const double mCross12 = (k1p + k2m).M(); // K+ from phi1 + K- from phi2 - const double mCross21 = (k2p + k1m).M(); // K+ from phi2 + K- from phi1 - - const bool cross12IsPhiLike = (mCross12 > cfgCrossPhiLow && mCross12 < cfgCrossPhiHigh); - const bool cross21IsPhiLike = (mCross21 > cfgCrossPhiLow && mCross21 < cfgCrossPhiHigh); - if (cross12IsPhiLike || cross21IsPhiLike) { - LOGF(info, - "Best-pairing rejected: mPhi1 = %3.4f, mPhi2 = %3.4f, mCross12 = %3.4f, mCross21 = %3.4f", - phi1.M(), - phi2.M(), - mCross12, - mCross21, - pair.Pt(), - pair.M()); - continue; + auto cross12 = k1p + k2m; + auto cross21 = k2p + k1m; + bool alternativePairValid = cross12.M() > cfgCrossPhiLow && cross12.M() < cfgCrossPhiHigh && cross21.M() > cfgCrossPhiLow && cross21.M() < cfgCrossPhiHigh; + if (alternativePairValid) { + float scoreOriginal = deltaMPhiNominal(phi1.M(), phi2.M()); + float scoreCross = deltaMPhiNominal(cross12.M(), cross21.M()); + if (useCrossPairRejection && (scoreCross < scoreOriginal)) { + LOGF(info, "Best-pairing rejected: original score = %3.4f, cross scoremPhi2 = %3.4f", scoreOriginal, scoreCross); + continue; // another pairing of these four tracks is better + } } histos.fill(HIST("hPhiMass"), phi1.M(), phi2.M(), pair.Pt()); histos.fill(HIST("hPhiMassNormalized"), getNormalizedMPhi(phi1.M(), phi1.Pt()), getNormalizedMPhi(phi2.M(), phi2.Pt()), pair.Pt()); @@ -1808,7 +1803,7 @@ struct doublephimeson { } if (pairPt > minExoticPt) { histos.fill(HIST("hPtCorrelation"), pairPt, ptcorr); - histos.fill(HIST("hMassCent"), p1.M(), p2.M(), collision.centrality()); + // histos.fill(HIST("hMassCent"), p1.M(), p2.M(), collision.centrality()); histos.fill(HIST("SEMassUnlike_AllVars"), M, pairPt, diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 76b9d788c63..d923663c5ce 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -47,8 +47,6 @@ #include #include #include -#include -#include #include #include #include @@ -221,19 +219,24 @@ struct lambdaspincorrderived { Configurable useNUA{"useNUA", false, "Apply single-candidate NUA weight in (phi,eta)"}; Configurable ConfNUAPathLambda{"ConfNUAPathLambda", "", "CCDB path for Lambda NUA TH2D(phi,eta)"}; Configurable ConfNUAPathAntiLambda{"ConfNUAPathAntiLambda", "", "CCDB path for AntiLambda NUA TH2D(phi,eta)"}; + + // Mixing selection and binning Configurable> massMixEdges{"massMixEdges", {1.09f, 1.108f, 1.122f, 1.14f}, "Mass-mixing region edges: [SB low | signal | SB high]"}; Configurable cfgMixLegMode{"cfgMixLegMode", 0, "0=replace leg-1 only, 1=replace leg-2 only, 2=do both one-leg replacements"}; Configurable cfgV5MassBins{"cfgV5MassBins", 5, "Number of fixed mass bins for V5 mixing"}; Configurable cfgV5NeighborPt{"cfgV5NeighborPt", 0, "v5: neighbor bins in pT (use symmetric ±N, edge-safe)"}; Configurable cfgV5NeighborEta{"cfgV5NeighborEta", 0, "v5: neighbor bins in eta (use symmetric ±N, edge-safe)"}; Configurable cfgV5NeighborPhi{"cfgV5NeighborPhi", 0, "v5: neighbor bins in phi (use symmetric ±N, periodic wrap)"}; - Configurable usePairKineMatch{"usePairKineMatch", true, "Require pair-level matching between (A,B) and (C,B)"}; Configurable cfgV5MaxMatches{"cfgV5MaxMatches", 50, "v5: max ME replacements per SE pair (after all cuts)"}; Configurable cfgMixSeed{"cfgMixSeed", 0xdecafbadULL, "RNG seed for downsampling matches (deterministic)"}; + Configurable cfgV6CarryUnmatched{"cfgV6CarryUnmatched", false, "Carry data and MC V6 replacement branches with zero matches to later data frames"}; + Configurable cfgV6MaxPendingBranches{"cfgV6MaxPendingBranches", 0, "Maximum retained V6 branches per data or MC buffer; <=0 means unlimited"}; + Configurable cfgV6MaxPendingAge{"cfgV6MaxPendingAge", 0, "Maximum later data frames searched by a retained V6 branch; <=0 means unlimited"}; + Configurable cfgV6LogPending{"cfgV6LogPending", false, "Print one V6 pending-branch summary per data frame"}; + + // Event and candidate mixing cuts Configurable centMin{"centMin", 0, "Minimum Centrality"}; Configurable centMax{"centMax", 80, "Maximum Centrality"}; - Configurable rngSeed{"rngSeed", 12345, "Seed for random mixing (reproducible)"}; - std::mt19937 rng{12345}; Configurable nEvtMixing{"nEvtMixing", 10, "Number of events to mix"}; ConfigurableAxis CfgVtxBins{"CfgVtxBins", {VARIABLE_WIDTH, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10}, "Mixing bins - z-vertex"}; ConfigurableAxis CfgMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0, 110}, "Mixing bins - centrality"}; @@ -478,7 +481,6 @@ struct lambdaspincorrderived { histos.add("hSparsePairMassAntiLambdaLambdaMixed", "hSparsePairMassAntiLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisPairMass}, true); histos.add("hSparsePairMassAntiLambdaAntiLambdaMixed", "hSparsePairMassAntiLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisPairMass}, true); } - rng.seed(static_cast(rngSeed.value)); ccdb->setURL(cfgCcdbParam.cfgURL); ccdbApi.init("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); @@ -1587,8 +1589,6 @@ struct lambdaspincorrderived { 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_); } - inline int radiusBin(float r) const { if (!std::isfinite(r) || nR_ <= 0) { @@ -1609,6 +1609,210 @@ struct lambdaspincorrderived { uint16_t ptBin, etaBin, phiBin, mBin, rBin; }; + struct StoredV6Candidate { + int64_t collisionIdx = -1; + int64_t globalIdx = -1; + int status = -1; + bool isDouble = false; + float cospa = 0.f; + float radius = 0.f; + float dcaPos = 0.f; + float dcaNeg = 0.f; + float dcaDau = 0.f; + float lPt = 0.f; + float lEta = 0.f; + float lPhi = 0.f; + float lMass = 0.f; + float pPt = 0.f; + float pEta = 0.f; + float pPhi = 0.f; + int64_t pIndex = -1; + int64_t piIndex = -1; + + int v0Status() const { return status; } + bool doubleStatus() const { return isDouble; } + float v0Cospa() const { return cospa; } + float v0Radius() const { return radius; } + float dcaPositive() const { return dcaPos; } + float dcaNegative() const { return dcaNeg; } + float dcaBetweenDaughter() const { return dcaDau; } + float lambdaPt() const { return lPt; } + float lambdaEta() const { return lEta; } + float lambdaPhi() const { return lPhi; } + float lambdaMass() const { return lMass; } + float protonPt() const { return pPt; } + float protonEta() const { return pEta; } + float protonPhi() const { return pPhi; } + int64_t protonIndex() const { return pIndex; } + int64_t pionIndex() const { return piIndex; } + int64_t globalIndex() const { return globalIdx; } + }; + + struct PendingV6Branch { + StoredV6Candidate target; + StoredV6Candidate fixed; + int colBin = -1; + int replacedLeg = 1; + int age = 0; + uint64_t seed = 0; + }; + + template + StoredV6Candidate storeV6Candidate(T const& t, int64_t collisionIdx) const + { + return {collisionIdx, static_cast(t.globalIndex()), static_cast(t.v0Status()), static_cast(t.doubleStatus()), + t.v0Cospa(), t.v0Radius(), t.dcaPositive(), t.dcaNegative(), t.dcaBetweenDaughter(), + t.lambdaPt(), t.lambdaEta(), t.lambdaPhi(), t.lambdaMass(), + t.protonPt(), t.protonEta(), t.protonPhi(), + static_cast(t.protonIndex()), static_cast(t.pionIndex())}; + } + + template + void fillV6MixedBranch(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) + { + const auto repProton = ROOT::Math::PtEtaPhiMVector(replacement.protonPt(), replacement.protonEta(), replacement.protonPhi(), o2::constants::physics::MassProton); + const auto repLambda = ROOT::Math::PtEtaPhiMVector(replacement.lambdaPt(), replacement.lambdaEta(), replacement.lambdaPhi(), replacement.lambdaMass()); + const auto fixedProton = ROOT::Math::PtEtaPhiMVector(fixed.protonPt(), fixed.protonEta(), fixed.protonPhi(), o2::constants::physics::MassProton); + const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(fixed.lambdaPt(), fixed.lambdaEta(), fixed.lambdaPhi(), fixed.lambdaMass()); + const int repStatus = replacement.v0Status(); + const int fixedStatus = fixed.v0Status(); + + if (replacedLeg == 1) { + fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); + fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); + if ((repStatus == 0 && fixedStatus == 1) || (repStatus == 1 && fixedStatus == 0)) { + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); + } + } + if (repStatus == 0 && fixedStatus == 1) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } else if (repStatus == 1 && fixedStatus == 0) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); + } else { + fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } + return; + } + + fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); + fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); + } + if (fixedStatus == 0 && repStatus == 1) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } else if (fixedStatus == 1 && repStatus == 0) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); + } else { + fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } + } + + struct StoredV6CandidateMC { + int64_t globalIdx = -1; + int status = -1; + bool isDouble = false; + float cospa = 0.f; + float radius = 0.f; + float dcaPos = 0.f; + float dcaNeg = 0.f; + float dcaDau = 0.f; + float lPt = 0.f; + float lEta = 0.f; + float lPhi = 0.f; + float lMass = 0.f; + float pPt = 0.f; + float pEta = 0.f; + float pPhi = 0.f; + int64_t pIndex = -1; + int64_t piIndex = -1; + + int v0Statusmc() const { return status; } + bool doubleStatusmc() const { return isDouble; } + float v0Cospamc() const { return cospa; } + float v0Radiusmc() const { return radius; } + float dcaPositivemc() const { return dcaPos; } + float dcaNegativemc() const { return dcaNeg; } + float dcaBetweenDaughtermc() const { return dcaDau; } + float lambdaPtmc() const { return lPt; } + float lambdaEtamc() const { return lEta; } + float lambdaPhimc() const { return lPhi; } + float lambdaMassmc() const { return lMass; } + float protonPtmc() const { return pPt; } + float protonEtamc() const { return pEta; } + float protonPhimc() const { return pPhi; } + int64_t protonIndexmc() const { return pIndex; } + int64_t pionIndexmc() const { return piIndex; } + int64_t globalIndex() const { return globalIdx; } + }; + + struct PendingV6BranchMC { + StoredV6CandidateMC target; + StoredV6CandidateMC fixed; + int colBin = -1; + int replacedLeg = 1; + int age = 0; + uint64_t seed = 0; + }; + + struct V6PendingState { + std::deque data; + std::deque mc; + }; + + V6PendingState v6Pending; + + template + StoredV6CandidateMC storeV6CandidateMC(T const& t) const + { + return {static_cast(t.globalIndex()), mcacc::v0Status(t), mcacc::doubleStatus(t), + mcacc::v0CosPA(t), mcacc::v0Radius(t), mcacc::dcaPos(t), mcacc::dcaNeg(t), mcacc::dcaDau(t), + mcacc::lamPt(t), mcacc::lamEta(t), mcacc::lamPhi(t), mcacc::lamMass(t), + mcacc::prPt(t), mcacc::prEta(t), mcacc::prPhi(t), + static_cast(mcacc::prIdx(t)), static_cast(mcacc::piIdx(t))}; + } + + template + void fillV6MixedBranchMC(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) + { + const auto repProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(replacement), mcacc::prEta(replacement), mcacc::prPhi(replacement), o2::constants::physics::MassProton); + const auto repLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(replacement), mcacc::lamEta(replacement), mcacc::lamPhi(replacement), mcacc::lamMass(replacement)); + const auto fixedProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(fixed), mcacc::prEta(fixed), mcacc::prPhi(fixed), o2::constants::physics::MassProton); + const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(fixed), mcacc::lamEta(fixed), mcacc::lamPhi(fixed), mcacc::lamMass(fixed)); + const int repStatus = mcacc::v0Status(replacement); + const int fixedStatus = mcacc::v0Status(fixed); + + if (replacedLeg == 1) { + fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); + fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); + } + if (repStatus == 0 && fixedStatus == 1) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } else if (repStatus == 1 && fixedStatus == 0) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); + } else { + fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } + return; + } + + fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); + fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); + } + if (fixedStatus == 0 && repStatus == 1) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } else if (fixedStatus == 1 && repStatus == 0) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); + } else { + fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } + } + static inline size_t linearKeyR(int colBin, int statBin, int ptBin, int etaBin, int phiBin, int mBin, int rBin, int nStatus, int nPt, int nEta, int nPhi, int nM, int nR) @@ -2129,6 +2333,68 @@ struct lambdaspincorrderived { } }; + const size_t pendingAtStart = v6Pending.data.size(); + size_t pendingMatched = 0; + size_t pendingExpired = 0; + size_t pendingAdded = 0; + if (!cfgV6CarryUnmatched) { + v6Pending.data.clear(); + } else { + for (auto it = v6Pending.data.begin(); it != v6Pending.data.end();) { + auto& pending = *it; + ++pending.age; + if (cfgV6MaxPendingAge.value > 0 && pending.age > cfgV6MaxPendingAge.value) { + ++pendingExpired; + it = v6Pending.data.erase(it); + continue; + } + + auto& matches = pending.replacedLeg == 1 ? matches1 : matches2; + collectMatchesForReplacedLeg(pending.target, pending.fixed, pending.colBin, -1, matches); + limitMatchesToNEvents(matches, nEvtMixing.value); + downsampleMatches(matches, pending.seed ^ splitmix64(static_cast(pending.age))); + + int nAccepted = 0; + for (auto const& m : matches) { + auto replacement = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(replacement) || !checkKinematics(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughters(replacement, pending.target) || hasSharedDaughters(replacement, pending.fixed)) { + continue; + } + ++nAccepted; + } + + if (nAccepted == 0) { + ++it; + continue; + } + + const float controlWeight = 1.0f / static_cast(nAccepted); + const float branchNorm = cfgMixLegMode.value == 2 ? 0.5f : 1.0f; + const float mixWeight = branchNorm * controlWeight; + for (auto const& m : matches) { + auto replacement = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(replacement) || !checkKinematics(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughters(replacement, pending.target) || hasSharedDaughters(replacement, pending.fixed)) { + continue; + } + fillV6MixedBranch(replacement, pending.fixed, pending.replacedLeg, controlWeight, mixWeight); + } + ++pendingMatched; + it = v6Pending.data.erase(it); + } + } + // -------- PASS 2: configurable one-leg / two-leg mixing -------- for (auto const& col1 : collisions) { const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); @@ -2235,6 +2501,22 @@ struct lambdaspincorrderived { } } + if (cfgV6CarryUnmatched) { + const auto hasPendingSpace = [&]() { + return cfgV6MaxPendingBranches.value <= 0 || static_cast(v6Pending.data.size()) < cfgV6MaxPendingBranches.value; + }; + if (doMixLeg1 && nFill1 == 0 && hasPendingSpace()) { + v6Pending.data.push_back({storeV6Candidate(t1, curColIdx), storeV6Candidate(t2, curColIdx), colBin, 1, 0, + static_cast(t1.globalIndex()) ^ splitmix64(static_cast(t2.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAdded; + } + if (doMixLeg2 && nFill2 == 0 && hasPendingSpace()) { + v6Pending.data.push_back({storeV6Candidate(t2, curColIdx), storeV6Candidate(t1, curColIdx), colBin, 2, 0, + static_cast(t2.globalIndex()) ^ splitmix64(static_cast(t1.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAdded; + } + } + if (nFill1 <= 0 && nFill2 <= 0) { continue; } @@ -2246,9 +2528,10 @@ struct lambdaspincorrderived { const int nActiveMixBranches = ((doMixLeg1 && nFill1 > 0) ? 1 : 0) + ((doMixLeg2 && nFill2 > 0) ? 1 : 0); - const float branchNorm = (cfgMixLegMode.value == 2 && nActiveMixBranches > 0) - ? 1.0f / static_cast(nActiveMixBranches) - : 1.0f; + float branchNorm = 1.0f; + if (cfgMixLegMode.value == 2) { + branchNorm = cfgV6CarryUnmatched ? 0.5f : 1.0f / static_cast(nActiveMixBranches); + } const float finalMixWeightLeg1 = branchNorm * wSELeg1; const float finalMixWeightLeg2 = branchNorm * wSELeg2; @@ -2356,6 +2639,10 @@ struct lambdaspincorrderived { } } } + if (cfgV6LogPending) { + LOGF(info, "MEV6 data pending branches: carriedIn=%zu matched=%zu expired=%zu newlyPropagated=%zu carriedToNext=%zu", + pendingAtStart, pendingMatched, pendingExpired, pendingAdded, v6Pending.data.size()); + } } PROCESS_SWITCH(lambdaspincorrderived, processMEV6, "Process data ME v6 with radius buffer", false); @@ -2592,6 +2879,68 @@ struct lambdaspincorrderived { } }; + const size_t pendingAtStartMC = v6Pending.mc.size(); + size_t pendingMatchedMC = 0; + size_t pendingExpiredMC = 0; + size_t pendingAddedMC = 0; + if (!cfgV6CarryUnmatched) { + v6Pending.mc.clear(); + } else { + for (auto it = v6Pending.mc.begin(); it != v6Pending.mc.end();) { + auto& pending = *it; + ++pending.age; + if (cfgV6MaxPendingAge.value > 0 && pending.age > cfgV6MaxPendingAge.value) { + ++pendingExpiredMC; + it = v6Pending.mc.erase(it); + continue; + } + + auto& matches = pending.replacedLeg == 1 ? matches1 : matches2; + collectMatchesForReplacedLeg(pending.target, pending.fixed, pending.colBin, -1, matches); + limitMatchesToNEvents(matches, nEvtMixing.value); + downsampleMatches(matches, pending.seed ^ splitmix64(static_cast(pending.age))); + + int nAccepted = 0; + for (auto const& m : matches) { + auto replacement = V0sMC.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0MC(replacement) || !checkKinematicsMC(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughtersMC(replacement, pending.target) || hasSharedDaughtersMC(replacement, pending.fixed)) { + continue; + } + ++nAccepted; + } + + if (nAccepted == 0) { + ++it; + continue; + } + + const float controlWeight = 1.0f / static_cast(nAccepted); + const float branchNorm = cfgMixLegMode.value == 2 ? 0.5f : 1.0f; + const float mixWeight = branchNorm * controlWeight; + for (auto const& m : matches) { + auto replacement = V0sMC.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0MC(replacement) || !checkKinematicsMC(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughtersMC(replacement, pending.target) || hasSharedDaughtersMC(replacement, pending.fixed)) { + continue; + } + fillV6MixedBranchMC(replacement, pending.fixed, pending.replacedLeg, controlWeight, mixWeight); + } + ++pendingMatchedMC; + it = v6Pending.mc.erase(it); + } + } + // -------- PASS 2: configurable one-leg / two-leg mixing -------- for (auto const& col1 : collisions) { const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col1), mcacc::cent(col1))); @@ -2699,6 +3048,22 @@ struct lambdaspincorrderived { } } + if (cfgV6CarryUnmatched) { + const auto hasPendingSpace = [&]() { + return cfgV6MaxPendingBranches.value <= 0 || static_cast(v6Pending.mc.size()) < cfgV6MaxPendingBranches.value; + }; + if (doMixLeg1 && nFill1 == 0 && hasPendingSpace()) { + v6Pending.mc.push_back({storeV6CandidateMC(t1), storeV6CandidateMC(t2), colBin, 1, 0, + static_cast(t1.globalIndex()) ^ splitmix64(static_cast(t2.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAddedMC; + } + if (doMixLeg2 && nFill2 == 0 && hasPendingSpace()) { + v6Pending.mc.push_back({storeV6CandidateMC(t2), storeV6CandidateMC(t1), colBin, 2, 0, + static_cast(t2.globalIndex()) ^ splitmix64(static_cast(t1.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAddedMC; + } + } + if (nFill1 <= 0 && nFill2 <= 0) { continue; } @@ -2708,9 +3073,10 @@ struct lambdaspincorrderived { const int nActiveMixBranches = ((doMixLeg1 && nFill1 > 0) ? 1 : 0) + ((doMixLeg2 && nFill2 > 0) ? 1 : 0); - const float branchNorm = (cfgMixLegMode.value == 2 && nActiveMixBranches > 0) - ? 1.0f / static_cast(nActiveMixBranches) - : 1.0f; + float branchNorm = 1.0f; + if (cfgMixLegMode.value == 2) { + branchNorm = cfgV6CarryUnmatched ? 0.5f : 1.0f / static_cast(nActiveMixBranches); + } const float finalMixWeightLeg1 = branchNorm * wSELeg1; const float finalMixWeightLeg2 = branchNorm * wSELeg2; @@ -2826,8 +3192,21 @@ struct lambdaspincorrderived { } } } + if (cfgV6LogPending) { + LOGF(info, "MEV6 MC pending branches: carriedIn=%zu matched=%zu expired=%zu newlyPropagated=%zu carriedToNext=%zu", + pendingAtStartMC, pendingMatchedMC, pendingExpiredMC, pendingAddedMC, v6Pending.mc.size()); + } } PROCESS_SWITCH(lambdaspincorrderived, processMCMEV6, "Process MC ME v6 with radius buffer", false); + + void endOfStream(EndOfStreamContext&) + { + if (cfgV6LogPending) { + LOGF(info, "MEV6 end of stream: unmatchedData=%zu unmatchedMC=%zu", v6Pending.data.size(), v6Pending.mc.size()); + } + v6Pending.data.clear(); + v6Pending.mc.clear(); + } }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {