Add tests against analytic expressions for three-body decays

mmikhasenko commented 7 months ago

alpha, beta, gamma = wigner_angle(
   list_of_vectors = ([0,0,p1z,E1],[p2x,0,p2z,E2],[p3x,0,p3z,E3]),
   topology=[[1,2],3],
   ref_topology=[[2,3],1],
   particle=2)

assert alpha == 0
assert gamma == 0
assert beta == zeta_31_for2 # zeta is the angle from DPD

mmikhasenko commented 7 months ago

Here is the expression ref: https://arxiv.org/abs/1910.04566

Here is a code for a cross-check

function cosζ(wr::Arg3WignerRotation, σs, msq)
    i, j, k = ijk(wr)
    # 
    msq[k] ≈ 0 && return one(σs[1])
    # 
    s = msq[4]
    EE4m1sq = (σs[i] - msq[j] - msq[k]) * (σs[j] - msq[k] - msq[i])
    pp4m1sq = sqrt(Kallen(σs[i], msq[j], msq[k]) * Kallen(σs[j], msq[k], msq[i]))
    rest = msq[i] + msq[j] - σs[k]
    return (2msq[k] * rest + EE4m1sq) / pp4m1sq
end

from TBDs.jl

for zeta_31_for2, the i, j, k is equal to 3,1,2, respectively

mmikhasenko commented 7 months ago

To help out further,

let me give you expressions for components,

[0,0,p1z,E1], [p2x,0,p2z,E2], [p3x,0,p3z,E3]

from math import sqrt

# Given values
m1, m2, m3, m0 = 1.1, 1.2, 1.3, 7
m12, m23 = 2.8, 3.3

# Squared masses
m0sq, m1sq, m2sq, m3sq, m12sq, m23sq = [x**2 for x in [m0, m1, m2, m3, m12, m23]]

# Källén function
def Kallen(x, y, z):
    return x**2 + y**2 + z**2 - 2*(x*y + x*z + y*z)

# Calculating missing mass squared using momentum conservation
m31sq = m0sq + m1sq + m2sq + m3sq - m12sq - m23sq

# Momenta magnitudes
p1a = sqrt(Kallen(m23sq, m1sq, m0sq)) / (2*m0)
p2a = sqrt(Kallen(m31sq, m2sq, m0sq)) / (2*m0)
p3a = sqrt(Kallen(m12sq, m3sq, m0sq)) / (2*m0)

# Directions and components
cos_zeta_12_for0 = 1.0 # Placeholder for actual expression
p1z = -p1a
p2z = p2a * cos_zeta_12_for0
p2x = sqrt(p2a**2 - p2z**2)
p3z = -p2z - p1z
p3x = -p2x

# Energy calculations based on the relativistic energy-momentum relation
E1 = sqrt(p1z**2 + m1sq)
E2 = sqrt(p2z**2 + p2x**2 + m2sq)
E3 = sqrt(p3z**2 + p3x**2 + m3sq)

# Vectors
p1 = [0, 0, p1z, E1]
p2 = [p2x, 0, p2z, E2]
p3 = [p3x, 0, p3z, E3]

# Mass squared function for a four-vector
def masssq(p):
    return p[3]**2 - (p[0]**2 + p[1]**2 + p[2]**2)

# Assertions to check calculations
assert round(masssq([p1[i]+p2[i] for i in range(4)]), 5) == round(m12sq, 5)
assert round(masssq([p2[i]+p3[i] for i in range(4)]), 5) == round(m23sq, 5)
assert round(masssq([p3[i]+p1[i] for i in range(4)]), 5) == round(m31sq, 5)

mmikhasenko commented 7 months ago

Happy coding!

from math import sqrt

def ijk(k):
    """Returns the indices based on k, adapted for Python's 0-based indexing."""
    return [(k + 1) % 3, (k + 2) % 3, k]

def cos_zeta_31_for2(msq, sigmas):
    """
    Calculate the cosine of the ζ angle for the case where k=2.

    :param msq: List containing squared masses, with msq[0] being m0ˆ2
    :param sigmas: List containing σ values, adjusted for Python.
    """
    # Adjusted indices for k=2, directly applicable without further modification
    i, j, k = ijk(2)

    s = msq[0]  # s is the first element, acting as a placeholder in this context
    EE4m1sq = (sigmas[i] - msq[j] - msq[k]) * (sigmas[j] - msq[k] - msq[i])
    pp4m1sq = sqrt(Kallen(sigmas[i], msq[j], msq[k]) * Kallen(sigmas[j], msq[k], msq[i]))
    rest = msq[i] + msq[j] - sigmas[k]

    return (2*msq[k] * rest + EE4m1sq) / pp4m1sq

# Example usage:
msq = [7**2, 1.1**2, 1.2**2, 1.3**2]  # m0^2, m1^2, m2^2, m3^2
sigmas = [0.0, 2.8**2, 3.3**2]  # Placeholder for σs[0], σs[1] = m23^2, σs[2] = m31^2

cos_zeta = cos_zeta_31_for2(msq, sigmas)
print(cos_zeta)