[PWGLF] Add new definition of spin correlation same as STAR to compare

PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx

@@ -30,7 +30,7 @@
 #include <Math/GenVector/Boost.h>
 #include <Math/Vector3D.h>
 #include <Math/Vector4D.h>
 #include <TLorentzVector.h>
 #include <TMath.h>
 
 #include <fairlogger/Logger.h>
@@ -38,7 +38,7 @@
 #include <algorithm>
 #include <cmath> // for std::fabs
 #include <deque>
 #include <iostream>
 #include <iterator>
 #include <set> // <<< CHANGED: for dedup sets
 #include <string>
@@ -112,7 +112,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&)
@@ -120,15 +120,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);
@@ -273,107 +273,129 @@
     auto proton1LambdaRF = boostLambda1ToCM(proton1pairCM);
     auto proton2LambdaRF = boostLambda2ToCM(proton2pairCM);
 
-    TVector3 zhat(lambda1CM.Px(), lambda1CM.Py(), lambda1CM.Pz());
-    if (zhat.Mag2() == 0)
-      zhat = TVector3(0, 0, 1);
-    zhat = zhat.Unit();
-
-    // pick a reference not collinear with ẑ (beam z works)
-    TVector3 ref(0., 0., 1.);
-    if (std::abs(zhat.Dot(ref)) > 0.999)
-      ref = TVector3(1., 0., 0.);
-
-    // build transverse axes
-    TVector3 xhat = (ref - (ref.Dot(zhat)) * zhat).Unit();
-    TVector3 yhat = (zhat.Cross(xhat)).Unit();
-
-    // proton unit vectors in their Λ rest frames
-    TVector3 n1(proton1LambdaRF.Px(), proton1LambdaRF.Py(), proton1LambdaRF.Pz());
-    n1 = n1.Unit();
-    TVector3 n2(proton2LambdaRF.Px(), proton2LambdaRF.Py(), proton2LambdaRF.Pz());
-    n2 = n2.Unit();
-
-    // cosθ* (Λ2 measured along −ẑ)
-    double c1 = n1.Dot(zhat);
-    double c2 = -n2.Dot(zhat);
-
-    // sinθ*
-    double s1 = std::sqrt(std::max(0.0, 1.0 - c1 * c1));
-    double s2 = std::sqrt(std::max(0.0, 1.0 - c2 * c2));
-
-    // φ1*, φ2* in common convention (flip axes for Λ2)
-    double phi1 = std::atan2(n1.Dot(yhat), n1.Dot(xhat));
-    double phi2 = std::atan2(n2.Dot(-yhat), n2.Dot(-xhat));
-
-    // helicity spin-correlation
-    double cosDelta = c1 * c2 + s1 * s2 * std::cos(phi1 - phi2);
-    // clamp for safety
-    if (cosDelta > 1.0)
-      cosDelta = 1.0;
-    if (cosDelta < -1.0)
-      cosDelta = -1.0;
-
-    //
-    // --- (B) Beam-based correlation (common z_B = beam; axes transported into each Λ rest frame) ---
-    //
-    TVector3 zB(0., 0., 1.); // beam axis in pair–CM coordinates
-
-    // x_B: projection of Λ1 direction onto plane ⟂ z_B; y_B completes RH basis
+    // =================== Opening-angle correlator: cos(Δθ) for helicity-z and beam-z ===================
+
+    // Proton unit directions in Λ rest frames
+    TVector3 k1(proton1LambdaRF.Px(), proton1LambdaRF.Py(), proton1LambdaRF.Pz());
+    k1 = k1.Unit();
+    TVector3 k2(proton2LambdaRF.Px(), proton2LambdaRF.Py(), proton2LambdaRF.Pz());
+    k2 = k2.Unit();
+
+    // Helper: boost a spacelike axis (t=0) from PRF into a Λ rest frame
+    auto transport = [](const TVector3& v, const ROOT::Math::Boost& B) -> TVector3 {
+      ROOT::Math::PxPyPzEVector a(v.X(), v.Y(), v.Z(), 0.0);
+      auto ar = B(a);
+      TVector3 out(ar.Px(), ar.Py(), ar.Pz());
+      return (out.Mag2() > 0) ? out.Unit() : out;
+    };
+
+    // ----------------------------- (1) Helicity-z construction -----------------------------
+    // z along Λ1 in PRF
+    TVector3 zPRF(lambda1CM.Px(), lambda1CM.Py(), lambda1CM.Pz());
+    if (zPRF.Mag2() == 0)
+      zPRF = TVector3(0, 0, 1);
+    zPRF = zPRF.Unit();
+
+    // transverse axes in PRF
+    TVector3 ref(0, 0, 1);
+    if (std::abs(zPRF.Dot(ref)) > 0.999)
+      ref = TVector3(1, 0, 0);
+    TVector3 xPRF = (ref - (ref.Dot(zPRF)) * zPRF).Unit();
+    TVector3 yPRF = (zPRF.Cross(xPRF)).Unit();
+
+    // carry PRF triad to Λ rest frames (flip triad for Λ2 to keep same PRF-handedness)
+    TVector3 z1_h = transport(zPRF, boostLambda1ToCM);
+    TVector3 x1_h = transport(xPRF, boostLambda1ToCM);
+    TVector3 y1_h = transport(yPRF, boostLambda1ToCM);
+
+    TVector3 z2_h = transport(-zPRF, boostLambda2ToCM);
+    TVector3 x2_h = transport(-xPRF, boostLambda2ToCM);
+    TVector3 y2_h = transport(-yPRF, boostLambda2ToCM);
+
+    // angles and cosΔθ (helicity)
+    double c1_h = k1.Dot(z1_h);
+    double s1_h = std::sqrt(std::max(0.0, 1.0 - c1_h * c1_h));
+    double phi1_h = std::atan2(k1.Dot(y1_h), k1.Dot(x1_h));
+
+    double c2_h = k2.Dot(z2_h);
+    double s2_h = std::sqrt(std::max(0.0, 1.0 - c2_h * c2_h));
+    double phi2_h = std::atan2(k2.Dot(y2_h), k2.Dot(x2_h));
+
+    double cosDeltaTheta_hel = c1_h * c2_h + s1_h * s2_h * std::cos(phi1_h - phi2_h);
+    if (cosDeltaTheta_hel > 1.0)
+      cosDeltaTheta_hel = 1.0;
+    if (cosDeltaTheta_hel < -1.0)
+      cosDeltaTheta_hel = -1.0;
+
+    // ------------------------------- (2) Beam-z construction -------------------------------
+    // z along beam in PRF; choose x by projecting Λ1 onto the ⟂ plane to fix azimuth zero
+    TVector3 zB(0, 0, 1);
     TVector3 L1dir(lambda1CM.Px(), lambda1CM.Py(), lambda1CM.Pz());
     L1dir = L1dir.Unit();
     TVector3 xB = L1dir - (L1dir.Dot(zB)) * zB;
     if (xB.Mag2() < 1e-12)
-      xB = TVector3(1., 0., 0.); // fallback if Λ1 ∥ beam
+      xB = TVector3(1, 0, 0);
     xB = xB.Unit();
     TVector3 yB = (zB.Cross(xB)).Unit();
 
-    // helper to transport an axis into a Λ rest frame (from pair–CM)
-    auto transport = [](const TVector3& v, const ROOT::Math::Boost& B) -> TVector3 {
-      ROOT::Math::PxPyPzEVector a(v.X(), v.Y(), v.Z(), 0.0); // spacelike 4-vector (t=0)
-      ROOT::Math::PxPyPzEVector ar = B(a);                   // boost into that Λ rest frame
-      TVector3 out(ar.Px(), ar.Py(), ar.Pz());
-      if (out.Mag2() == 0)
-        return out;
-      return out.Unit();
-    };
+    // carry beam triad to Λ rest frames (no flip for a common external axis)
+    TVector3 z1_b = transport(zB, boostLambda1ToCM);
+    TVector3 x1_b = transport(xB, boostLambda1ToCM);
+    TVector3 y1_b = transport(yB, boostLambda1ToCM);
+
+    TVector3 z2_b = transport(zB, boostLambda2ToCM);
+    TVector3 x2_b = transport(xB, boostLambda2ToCM);
+    TVector3 y2_b = transport(yB, boostLambda2ToCM);
+
+    // angles and cosΔθ (beam)
+    double c1_b = k1.Dot(z1_b);
+    double s1_b = std::sqrt(std::max(0.0, 1.0 - c1_b * c1_b));
+    double phi1_b = std::atan2(k1.Dot(y1_b), k1.Dot(x1_b));
+
+    double c2_b = k2.Dot(z2_b);
+    double s2_b = std::sqrt(std::max(0.0, 1.0 - c2_b * c2_b));
+    double phi2_b = std::atan2(k2.Dot(y2_b), k2.Dot(x2_b));
+
+    double cosDeltaTheta_beam = c1_b * c2_b + s1_b * s2_b * std::cos(phi1_b - phi2_b);
+    if (cosDeltaTheta_beam > 1.0)
+      cosDeltaTheta_beam = 1.0;
+    if (cosDeltaTheta_beam < -1.0)
+      cosDeltaTheta_beam = -1.0;
+
+    // --- STAR-style Δθ (as written: dot product of proton directions in their own Λ RFs) ---
+
+    // Boost each proton into its parent's rest frame
+    ROOT::Math::Boost boostL1_LabToRF{particle1Dummy.BoostToCM()}; // Λ1 velocity in lab
+    ROOT::Math::Boost boostL2_LabToRF{particle2Dummy.BoostToCM()}; // Λ2 velocity in lab
+
+    auto p1_LRF = boostL1_LabToRF(daughpart1);
+    auto p2_LRF = boostL2_LabToRF(daughpart2);
+
+    // Unit 3-vectors (in different rest frames!)
+    TVector3 u1 = TVector3(p1_LRF.Px(), p1_LRF.Py(), p1_LRF.Pz()).Unit();
+    TVector3 u2 = TVector3(p2_LRF.Px(), p2_LRF.Py(), p2_LRF.Pz()).Unit();
 
-    // transport beam triad to each Λ rest frame
-    TVector3 zB1 = transport(zB, boostLambda1ToCM);
-    TVector3 xB1 = transport(xB, boostLambda1ToCM);
-    TVector3 yB1 = transport(yB, boostLambda1ToCM);
-
-    TVector3 zB2 = transport(zB, boostLambda2ToCM);
-    TVector3 xB2 = transport(xB, boostLambda2ToCM);
-    TVector3 yB2 = transport(yB, boostLambda2ToCM);
-
-    // angles w.r.t. beam triad (same z_B convention for both Λ’s; no flips)
-    double c1B = n1.Dot(zB1), c2B = n2.Dot(zB2);
-    double s1B = std::sqrt(std::max(0.0, 1.0 - c1B * c1B));
-    double s2B = std::sqrt(std::max(0.0, 1.0 - c2B * c2B));
-    double phi1B = std::atan2(n1.Dot(yB1), n1.Dot(xB1));
-    double phi2B = std::atan2(n2.Dot(yB2), n2.Dot(xB2));
-
-    // beam correlation
-    double cosDeltaBeam = c1B * c2B + s1B * s2B * std::cos(phi1B - phi2B);
-    if (cosDeltaBeam > 1.0)
-      cosDeltaBeam = 1.0;
-    if (cosDeltaBeam < -1.0)
-      cosDeltaBeam = -1.0;
+    // STAR-style cosΔθ definition
+    double cosDeltaTheta_STAR_naive = u1.Dot(u2);
+    if (cosDeltaTheta_STAR_naive > 1.0)
+      cosDeltaTheta_STAR_naive = 1.0;
+    if (cosDeltaTheta_STAR_naive < -1.0)
+      cosDeltaTheta_STAR_naive = -1.0;
 
     auto cosThetaDiff = -999.0;
     auto costhetaz1costhetaz2 = -999.0;
     if (cosDef == 0) {
-      cosThetaDiff = proton1LambdaRF.Vect().Unit().Dot(proton2LambdaRF.Vect().Unit());
+      cosThetaDiff = cosDeltaTheta_STAR_naive;
       costhetaz1costhetaz2 = (proton1LambdaRF.Pz() * proton2LambdaRF.Pz()) / (proton1LambdaRF.P() * proton2LambdaRF.P());
     } else {
-      cosThetaDiff = cosDelta;
-      costhetaz1costhetaz2 = cosDeltaBeam;
+      cosThetaDiff = cosDeltaTheta_hel;
+      costhetaz1costhetaz2 = cosDeltaTheta_beam;
     }
 
     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;