diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 7382a32c860..daeefca4adc 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -273,101 +273,123 @@ struct lambdaspincorrderived { 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));