[PWGLF] Add new event-mixing process function

PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx

@@ -30,14 +30,14 @@
 #include <Math/GenVector/Boost.h>
 #include <Math/Vector3D.h>
 #include <Math/Vector4D.h>
 #include <TLorentzVector.h>
 #include <TMath.h>
 
 #include <fairlogger/Logger.h>
 
 #include <cmath> // for std::fabs
 #include <deque>
 #include <iostream>
 #include <iterator>
 #include <set> // <<< CHANGED: for dedup sets
 #include <string>
@@ -110,7 +110,7 @@
   ConfigurableAxis configThnAxisPol{"configThnAxisPol", {80, 0.0, 8.0}, "cos#it{#theta *}"};
   ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {8, 0.0, 80.0}, "Centrality"};
   ConfigurableAxis configThnAxisRapidity{"configThnAxisRapidity", {5, 0.0, 1.0}, "Rapidity"};
   ConfigurableAxis configThnAxisPhi{"configThnAxisPhi", {18, 0.0, 2.0 * TMath::Pi()}, "Phi"};
   HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
   void init(o2::framework::InitContext&)
@@ -118,15 +118,15 @@
     histos.add("hPtYSame", "hPtYSame", kTH2F, {{100, 0.0, 10.0}, {200, -1.0, 1.0}});
     histos.add("hPtYMix", "hPtYMix", kTH2F, {{100, 0.0, 10.0}, {200, -1.0, 1.0}});
     histos.add("hCentrality", "Centrality distribution", kTH1F, {{configThnAxisCentrality}});
     histos.add("deltaPhiSame", "deltaPhiSame", HistType::kTH1D, {{72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("deltaPhiMix", "deltaPhiMix", HistType::kTH1D, {{72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("ptCent", "ptCent", HistType::kTH2D, {{100, 0.0, 10.0}, {8, 0.0, 80.0}}, true);
     histos.add("etaCent", "etaCent", HistType::kTH2D, {{32, -0.8, 0.8}, {8, 0.0, 80.0}}, true);
 
     histos.add("hLambdaSameForLL", "hLambdaSameForLL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("hLambdaSameForLAL", "hLambdaSameForLAL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("hLambdaSameForALL", "hLambdaSameForALL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("hAntiLambdaSameForALAL", "hAntiLambdaSameForALAL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
 
     histos.add("hLambdaMixForLL", "hLambdaMixForLL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
     histos.add("hLambdaMixForLAL", "hLambdaMixForLAL", HistType::kTH3D, {{50, 0.0, 5.0}, {32, -0.8, 0.8}, {72, 0.0, 2.0 * TMath::Pi()}}, true);
@@ -138,10 +138,10 @@
     histos.add("hSparseAntiLambdaLambda", "hSparseAntiLambdLambda", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
     histos.add("hSparseAntiLambdaAntiLambda", "hSparseAntiLambdaAntiLambda", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
 
-    histos.add("hSparseLambdaLambdaMixed", "hSparseLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity, configThnAxisR}, true);
-    histos.add("hSparseLambdaAntiLambdaMixed", "hSparseLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity, configThnAxisR}, true);
-    histos.add("hSparseAntiLambdaLambdaMixed", "hSparseAntiLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity, configThnAxisR}, true);
-    histos.add("hSparseAntiLambdaAntiLambdaMixed", "hSparseAntiLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity, configThnAxisR}, true);
+    histos.add("hSparseLambdaLambdaMixed", "hSparseLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
+    histos.add("hSparseLambdaAntiLambdaMixed", "hSparseLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
+    histos.add("hSparseAntiLambdaLambdaMixed", "hSparseAntiLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
+    histos.add("hSparseAntiLambdaAntiLambdaMixed", "hSparseAntiLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true);
 
     histos.add("hSparseRapLambdaLambda", "hSparseRapLambdaLambda", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity}, true);
     histos.add("hSparseRapLambdaAntiLambda", "hSparseRapLambdaAntiLambda", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisRapidity}, true);
@@ -279,7 +279,7 @@
     double deltaPhi = std::abs(RecoDecay::constrainAngle(particle1Dummy.Phi(), 0.0F, harmonic) - RecoDecay::constrainAngle(particle2Dummy.Phi(), 0.0F, harmonic));
     double deltaEta = particle1Dummy.Eta() - particle2Dummy.Eta();
     double deltaRap = std::abs(particle1Dummy.Rapidity() - particle2Dummy.Rapidity());
     double deltaR = TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
 
     if (datatype == 0) {
       mixpairweight = 1.0;
@@ -316,24 +316,27 @@
         weight3 = mixpairweight * hweight3->GetBinContent(hweight3->FindBin(particle1.Pt(), particle1.Eta(), RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic)));
         weight4 = mixpairweight * hweight4->GetBinContent(hweight4->FindBin(particle1.Pt(), particle1.Eta(), RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic)));
       }
-      histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity());
+
       if (tag1 == 0 && tag2 == 0) {
+        histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity(), weight1);
         histos.fill(HIST("hSparseLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight1);
         histos.fill(HIST("hSparseRapLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight1);
         histos.fill(HIST("hSparsePhiLambdaLambdaMixed"), particle1.M(), particle2.M(), costhetaz1costhetaz2, deltaPhi, weight1);
         histos.fill(HIST("hLambdaMixForLL"), particle1.Pt(), particle1.Eta(), RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic), weight1);
       } else if (tag1 == 0 && tag2 == 1) {
+        histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity(), weight2);
         histos.fill(HIST("hSparseLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight2);
         histos.fill(HIST("hSparseRapLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight2);
         histos.fill(HIST("hSparsePhiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), costhetaz1costhetaz2, deltaPhi, weight2);
         histos.fill(HIST("hLambdaMixForLAL"), particle1.Pt(), particle1.Eta(), RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic), weight2);
       } else if (tag1 == 1 && tag2 == 0) {
+        histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity(), weight3);
         histos.fill(HIST("hSparseAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight3);
         histos.fill(HIST("hSparseRapAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight3);
         histos.fill(HIST("hSparsePhiAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), costhetaz1costhetaz2, deltaPhi, weight3);
-
         histos.fill(HIST("hLambdaMixForALL"), particle1.Pt(), particle1.Eta(), RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic), weight3);
       } else if (tag1 == 1 && tag2 == 1) {
+        histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity(), weight4);
         histos.fill(HIST("hSparseAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight4);
         histos.fill(HIST("hSparseRapAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight4);
         histos.fill(HIST("hSparsePhiAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), costhetaz1costhetaz2, deltaPhi, weight4);
@@ -554,6 +557,105 @@
     } // end primary-event loop
   }
   PROCESS_SWITCH(lambdaspincorrderived, processMEV2, "Process data ME", false);
+
+  void processMEV3(EventCandidates const& collisions, AllTrackCandidates const& V0s)
+  {
+    // one pool (deque) per mixing bin; each entry holds (collision index, slice of its V0s)
+    auto nBins = colBinning.getAllBinsCount();
+    std::vector<std::deque<std::pair<int, AllTrackCandidates>>> eventPools(nBins);
+
+    for (auto& collision1 : collisions) {
+      // select mixing bin for this event
+      const int bin = colBinning.getBin(std::make_tuple(collision1.posz(), collision1.cent()));
+
+      // all V0s from the current event
+      auto poolA = V0s.sliceBy(tracksPerCollisionV0, collision1.index());
+
+      // loop over same-event candidate pairs (t1,t2)
+      for (auto& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) {
+        if (!selectionV0(t1) || !selectionV0(t2))
+          continue;
+        if (t2.index() <= t1.index())
+          continue; // unique unordered pairs
+        if (t1.protonIndex() == t2.protonIndex())
+          continue; // no shared daughter
+        if (t1.pionIndex() == t2.pionIndex())
+          continue;
+
+        // --- First pass over previous events in this bin: count replacements and build a list of usable pools
+        int mixes = 0;
+        struct PoolView {
+          AllTrackCandidates* pool;
+          int nRepl;
+          int collIdx;
+        };
+        std::vector<PoolView> usable;
+        int totalRepl = 0;
+
+        for (auto it = eventPools[bin].rbegin();
+             it != eventPools[bin].rend() && mixes < nEvtMixing; ++it, ++mixes) {
+          const int collision2idx = it->first;
+          auto& poolB = it->second;
+
+          // (defensive; shouldn't happen because we push the current event after mixing)
+          if (collision2idx == collision1.index())
+            continue;
+
+          int nRepl = 0;
+          for (auto& t3 : poolB) {
+            if (selectionV0(t3) && checkKinematics(t1, t3))
+              ++nRepl;
+          }
+          if (nRepl > 0) {
+            usable.push_back(PoolView{&poolB, nRepl, collision2idx});
+            totalRepl += nRepl;
+          }
+        }
+
+        if (totalRepl == 0)
+          continue;
+        const float w = 1.0f / static_cast<float>(totalRepl); // global normalization: sum of weights over all replacements = 1
+
+        // --- Second pass: fill with normalized weight w
+        for (auto& pv : usable) {
+          auto& poolB = *pv.pool;
+          for (auto& t3 : poolB) {
+            if (!(selectionV0(t3) && checkKinematics(t1, t3)))
+              continue;
+
+            // build 4-vectors for the mixed pair (t3 from prior event replaces t1; t2 stays from current event)
+            proton = ROOT::Math::PtEtaPhiMVector(t3.protonPt(), t3.protonEta(), t3.protonPhi(), o2::constants::physics::MassProton);
+            lambda = ROOT::Math::PtEtaPhiMVector(t3.lambdaPt(), t3.lambdaEta(), t3.lambdaPhi(), t3.lambdaMass());
+            proton2 = ROOT::Math::PtEtaPhiMVector(t2.protonPt(), t2.protonEta(), t2.protonPhi(), o2::constants::physics::MassProton);
+            lambda2 = ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass());
+
+            float dPhi = std::fabs(
+              RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) -
+              RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic));
+            histos.fill(HIST("deltaPhiMix"), dPhi, w);
+
+            if (t3.v0Status() == 0 && t2.v0Status() == 0) {
+              fillHistograms(0, 0, lambda, lambda2, proton, proton2, 1, w);
+            } else if (t3.v0Status() == 0 && t2.v0Status() == 1) {
+              fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, w);
+            } else if (t3.v0Status() == 1 && t2.v0Status() == 0) {
+              fillHistograms(1, 0, lambda, lambda2, proton, proton2, 1, w);
+            } else if (t3.v0Status() == 1 && t2.v0Status() == 1) {
+              fillHistograms(1, 1, lambda, lambda2, proton, proton2, 1, w);
+            }
+          }
+        }
+      } // end same-event pair loop
+
+      // after mixing with prior events, push current event into the pool
+      auto sliced = V0s.sliceBy(tracksPerCollisionV0, collision1.index());
+      eventPools[bin].emplace_back(collision1.index(), std::move(sliced));
+      if (static_cast<int>(eventPools[bin].size()) > nEvtMixing) {
+        eventPools[bin].pop_front();
+      }
+    } // end primary-event loop
+  }
+  PROCESS_SWITCH(lambdaspincorrderived, processMEV3, "Process data ME", false);
 };
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {