Construct three-body decay

\(\Lambda_c^+ \to p K^- \pi^+\)

Model definition: lc2ppik-lhcb-2683025.json.

This notebooks illustrates the use of the ampform_dpd.io.serialization module for the decay \(\Lambda_c^+ \to p K^- \pi^+\). The corresponding model was optimized to a data sample of roughly half a million \(\Lambda_c^{\pm}\) decay candidates by the LHCb collaboration, INSPIRE-HEP 2683025.

Warning

The ampform_dpd.io.serialization module is a preview feature. This notebook illustrates the deserialization of the amplitude model JSON file to symbolic expressions. Keep an eye on ComPWA/ampform-dpd#133 for a list of tracked issues.

Import model

from pathlib import Path
import json
import logging
import os

import jax.numpy as jnp
import matplotlib.pyplot as plt
import pandas as pd
import sympy as sp
from ampform.dynamics.form_factor import FormFactor
from ampform.sympy import perform_cached_doit
from ampform_dpd.decay import FinalStateID, State, ThreeBodyDecay
from ampform_dpd.dynamics import (
    BreitWigner,
    ChannelArguments,
    MultichannelBreitWigner,
)
from ampform_dpd.io import (
    aslatex,
    cached,
    perform_cached_lambdify,
    simplify_latex_rendering,
    unfold_definitions,
)
from ampform_dpd.io.serialization.amplitude import (
    HelicityRecoupling,
    LSRecoupling,
    ParityRecoupling,
    formulate,
    formulate_aligned_amplitude,
    formulate_chain_amplitude,
    formulate_recoupling,
)
from ampform_dpd.io.serialization.decay import to_decay
from ampform_dpd.io.serialization.dynamics import (
    formulate_breit_wigner,
    formulate_dynamics,
    formulate_form_factor,
    formulate_multichannel_breit_wigner,
)
from ampform_dpd.io.serialization.format import (
    get_decay_chains,
    get_function_definition,
)
from IPython.display import JSON, Math
from matplotlib_inline.backend_inline import set_matplotlib_formats
from mpl_toolkits.axes_grid1 import make_axes_locatable  # cspell:ignore mpl_toolkits
from tqdm.auto import tqdm

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

with open(THIS_DIR.parent.parent / "models" / "lc2ppik-lhcb-2683025.json") as f:
    MODEL_DEFINITION = json.load(f)

def to_latex(name: str) -> str:
    latex = {
        "Lc": R"\Lambda_c^+",
        "pi": R"\pi^+",
        "K": "K^-",
        "p": "p",
    }.get(name)
    if latex is not None:
        return latex
    mass_str = name[1:].strip("(").strip(")")
    subsystem_letter = name[0]
    subsystem = {"D": "D", "K": "K", "L": R"\Lambda"}.get(subsystem_letter)
    if subsystem is None:
        return name
    return f"{subsystem}({mass_str})"

DECAY = to_decay(MODEL_DEFINITION, to_latex=to_latex)
Math(aslatex(DECAY, with_jp=True))

\(\displaystyle \begin{array}{c} \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1232)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1600)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1700)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(1430)\left[J=0\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(700)\left[J=0\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(892)\left[J=1\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1405)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1520)\left[J=\frac{3}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1600)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1670)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1690)\left[J=\frac{3}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(2000)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \end{array}\)

Dynamics

See also RUB-EP1/amplitude-serialization#22 about serialization of custom lineshapes.

CHAIN_DEFS = get_decay_chains(MODEL_DEFINITION)

Vertices

Blatt-Weisskopf form factor

s, m1, m2, L, d = sp.symbols("s m1 m2 L R", nonnegative=True)
expr = FormFactor(s, m1, m2, L, d)
Math(aslatex(unfold_definitions(expr)))

\(\displaystyle \begin{aligned} \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right) \;&=\; \sqrt{B_{L}^2\left(R^{2} q^2\left(s\right)\right)} \\ B_{L}^2\left(R^{2} q^2\left(s\right)\right) \;&=\; \frac{\left|{h_{L}^{(1)}\left(1\right)}\right|^{2}}{R^{2} \left|{h_{L}^{(1)}\left(R \sqrt{q^2\left(s\right)}\right)}\right|^{2} q^2\left(s\right)} \\ q^2\left(s\right) \;&=\; \frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{4 s} \\ h_{L}^{(1)}\left(z\right) \;&=\; \frac{\left(- i\right)^{L + 1} e^{i z} \sum_{k=0}^{L} \frac{\left(\frac{i}{2 z}\right)^{k} \left(k + L\right)!}{k! \left(- k + L\right)!}}{z} \\ \end{aligned}\)

ff_L1520 = formulate_form_factor(
    vertex=CHAIN_DEFS[2]["vertices"][0],
    model=MODEL_DEFINITION,
)
Math(aslatex(ff_L1520))

\(\displaystyle \begin{array}{rcl} \frac{\sqrt{2} \mathcal{F}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right)}{2} \\ R_{Lc} &=& 5.0 \\ \end{array}\)

Propagators

Breit-Wigner

s, m0, Γ0, m1, m2, L, d = sp.symbols("s m0 Gamma0 m1 m2 L R", nonnegative=True)
expr = BreitWigner(s, m0, Γ0, m1, m2, L, d)
Math(aslatex(unfold_definitions(expr)))

\(\displaystyle \begin{aligned} \mathcal{R}^\mathrm{BW}_{L}\left(s; m_{0}, \Gamma_{0}\right) \;&=\; \mathcal{R}^\mathrm{BW}\left(s; m_{0}, \Gamma_{0}\left(s\right)\right) \\ \mathcal{R}^\mathrm{BW}\left(s; m_{0}, \Gamma_{0}\left(s\right)\right) \;&=\; \frac{1}{m_{0}^{2} - i m_{0} \Gamma_{0}\left(s\right) - s} \\ \Gamma_{0}\left(s\right) \;&=\; \frac{\Gamma_{0} \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right)^{2} \rho\left(s\right)}{\mathcal{F}_{L}\left(m_{0}^{2}, m_{1}, m_{2}\right)^{2} \rho_{0}\left(m_{0}^{2}\right)} \\ \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right) \;&=\; \sqrt{B_{L}^2\left(R^{2} q^2\left(s\right)\right)} \\ \rho\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}}{s} \\ B_{L}^2\left(R^{2} q^2\left(s\right)\right) \;&=\; \frac{\left|{h_{L}^{(1)}\left(1\right)}\right|^{2}}{R^{2} \left|{h_{L}^{(1)}\left(R \sqrt{q^2\left(s\right)}\right)}\right|^{2} q^2\left(s\right)} \\ q^2\left(s\right) \;&=\; \frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{4 s} \\ h_{L}^{(1)}\left(z\right) \;&=\; \frac{\left(- i\right)^{L + 1} e^{i z} \sum_{k=0}^{L} \frac{\left(\frac{i}{2 z}\right)^{k} \left(k + L\right)!}{k! \left(- k + L\right)!}}{z} \\ \end{aligned}\)

K892_BW = formulate_breit_wigner(
    propagator=CHAIN_DEFS[20]["propagators"][0],
    resonance=to_latex(CHAIN_DEFS[20]["name"]),
    model=MODEL_DEFINITION,
)
Math(aslatex(K892_BW))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{K(892)}, \Gamma_{K(892)}\right) &=& \mathcal{R}^\mathrm{BW}\left(\sigma_{1}; m_{K(892)}, \Gamma_{892}\left(\sigma_{1}\right)\right) \\ m_{K(892)} &=& 0.8955 \\ \Gamma_{K(892)} &=& 0.047299999999999995 \\ m_{2} &=& 0.13957018 \\ m_{3} &=& 0.493677 \\ R_\mathrm{res} &=& 1.5 \\ \end{array}\)

Multi-channel Breit-Wigner

The gsq value that is serialized for each channel is the coupling squared, not an energy width. The channel term \(\Gamma^\text{ch}\) below follows the same convention as HadronicLineshapes.jl. See ComPWA/ampform-dpd#199.

s, m0, m1, m2, L, d = sp.symbols("s m0 m1 m2 L R", nonnegative=True)
channels = tuple(
    ChannelArguments(
        s,
        m0,
        coupling_squared=sp.Symbol(f"g_{{{i}}}^2", nonnegative=True),
        m1=sp.Symbol(f"m_{{a,{i}}}", nonnegative=True),
        m2=sp.Symbol(f"m_{{b,{i}}}", nonnegative=True),
        angular_momentum=sp.Symbol(f"L{i}", integer=True, nonnegative=True),
        meson_radius=d,
    )
    for i in [1, 2]
)
expr = MultichannelBreitWigner(s, m0, channels)
Math(aslatex(unfold_definitions(expr)))

\(\displaystyle \begin{aligned} \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(s; g_{1}^2, g_{2}^2\right) \;&=\; \mathcal{R}^\mathrm{BW}_{L=0}\left(s; m_{0}, \Gamma^\text{ch}\left(s; m_{0}, g_{1}^2\right) + \Gamma^\text{ch}\left(s; m_{0}, g_{2}^2\right)\right) \\ \mathcal{R}^\mathrm{BW}_{L=0}\left(s; m_{0}, \Gamma^\text{ch}\left(s; m_{0}, g_{1}^2\right) + \Gamma^\text{ch}\left(s; m_{0}, g_{2}^2\right)\right) \;&=\; \frac{1}{m_{0}^{2} - i m_{0} \left(\Gamma^\text{ch}\left(s; m_{0}, g_{1}^2\right) + \Gamma^\text{ch}\left(s; m_{0}, g_{2}^2\right)\right) - s} \\ \Gamma^\text{ch}\left(s; m_{0}, g_{1}^2\right) \;&=\; \frac{g_{1}^2 \mathcal{F}_{L_{1}}\left(s, m_{a,1}, m_{b,1}\right)^{2} \rho\left(s\right)}{m_{0}} \\ \mathcal{F}_{L_{1}}\left(s, m_{a,1}, m_{b,1}\right) \;&=\; \sqrt{B_{L_{1}}^2\left(R^{2} q^2\left(s\right)\right)} \\ \rho\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{a,1} - m_{b,1}\right)^{2}\right) \left(s - \left(m_{a,1} + m_{b,1}\right)^{2}\right)}}{s} \\ B_{L_{1}}^2\left(R^{2} q^2\left(s\right)\right) \;&=\; \frac{\left|{h_{L_{1}}^{(1)}\left(1\right)}\right|^{2}}{R^{2} \left|{h_{L_{1}}^{(1)}\left(R \sqrt{q^2\left(s\right)}\right)}\right|^{2} q^2\left(s\right)} \\ q^2\left(s\right) \;&=\; \frac{\left(s - \left(m_{a,1} - m_{b,1}\right)^{2}\right) \left(s - \left(m_{a,1} + m_{b,1}\right)^{2}\right)}{4 s} \\ h_{L_{1}}^{(1)}\left(z\right) \;&=\; \frac{\left(- i\right)^{L_{1} + 1} e^{i z} \sum_{k=0}^{L_{1}} \frac{\left(\frac{i}{2 z}\right)^{k} \left(k + L_{1}\right)!}{k! \left(- k + L_{1}\right)!}}{z} \\ \end{aligned}\)

L1405_Flatte = formulate_multichannel_breit_wigner(
    propagator=CHAIN_DEFS[0]["propagators"][0],
    resonance=to_latex(CHAIN_DEFS[0]["name"]),
    model=MODEL_DEFINITION,
)
Math(aslatex(L1405_Flatte))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(\sigma_{2}; \Gamma_{\Lambda(1405)}, \Gamma_{\Lambda(1405)}^\text{ch. 2}\right) &=& \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{\Lambda(1405)}, \Gamma^\text{ch}\left(\sigma_{2}; m_{\Lambda(1405)}, \Gamma_{\Lambda(1405)}\right) + \Gamma^\text{ch}\left(\sigma_{2}; m_{\Lambda(1405)}, \Gamma_{\Lambda(1405)}^\text{ch. 2}\right)\right) \\ m_{\Lambda(1405)} &=& 1.4051 \\ \Gamma_{\Lambda(1405)} &=& 0.328725260215546 \\ m_{3} &=& 0.938272046 \\ m_{1} &=& 0.493677 \\ R_{\Lambda(1405)} &=& 0 \\ m_{a,2} &=& 1.18937 \\ m_{b,2} &=& 0.13957018 \\ \Gamma_{\Lambda(1405)}^\text{ch. 2} &=& 0.328725260215546 \\ \end{array}\)

Breit-Wigner with exponential

The Bugg lineshapes are serialized as generic_function expression strings. The built-in dynamics builder parses these expressions and substitutes \(i\) and \(\sigma\) with the imaginary unit and the Mandelstam variable of the propagator node, respectively.

get_function_definition("K700_BuggBW", MODEL_DEFINITION)

{'name': 'K700_BuggBW',
 'type': 'generic_function',
 'expression': '1/(0.824^2 - σ - i * 0.824 * (σ - 0.23397706275638377) / (0.824^2 - 0.23397706275638377) * 0.478 * exp(-0.941060 * σ))'}

K700_BuggBW = formulate_dynamics(CHAIN_DEFS[18], MODEL_DEFINITION, to_latex)
Math(aslatex(K700_BuggBW))

\(\displaystyle \begin{array}{rcl} \frac{1}{- \sigma_{1} - 0.88510773180649964 i \left(\sigma_{1} - 0.23397706275638377\right) e^{- 0.94106 \sigma_{1}} + 0.678976} \\ \end{array}\)

Construct amplitude model

Unpolarized intensity

λ0, λ1, λ2, λ3 = sp.symbols("lambda(:4)", rational=True)
amplitude_expr, _ = formulate_aligned_amplitude(MODEL_DEFINITION, λ0, λ1, λ2, λ3)
amplitude_expr.cleanup()

\(\displaystyle \sum_{\lambda_0^{\prime}=-1/2}^{1/2} \sum_{\lambda_1^{\prime}=-1/2}^{1/2}{A^{1}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{1(1)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{1(1)}\right) + A^{2}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{2(1)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{2(1)}\right) + A^{3}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{3(1)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{3(1)}\right)}\)

Amplitude for the decay chain

Helicity recouplings

λa = sp.Symbol(R"\lambda_a", rational=True)
λb = sp.Symbol(R"\lambda_b", rational=True)
λa0 = sp.Symbol(R"\lambda_a^0", rational=True)
λb0 = sp.Symbol(R"\lambda_b^0", rational=True)
f = sp.Symbol("f", integer=True)
l = sp.Symbol("l", integer=True, nonnegative=True)
s = sp.Symbol("s", nonnegative=True, rational=True)
ja = sp.Symbol("j_a", nonnegative=True, rational=True)
jb = sp.Symbol("j_b", nonnegative=True, rational=True)
j = sp.Symbol("j", nonnegative=True, rational=True)
exprs = [
    HelicityRecoupling(λa, λb, λa0, λb0),
    ParityRecoupling(λa, λb, λa0, λb0, f),
    LSRecoupling(λa, λb, l, s, ja, jb, j),
]
Math(aslatex({e: e.doit(deep=False) for e in exprs}))

\(\displaystyle \begin{aligned} \mathcal{H}^\text{helicity}\left(\lambda_{a},\lambda_{b}\middle|\lambda^{0}_{a},\lambda^{0}_{b}\right) \;&=\; \delta_{\lambda_{a} \lambda^{0}_{a}} \delta_{\lambda_{b} \lambda^{0}_{b}} \\ \mathcal{H}^\text{parity}\left(\lambda_{a},\lambda_{b}\middle|\lambda^{0}_{a},\lambda^{0}_{b},f\right) \;&=\; f \delta_{\lambda_{a}, - \lambda^{0}_{a}} \delta_{\lambda_{b}, - \lambda^{0}_{b}} + \delta_{\lambda_{a} \lambda^{0}_{a}} \delta_{\lambda_{b} \lambda^{0}_{b}} \\ \mathcal{H}^\text{parity}\left(\lambda_{a},\lambda_{b}\middle|l,s,j_{a},j_{b},j\right) \;&=\; \frac{\sqrt{2 l + 1} C^{s,\lambda_{a} - \lambda_{b}}_{j_{a},\lambda_{a},j_{b},- \lambda_{b}} C^{j,\lambda_{a} - \lambda_{b}}_{l,0,s,\lambda_{a} - \lambda_{b}}}{\sqrt{2 j + 1}} \\ \end{aligned}\)

Recoupling deserialization

recouplings = [
    formulate_recoupling(MODEL_DEFINITION, chain_idx=0, vertex_idx=i) for i in range(2)
]
Math(aslatex({e: e.doit(deep=False) for e in recouplings}))

\(\displaystyle \begin{aligned} \mathcal{H}^\text{helicity}\left(\lambda_{R},\lambda_{2}\middle|\frac{1}{2},0\right) \;&=\; \delta_{0 \lambda_{2}} \delta_{\frac{1}{2} \lambda_{R}} \\ \mathcal{H}^\text{parity}\left(\lambda_{3},\lambda_{1}\middle|0,\frac{1}{2},1\right) \;&=\; \delta_{- \frac{1}{2} \lambda_{1}} \delta_{0 \lambda_{3}} + \delta_{0 \lambda_{3}} \delta_{\frac{1}{2} \lambda_{1}} \\ \end{aligned}\)

Chain amplitudes

definitions = formulate_chain_amplitude(λ0, λ1, λ2, λ3, MODEL_DEFINITION, chain_idx=0)
Math(aslatex(definitions))

\(\displaystyle \begin{aligned} A^{2}_{\lambda_{0}, \lambda_{1}, \lambda_{2}, \lambda_{3}} \;&=\; \sum_{\lambda_{R}=-1/2}^{1/2}{\left(-1\right)^{- \lambda_{2}} \left(-1\right)^{\frac{1}{2} - \lambda_{1}} \sqrt{2} c^{L1405[1/2]}_{\frac{1}{2}, 0, 0} \delta_{\lambda_{0}, \lambda_{R} - \lambda_{2}} \mathcal{H}^\text{helicity}\left(\lambda_{R},\lambda_{2}\middle|\frac{1}{2},0\right) \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(\sigma_{2}; \Gamma_{L1405}, \Gamma_{L1405}^\text{ch. 2}\right) \mathcal{H}^\text{parity}\left(\lambda_{3},\lambda_{1}\middle|0,\frac{1}{2},1\right) d^{\frac{1}{2}}_{\lambda_{R},- \lambda_{1} + \lambda_{3}}\left(\theta_{31}\right)} \\ c^{L1405[1/2]}_{\frac{1}{2}, 0, 0} \;&=\; 7.38649400481717+1.971018433257411i \\ m_{L1405} \;&=\; 1.4051 \\ \Gamma_{L1405} \;&=\; 0.328725260215546 \\ m_{3} \;&=\; 0.938272046 \\ m_{1} \;&=\; 0.493677 \\ R_{L1405} \;&=\; 0 \\ m_{a,2} \;&=\; 1.18937 \\ m_{b,2} \;&=\; 0.13957018 \\ \Gamma_{L1405}^\text{ch. 2} \;&=\; 0.328725260215546 \\ \theta_{31} \;&=\; \operatorname{acos}{\left(\frac{2 \sigma_{2} \left(- m_{2}^{2} - m_{3}^{2} + \sigma_{1}\right) - \left(m_{0}^{2} - m_{2}^{2} - \sigma_{2}\right) \left(- m_{1}^{2} + m_{3}^{2} + \sigma_{2}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{2}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{3}^{2}, m_{1}^{2}\right)}} \right)} \\ \end{aligned}\)

Validation

checksum_points = {
    point["name"]: {par["name"]: par["value"] for par in point["parameters"]}
    for point in MODEL_DEFINITION["parameter_points"]
}

The table lists the validation checks and their status. The marks “🟢”, “🟡”, and “🔴” indicate an accuracy of \(<10^{-10}\), \(<10^{-2}\), or \(\ge10^{-2}\), respectively, for the difference between the reference and computed values.

def to_complex(value: float | str) -> complex:
    if isinstance(value, str):
        return complex(value.replace(" ", "").replace("i", "j"))
    return complex(value)


def label_diff(difference: complex) -> str:
    absolute_difference = abs(difference)
    if absolute_difference < 1e-10:
        return "🟢"
    if absolute_difference < 1e-2:
        return "🟡"
    return "🔴"


chains_by_propagator = {
    propagator["parametrization"]: chain
    for chain in CHAIN_DEFS
    for propagator in chain["propagators"]
}
validation_results = []
for checksum in MODEL_DEFINITION["misc"]["amplitude_model_checksums"]:
    distribution_name = checksum["distribution"]
    chain = chains_by_propagator.get(distribution_name)
    if chain is None:
        continue
    dynamics = formulate_dynamics(chain, MODEL_DEFINITION)
    function = perform_cached_lambdify(
        dynamics.expression.doit(),
        parameters=dynamics.parameters,
    )
    point = checksum_points[checksum["point"]]
    variable = next(iter(dynamics.expression.free_symbols - dynamics.parameters.keys()))
    computed_value = function({str(variable): next(iter(point.values()))})
    reference_value = to_complex(checksum["value"])
    difference = abs(reference_value - computed_value)
    assert difference < 1e-10, (
        f"{distribution_name} at {checksum['point']} differs by {difference}"
    )
    validation_results.append({
        "Distribution": distribution_name,
        "Point": checksum["point"],
        "Status": label_diff(reference_value - computed_value),
    })

def evaluate_intensity(model, point: dict[str, float]) -> float:
    intensity_function = create_intensity_function(model, subsystem=3)
    sigma2_value = point["m_31"] ** 2
    sigma1_value = solve_sigma1_from_cos_theta31(
        model,
        cos_theta31=point["cos_theta_31"],
        sigma2_value=sigma2_value,
    )
    return float(
        intensity_function({
            "sigma1": sigma1_value,
            "sigma2": sigma2_value,
        })
    )


def create_intensity_function(
    model,
    subsystem: int,
    resonance_name: str | None = None,
):
    invariant = next(s for s in model.invariants if str(s) == f"sigma{subsystem}")
    intensity_expr = cached.xreplace(cached.unfold(model), model.variables)
    parameter_defaults = select_resonance_parameters(model, resonance_name)
    intensity_expr = cached.xreplace(intensity_expr, parameter_defaults)
    invariant_expr = model.invariants[invariant].xreplace(model.masses).doit()
    intensity_expr = cached.doit(intensity_expr.xreplace({invariant: invariant_expr}))
    return cached.lambdify(intensity_expr, backend="jax")


def select_resonance_parameters(model, resonance_name: str | None):
    if resonance_name is None:
        return model.parameter_defaults
    return {
        symbol: value
        if not str(symbol).startswith("c^") or resonance_name in str(symbol)
        else 0
        for symbol, value in model.parameter_defaults.items()
    }


def solve_sigma1_from_cos_theta31(
    model,
    cos_theta31: float,
    sigma2_value: float,
) -> float:
    sigma1, sigma2, _ = sorted(model.invariants, key=str)
    theta31_expr = next(
        expression
        for angle, expression in model.variables.items()
        if str(angle) == "theta_31"
    )
    solutions = sp.solve(
        sp.Eq(
            sp.cos(theta31_expr).xreplace(model.parameter_defaults),
            cos_theta31,
        ),
        sigma1,
    )
    return float(solutions[0].xreplace({sigma2: sigma2_value}))


intensity_checksum = next(
    checksum
    for checksum in MODEL_DEFINITION["misc"]["amplitude_model_checksums"]
    if checksum["distribution"] == MODEL_DEFINITION["distributions"][0]["name"]
)
model = formulate(MODEL_DEFINITION, cleanup_summations=True, to_latex=to_latex)
point = checksum_points[intensity_checksum["point"]]
computed_value = evaluate_intensity(model, point)
reference_value = to_complex(intensity_checksum["value"])
difference = abs(reference_value - computed_value)
assert difference < 1e-10, (
    f"default_model at {intensity_checksum['point']} differs by {difference}"
)
validation_results.append({
    "Distribution": intensity_checksum["distribution"],
    "Point": intensity_checksum["point"],
    "Status": label_diff(reference_value - computed_value),
})
pd.DataFrame(validation_results)

	Distribution	Point	Status
0	L1405_Flatte	validation_point_m31sq	🟢
1	L1690_BW	validation_point_m31sq	🟢
2	D1232_BW	validation_point_m12sq	🟢
3	L1520_BW	validation_point_m31sq	🟢
4	L1600_BW	validation_point_m31sq	🟢
5	L2000_BW	validation_point_m31sq	🟢
6	D1600_BW	validation_point_m12sq	🟢
7	D1700_BW	validation_point_m12sq	🟢
8	K892_BW	validation_point_m23sq	🟢
9	L1670_BW	validation_point_m31sq	🟢
10	default_model	validation_point	🟢

Visualization

i, j = 2, 1
k, *_ = {1, 2, 3} - {i, j}
resolution = 1_000
masses = sorted(model.masses, key=str)
x_min = float(((masses[j] + masses[k]) ** 2).xreplace(model.masses))
x_max = float(((masses[0] - masses[i]) ** 2).xreplace(model.masses))
y_min = float(((masses[i] + masses[k]) ** 2).xreplace(model.masses))
y_max = float(((masses[0] - masses[j]) ** 2).xreplace(model.masses))
x_margin = 0.05 * (x_max - x_min)
y_margin = 0.05 * (y_max - y_min)
X, Y = jnp.meshgrid(
    jnp.linspace(x_min - x_margin, x_max + x_margin, num=resolution),
    jnp.linspace(y_min - y_margin, y_max + y_margin, num=resolution),
)
intensity_function = create_intensity_function(model, subsystem=k)
intensities = intensity_function({f"sigma{i}": X, f"sigma{j}": Y})
normalized_intensities = intensities / jnp.nansum(intensities)


def get_decay_products(
    decay: ThreeBodyDecay,
    subsystem: FinalStateID,
) -> tuple[State, State]:
    return tuple(s for s in decay.final_state.values() if s.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 savefig
# cspell:ignore 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.set_aspect("equal")
ax.set_xlabel(sigma_labels[i])
ax.set_ylabel(sigma_labels[j])
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()

resonance_names = sorted({chain["name"].split("[")[0] for chain in CHAIN_DEFS})


def project(intensity_values, subsystem: int):
    if subsystem == i:
        mass_values = jnp.sqrt(X[0])
        projection = jnp.nansum(intensity_values, axis=0)
    else:
        mass_values = jnp.sqrt(Y[:, 0])
        projection = jnp.nansum(intensity_values, axis=1)
    return mass_values, mass_values * projection


projections = {}
normalizations = {}
for subsystem in (i, j):
    mass_values, total_projection = project(intensities, subsystem)
    projections[subsystem] = (mass_values, total_projection)
    normalizations[subsystem] = jnp.trapezoid(total_projection, mass_values)

fig, axes = plt.subplots(1, 2, figsize=(14, 6), sharey=True)
for axis, subsystem in zip(axes, (i, j), strict=True):
    mass_values, total_projection = projections[subsystem]
    axis.plot(
        mass_values,
        total_projection / normalizations[subsystem],
        color="#808080",
        lw=3,
        label="Total",
    )
for resonance_name in resonance_names:
    component_function = create_intensity_function(
        model,
        subsystem=k,
        resonance_name=to_latex(resonance_name),
    )
    component = component_function({"sigma1": Y, "sigma2": X})
    component = jnp.where(jnp.isnan(intensities), jnp.nan, component)
    for axis, subsystem in zip(axes, (i, j), strict=True):
        mass_values, component_projection = project(component, subsystem)
        axis.plot(
            mass_values,
            component_projection / normalizations[subsystem],
            label=resonance_name,
        )
axes[0].set_xlabel(R"$m_{13}$ [GeV]")
axes[1].set_xlabel(R"$m_{23}$ [GeV]")
axes[0].set_ylabel("Normalized intensity (a.u.)")
axes[0].set_ylim(bottom=0)
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(
    handles,
    labels,
    bbox_to_anchor=(0.5, 0.97),
    fontsize="small",
    loc="upper center",
    ncol=7,
)
fig.subplots_adjust(top=0.8, wspace=0.05)
plt.show()