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 reduce_svg_size, 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{split}\displaystyle \begin{cases} 1 & \text{for}\: \phi\left(\sigma_{i}, \sigma_{j}\right) \leq 0 \\\text{NaN} & \text{otherwise} \end{cases}\end{split}\]

Show code cell source

Hide code cell source

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

\[\begin{split}\displaystyle \begin{array}{rcl} \phi\left(\sigma_{i}, \sigma_{j}\right) &=& \lambda\left(\lambda\left(\sigma_{j}, m_{j}^{2}, m_{0}^{2}\right), \lambda\left(\sigma_{k}, m_{k}^{2}, m_{0}^{2}\right), \lambda\left(\sigma_{i}, m_{i}^{2}, m_{0}^{2}\right)\right) \\ \end{array}\end{split}\]

Show code cell source

Hide code cell source

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

\[\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}\]

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

\[\begin{split}\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{split}\]

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{split}\displaystyle \begin{array}{rcl} m_{0} &=& 2.28646 \\ m_{1} &=& 0.938272046 \\ m_{2} &=& 0.13957018 \\ m_{3} &=& 0.49367700000000003 \\ \end{array}\end{split}\]

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
    )
    xyz = (
        phsp[f"sigma{x_mandelstam}"],
        phsp[f"sigma{y_mandelstam}"],
        phsp_filter(phsp),
    )
    ax.contourf(*xyz, colors="lightgray", levels=[0.5, 1])
    ax.contour(*xyz, colors="black", levels=[0.5, 1])


resolution = 500
use_mpl_latex_fonts()
plt.rc("font", size=18)
fig, ax = plt.subplots(figsize=(6, 5), tight_layout=True)
fig.patch.set_facecolor("none")
ax.patch.set_color("none")
plot_phsp_boundary(ax, 1, 2)
output_path = "../_images/phase-space-boundary.svg"
fig.savefig(output_path, bbox_inches="tight")
reduce_svg_size(output_path)
plt.show(fig)

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(),
        rasterized=True,
    )


phsp = generate_phasespace_sample(decay, n_events=10_000_000, seed=0)
fig, axes = plt.subplots(dpi=200, figsize=(15, 4.9), ncols=3, tight_layout=True)
fig.patch.set_facecolor("none")
for ax in axes.flatten():
    ax.patch.set_color("none")
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(fig)