Construct three-body decay

\(\Lambda_b^0 \to p K^- \gamma\)

Model definition: lb2pkg-lhcb-2765817.json.

This page demonstrates deserialization and evaluation of an amplitude model for the decay \(\Lambda_b^0 \to p K^- \gamma\). The resonant structure was studied using proton-proton collision data recorded at centre-of-mass energies of \(7\), \(8\), and \(13\) TeV with the LHCb detector, INSPIRE-HEP 2765817.

import logging
from pathlib import Path

import jax.numpy as jnp
import matplotlib.pyplot as plt
import pandas as pd
import sympy as sp
from ampform_dpd.decay import FinalStateID, State, ThreeBodyDecay
from ampform_dpd.io import aslatex
from ampform_dpd.io.serialization import (
    compile_workspace,
    load_workspace,
    validate_checksums,
)
from ampform_dpd.io.serialization.format import get_decay_chains
from IPython.display import Math
from matplotlib_inline.backend_inline import set_matplotlib_formats
from mpl_toolkits.axes_grid1 import make_axes_locatable  # cspell:ignore mpl_toolkits

this_dir = Path(".").absolute()
logging.getLogger("ampform.sympy").setLevel(logging.ERROR)
set_matplotlib_formats("svg")

model_path = this_dir.parent.parent / "models" / "lb2pkg-lhcb-2765817.json"

def to_latex(name: str) -> str:
    latex = {
        "Lb": R"\Lambda_b^0",
        "g": R"\gamma",
        "K": "K^-",
        "p": "p",
    }.get(name)
    if latex is not None:
        return latex
    if name.startswith("LNR"):
        return R"\Lambda_\mathrm{NR}"
    if name.startswith("L"):
        return Rf"\Lambda({name[1:]})"
    return name

workspace = load_workspace(model_path, to_latex=to_latex)
model_definition = workspace.definition
((intensity_name, decay),) = workspace.decays.items()
Math(aslatex(decay, with_jp=True))

\(\displaystyle \begin{array}{c} \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1405)\left[J=\frac{1}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1520)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1600)\left[J=\frac{1}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1670)\left[J=\frac{1}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1690)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1800)\left[J=\frac{1}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1810)\left[J=\frac{1}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1820)\left[J=\frac{5}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1830)\left[J=\frac{5}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(1890)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(2100)\left[J=\frac{7}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(2110)\left[J=\frac{5}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda(2350)\left[J=\frac{9}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \Lambda_b^0\left[J=\frac{1}{2}\right] \to \left(\Lambda_\mathrm{NR}\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] K^-\left[J=0\right]\right) \gamma\left[J=1\right] \\ \end{array}\)

Dynamics

The model contains Breit-Wigner and multichannel Breit-Wigner propagators, as well as a constant non-resonant term. The latter is serialized as a zeroth-order polynomial.

chain_defs = get_decay_chains(model_definition)

l1520_bw = workspace.functions["L1520_BW"]
Math(aslatex(l1520_bw))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_{L=2}\left(\sigma_{3}; m_{\Lambda(1520)}, \Gamma_{\Lambda(1520)}\right) &=& \mathcal{R}^\mathrm{BW}\left(\sigma_{3}; m_{\Lambda(1520)}, \Gamma_{1520}\left(\sigma_{3}\right)\right) \\ m_{\Lambda(1520)} &=& 1.519 \\ \Gamma_{\Lambda(1520)} &=& 0.016 \\ m_{1} &=& 0.938 \\ m_{2} &=& 0.493 \\ R_\mathrm{res} &=& 1.5 \\ \end{array}\)

l1405_flatte = workspace.functions["L1405_Flatte"]
Math(aslatex(l1405_flatte))

\(\displaystyle \begin{array}{rcl} \frac{1}{m_{\Lambda(1405)}^{2} - \sigma_{3} - i \left(\Gamma_{\Lambda(1405)} \mathcal{F}_{0}\left(\sigma_{3}, m_{1}, m_{2}\right)^{2} \rho^\mathrm{c}_{3}\left(\sigma_{3}\right) + \Gamma_{\Lambda(1405)}^\text{ch. 2} \mathcal{F}_{0}\left(\sigma_{3}, m_{a,2}, m_{b,2}\right)^{2} \rho^\mathrm{c}_{3}\left(\sigma_{3}\right)\right)} \\ m_{\Lambda(1405)} &=& 1.405 \\ \Gamma_{\Lambda(1405)} &=& 0.3504277764790773 \\ m_{1} &=& 0.938 \\ m_{2} &=& 0.493 \\ R_{\Lambda(1405)} &=& 0 \\ m_{a,2} &=& 1.197 \\ m_{b,2} &=& 0.14 \\ \Gamma_{\Lambda(1405)}^\text{ch. 2} &=& 0.3504277764790773 \\ \end{array}\)

Construct amplitude model

load_workspace formulates LS-coupled models and their serialized dynamics without notebook-specific compatibility adapters.

model = workspace.distributions[intensity_name]

Validation

validate_checksums compiles each serialized target independently and records numeric differences or diagnostics without aborting the remaining checks.

The status marks 🟢, 🟡, and 🔴 indicate absolute differences of \(<10^{-10}\), \(<10^{-2}\), and \(\ge10^{-2}\), respectively. Failed evaluations are marked 🔴; the diagnostic column provides details.

checksum_results = validate_checksums(workspace, backend="jax")
pd.DataFrame({
    "Distribution": [result.target for result in checksum_results],
    "Point": [result.point_name for result in checksum_results],
    "Reference": [result.reference for result in checksum_results],
    "Computed": [result.value for result in checksum_results],
    "Difference": [result.difference for result in checksum_results],
    "Status": [
        "🟢" if result.difference < 1e-10 else "🟡" if result.difference < 1e-2 else "🔴"
        for result in checksum_results
    ],
    "Diagnostic": [result.diagnostic for result in checksum_results],
})

	Distribution	Point	Reference	Computed	Difference	Status	Diagnostic
0	default_model	validation_point1	645.687778+  0.000000j	645.687778+  0.000000j	4.547474e-13	🟢	None
1	default_model	validation_point2	33.313132+  0.000000j	33.313132+  0.000000j	2.842171e-14	🟢	None
2	default_model	validation_point3	38.286729+  0.000000j	38.286729+  0.000000j	1.492140e-13	🟢	None
3	L1405_Flatte	validation_point_m12sq	-0.740198+  0.236368j	-0.740198+  0.236368j	2.237726e-16	🟢	None
4	L1890_BW	validation_point_m12sq	2.198693+  1.036641j	2.198693+  1.036641j	2.220446e-16	🟢	None
5	L1800_BW	validation_point_m12sq	0.311734+  2.774200j	0.311734+  2.774200j	2.275957e-15	🟢	None
6	L1520_BW	validation_point_m12sq	-0.643171+  0.553948j	-0.643171+  0.553948j	1.831027e-15	🟢	None
7	L2110_BW	validation_point_m12sq	0.854068+  0.041564j	0.854068+  0.041564j	6.722137e-16	🟢	None
8	L1680_BW	validation_point_m12sq	-2.461330+  0.361098j	-2.461330+  0.361098j	1.042960e-15	🟢	None
9	L2350_BW	validation_point_m12sq	0.430570+  0.000197j	0.430570+  0.000197j	1.665340e-16	🟢	None
10	L1600_BW	validation_point_m12sq	-0.629331+  0.766337j	-0.629331+  0.766337j	6.280370e-16	🟢	None
11	L1810_BW	validation_point_m12sq	0.106436+  5.093986j	0.106436+  5.093986j	3.555992e-15	🟢	None
12	L1820_BW	validation_point_m12sq	4.647408+  4.443945j	4.647408+  4.443945j	1.776357e-14	🟢	None
13	L1830_BW	validation_point_m12sq	3.824528+  3.827750j	3.824528+  3.827750j	4.528839e-15	🟢	None
14	L2100_BW	validation_point_m12sq	0.826225+  0.013532j	0.826225+  0.013532j	1.040834e-17	🟢	None
15	L1690_BW	validation_point_m12sq	-1.859729+  1.396126j	-1.859729+  1.396126j	3.884194e-15	🟢	None
16	LNR30_NR	validation_point_m12sq	1.000000+  0.000000j	1.000000+  0.000000j	0.000000e+00	🟢	None

Visualization

Dalitz plot

The Dalitz plot shows the intensity up to \(m_{pK^-}=2.5\ \mathrm{GeV}\). The plotting range includes a margin below the physical threshold so that the left phase-space boundary remains visible, as in the Julia implementation.

x_subsystem, y_subsystem = 3, 2
resolution = 250
x_min, x_max = 1.9, 2.5**2
y_min, y_max = 5, 27
x, y = jnp.meshgrid(
    jnp.linspace(x_min, x_max, resolution),
    jnp.linspace(y_min, y_max, resolution),
)
dalitz_data = {
    f"sigma{x_subsystem}": x,
    f"sigma{y_subsystem}": y,
}
compiled_workspace = compile_workspace(
    workspace,
    backend="jax",
    targets=[intensity_name],
)
intensity_function = compiled_workspace.functions[intensity_name]
intensities = jnp.real(intensity_function(dalitz_data))
normalized_intensities = intensities / jnp.nansum(intensities)


def get_decay_products(
    decay: ThreeBodyDecay,
    subsystem: FinalStateID,
) -> tuple[State, State]:
    return tuple(state for state in decay.final_state.values() if state.index != subsystem)


sigma_labels = {
    subsystem: Rf"$\sigma_{subsystem} = M^2\left({' '.join(p.latex for p in get_decay_products(decay, subsystem))}\right)$"
    for subsystem in (1, 2, 3)
}
# cspell:ignore bbox edgecolor facecolor fontset fontsize labelcolor mathtext ncol
# cspell:ignore savefig sharey startswith STIX wspace xtick ytick
plt.rcParams.update({
    "axes.edgecolor": "#808080",
    "axes.facecolor": "none",
    "axes.labelcolor": "#808080",
    "figure.facecolor": "none",
    "font.family": "serif",
    "font.serif": ["STIXGeneral", "DejaVu Serif"],
    "font.size": 18,
    "mathtext.fontset": "stix",
    "savefig.transparent": True,
    "text.color": "#808080",
    "xtick.color": "#808080",
    "ytick.color": "#808080",
})

fig, ax = plt.subplots(figsize=(10, 8))
mesh = ax.pcolormesh(x, y, normalized_intensities, rasterized=True)
ax.contour(
    x,
    y,
    jnp.isfinite(intensities),
    levels=[0.5],
    colors="black",
    linewidths=3,
)
ax.set_xlabel(sigma_labels[x_subsystem])
ax.set_ylabel(sigma_labels[y_subsystem])
divider = make_axes_locatable(ax)
colorbar_axes = divider.append_axes("right", size="5%", pad=0.1)
colorbar = fig.colorbar(mesh, cax=colorbar_axes)
colorbar.ax.set_ylabel("Normalized intensity (a.u.)")
plt.show()

Mass projection

The projection onto the \(pK^-\) mass is shown as a black line. The colored curves show the incoherent contribution from each \(\Lambda\) resonance.

resonance_names = sorted({chain["name"] for chain in chain_defs})
mass_values = jnp.sqrt(x[0])
minimum_mass = float(
    (
        sp.Symbol("m1", nonnegative=True) + sp.Symbol("m2", nonnegative=True)
    ).xreplace(model.masses)
)
physical_region = mass_values >= minimum_mass
mass_values = mass_values[physical_region]


def project_to_pk(intensity_values):
    projection = jnp.nansum(intensity_values, axis=0)
    return mass_values * projection[physical_region]


total_projection = project_to_pk(intensities)
normalization = jnp.trapezoid(total_projection, mass_values)

fig, ax = plt.subplots(figsize=(12, 7))
ax.plot(
    mass_values,
    total_projection / normalization,
    color="#202020",
    lw=3,
    label="Total",
)
for resonance_name in resonance_names:
    resonance_latex = to_latex(resonance_name)
    coupling_overrides = {
        symbol: value if resonance_latex in str(symbol) else 0
        for symbol, value in model.parameter_defaults.items()
        if str(symbol).startswith("c^")
    }
    component_workspace = compile_workspace(
        workspace,
        backend="jax",
        targets=[intensity_name],
        parameter_overrides=coupling_overrides,
    )
    component = jnp.real(component_workspace.functions[intensity_name](dalitz_data))
    component = jnp.where(jnp.isnan(intensities), jnp.nan, component)
    label = Rf"${to_latex(resonance_name)}$"
    ax.plot(
        mass_values,
        project_to_pk(component) / normalization,
        label=label,
    )
ax.set_xlabel(R"$m_{pK^-}$ [GeV]")
ax.set_ylabel("Normalized intensity (a.u.)")
ax.set_xlim(minimum_mass, 2.5)
ax.set_ylim(bottom=0)
ax.legend(fontsize="small", loc="upper right", ncol=3)
plt.show()