7.3. Phase space sample

Import Python libraries Hide code cell content

import matplotlib.pyplot as plt
import sympy as sp
from ampform.kinematics.phasespace import (
    Kallen,
    Kibble,
    compute_third_mandelstam,
    is_within_phasespace,
)
from matplotlib.colors import LogNorm

from polarimetry.data import (
    create_mass_symbol_mapping,
    create_phase_space_filter,
    generate_meshgrid_sample,
    generate_phasespace_sample,
)
from polarimetry.io import display_doit, display_latex, mute_jax_warnings
from polarimetry.lhcb import load_model_builder
from polarimetry.lhcb.particle import load_particles
from polarimetry.plot import use_mpl_latex_fonts

mute_jax_warnings()

7.3.1. Definition

See also

AmpForm’s Kinematics page.

Show code cell source Hide code cell source

m0, mi, mj, mk = sp.symbols("m0 m_(i:k)", nonnegative=True)
σi, σj, σk = sp.symbols("sigma_(i:k)", nonnegative=True)
is_within_phasespace(σi, σj, m0, mi, mj, mk)

$$\begin{array}{r}{\displaystyle \{\begin{array}{ll}1& \text{for}\varphi ({\sigma }_i,{\sigma }_j)\le 0\\ \text{NaN}& \text{otherwise}\end{array}}\end{array}$$

Show code cell source Hide code cell source

display_doit(Kibble(σi, σj, σk, m0, mi, mj, mk))

$$\begin{array}{r}{\displaystyle \begin{array}{rcl}\varphi ({\sigma }_i,{\sigma }_j)& =& \lambda (\lambda ({\sigma }_j,m_j^2,m_0^2),\lambda ({\sigma }_k,m_k^2,m_0^2),\lambda ({\sigma }_i,m_i^2,m_0^2))\end{array}}\end{array}$$

Show code cell source Hide code cell source

display_doit(Kallen(*sp.symbols("x:z")))

$$\begin{array}{r}{\displaystyle \begin{array}{rcl}\lambda (x,y,z)& =& x^2-2xy-2xz+y^2-2yz+z^2\end{array}}\end{array}$$

m1, m2, m3 = sp.symbols("m1:4")
display_latex({σk: compute_third_mandelstam(σi, σj, m0, m1, m2, m3)})

$$\begin{array}{r}{\displaystyle \begin{array}{rcl}{\sigma }_k& =& m_0^2+m_1^2+m_2^2+m_3^2-{\sigma }_i-{\sigma }_j\end{array}}\end{array}$$

7.3.2. Visualization

Show code cell source Hide code cell source

model_choice = 0
model_file = "../../data/model-definitions.yaml"
particles = load_particles("../../data/particle-definitions.yaml")
amplitude_builder = load_model_builder(model_file, particles, model_choice)
decay = amplitude_builder.decay
display_latex(create_mass_symbol_mapping(decay))

$$\begin{array}{r}{\displaystyle \begin{array}{rcl}{m}_0& =& 2.28646\\ m_1& =& 0.938272046\\ m_2& =& 0.13957018\\ m_3& =& 0.49367700000000003\end{array}}\end{array}$$

Show code cell source Hide code cell source

def plot_phsp_boundary(ax, x_mandelstam: int, y_mandelstam: int):
    ax.set_xticks([])
    ax.set_xlabel(Rf"$\sigma_{x_mandelstam}$ [GeV$^2$]")
    ax.set_ylabel(Rf"$\sigma_{y_mandelstam}$ [GeV$^2$]")
    phsp = generate_meshgrid_sample(decay, resolution, x_mandelstam, y_mandelstam)
    phsp_filter = create_phase_space_filter(
        decay, x_mandelstam, y_mandelstam, outside_value=0
    )
    mesh = ax.contour(
        phsp[f"sigma{x_mandelstam}"],
        phsp[f"sigma{y_mandelstam}"],
        phsp_filter(phsp),
        colors="black",
    )
    contour = mesh.collections[0]
    contour.set_facecolor("lightgray")


%config InlineBackend.figure_formats = ['svg']
resolution = 500
use_mpl_latex_fonts()
plt.rc("font", size=18)
fig, ax = plt.subplots(figsize=(5, 5), tight_layout=True)
plot_phsp_boundary(ax, 1, 2)
fig.savefig("../_images/phase-space-boundary.svg")
plt.show()

Show code cell source Hide code cell source

def plot_phsp_distribution(ax, x_mandelstam: int, y_mandelstam: int):
    ax.set_xticks([])
    ax.set_xlabel(Rf"$\sigma_{x_mandelstam}$")
    ax.set_ylabel(Rf"$\sigma_{y_mandelstam}$")
    ax.hist2d(
        phsp[f"sigma{x_mandelstam}"],
        phsp[f"sigma{y_mandelstam}"],
        bins=500,
        norm=LogNorm(),
    )


%config InlineBackend.figure_formats = ['png']
phsp = generate_phasespace_sample(decay, n_events=10_000_000, seed=0)
fig, axes = plt.subplots(ncols=3, dpi=200, figsize=(15, 4.9), tight_layout=True)
plot_phsp_distribution(axes[0], 1, 2)
plot_phsp_distribution(axes[1], 2, 3)
plot_phsp_distribution(axes[2], 3, 1)
fig.tight_layout()
plt.show()