Example for Numeric Amplitude

This example give a full formula and code for numeric results of a total amplitude. Here we use \(\Lambda_c^{+} \rightarrow \Lambda (\rightarrow p \pi^{-}) \pi^{+} \pi^{0}\) as an example to show how its amplitudes are calculated.

0. The inputs

A total amplitude requires two inputs, the ampulitude parameters and data. In TFPWA, amplitude parameters include the coupling, mass, widths and so on. And the data include the 4-momenta \(p^{\mu}\), and some addtional datas, such as charge and weight. In the simple example, we only consider \(p^{\mu}\).

In this example, we will calculate the amplitude with the following data

import cmath
import itertools
import math
from pprint import pprint

m_lc, m_l, m_p, m_pi, m_pi0 = (
    2.28646,
    1.115683,
    0.93827208816,
    0.13957039,
    0.1349768,
)
m_sigma_1660, g_sigma_1660 = 1.665, 0.3
m_rho770, g_rho770 = 0.77511, 0.1491

# the example momentum is generated by the follwing code
if False:
    from tf_pwa import set_random_seed
    from tf_pwa.phasespace import generate_phsp

    set_random_seed(2025)
    ((p_p_tf, p_pim_tf), p_pip_tf, p_pi0_tf), _ = generate_phsp(
        4.6, [(m_lc, [(m_l, [m_p, m_pi]), m_pi, m_pi0]), m_lc], N=1
    )
    [p_p, p_pim, p_pip, p_pi0] = [
        i[0].numpy().tolist() for i in [p_p_tf, p_pim_tf, p_pip_tf, p_pi0_tf]
    ]

p_p = [
    0.9798045738146876,
    0.023972732547648895,
    0.16664847395929822,
    -0.22653054025735844,
]
p_pim = [
    0.19300790258521403,
    0.12035810723629733,
    0.004987489309569788,
    -0.057106984410615395,
]
p_pip = [
    0.8227454084228714,
    -0.06451285582391393,
    -0.5452445877195563,
    0.5966376993721896,
]
p_pi0 = [
    0.3044421333949209,
    -0.11013601413926419,
    0.15645293575957042,
    -0.19457341372741113,
]

print("p", p_p)
print("pim", p_pim)
print("pip", p_pip)
print("pi0", p_pi0)

p [0.9798045738146876, 0.023972732547648895, 0.16664847395929822, -0.22653054025735844]
pim [0.19300790258521403, 0.12035810723629733, 0.004987489309569788, -0.057106984410615395]
pip [0.8227454084228714, -0.06451285582391393, -0.5452445877195563, 0.5966376993721896]
pi0 [0.3044421333949209, -0.11013601413926419, 0.15645293575957042, -0.19457341372741113]

The example using the following config and parameters

config_str = f"""
data:
   dat_order: [p, pim, pip, pi0]
   random_z: False
decay:
   Lc:
   - [Sigma, pi0, p_break: True, barrier_factor_norm: True]
   - [rho, Lambda, p_break: True, barrier_factor_norm: True]
   Sigma: [Lambda, pip, barrier_factor_norm: True]
   rho: [pip, pi0, barrier_factor_norm: True]
   Lambda: [p, pim, p_break: True, barrier_factor_norm: True]

particle:
    $top:
        Lc:
            J: 1/2
            P: +1
            mass: {m_lc}
    $finals:
        p:
            J: 1/2
            P: +1
            mass: {m_p}
        pip:
            J: 0
            P: -1
            mass: {m_pi}
        pim:
            J: 0
            P: -1
            mass: {m_pi}
        pi0:
            J: 0
            P: -1
            mass: {m_pi0}
    Lambda:
        J: 1/2
        P: +1
        mass: {m_l}
        model: one
    Sigma:
        J: 1/2
        P: +1
        mass: {m_sigma_1660}
        width: {g_sigma_1660}
        model: BW
    rho:
        J: 1
        P: -1
        mass: {m_rho770}
        width: {g_rho770}
        model: BW
"""

params = {
    "Lc->Sigma.pi0Sigma->Lambda.pipLambda->p.pim_total_0r": 1.0,
    "Lc->Sigma.pi0Sigma->Lambda.pipLambda->p.pim_total_0i": 0.0,
    "Lc->Sigma.pi0_g_ls_0r": 1.0,
    "Lc->Sigma.pi0_g_ls_0i": 0.0,
    "Lc->Sigma.pi0_g_ls_1r": 1.0,
    "Lc->Sigma.pi0_g_ls_1i": -1.0,
    "Sigma->Lambda.pip_g_ls_0r": 1.0,
    "Sigma->Lambda.pip_g_ls_0i": 0.0,
    "Lambda->p.pim_g_ls_0r": 1.0,
    "Lambda->p.pim_g_ls_0i": 0.0,
    "Lambda->p.pim_g_ls_1r": 0.435376,
    "Lambda->p.pim_g_ls_1i": 0.0,
    "Lc->rho.Lambdarho->pip.pi0Lambda->p.pim_total_0r": 2.0,
    "Lc->rho.Lambdarho->pip.pi0Lambda->p.pim_total_0i": 1.0,
    "Lc->rho.Lambda_g_ls_0r": 1.0,
    "Lc->rho.Lambda_g_ls_0i": 0.0,
    "Lc->rho.Lambda_g_ls_1r": 1.0,
    "Lc->rho.Lambda_g_ls_1i": 0.0,
    "Lc->rho.Lambda_g_ls_2r": 1.0,
    "Lc->rho.Lambda_g_ls_2i": 0.0,
    "Lc->rho.Lambda_g_ls_3r": 1.0,
    "Lc->rho.Lambda_g_ls_3i": 0.0,
    "rho->pip.pi0_g_ls_0r": 1.0,
    "rho->pip.pi0_g_ls_0i": 0.0,
}

Here are some basic functions for momentum and angles

def M(p):
    return math.sqrt(p[0] * p[0] - p[1] * p[1] - p[2] * p[2] - p[3] * p[3])


def dot(p1, p2):
    return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]


def boost(p, v):
    E, px, py, pz = p
    vx, vy, vz = v
    beta2 = vx * vx + vy * vy + vz * vz
    gamma = 1.0 / math.sqrt(1 - beta2)
    bp = dot([px, py, pz], [vx, vy, vz])
    gamma2 = 1 / (math.sqrt(1 - beta2) + 1 - beta2)  # (gamma - 1) / beta2
    r = gamma2 * bp + gamma * E
    px_prime = px + r * vx
    py_prime = py + r * vy
    pz_prime = pz + r * vz
    E_prime = gamma * (E + bp)
    return [E_prime, px_prime, py_prime, pz_prime]


def restframe(base, other):
    E, px, py, pz = base
    p = -px / E, -py / E, -pz / E
    return boost(other, p)


def add(*args):  # (E1+E2, x1+x2, y1+y2, z1+z2)
    return [sum(i) for i in zip(*args)]


def unit(p):
    n = math.sqrt(dot(p, p))
    return [i / n for i in p]


def cross_unit(p1, p2):  # p1 x p2 = (y1z2 - y2z1, z1x2 - z2x1, x1y2 - x2y1)
    px = p1[1] * p2[2] - p2[1] * p1[2]
    py = p1[2] * p2[0] - p2[2] * p1[0]
    pz = p1[0] * p2[1] - p2[0] * p1[1]
    return unit([px, py, pz])


def angle_phi(xr, x, y):  # arg(xr.x +i xr.y)
    costheta = dot(xr, x)
    sintheta = dot(xr, y)
    return math.atan2(sintheta, costheta)


def angle_theta(p1, p2):
    costheta = dot(p1, p2)
    return math.acos(costheta)


def angle_zx_p_getzx(z0, x0, p):
    z1 = unit(p[1:])
    y1 = cross_unit(z0, z1)
    y0 = cross_unit(z0, x0)
    x1 = cross_unit(y1, z1)
    phi = angle_phi(y1, y0, [-i for i in x0])
    theta = angle_theta(z0, z1)
    return phi, theta, z1, x1


def angle_zx_zx(z0, x0, z1, x1):
    yr = cross_unit(z0, z1)
    y0 = cross_unit(z0, x0)
    y1 = cross_unit(z1, x1)
    phi = angle_phi(yr, y0, [-i for i in x0])
    theta = angle_theta(z0, z1)
    psi = -angle_phi(yr, y1, [-i for i in x1])
    return phi, theta, psi

1. Calculate the angle

Before calculate the amplitude, we calculate the angle in the data based on the 4-momenta and decay structure.

1.1 First Decay chain

The first decay chain is \(\Lambda_c^{+} \rightarrow \Sigma^{+} \pi^{0},\Sigma^{+} \rightarrow \Lambda \pi^{+}, \Lambda \rightarrow p \pi^{-}\).

Firstly we boost all the momentum into the rest frame of \(\Lambda_c^{+}\).

p_total = add(p_p, p_pim, p_pip, p_pi0)
p1_p, p1_pim, p1_pip, p1_pi0 = [
    restframe(p_total, p) for p in [p_p, p_pim, p_pip, p_pi0]
]
pprint([p1_p, p1_pim, p1_pip, p1_pi0])

[[1.013485155460205,
  0.037149042316438785,
  0.2610250084944836,
  -0.2779991771813273],
 [0.19917831654329266,
  0.12295058891840963,
  0.02355637857348515,
  -0.06723360565511058],
 [0.7440749913533434,
  -0.054155650528414465,
  -0.4710601491012403,
  0.556180906214235],
 [0.32972155378779994,
  -0.10594398070643395,
  0.18647876203327152,
  -0.2109481233777971]]

Then we choose the initial coordinate system as \(\vec{x}_{0}=(1,0,0)\), \(\vec{z}_{0}=(0,0,1)\), \(\vec{y}_{0} = \vec{z}_{0}\times \vec{x}_{0}=(0,1,0)\), all the angle will be calculated based on it.

z0 = [0.0, 0.0, 1.0]
y0 = [0.0, 1.0, 0.0]
x0 = [1.0, 0.0, 0.0]

In the first decay \(\Lambda_c^{+} \rightarrow \Sigma^{+} \pi^{0}\), we choose the direction of \(\Sigma^{+}\) as the \(z\)-axis of new coordinate system (\(\vec{z}_1 = \frac{\vec{p}_{\Sigma}}{|\vec{p}_{\Sigma}|}\)). In the convention we used, , the third Euler angle is zero in the formula. We have to choose \(\vec{y} = \frac{\vec{z}_{0} \times \vec{z}_{1}}{|\vec{z}_{0} \times \vec{z}_{1}|}\), which is the rotation axis of second rotaion that rotate \(\vec{z}_{0}\) to \(\vec{z}_{1}\). The angle can be calulated by the rotation of axis.

p1_sigma = add(p1_p, p1_pim, p1_pip)
z1 = unit(p1_sigma[1:])  # using p1 as z axis
y1 = cross_unit(z0, z1)  # the rotation axis is z0 x z1
x1 = cross_unit(y1, z1)  # the x-axis of new coordinate system
xr = cross_unit(y1, z0)  # the x-axis after the first rotation

phi1_test = angle_phi(
    xr, x0, y0
)  # the angle is calculate based on new x from axis (x, y)
# or without xr
phi1 = angle_phi(
    y1, y0, [-i for i in x0]
)  # the angle is calculate based on new y from axis (y, -x)

theta1 = angle_theta(z0, z1)  # theta is arccos( z0 . z1 )

assert phi1_test == phi1
print(phi1, theta1)

-1.0541411652888346 0.7936824390985183

In the second decay \(\Sigma^{+} \rightarrow \Lambda \pi^{+}\), the initial coordinate system inherits from the above decay that produce \(\Sigma^{+}\). We choose the direction of \(\Lambda^{+}\) as the \(z\)-axis of new coordinate system (\(\vec{z}_2 = \frac{\vec{p}_{\Lambda}}{|\vec{p}_{\Lambda}|}\)). And also \(\vec{y} = \frac{\vec{z}_{1} \times \vec{z}_{2}}{|\vec{z}_{1} \times \vec{z}_{2}|}\). The angle calulation is similar, we can use a common function angle_zx_p_getzx.

p2_p, p2_pim, p2_pip = [restframe(p1_sigma, p) for p in [p1_p, p1_pim, p1_pip]]
p2_lambda = add(p2_p, p2_pim)

phi2, theta2, z2, x2 = angle_zx_p_getzx(z1, x1, p2_lambda)

print(phi2, theta2)

1.141488741005533 2.635475098581215

In the last decay \(\Lambda^{+} \rightarrow p \pi^{-}\), the initial coordinate system inherits from the decay that produce \(\Lambda^{+}\). We choose the direction of \(p\) as the \(z\)-axis of new coordinate system (\(\vec{z}_3 = \frac{\vec{p}_{p}}{|\vec{p}_{p}|}\)). And also \(\vec{y} = \frac{\vec{z}_{3} \times \vec{z}_{3}}{|\vec{z}_{3} \times \vec{z}_{3}|}\).

p3_p, p3_pim = [restframe(p2_lambda, p) for p in [p2_p, p2_pim]]
phi3, theta3, z3, x3 = angle_zx_p_getzx(z2, x2, p3_p)

print(phi3, theta3)

0.2005695899392871 1.6849581714082897

1.2 Second decay Chain

In the second decay chain \(\Lambda_c^{+} \rightarrow \rho \Lambda, \rho\rightarrow \pi^{+}\pi^{0}, \Lambda\rightarrow p \pi^{-}\), the procedure is similar. In the first decay \(\Lambda_c^{+} \rightarrow \rho \Lambda\), we choose the direction of \(\rho\) as the \(z\)-axis of new coordinate system (\(\vec{z}_1' = \frac{\vec{p}_{\rho}}{|\vec{p}_{\rho}|}\)), and also \(\vec{y} = \frac{\vec{z}_{0} \times \vec{z}_{1}}{|\vec{z}_{0} \times \vec{z}_{1}|}\).

p_total = add(p_p, p_pim, p_pip, p_pi0)
p1_p, p1_pim, p1_pip, p1_pi0 = [
    restframe(p_total, p) for p in [p_p, p_pim, p_pip, p_pi0]
]
p1_rho_2 = add(p1_pip, p1_pi0)

phi1_2, theta1_2, z1_2, x1_2 = angle_zx_p_getzx(z0, x0, p1_rho_2)

print(phi1_2, theta1_2)

-2.083246114791286 0.7575560832474785

In the second decay \(\rho^{+} \rightarrow \pi^{+} \pi^{0}\), the initial coordinate system inherits from the above decay that produce \(\rho^{+}\). We choose the direction of \(\pi^{+}\) as the \(z\)-axis of new coordinate system (\(\vec{z}_2' = \frac{\vec{p}_{\pi^{+}} }{|\vec{p}_{\pi^{+}}|}\)). And also \(\vec{y}_1' = \frac{\vec{z}_{1}' \times \vec{z}_{2}'}{|\vec{z}_{1}' \times \vec{z}_{2}'|}\).

p2_pip_2, p2_pi0_2 = [restframe(p1_rho_2, p) for p in [p1_pip, p1_pi0]]
phi2_2, theta2_2, z2_2, x2_2 = angle_zx_p_getzx(z1_2, x1_2, p2_pip_2)
print(phi2_2, theta2_2)

1.9089059699896964 0.43577968070080797

There are some differences of \(\Lambda \rightarrow p \pi^{-}\). Firstly the initial coordinate system inherits from the decay that produce \(\Lambda\), which is the first decay, no the second decay. And then the direction of \(\Lambda\) is the opposite of \(\rho^{+}\). We need to consider the transform to the opposite aixs. Here we simplely choose \((\vec{x}_{1}', -\vec{y}_{1}',-\vec{z}_{1}')\) as the opposite of \((\vec{x}_{1}', \vec{y}_{1}',\vec{z}_{1}')\). It means a \(\pi\) rotation around x-axis. Then we can calculate the angle using the same procedure.

p1_lambda = add(p1_p, p1_pim)
p3_p_2, p3_pim_2 = [restframe(p1_lambda, p) for p in [p1_p, p1_pim]]

#   x1', -y1', -z1'
o_z1_2 = [-i for i in z1_2]
o_x1_2 = x1_2

phi3_2, theta3_2, z3_2, x3_2 = angle_zx_p_getzx(o_z1_2, o_x1_2, p3_p_2)
print(phi3_2, theta3_2)

1.4406467063875559 1.8593100836742338

1.3 Alignment angle

From the above processes, we can see that the final coordinate system of particle \(p\) (and also other final particles) is not the same.

print("x")
print(x3_2)
print(x3)
print("z")
print(z3_2)
print(z3)

x
[-0.06395677311665698, -0.6890741879517143, 0.7218630719680776]
[-0.034929150029978914, -0.6644561028949505, 0.7465105771546967]
z
[-0.9686396455654857, 0.21689353360891328, 0.12122059279433339]
[-0.9694191132638728, 0.20408496607719684, 0.13629346814853677]

We need to add additonal rotation to align them. In addtion, we also need to consider the effect of boost sequence in it. The sequence boost will introduce an additional rotation. In TFPWA, we record all the boost and rotation in the decay chain, and get a total Lorentz transform for a particle. The Lorentz transform using the following 2-D matrix form.

\[\begin{split}R_z(\alpha) = \begin{pmatrix} e^{-i\frac{\alpha}{2}} & 0 \\ 0 & e^{i\frac{\alpha}{2}} \end{pmatrix}, R_y(\beta) = \begin{pmatrix} \cos\frac{\beta}{2} & -\sin\frac{\beta}{2} \\ \sin\frac{\beta}{2} & \cos\frac{\beta}{2} \end{pmatrix}, B_z(\omega) = \begin{pmatrix} e^{-\frac{\omega}{2}} & 0 \\ 0 & e^{\frac{\omega}{2}} \end{pmatrix}\end{split}\]

Here, \(\omega=\tanh^{-1}(|p|/E)\). Base on the 2-D matrix form, we can solve the angle from \(L=R_z(\gamma)R_y(\beta)R_z(\alpha)B_z(\omega)\) as

\[\begin{split}\begin{cases} \cos\beta = & L_{11} L_{22} + L_{12} L_{21} \\ \alpha + \gamma =& 2\arg(L_{22}) \\ \alpha - \gamma =& 2\arg(L_{21}) \\ \end{cases}\end{split}\]

Here is some code for this form.

class LorentzGroup:
    def __init__(self, x):
        self.x = x

    def __mul__(self, other):  # matrix multiply
        x, y = self.x, other.x
        aa = x[0][0] * y[0][0] + x[0][1] * y[1][0]
        ab = x[0][0] * y[0][1] + x[0][1] * y[1][1]
        ba = x[1][0] * y[0][0] + x[1][1] * y[1][0]
        bb = x[1][0] * y[0][1] + x[1][1] * y[1][1]
        return LorentzGroup([[aa, ab], [ba, bb]])

    def inv(self):  # matrix inverse, it has det(M) =1
        aa, ab = self.x[0]
        ba, bb = self.x[1]
        return LorentzGroup([[bb, -ab], [-ba, aa]])

    def get_euler_angle(self):  # solve euler angle from the matrix
        x = self.x
        cosbeta = (x[0][0] * x[1][1] + x[0][1] * x[1][0]).real
        cosbeta = min(
            1.0, cosbeta
        )  # due to the pricision it could be 1 + epsilon
        beta = math.acos(cosbeta)
        alpha_p_gamma = math.atan2(x[1][1].imag, x[1][1].real)
        alpha_m_gamma = math.atan2(x[1][0].imag, x[1][0].real)
        alpha = alpha_p_gamma + alpha_m_gamma
        gamma = alpha_p_gamma - alpha_m_gamma
        return alpha, beta, gamma

    def __repr__(self):
        return str(self.x)


def Bz(omega):
    a = math.exp(omega / 2)
    return LorentzGroup([[1 / a, 0.0], [0.0, a]])


def Bz_from_p(p4):
    E, px, py, pz = p4
    p = math.sqrt(px * px + py * py + pz * pz)
    omega = math.atanh(p / E)
    return Bz(omega)


def Rz(alpha):
    a = cmath.exp(1j * alpha / 2)
    return LorentzGroup([[1 / a, 0.0], [0.0, a]])


def Ry(beta):
    s = math.sin(beta / 2)
    c = math.cos(beta / 2)
    return LorentzGroup([[c, -s], [s, c]])

For the first decay chain, \(\Lambda_c^{+} \rightarrow \Sigma^{+} \pi^{0},\Sigma^{+} \rightarrow \Lambda \pi^{+}, \Lambda \rightarrow p \pi^{-}\), the total Lorentz transform of \(p\) is

\[L_{\Sigma} = R_y(\theta_3)R_z(\phi_3) B_z(\omega_{3}) R_y(\theta_2)R_z(\phi_2) B_z(\omega_{2}) R_y(\theta_1)R_z(\phi_1) B_z(\omega_{1})\]

we can remove the first \(B_z(\omega_{1})\), since it is same for all decay chains.

L1 = Ry(theta1) * Rz(phi1) * Bz_from_p(p_total)
L2 = Ry(theta2) * Rz(phi2) * Bz_from_p(p1_sigma)
L3 = Ry(theta3) * Rz(phi3) * Bz_from_p(p2_lambda)
L_Sigma = L3 * L2 * L1

For the second decay chain, \(\Lambda_c^{+} \rightarrow \rho \Lambda, \rho\rightarrow \pi^{+}\pi^{0}, \Lambda\rightarrow p \pi^{-}\), the total Lorentz transform of \(p\) is

\[L_{\rho} = R_y(\theta_{3}')R_z(\phi_{3}') B_z(\omega_{3}') {\color{ead}R_{x}(\pi)} R_y(\theta_1') R_z(\phi_1') B_z(\omega_{1})\]

There is an additional retation \({\color{ead}R_{x}(\pi)}\), since \(\Lambda\) is in the opposite direction of \(\rho\). Here we choose \(R_{x}(\pi) R_y(\theta_1') R_z(\phi_1') = R_y(\pi-\theta_1') R_z(\phi_1'{\color{red}-\pi})\), so that

\[L_{\rho} = R_y(\theta_{3}')R_z(\phi_{3}') B_z(\omega_{3}') R_y(\pi-\theta_1') R_z(\phi_1'-\pi) B_z(\omega_{1}')\]

The \({\color{red}-\pi}\) can also be \({\color{red}+\pi}\). But for baryon, rotation of \(-\pi\) and \(+\pi\) have different factor \(i\) and \(-i\), we need to choose one of them in the conversion.

L1_2 = Ry(math.pi - theta1_2) * Rz(phi1_2 - math.pi) * Bz_from_p(p_total)
L3_2 = Ry(theta3_2) * Rz(phi3_2) * Bz_from_p(p1_lambda)
L_rho = L3_2 * L1_2

The relative Lorentz transform between the two coordinate systems is

\[L_{\text{align}} = L_{\Sigma} L^{-1}_{\rho}\]

We can solve the angle form it

L_p = L_Sigma * L_rho.inv()

ang = L_p.get_euler_angle()
print(ang)
print((ang[0] + ang[2]) % (math.pi * 2))

(-4.230335631198603, 0.0, -2.016038195610422)
0.036811480370561256

To inclue the effect \({\color{red}-\pi}\), the range of \(\alpha\) and \(\gamma\) is \((-2\pi,2\pi)\). Here is the same rest frame of \(p\), \(\beta\) is always zero. But it has a non-zero phase \(\alpha + \gamma\). If we do not consider the effect of Lorentz boost, we can use the coordinate system to calculate the angle directly. The value would be a little different, and do not have the property that \(\beta=0\).

ang_prime = angle_zx_zx(z3_2, x3_2, z3, x3)
print(ang_prime)

(0.061965116196543796, 0.019795740458797503, -0.021156236313810133)

For \(\pi^{\pm,0}\), we can also calculate the angle for alignment. Each particle require independent alignment angles. It could have no-zere \(\beta\). But their spins is 0, no effect will appear. We can skip those steps.

2. Calculate the amplitude

2.1 First decay chain

The total amplitude of decay chain include \(\Sigma(1660)^{+}\) is expanded by the three decays in the decay chain The first one is \(\Lambda_c^{+}\rightarrow \Sigma(1660)^{+} \pi^0\). We start from the LS amplitude \(G_{ls} = g_{ls} q^{l}B_l'(q,q0,d)\). \(g_{ls}\) is the fit parameters.

gls0 = params["Lc->Sigma.pi0_g_ls_0r"] * cmath.exp(
    1j * params["Lc->Sigma.pi0_g_ls_0i"]
)
gls1 = params["Lc->Sigma.pi0_g_ls_1r"] * cmath.exp(
    1j * params["Lc->Sigma.pi0_g_ls_1i"]
)

\(q\) is the momentum in the restframe

\[q = \frac{\sqrt{(m^2 - (m_1+m_2)^2 )(m^2 - (m_1-m_2)^2 )}}{2m}\]

def get_relative_p(m0, m1, m2):
    return (
        math.sqrt((m0**2 - (m1 + m2) ** 2) * (m0**2 - (m1 - m2) ** 2)) / 2 / m0
    )

\(q\) is the value in data. And \(q_0\) is the value in model, as a normalize factor.

p_lambda = add(p_p, p_pim)
p_sigma = add(p_lambda, p_pip)
p_rho = add(p_pip, p_pi0)

q = get_relative_p(M(p_total), M(p_sigma), M(p_pi0))
q0 = get_relative_p(m_lc, m_sigma_1660, m_pi0)

and \(B_l\) is the Blatt-Weisskopf barrier factors.

\[\begin{split}B_l'(q,q0,d) = \begin{cases} 1 & \text{ for } l=0 \\ \sqrt{\frac{1+(q_0 d)^2}{1+(qd)^2}} & \text{ for } l= 1 \\ \sqrt{\frac{9+3(q_0 d)^2+(q_0 d)^4}{9+3(qd)^2+(q d)^4}} & \text{ for } l= 2 \\ \end{cases}\end{split}\]

def Bprime(l, q, q0, d=3.0):
    if l == 0:
        return 1.0
    if l == 1:
        return math.sqrt((1 + (q0 * d) ** 2) / (1 + (q * d) ** 2))
    if l == 2:
        return math.sqrt(
            (9 + 3 * (q0 * d) ** 2 + (q0 * d) ** 4)
            / (9 + 3 * (q * d) ** 2 + (q * d) ** 4)
        )

For this process, we add barrier_factor_norm: True option in config, it becomes \(G_{ls} = g_{ls}\left(\frac{q}{q_0}\right)^{l}B_l'(q,q0,d)\).

G0 = gls0 * (q / q0) ** 0 * Bprime(0, q, q0)
G1 = gls1 * (q / q0) ** 1 * Bprime(1, q, q0)

Next, we calculate the helicity amplitude, it is

\[H_{\lambda_{1},\lambda_{2} } = \sum_{l,s} G_{l,s} \sqrt{\frac{2l+1}{2J_0+1}} \langle J_{1}, \lambda_{1}; J_{2}, -\lambda_{2}| s, \delta \rangle \langle l,0;s,\delta | J_{0}, \delta\rangle\]

\(\delta=\lambda_{1}-\lambda_{2}\). \(l,s\) is all possible in the range

\[|J_{1} - J_{2} | \leq s \leq J_{1} + J_{2}, |l -s | \leq J_{0} \leq l + s\]

For the first decay, \(\Lambda_c^{+}\rightarrow \Sigma(1665)^{+} \pi^0\). \(s\) in range \([|\frac{1}{2}-0|, \frac{1}{2}+0]\), so it is only \(\frac{1}{2}\). And \(l\) in range \([|\frac{1}{2}-\frac{1}{2}|, \frac{1}{2}+ \frac{1}{2}]\), could be \(0\) or \(1\). All possible \((l,s)\) is \((0, \frac{1}{2})\) and \((1, \frac{1}{2})\).

For \((l,s)=(0, \frac{1}{2})\), we have the following Clebsch–Gordan coefficients,

\[\langle \frac{1}{2}, -\frac{1}{2}; 0, 0| \frac{1}{2}, -\frac{1}{2}-0 \rangle=1, \langle 0,0;\frac{1}{2}, - \frac{1}{2}-0 | \frac{1}{2}, -\frac{1}{2} - 0\rangle = 1\]

\[\langle \frac{1}{2}, \frac{1}{2}; 0, 0| \frac{1}{2}, \frac{1}{2}-0 \rangle=1, \langle 0,0;\frac{1}{2}, \frac{1}{2}-0 | \frac{1}{2}, \frac{1}{2} - 0\rangle = 1\]

For \((l,s)=(1, \frac{1}{2})\), we have the following Clebsch–Gordan coefficients

\[\langle \frac{1}{2}, -\frac{1}{2}; 0, 0| \frac{1}{2}, -\frac{1}{2}-0 \rangle=1, \langle 1,0;\frac{1}{2}, - \frac{1}{2}-0 | \frac{1}{2}, -\frac{1}{2} - 0\rangle = \frac{\sqrt{3}}{3}\]

\[\langle \frac{1}{2}, \frac{1}{2}; 0, 0| \frac{1}{2}, \frac{1}{2}-0 \rangle=1, \langle 1,0;\frac{1}{2}, \frac{1}{2}-0 | \frac{1}{2}, \frac{1}{2} - 0\rangle = -\frac{\sqrt{3}}{3}\]

Combine with the fractor \(\sqrt{\frac{2l+1}{2J_0 + 1}}\), we can get that

\[H_{-\frac{1}{2},0} = \frac{1}{\sqrt{2}}(G_{0,\frac{1}{2}} + G_{1,\frac{1}{2}}), H_{\frac{1}{2},0} = \frac{1}{\sqrt{2}}( G_{0,\frac{1}{2}} - G_{1,\frac{1}{2}})\]

H1 = {
    (-1 / 2, 0): (G0 + G1) / math.sqrt(2),
    (1 / 2, 0): (G0 - G1) / math.sqrt(2),
}

The amplitude of this decay is also related to function

\[D_{m,n}^{J_0}(\alpha,\beta,\gamma) = e^{-im \alpha} d_{m,n}^{J}(\beta) e^{-in\gamma}.\]

The small d matrix \(d^{J}_{m,n}(\beta)\) is \(d^{0}_{m,n}(\beta)=1\) and

\[\begin{split}d_{m,n}^{1/2}(\beta) = \begin{pmatrix} \cos(\frac{\beta}{2}) & - \sin(\frac{\beta}{2}) \\ \sin(\frac{\beta}{2}) & \cos(\frac{\beta}{2}) \end{pmatrix}\end{split}\]

\[\begin{split}d_{m,n}^{1}(\beta) = \begin{pmatrix} \frac{1}{2}[1+\cos(\beta)] & \frac{1}{\sqrt{2}}\sin(\beta) & \frac{1}{2}[1-\cos(\beta)] \\ -\frac{1}{\sqrt{2}}\sin(\beta) & \cos(\beta) & \frac{1}{\sqrt{2}}\sin(\beta) \\ \frac{1}{2}[1-\cos(\beta)] & -\frac{1}{\sqrt{2}}\sin(\beta) & \frac{1}{2}[1+\cos(\beta)] \\ \end{pmatrix}\end{split}\]

here \(m,n = [-J, -J+1,\cdots, J]\). And we use the complex conjugation version.

def small_d_matrix(beta, j):
    if j == 0:
        return [[1.0]]
    if j == 1:
        return [
            [math.cos(beta / 2), math.sin(beta / 2)],
            [-math.sin(beta / 2), math.cos(beta / 2)],
        ]
    if j == 2:
        c = math.cos(beta)
        s = math.sin(beta)
        k = 1 / math.sqrt(2)
        return [
            [(1 + c) / 2, k * s, (1 - c) / 2],
            [-k * s, c, k * s],
            [(1 - c) / 2, -k * s, (1 + c) / 2],
        ]


def D_matrix_conj(alpha, beta, gamma, j):
    d = small_d_matrix(beta, j)
    ret = []
    for i in range(int(j + 1)):
        tmp = []
        for k in range(int(j + 1)):
            tmp.append(
                cmath.exp(1.0j * (i - j / 2) * alpha)
                * d[i][k]
                * cmath.exp(1.0j * (k - j / 2) * gamma)
            )
        ret.append(tmp)
    return ret


j0 = 1 / 2  # spin of Lc+
Dm = D_matrix_conj(phi1, theta1, 0.0, 2 * j0)

With helicity amplitude and D-matrix we can calculate the amplitude of this decay as

\[A_{\lambda_0,\lambda_1,\lambda_2} = H_{\lambda_1,\lambda_2} D_{\lambda_0,\lambda_1-\lambda_2}^{J_0*}(\phi, \theta,0)\]

def sum_decay_amplitude(H, D, j0, l0, l1, l2):
    ret = {}
    for i in l0:
        for j in l1:
            for k in l2:
                if abs(j - k) <= j0:
                    idx1 = l0.index(i)
                    idx2 = l0.index(j - k)
                    ret[(i, j, k)] = H[(j, k)] * Dm[idx1][idx2]
                else:
                    ret[(i, j, k)] = 0.0j
    return ret


A1 = sum_decay_amplitude(H1, Dm, 1 / 2, [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0])

pprint(A1)

{(-0.5, -0.5, 0): (1.0246845284339865+0.09214359867084504j),
 (-0.5, 0.5, 0): (0.043001323998895086+0.23632845496901808j),
 (0.5, -0.5, 0): (-0.24569717760383605+0.3542951434609425j),
 (0.5, 0.5, 0): (0.5410098258810769+0.18934959397575796j)}

For the second decay, \(\Sigma^{+}\rightarrow \Lambda \pi^{+}\), the procedure is similar. But it is strong decay (as we do not use p_break: True), so we required that \((-1)^{l} = P_0 P_1 P_2 = 1 \times 1 \times -1\), only \(l=1\) remain.

q = get_relative_p(M(p_sigma), M(p_lambda), M(p_pip))
q0 = get_relative_p(m_sigma_1660, m_l, m_pi)

gls0 = params["Sigma->Lambda.pip_g_ls_0r"] * cmath.exp(
    1j * params["Sigma->Lambda.pip_g_ls_0i"]
)

G1 = gls0 * (q / q0) ** 1 * Bprime(1, q, q0)

H2 = {(-1 / 2, 0): G1 / math.sqrt(2), (1 / 2, 0): -G1 / (math.sqrt(2))}

j0 = 1 / 2  # Sigma+
Dm = D_matrix_conj(phi2, theta2, 0.0, 2 * j0)

A2 = sum_decay_amplitude(H2, Dm, 1 / 2, [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0])

pprint(A2)

{(-0.5, -0.5, 0): (0.16504850460193102-0.10596431571836415j),
 (-0.5, 0.5, 0): (-0.638231980063331+0.4097572116760412j),
 (0.5, -0.5, 0): (-0.638231980063331-0.4097572116760412j),
 (0.5, 0.5, 0): (-0.16504850460193102-0.10596431571836415j)}

For the decay, \(\Lambda \rightarrow p \pi^{-}\), the process is similar as \(\Lambda_c^{+}\rightarrow \Sigma(1660)^{+} \pi^0\).

q = get_relative_p(M(p_lambda), M(p_p), M(p_pim))
q0 = get_relative_p(m_l, m_p, m_pi)
gls0 = params["Lambda->p.pim_g_ls_0r"] * cmath.exp(
    1j * params["Lambda->p.pim_g_ls_0i"]
)
gls1 = params["Lambda->p.pim_g_ls_1r"] * cmath.exp(
    1j * params["Lambda->p.pim_g_ls_1i"]
)


G0 = gls0 * (q / q0) ** 0 * Bprime(0, q, q0)
G1 = gls1 * (q / q0) ** 1 * Bprime(1, q, q0)

H3 = {
    (-1 / 2, 0): (G0 + G1) / math.sqrt(2),
    (1 / 2, 0): (G0 - G1) / math.sqrt(2),
}


Dm = D_matrix_conj(phi3, theta3, 0.0, 2 * 1 / 2)

A3 = sum_decay_amplitude(H3, Dm, 1 / 2, [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0])

pprint(A3)

{(-0.5, -0.5, 0): (0.6721807828658565-0.06763640483927624j),
 (-0.5, 0.5, 0): (0.2964610424603069-0.029830604501111965j),
 (0.5, -0.5, 0): (-0.7536574536506793-0.07583477829866517j),
 (0.5, 0.5, 0): (0.264411125564959+0.02660566679171078j)}

Then the total amplitude of the decay chain is

\[A_{ \lambda_{\Lambda_c^{+}},\lambda_{p},\lambda_{\pi^{-}},\lambda_{\pi^{+}},\lambda_{\pi^{0}} } = a_{total} \sum_{\lambda_{\Sigma^{+}},\lambda_{\Lambda}} A_{\lambda_{\Lambda_c^{+}},\lambda_{\Sigma^{+}},\lambda_{\pi^{0}} } R_{\Sigma^{+}}(m_{\Sigma^{+}}) A_{\lambda_{\Sigma^{+}},\lambda_{\Lambda},\lambda_{\pi^{+}} } R_{\Lambda}(m_{\Lambda}) A_{\lambda_{\Lambda},\lambda_{p},\lambda_{\pi^{-}} }\]

\(a_{total}\) is a factor in fit parameters. The amplitude sum of all possible helicities of intermidiate states. To be simple, we use the simplest BW shape \(\Sigma(1660)^{+}\), as

\[R(m) = \frac{1}{m_0^2 - m^2 - i m_0 \Gamma_0 }.\]

and we do not consider the shape of \(\Lambda\)

def bw(m, m0, g0):
    return 1 / (m0 * m0 - m * m - 1.0j * m0 * g0)


r_sigma = bw(M(p_sigma), m_sigma_1660, g_sigma_1660)
r_lambda = 1.0

total = params[
    "Lc->Sigma.pi0Sigma->Lambda.pipLambda->p.pim_total_0r"
] * cmath.exp(
    1j * params["Lc->Sigma.pi0Sigma->Lambda.pipLambda->p.pim_total_0i"]
)

A = {}
for l_lc, l_p, l_pim, l_pip, l_pi0 in itertools.product(
    [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0], [0], [0]
):
    tmp = 0.0
    for l_sigma in [-1 / 2, 1 / 2]:
        for l_l in [-1 / 2, 1 / 2]:
            tmp = (
                tmp
                + A1[(l_lc, l_sigma, l_pi0)]
                * A2[(l_sigma, l_l, l_pip)]
                * A3[(l_l, l_p, l_pim)]
            )
    A[(l_lc, l_p, l_pim, l_pip, l_pi0)] = total * tmp * r_sigma * r_lambda

pprint(A)

{(-0.5, -0.5, 0, 0, 0): (-0.3953877755784625+0.5987680982164796j),
 (-0.5, 0.5, 0, 0, 0): (0.10669876852931258-0.03549608323128384j),
 (0.5, -0.5, 0, 0, 0): (0.08461940080832495-0.21173858458097958j),
 (0.5, 0.5, 0, 0, 0): (0.15116503907073653+0.107932131617686j)}

2.1 Second decay chain

For the first decay, \(\Lambda_c^{+}\rightarrow \rho^{+} \Lambda\), \(s\) in range \([|1-\frac{1}{2}|, 1+\frac{1}{2}]\), so it could be \(\frac{1}{2}\) or \(\frac{3}{2}\). when \(s=\frac{1}{2}\), \(l\) in range \([|\frac{1}{2}-\frac{1}{2}|, \frac{1}{2}+ \frac{1}{2}]\), could be \(0\) or \(1\). when \(s=\frac{3}{2}\), \(l\) in range \([|\frac{3}{2}-\frac{1}{2}|, \frac{3}{2}+ \frac{1}{2}]\), could be \(1\) or \(2\). All possible \((l,s)\) is \((0, \frac{1}{2})\), \((1, \frac{1}{2})\), \((1, \frac{3}{2})\) and \((2, \frac{3}{2})\).

we can find the relation as

\[H_{-1, -\frac{1}{2}} = -\frac{ g_{0,\frac{1}{2}} }{\sqrt{3}} - \frac{ g_{1,\frac{1}{2}} }{\sqrt{3}} - \frac{ g_{1,\frac{3}{2}} }{\sqrt{6}} - \frac{ g_{2,\frac{3}{2}} }{\sqrt{6}}\]

\[H_{0, -\frac{1}{2}} = -\frac{ g_{0,\frac{1}{2}} }{\sqrt{6}} + \frac{ g_{1,\frac{1}{2}} }{\sqrt{6}} - \frac{ g_{1,\frac{3}{2}} }{\sqrt{3}} + \frac{ g_{2,\frac{3}{2}} }{\sqrt{3}}\]

\[H_{0, \frac{1}{2}} = \frac{ g_{0,\frac{1}{2}} }{\sqrt{6}} + \frac{ g_{1,\frac{1}{2}} }{\sqrt{6}} - \frac{ g_{1,\frac{3}{2}} }{\sqrt{3}} - \frac{ g_{2,\frac{3}{2}} }{\sqrt{3}}\]

\[H_{\frac{1}{2},1} = \frac{ g_{0,\frac{1}{2}} }{\sqrt{3}} - \frac{ g_{1,\frac{1}{2}} }{\sqrt{3}} - \frac{ g_{1,\frac{3}{2}} }{\sqrt{6}} + \frac{ g_{2,\frac{3}{2}} }{\sqrt{6}}\]

m0, m1, m2 = M(p_total), M(p_rho), M(p_lambda)

q = get_relative_p(m0, m1, m2)
q0 = get_relative_p(m_lc, m_rho770, m_l)

gls01 = params["Lc->rho.Lambda_g_ls_0r"] * cmath.exp(
    1j * params["Lc->rho.Lambda_g_ls_0i"]
)
gls11 = params["Lc->rho.Lambda_g_ls_1r"] * cmath.exp(
    1j * params["Lc->rho.Lambda_g_ls_1i"]
)
gls13 = params["Lc->rho.Lambda_g_ls_2r"] * cmath.exp(
    1j * params["Lc->rho.Lambda_g_ls_2i"]
)
gls23 = params["Lc->rho.Lambda_g_ls_3r"] * cmath.exp(
    1j * params["Lc->rho.Lambda_g_ls_3i"]
)


G01 = gls01 * (q / q0) ** 0 * Bprime(0, q, q0)
G11 = gls11 * (q / q0) ** 1 * Bprime(1, q, q0)
G13 = gls13 * (q / q0) ** 1 * Bprime(1, q, q0)
G23 = gls23 * (q / q0) ** 2 * Bprime(2, q, q0)


H1_prime = {
    (-1, -1 / 2): 1
    / math.sqrt(6)
    * (math.sqrt(2) * (-G01 - G11) + (-G13 - G23)),
    (0, -1 / 2): 1
    / math.sqrt(6)
    * ((-G01 + G11) + math.sqrt(2) * (-G13 + G23)),
    (0, 1 / 2): 1 / math.sqrt(6) * ((G01 + G11) + math.sqrt(2) * (-G13 - G23)),
    (1, 1 / 2): 1 / math.sqrt(6) * (math.sqrt(2) * (G01 - G11) + (-G13 + G23)),
}

Dm = D_matrix_conj(phi1_2, theta1_2, 0, 2 * 1 / 2)

A1_prime = sum_decay_amplitude(
    H1_prime, Dm, 1 / 2, [-1 / 2, 1 / 2], [-1, 0, 1], [-1 / 2, 1 / 2]
)
pprint(A1_prime)

{(-0.5, -1, -0.5): (-0.8386137946350432-1.4340009423866682j),
 (-0.5, -1, 0.5): 0j,
 (-0.5, 0, -0.5): (-0.02632096486978926-0.04500795082226156j),
 (-0.5, 0, 0.5): (-0.08050457224207319-0.1376600685561207j),
 (-0.5, 1, -0.5): 0j,
 (-0.5, 1, 0.5): (-0.00639040603033944-0.01092737602024329j),
 (0.5, -1, -0.5): (0.33376536678723034-0.5707273760232039j),
 (0.5, -1, 0.5): 0j,
 (0.5, 0, -0.5): (-0.06613365682719523+0.11308629409673149j),
 (0.5, 0, 0.5): (0.03204053910670293-0.054788227390874515j),
 (0.5, 1, -0.5): 0j,
 (0.5, 1, 0.5): (-0.016056437197026326+0.027455959130660786j)}

For the second decay, \(\rho^{+}\rightarrow \pi^{+} \pi^{0}\), \(s\) in range \([|0-0|, 0+0]\), so it is only \(0\). \(l\) in range \([|1-0|, 1+0]\), so it is only \(1\) . All possible \((l,s)\) is \((1, 0)\). It also satisfies the requirement that \((-1)^l = P_0 P_1 P_2 = -1\). we can find the relation as

\[H_{0,0} = g_{1,0}\]

m0, m1, m2 = M(p_rho), M(p_pip), M(p_pi0)

q = get_relative_p(m0, m1, m2)
q0 = get_relative_p(m_rho770, m_pi, m_pi0)

gls = params["rho->pip.pi0_g_ls_0r"] * cmath.exp(
    1j * params["rho->pip.pi0_g_ls_0i"]
)

G = gls * (q / q0) ** 1 * Bprime(1, q, q0)

H2_prime = {(0, 0): G}

Dm = D_matrix_conj(phi2_2, theta2_2, 0.0, 2 * 1)

A2_prime = sum_decay_amplitude(H2_prime, Dm, 1, [-1, 0, 1], [0], [0])

pprint(A2_prime)

{(-1, 0, 0): (-0.10904456128821746-0.3101280597575634j),
 (0, 0, 0): (0.9984404935477759+0j),
 (1, 0, 0): (0.10904456128821746-0.3101280597575634j)}

For the last decay, \(\Lambda \rightarrow p \pi^{-}\), the process is similar as the deacy chain of \(\Sigma^{+}\), but it has different angles.

H3_prime = H3

Dm = D_matrix_conj(phi3_2, theta3_2, 0.0, 2 * 1 / 2)
A3_prime = sum_decay_amplitude(
    H3_prime, Dm, 1 / 2, [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0]
)

pprint(A3_prime)

{(-0.5, -0.5, 0): (0.45626222979386544-0.4004340537172614j),
 (-0.5, 0.5, 0): (0.24048250582071312-0.21105710349370635j),
 (0.5, -0.5, 0): (-0.6113499145798885-0.536545233309202j),
 (0.5, 0.5, 0): (0.17947673127193456+0.1575159861109122j)}

The total amplitude of this decay chain is

r_rho = bw(M(p_rho), m_rho770, g_rho770)

total = params["Lc->rho.Lambdarho->pip.pi0Lambda->p.pim_total_0r"] * cmath.exp(
    1j * params["Lc->rho.Lambdarho->pip.pi0Lambda->p.pim_total_0i"]
)


A_prime = {}
for l_lc, l_p, l_pim, l_pip, l_pi0 in itertools.product(
    [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0], [0], [0]
):
    tmp = 0.0
    for l_rho in [-1, 0, 1]:
        for l_l in [-1 / 2, 1 / 2]:
            tmp = (
                tmp
                + A1_prime[(l_lc, l_rho, l_l)]
                * A2_prime[(l_rho, l_pip, l_pi0)]
                * A3_prime[(l_l, l_p, l_pim)]
            )
    A_prime[(l_lc, l_p, l_pim, l_pip, l_pi0)] = total * tmp * r_rho * r_lambda

pprint(A_prime)

{(-0.5, -0.5, 0, 0, 0): (1.8045559263449882-1.8920214364097845j),
 (-0.5, 0.5, 0, 0, 0): (0.4955793141170991-0.5776277297411841j),
 (0.5, -0.5, 0, 0, 0): (1.2127312641251653-0.17286687566831907j),
 (0.5, 0.5, 0, 0, 0): (0.4339860510630706-0.20365827206528628j)}

2.3 Total amplitude

Before sum of all decay chain, we align the coordinate system of all the finnale particles. For \(\pi^{\pm,0}\), \(J=0\) so \(D^{J}(\alpha,\beta,\gamma)=1\), we can remove them. But we need to align \(p\) in the decay chain of \(\rho(770)^{+}\).

\[A^{aligned}_{\lambda_{\Lambda_c^+},\lambda_{p},\lambda_{\pi^-},\lambda_{\pi^+},\lambda_{\pi^0}} = \sum_{\lambda_{p}'} A_{\lambda_{\Lambda_c^+},\lambda_{p}',\lambda_{\pi^-},\lambda_{\pi^+},\lambda_{\pi^0}} D^{J_p*}_{\lambda_{p}',\lambda_{p}}(\alpha_p,\beta_p,\gamma_p)\]

Dm = D_matrix_conj(ang[0], ang[1], ang[2], 2 * 1 / 2)


A_prime_aligned = {}
for l_lc, l_p, l_pim, l_pip, l_pi0 in itertools.product(
    [-1 / 2, 1 / 2], [-1 / 2, 1 / 2], [0], [0], [0]
):
    tmp = 0.0
    for l_p_a in [-1 / 2, 1 / 2]:
        idx1 = [-1 / 2, 1 / 2].index(l_p_a)
        idx2 = [-1 / 2, 1 / 2].index(l_p)
        tmp = (
            tmp + A_prime[(l_lc, l_p_a, l_pim, l_pip, l_pi0)] * Dm[idx1][idx2]
        )
    A_prime_aligned[(l_lc, l_p, l_pim, l_pip, l_pi0)] = tmp

pprint(A_prime_aligned)

{(-0.5, -0.5, 0, 0, 0): (-1.7694281803357732+1.9249132764295112j),
 (-0.5, 0.5, 0, 0, 0): (-0.5061264381134051+0.5684089015747082j),
 (0.5, -0.5, 0, 0, 0): (-1.2093442875118885+0.1951575516148681j),
 (0.5, 0.5, 0, 0, 0): (-0.4376608117346166+0.19563639275383493j)}

And then we can sum over all decay chains to get the total amplutude of the process

A_total = {}
for k in A:
    tmp = 0.0
    for chain_amps in [A, A_prime_aligned]:
        tmp = tmp + chain_amps[k]
    A_total[k] = tmp

pprint(A_total)

{(-0.5, -0.5, 0, 0, 0): (-2.1648159559142357+2.5236813746459905j),
 (-0.5, 0.5, 0, 0, 0): (-0.3994276695840925+0.5329128183434244j),
 (0.5, -0.5, 0, 0, 0): (-1.1247248867035635-0.016581032966111464j),
 (0.5, 0.5, 0, 0, 0): (-0.28649577266388004+0.30356852437152093j)}

2.4 Prbability density

If no polarization in the process, the probability density function (not normalized) is

\[P \propto \sum_{\lambda_{\Lambda_c^+},\lambda_{p},\lambda_{\pi^-},\lambda_{\pi^+},\lambda_{\pi^0}}|A_{\lambda_{\Lambda_c^+},\lambda_{p},\lambda_{\pi^-},\lambda_{\pi^+},\lambda_{\pi^0}}|^2\]

prob = 0.0
for i in A_total.values():
    prob = prob + abs(i) ** 2

print("prob", prob)

prob 12.938449017067997

3. Validation

we can quackly check the results with

import tensorflow as tf
import yaml

from tf_pwa.config_loader import ConfigLoader

config = ConfigLoader(yaml.full_load(config_str))

config.set_params(params)

angle = config.data.cal_angle(
    tf.unstack(
        tf.cast(tf.stack([[p_p], [p_pim], [p_pip], [p_pi0]]), dtype=tf.float64)
    )
)

amp = config.get_amplitude()

prob2 = amp(angle)

print("prob ", prob2)

assert abs(prob2 - prob) / prob < 1e-6

prob  tf.Tensor([12.93845315], shape=(1,), dtype=float64)

We can get the same prob values, (would be some difference due to the precision).

The following code is used for an icon in Sphinx-Gallery.

import matplotlib.pyplot as plt

plt.axis("off")
plt.text(
    0, 0, "$P\\propto \\sum_{\\lambda}|A_{\\lambda}(p^{\mu})|^2$", fontsize=60
)
plt.tight_layout()