[PWGLF] Updated mixing

PWGLF/Tasks/Resonances/phiflowder.cxx

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -17,7 +17,9 @@
 #include "PWGLF/DataModel/LFPhiFlowTables.h"
 
 #include <CommonConstants/PhysicsConstants.h>
+#include <Framework/ASoAHelpers.h>
 #include <Framework/AnalysisTask.h>
+#include <Framework/BinningPolicy.h>
 #include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
 #include <Framework/HistogramSpec.h>
@@ -27,9 +29,11 @@
 
 #include <Math/Vector4D.h>
 
+#include <cmath>
 #include <cstdint>
 #include <map>
 #include <string>
+#include <tuple>
 
 using namespace o2;
 using namespace o2::framework;
@@ -63,6 +67,7 @@
   /*Configurable<float> ptMix{"ptMix", 0.2f, "ME: pT bin width"};
   Configurable<float> etaMix{"etaMix", 0.2f, "ME: eta bin width"};
   Configurable<float> phiMix{"phiMix", 0.3f, "ME: phi bin width"};*/
+  ConfigurableAxis cfgSPAngleBins{"cfgSPAngleBins", {12, 0.0, 2.0 * TMath::Pi()}, "Mixing bins - odd spectator-plane angle"};
 
   HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
@@ -84,6 +89,10 @@
     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);
+    histos.add("hSparseSameUx", "Same-event #Sigma u_{x};M;pT;centrality;eta", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}, true);
+    histos.add("hSparseSameUy", "Same-event #Sigma u_{y};M;pT;centrality;eta", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}, true);
+    histos.add("hSparseSameUQA", "Same-event #Sigma u#upoint Q_{A};M;pT;centrality;eta", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}, true);
+    histos.add("hSparseSameUQC", "Same-event #Sigma u#upoint Q_{C};M;pT;centrality;eta", kTHnSparseF, {axisInvMass, axisPhiPt, axisCent, axisEta, axisV1}, true);
   }
 
   uint8_t getRequiredPidBit() const
@@ -184,7 +193,7 @@
         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);
@@ -282,7 +291,7 @@
         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);
@@ -291,12 +300,18 @@
         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);
+        auto uQA = ux * qxZDCA + uy * qyZDCA;
+        auto uQC = ux * qxZDCC + uy * qyZDCC;
 
         if (phiMass < massMin || phiMass > massMax) {
           continue;
         }
 
         histos.fill(HIST("hSparseSame"), phiMass, phiPt, centrality, etaCand, oddv1);
+        histos.fill(HIST("hSparseSameUx"), phiMass, phiPt, centrality, etaCand, ux);
+        histos.fill(HIST("hSparseSameUy"), phiMass, phiPt, centrality, etaCand, uy);
+        histos.fill(HIST("hSparseSameUQA"), phiMass, phiPt, centrality, etaCand, uQA);
+        histos.fill(HIST("hSparseSameUQC"), phiMass, phiPt, centrality, etaCand, uQC);
         // nPhiSame++;
       }
     }
@@ -466,6 +481,127 @@
   }
 
   PROCESS_SWITCH(phiflowder, processMixedData, "Process mixed-event K+K- pairs", true);
+
+  // Processing Event Mixing
+  void processMixedData2(EventCandidates const& collisions,
+                         aod::KaonTracks const& /*kaontracks*/)
+  {
+    // Calculate the odd spectator-plane angle directly
+    // from the Q vectors already stored in KaonEvents.
+    auto getSPAngle =
+      [](EventCandidates::iterator const& collision) -> float {
+      const float qxOdd =
+        collision.qxA() - collision.qxC();
+
+      const float qyOdd =
+        collision.qyA() - collision.qyC();
+
+      float psiSP =
+        std::atan2(qyOdd, qxOdd);
+
+      // Put first-harmonic angle in [0, 2pi).
+      if (psiSP < 0.f) {
+        psiSP += 2.f * TMath::Pi();
+      }
+
+      return psiSP;
+    };
+
+    using BinningTypeSP =
+      FlexibleBinningPolicy<
+        std::tuple<decltype(getSPAngle)>,
+        aod::kaonevent::Posz,
+        aod::kaonevent::Cent,
+        decltype(getSPAngle)>;
+
+    BinningTypeSP binningOnSPAngle{
+      {getSPAngle},
+      {cfgVtxBins,
+       cfgCentBins,
+       cfgSPAngleBins},
+      true};
+
+    for (const auto& [collision1, collision2] :
+         selfCombinations(binningOnSPAngle,
+                          nEvtMixing.value,
+                          -1,
+                          collisions,
+                          collisions)) {
+
+      if (collision1.globalIndex() ==
+          collision2.globalIndex()) {
+        continue;
+      }
+
+      const float centrality =
+        collision1.cent();
+
+      const float qxZDCA =
+        collision1.qxA();
+
+      const float qyZDCA =
+        collision1.qyA();
+
+      const float qxZDCC =
+        collision1.qxC();
+
+      const float 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;
+
+      // Both combinations use centrality and Q vectors
+      // from collision1.
+      nMixedPairs +=
+        fillMixedPairs(posGroup1,
+                       negGroup2,
+                       centrality,
+                       qxZDCA,
+                       qxZDCC,
+                       qyZDCA,
+                       qyZDCC);
+
+      nMixedPairs +=
+        fillMixedPairs(posGroup2,
+                       negGroup1,
+                       centrality,
+                       qxZDCA,
+                       qxZDCC,
+                       qyZDCA,
+                       qyZDCC);
+
+      histos.fill(HIST("hMixpairs"),
+                  nMixedPairs);
+    }
+  }
+
+  PROCESS_SWITCH(phiflowder, processMixedData2, "Process mixed-event K+K- pairs in SP-angle bins", true);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)