Kinematics

Kinematics#

Hide code cell content
import sympy as sp
from IPython.display import Math

from ampform.io import aslatex

Lorentz vectors#

AmpForm provides classes for formulating symbolic expressions for boosting and rotating Lorentz vectors. Usually, when building an amplitude model, you don’t have to use these classes, but sometimes you want to boost some four-momenta yourself (for instance to boost into the center-of-mass frame of your experiment. Here, we boost a four-momentum \(q\) from the lab frame of the BESIII detector into the center-of-mass frame \(p\) of the \(e^-e^+\) collision. Symbolically, this looks like this:

from ampform.kinematics.lorentz import (
    ArrayMultiplication,
    ArraySize,
    BoostZMatrix,
    Energy,
    FourMomentumSymbol,
    three_momentum_norm,
)

p = FourMomentumSymbol("p", shape=[])
q = FourMomentumSymbol("q", shape=[])
beta = three_momentum_norm(p) / Energy(p)
Bz = BoostZMatrix(beta, n_events=ArraySize(beta))
Bz_expr = ArrayMultiplication(Bz, q)
Bz_expr
\[\displaystyle \boldsymbol{B_z}\left(\frac{\left|\vec{p}\right|}{E\left(p\right)}\right) q\]

We now use SymPy to create a numerical function. (Of course, you can use TensorWaves instead to use other numerical backends.)

Bz_func = sp.lambdify([p, q], Bz_expr.doit(), cse=True)
Hide code cell content
import inspect

from black import FileMode, format_str

src = inspect.getsource(Bz_func)
src = format_str(src, mode=FileMode())
print(src)
def _lambdifygenerated(p, q):
    x0 = p
    x1 = x0[:, 0]
    x2 = sum(x0[:, 1:] ** 2, axis=1)
    x3 = sqrt(x2) / x1
    x4 = 1 / sqrt(1 - x2 / x1**2)
    x5 = len(x3)
    return einsum(
        "...ij,...j->...i",
        array(
            [
                [x4, zeros(x5), zeros(x5), -x3 * x4],
                [zeros(x5), ones(x5), zeros(x5), zeros(x5)],
                [zeros(x5), zeros(x5), ones(x5), zeros(x5)],
                [-x3 * x4, zeros(x5), zeros(x5), x4],
            ]
        ).transpose((2, 0, 1)),
        q,
    )

Finally, plugging in some numbers that represent data, we get the \(q\) in the rest frame of \(p\):

import numpy as np

pz_array = np.array([[3.0971, 0, 0, 30e-3]])  # J/psi in BESIII lab frame
q_array = np.array([
    [2.4, 0.3, -1.5, 0.02],
    [3.4, -0.045, 0.6, 1.1],
    # list of measured four-momenta q in lab frame
])
Bz_func(pz_array, q_array)
array([[ 2.39991886e+00,  3.00000000e-01, -1.50000000e+00,
        -3.24770653e-03],
       [ 3.38950389e+00, -4.50000000e-02,  6.00000000e-01,
         1.06711603e+00]])

Four-vector array format

Lambdified expressions that involve Lorentz vector computations, expect the format \(p = \left(E, p_x, p_y, p_z\right)\). In addition, the shape of input arrays should be (n, 4) with n the number of events.

As a cross-check, notice how boosting the original boost momentum into its own rest frame, results in \(B_z(p) p = \left(m_{J/\psi}, 0, 0, 0\right)\):

Bz_func(pz_array, pz_array)
array([[3.0969547, 0.       , 0.       , 0.       ]])

Note that in this case, boost vector \(p\) was in the \(z\) direction, so we were able to just boost with BoostZMatrix. In the more general case, we can use:

from ampform.kinematics.lorentz import BoostMatrix

B = BoostMatrix(p)
B_expr = ArrayMultiplication(B, q)
B_expr
\[\displaystyle \boldsymbol{B}\left(p\right) q\]
B_func = sp.lambdify([p, q], B_expr.doit(), cse=True)
px_array = np.array([[3.0971, 30e-3, 0, 0]])  # x direction!
B_func(px_array, q_array)
array([[ 2.39720652,  0.27676543, -1.5       ,  0.02      ],
       [ 3.40059543, -0.07793769,  0.6       ,  1.1       ]])

And again, \(B(p) p = \left(m_{J/\psi}, 0, 0, 0\right)\):

B_func(px_array, px_array)
array([[3.0969547, 0.       , 0.       , 0.       ]])

Phase space#

Kinematics for a three-body decay \(0 \to 123\) can be fully described by two Mandelstam variables \(\sigma_1, \sigma_2\), because the third variable \(\sigma_3\) can be expressed in terms \(\sigma_1, \sigma_2\), the mass \(m_0\) of the initial state, and the masses \(m_1, m_2, m_3\) of the final state. As can be seen, the roles of \(\sigma_1, \sigma_2, \sigma_3\) are interchangeable.

Hide code cell source
from ampform.kinematics.phasespace import compute_third_mandelstam

m0, m1, m2, m3 = sp.symbols("m:4")
s1, s2, s3 = sp.symbols("sigma1:4")
s3_expr = compute_third_mandelstam(s1, s2, m0, m1, m2, m3)

latex = aslatex({s3: s3_expr})
Math(latex)
\[\begin{split}\displaystyle \begin{array}{rcl} \sigma_{3} &=& m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \sigma_{1} - \sigma_{2} \\ \end{array}\end{split}\]

The phase space is defined by the closed area that satisfies the condition \(\phi(\sigma_1,\sigma_2) \leq 0\), where \(\phi\) is a Kibble function:

Hide code cell source
from ampform.kinematics.phasespace import Kibble

kibble = Kibble(s1, s2, s3, m0, m1, m2, m3)

latex = aslatex({kibble: kibble.evaluate()})
Math(latex)
\[\begin{split}\displaystyle \begin{array}{rcl} \phi\left(\sigma_{1}, \sigma_{2}\right) &=& \lambda\left(\lambda\left(\sigma_{2}, m_{2}^{2}, m_{0}^{2}\right), \lambda\left(\sigma_{3}, m_{3}^{2}, m_{0}^{2}\right), \lambda\left(\sigma_{1}, m_{1}^{2}, m_{0}^{2}\right)\right) \\ \end{array}\end{split}\]

and \(\lambda\) is the Källén function:

Hide code cell source
from ampform.kinematics.phasespace import Kallen

x, y, z = sp.symbols("x:z")
kallen = Kallen(x, y, z)

latex = aslatex({kallen: kallen.evaluate()})
Math(latex)
\[\begin{split}\displaystyle \begin{array}{rcl} \lambda\left(x, y, z\right) &=& x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \\ \end{array}\end{split}\]

Any distribution over the phase space can now be defined using a two-dimensional grid over a Mandelstam pair \(\sigma_1,\sigma_2\) of choice, with the condition \(\phi(\sigma_1,\sigma_2)<0\) selecting the values that are physically allowed.

Hide code cell source
from ampform.kinematics.phasespace import is_within_phasespace

is_within_phasespace(s1, s2, m0, m1, m2, m3)
\[\begin{split}\displaystyle \begin{cases} 1 & \text{for}\: \phi\left(\sigma_{1}, \sigma_{2}\right) \leq 0 \\\text{NaN} & \text{otherwise} \end{cases}\end{split}\]

See Phase space for a three-body decay for an interactive visualization of the phase space region and an analytic expression for the phase space boundary.