Amplitude model with ampform-dpd

PWA study on \(p \gamma \to \Lambda K^+ \pi^0\). We formulate the helicity amplitude model symbolically using AmpForm-DPD here.

Import Python libraries

Hide code cell source

import logging
import os
import warnings
from collections import defaultdict
from fractions import Fraction
from textwrap import dedent

import ampform
import graphviz
import ipywidgets as w
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import qrules
import sympy as sp
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from ampform.io import improve_latex_rendering
from ampform_dpd import DalitzPlotDecompositionBuilder
from ampform_dpd.adapter.qrules import normalize_state_ids, to_three_body_decay
from ampform_dpd.decay import DecayNode, ThreeBodyDecayChain
from ampform_dpd.dynamics.builder import create_mass_symbol, get_mandelstam_s
from ampform_dpd.io import aslatex
from IPython.display import SVG, Image, Markdown, Math, display
from qrules.particle import Particle, Spin, create_particle, load_pdg
from tensorwaves.data import SympyDataTransformer
from tensorwaves.function.sympy import create_parametrized_function

STATIC_PAGE = "EXECUTE_NB" in os.environ

os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
logging.disable(logging.WARNING)
warnings.filterwarnings("ignore")

improve_latex_rendering()
particle_db = load_pdg()

Decay definition

Particle definitions

Show code cell source

Hide code cell source

def generate_markdown_table(particles: list[str]):
    src = dedent(r"""
    | Particle | Name | PID | $J^{PC} (I^G)$ | $I_3$ | $M$ | $\Gamma$ | $Q$ | $S$ | $B$ |
    | :------- |------|-----|----------------|-------|-----|----------|-----|-----|-----|
    """)
    for name in particles:
        p = particle_db[name]
        src += f"| ${p.latex}$ | `{p.name}` | {p.pid} | {jpc_ig(p)} | {i_3(p)} |  {p.mass:.3g}| {p.width:g} | {p.charge} |{p.strangeness} | {p.baryon_number}|\n"
    return src


def jpc_ig(particle: Particle) -> str:
    j = format_fraction(particle.spin)
    p = format_parity(particle.parity)
    c = format_parity(particle.c_parity)
    if particle.isospin is None:
        return f"${j}^{{{p}{c}}}$"
    i = format_fraction(particle.isospin.magnitude)
    g = format_parity(particle.g_parity)
    return rf"${j}^{{{p}{c}}} \; ({i}^{{{g}}})$"


def i_3(particle: Particle) -> str:
    if particle.isospin is None:
        return "N/A"
    return f"${format_fraction(particle.isospin.projection)}$"


def format_fraction(value: float) -> str:
    value = Fraction(value)
    if value.denominator == 1:
        return str(value.numerator)
    return rf"\frac{{{value.numerator}}}{{{value.denominator}}}"


def format_parity(parity: int | None) -> str:
    if parity is None:
        return " "
    if parity == -1:
        return "-"
    if parity == 1:
        return "+"
    raise NotImplementedError


particles = ["Lambda", "K+", "pi0", "gamma", "p"]
src = generate_markdown_table(particles)
Markdown(src)

Particle	Name	PID	\(J^{PC} (I^G)\)	\(I_3\)	\(M\)	\(\Gamma\)	\(Q\)	\(S\)	\(B\)
\(\Lambda\)	Lambda	3122	\(\frac{1}{2}^{+ } \; (0^{ })\)	\(0\)	1.12	2.515e-15	0	-1	1
\(K^{+}\)	K+	321	\(0^{- } \; (\frac{1}{2}^{ })\)	\(\frac{1}{2}\)	0.494	5.317e-17	1	1	0
\(\pi^{0}\)	pi0	111	\(0^{-+} \; (1^{-})\)	\(0\)	0.135	7.81e-09	0	0	0
\(\gamma\)	gamma	22	\(1^{--}\)	N/A	0	0	0	0	0
\(p\)	p	2212	\(\frac{1}{2}^{+ } \; (\frac{1}{2}^{ })\)	\(\frac{1}{2}\)	0.938	0	1	0	1

In the table above, PID is the PDG ID from PDG particle numbering scheme, \(J\) is the spin, \(P\) is the parity, \(C\) is the C parity, \(I\) is the isospin (magnitude), \(G\) is the G parity. \(I_3\) is the isospin projection (or the 3rd component), \(M\) is the mass, \(\Gamma\) is the width, \(Q\) is the charge, \(S\) is the strangeness number, and \(B\) is the baryon number.

Initial state definition

Mass for \(p \gamma\) system

E_lab_gamma = 8.5
m_proton = 0.938
m_0 = np.sqrt(2 * E_lab_gamma * m_proton + m_proton**2)
m_eta = 0.548
m_pi = 0.135
m_0

4.101931740046389

Add custom particle \(p \gamma\)

pgamma1 = Particle(
    name="pgamma1",
    latex=r"p\gamma (s1/2)",
    spin=0.5,
    mass=m_0,
    charge=1,
    isospin=Spin(1 / 2, +1 / 2),
    baryon_number=1,
    parity=-1,
    pid=99990,
)
pgamma2 = create_particle(
    template_particle=pgamma1,
    name="pgamma2",
    latex=R"p\gamma (s3/2)",
    spin=1.5,
    pid=pgamma1.pid + 1,
)
particle_db.update([pgamma1, pgamma2])

Generate transitions

For simplicity, we use the initial state \(p \gamma\) (with spin-\(\frac{1}{2}\)), and set the allowed interaction type to be strong only, the formalism is selected to be helicity formalism instead of canonical.

See also

Helicity versus canonical

reaction = qrules.generate_transitions(
    initial_state=("pgamma1"),
    final_state=["Lambda", "K+", "pi0"],
    allowed_interaction_types=["strong"],
    formalism="helicity",
    particle_db=particle_db,
    max_angular_momentum=4,
    max_spin_magnitude=4,
    mass_conservation_factor=0,
)
reaction = normalize_state_ids(reaction)

Show code cell source

Hide code cell source

dot = qrules.io.asdot(reaction, collapse_graphs=True)
graphviz.Source(dot)

decay = to_three_body_decay(reaction.transitions)
Math(aslatex(decay, with_jp=True))

\[\begin{split}\displaystyle \begin{array}{c} p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{0}^{*}(1430)^{+}\left[0^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{0}^{*}(1950)^{+}\left[0^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{0}^{*}(700)^{+}\left[0^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{2}^{*}(1430)^{+}\left[2^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{2}^{*}(1980)^{+}\left[2^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{3}^{*}(1780)^{+}\left[3^-\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K_{4}^{*}(2045)^{+}\left[4^+\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K^{*}(1410)^{+}\left[1^-\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K^{*}(1680)^{+}\left[1^-\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(K^{*}(892)^{+}\left[1^-\right] \to K^{+}\left[0^-\right] \pi^{0}\left[0^-\right]\right) \Lambda\left[\frac{1}{2}^+\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1650)^{+}\left[\frac{1}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1675)^{+}\left[\frac{5}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1680)^{+}\left[\frac{5}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1700)^{+}\left[\frac{3}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1710)^{+}\left[\frac{1}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(1720)^{+}\left[\frac{3}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(N(2190)^{+}\left[\frac{7}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] K^{+}\left[0^-\right]\right) \pi^{0}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1385)^{0}\left[\frac{3}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1660)^{0}\left[\frac{1}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1670)^{0}\left[\frac{3}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1750)^{0}\left[\frac{1}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1775)^{0}\left[\frac{5}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1940)^{0}\left[\frac{3}{2}^-\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(1915)^{0}\left[\frac{5}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ p\gamma (s1/2)\left[\frac{1}{2}^-\right] \to \left(\Sigma(2030)^{0}\left[\frac{7}{2}^+\right] \to \Lambda\left[\frac{1}{2}^+\right] \pi^{0}\left[0^-\right]\right) K^{+}\left[0^-\right] \\ \end{array}\end{split}\]

Formulate amplitude model

model_builder = ampform.get_builder(reaction)
model_builder.config.scalar_initial_state_mass = True
model_builder.config.stable_final_state_ids = list(reaction.final_state)
bw_builder = RelativisticBreitWignerBuilder(
    energy_dependent_width=False,
    form_factor=False,
)
for name in reaction.get_intermediate_particles().names:
    model_builder.dynamics.assign(name, bw_builder)
model = model_builder.formulate()

model.intensity

\[\displaystyle \sum_{m_{0}=-1/2}^{1/2} \sum_{m_{1}=-1/2}^{1/2} \sum_{m_{2}=0} \sum_{m_{3}=0}{\left|{A^{12}_{m_{0}, m_{1}, m_{2}, m_{3}} + A^{13}_{m_{0}, m_{1}, m_{2}, m_{3}} + A^{23}_{m_{0}, m_{1}, m_{2}, m_{3}}}\right|^{2}}\]

def formulate_breit_wigner(
    decay: ThreeBodyDecayChain,
) -> tuple[sp.Expr, dict[sp.Symbol, complex | float]]:
    decay_node = decay.decay_node
    s = get_mandelstam_s(decay_node)
    parameter_defaults = {}
    breit_wigner, new_pars = _create_breit_wigner(s, decay_node)
    parameter_defaults.update(new_pars)
    return (
        breit_wigner,
        parameter_defaults,
    )


def _create_breit_wigner(
    s: sp.Symbol, isobar: DecayNode
) -> tuple[sp.Expr, dict[sp.Symbol, complex | float]]:
    mass = create_mass_symbol(isobar.parent)
    width = sp.Symbol(Rf"\Gamma_{{{isobar.parent.latex}}}", nonnegative=True)
    breit_wigner_expr = (mass * width) / (mass**2 - s - width * mass * sp.I)
    parameter_defaults: dict[sp.Symbol, complex | float] = {
        mass: isobar.parent.mass,
        width: isobar.parent.width,
    }
    return breit_wigner_expr, parameter_defaults

model_builder = DalitzPlotDecompositionBuilder(decay)
for chain in model_builder.decay.chains:
    model_builder.dynamics_choices.register_builder(chain, formulate_breit_wigner)
model = model_builder.formulate(reference_subsystem=1)

model.intensity

\[\displaystyle \sum_{\lambda_{0}=-1/2}^{1/2} \sum_{\lambda_{1}=-1/2}^{1/2} \sum_{\lambda_{2}=0} \sum_{\lambda_{3}=0}{\left|{\sum_{\lambda_0^{\prime}=-1/2}^{1/2} \sum_{\lambda_1^{\prime}=-1/2}^{1/2} \sum_{\lambda_2^{\prime}=0} \sum_{\lambda_3^{\prime}=0}{A^{1}_{\lambda_0^{\prime}, \lambda_1^{\prime}, \lambda_2^{\prime}, \lambda_3^{\prime}} 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}, \lambda_2^{\prime}, \lambda_3^{\prime}} 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}, \lambda_2^{\prime}, \lambda_3^{\prime}} 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)}}\right|^{2}}\]

The first term in the amplitude model:

Show code cell source

Hide code cell source

(symbol, expr), *_ = model.amplitudes.items()
Math(aslatex({symbol: expr}, terms_per_line=1))

\[\begin{split}\displaystyle \begin{array}{rcl} A^{1}_{- \frac{1}{2}, - \frac{1}{2}, 0, 0} &=& \sum_{\lambda_{R}=-1}^{1}{- \frac{\Gamma_{K^{*}(1410)^{+}} m_{K^{*}(1410)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(1410)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(1410)^{+}, \lambda_{R}, - \frac{1}{2}} d^{1}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K^{*}(1410)^{+}} m_{K^{*}(1410)^{+}} + \left(m_{K^{*}(1410)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-1}^{1}{- \frac{\Gamma_{K^{*}(1680)^{+}} m_{K^{*}(1680)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(1680)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(1680)^{+}, \lambda_{R}, - \frac{1}{2}} d^{1}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K^{*}(1680)^{+}} m_{K^{*}(1680)^{+}} + \left(m_{K^{*}(1680)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-1}^{1}{- \frac{\Gamma_{K^{*}(892)^{+}} m_{K^{*}(892)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{+}, \lambda_{R}, - \frac{1}{2}} d^{1}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K^{*}(892)^{+}} m_{K^{*}(892)^{+}} + \left(m_{K^{*}(892)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=0}{- \frac{\Gamma_{K_{0}^{*}(1430)^{+}} m_{K_{0}^{*}(1430)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(1430)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{+}, \lambda_{R}, - \frac{1}{2}} d^{0}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{0}^{*}(1430)^{+}} m_{K_{0}^{*}(1430)^{+}} + \left(m_{K_{0}^{*}(1430)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=0}{- \frac{\Gamma_{K_{0}^{*}(1950)^{+}} m_{K_{0}^{*}(1950)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(1950)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1950)^{+}, \lambda_{R}, - \frac{1}{2}} d^{0}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{0}^{*}(1950)^{+}} m_{K_{0}^{*}(1950)^{+}} + \left(m_{K_{0}^{*}(1950)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=0}{- \frac{\Gamma_{K_{0}^{*}(700)^{+}} m_{K_{0}^{*}(700)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(700)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{0}^{*}(700)^{+}, \lambda_{R}, - \frac{1}{2}} d^{0}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{0}^{*}(700)^{+}} m_{K_{0}^{*}(700)^{+}} + \left(m_{K_{0}^{*}(700)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-2}^{2}{- \frac{\Gamma_{K_{2}^{*}(1430)^{+}} m_{K_{2}^{*}(1430)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{2}^{*}(1430)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{+}, \lambda_{R}, - \frac{1}{2}} d^{2}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{2}^{*}(1430)^{+}} m_{K_{2}^{*}(1430)^{+}} + \left(m_{K_{2}^{*}(1430)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-2}^{2}{- \frac{\Gamma_{K_{2}^{*}(1980)^{+}} m_{K_{2}^{*}(1980)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{2}^{*}(1980)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1980)^{+}, \lambda_{R}, - \frac{1}{2}} d^{2}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{2}^{*}(1980)^{+}} m_{K_{2}^{*}(1980)^{+}} + \left(m_{K_{2}^{*}(1980)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-3}^{3}{- \frac{\Gamma_{K_{3}^{*}(1780)^{+}} m_{K_{3}^{*}(1780)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{3}^{*}(1780)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{3}^{*}(1780)^{+}, \lambda_{R}, - \frac{1}{2}} d^{3}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{3}^{*}(1780)^{+}} m_{K_{3}^{*}(1780)^{+}} + \left(m_{K_{3}^{*}(1780)^{+}}\right)^{2} - \sigma_{1}}} \\ &+& \sum_{\lambda_{R}=-4}^{4}{- \frac{\Gamma_{K_{4}^{*}(2045)^{+}} m_{K_{4}^{*}(2045)^{+}} \delta_{- \frac{1}{2}, \lambda_{R} + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K_{4}^{*}(2045)^{+}, 0, 0} \mathcal{H}^\mathrm{production}_{K_{4}^{*}(2045)^{+}, \lambda_{R}, - \frac{1}{2}} d^{4}_{\lambda_{R},0}\left(\theta_{23}\right)}{- i \Gamma_{K_{4}^{*}(2045)^{+}} m_{K_{4}^{*}(2045)^{+}} + \left(m_{K_{4}^{*}(2045)^{+}}\right)^{2} - \sigma_{1}}} \\ \end{array}\end{split}\]

Model parameters

Hide code cell source

if STATIC_PAGE:
    src = aslatex(model.parameter_defaults)
    display(Math(src))

\[\begin{split}\displaystyle \begin{array}{rcl} m_{K_{0}^{*}(1430)^{+}} &=& 1.43 \\ \Gamma_{K_{0}^{*}(1430)^{+}} &=& 0.27 \\ m_{K_{0}^{*}(1950)^{+}} &=& 1.957 \\ \Gamma_{K_{0}^{*}(1950)^{+}} &=& 0.17 \\ m_{K_{0}^{*}(700)^{+}} &=& 0.845 \\ \Gamma_{K_{0}^{*}(700)^{+}} &=& 0.468 \\ m_{K_{2}^{*}(1430)^{+}} &=& 1.4273 \\ \Gamma_{K_{2}^{*}(1430)^{+}} &=& 0.1 \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{2}^{*}(1430)^{+}, 0, 0} &=& 1 \\ m_{K_{2}^{*}(1980)^{+}} &=& 1.99 \\ \Gamma_{K_{2}^{*}(1980)^{+}} &=& 0.348 \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1980)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{2}^{*}(1980)^{+}, 0, 0} &=& 1 \\ m_{K_{3}^{*}(1780)^{+}} &=& 1.779 \\ \Gamma_{K_{3}^{*}(1780)^{+}} &=& 0.161 \\ \mathcal{H}^\mathrm{production}_{K_{3}^{*}(1780)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{3}^{*}(1780)^{+}, 0, 0} &=& 1 \\ m_{K_{4}^{*}(2045)^{+}} &=& 2.048 \\ \Gamma_{K_{4}^{*}(2045)^{+}} &=& 0.199 \\ \mathcal{H}^\mathrm{production}_{K_{4}^{*}(2045)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{4}^{*}(2045)^{+}, 0, 0} &=& 1 \\ m_{K^{*}(1410)^{+}} &=& 1.414 \\ \Gamma_{K^{*}(1410)^{+}} &=& 0.232 \\ \mathcal{H}^\mathrm{production}_{K^{*}(1410)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K^{*}(1410)^{+}, 0, 0} &=& 1 \\ m_{K^{*}(1680)^{+}} &=& 1.718 \\ \Gamma_{K^{*}(1680)^{+}} &=& 0.32 \\ \mathcal{H}^\mathrm{production}_{K^{*}(1680)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K^{*}(1680)^{+}, 0, 0} &=& 1 \\ m_{K^{*}(892)^{+}} &=& 0.89167 \\ \Gamma_{K^{*}(892)^{+}} &=& 0.0514 \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{+}, -1, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{+}, 0, 0} &=& 1 \\ m_{\Sigma(1385)^{0}} &=& 1.3837000000000002 \\ \Gamma_{\Sigma(1385)^{0}} &=& 0.036 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1385)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1385)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1660)^{0}} &=& 1.66 \\ \Gamma_{\Sigma(1660)^{0}} &=& 0.2 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1660)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1660)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1670)^{0}} &=& 1.675 \\ \Gamma_{\Sigma(1670)^{0}} &=& 0.07 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1670)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1670)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1750)^{0}} &=& 1.75 \\ \Gamma_{\Sigma(1750)^{0}} &=& 0.15 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1750)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1750)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1775)^{0}} &=& 1.775 \\ \Gamma_{\Sigma(1775)^{0}} &=& 0.12 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1775)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1775)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1940)^{0}} &=& 1.91 \\ \Gamma_{\Sigma(1940)^{0}} &=& 0.22 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1940)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1940)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(1915)^{0}} &=& 1.915 \\ \Gamma_{\Sigma(1915)^{0}} &=& 0.12 \\ \mathcal{H}^\mathrm{production}_{\Sigma(1915)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1915)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{\Sigma(2030)^{0}} &=& 2.03 \\ \Gamma_{\Sigma(2030)^{0}} &=& 0.18 \\ \mathcal{H}^\mathrm{production}_{\Sigma(2030)^{0}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(2030)^{0}, 0, - \frac{1}{2}} &=& 1 \\ m_{N(1650)^{+}} &=& 1.65 \\ \Gamma_{N(1650)^{+}} &=& 0.125 \\ \mathcal{H}^\mathrm{production}_{N(1650)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1650)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(1675)^{+}} &=& 1.675 \\ \Gamma_{N(1675)^{+}} &=& 0.145 \\ \mathcal{H}^\mathrm{production}_{N(1675)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1675)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(1680)^{+}} &=& 1.685 \\ \Gamma_{N(1680)^{+}} &=& 0.12 \\ \mathcal{H}^\mathrm{production}_{N(1680)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1680)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(1700)^{+}} &=& 1.72 \\ \Gamma_{N(1700)^{+}} &=& 0.2 \\ \mathcal{H}^\mathrm{production}_{N(1700)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1700)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(1710)^{+}} &=& 1.71 \\ \Gamma_{N(1710)^{+}} &=& 0.14 \\ \mathcal{H}^\mathrm{production}_{N(1710)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1710)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(1720)^{+}} &=& 1.72 \\ \Gamma_{N(1720)^{+}} &=& 0.25 \\ \mathcal{H}^\mathrm{production}_{N(1720)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(1720)^{+}, - \frac{1}{2}, 0} &=& 1 \\ m_{N(2190)^{+}} &=& 2.18 \\ \Gamma_{N(2190)^{+}} &=& 0.4 \\ \mathcal{H}^\mathrm{production}_{N(2190)^{+}, - \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{N(2190)^{+}, - \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(1430)^{+}, 0, 0} &=& 1 \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1950)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(1950)^{+}, 0, 0} &=& 1 \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(700)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{K_{0}^{*}(700)^{+}, 0, 0} &=& 1 \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1980)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{3}^{*}(1780)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{4}^{*}(2045)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1410)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1680)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{+}, 0, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1385)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1660)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1670)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1750)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1775)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1940)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(1915)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{\Sigma(2030)^{0}, 0, \frac{1}{2}} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1650)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1675)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1680)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1700)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1710)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(1720)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{decay}_{N(2190)^{+}, \frac{1}{2}, 0} &=& 1 \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1950)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(700)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1980)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{3}^{*}(1780)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{4}^{*}(2045)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1410)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1680)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{+}, 0, - \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1385)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1660)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1670)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1750)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1775)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1940)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(1915)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{\Sigma(2030)^{0}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1650)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1675)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1680)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1700)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1710)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(1720)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{N(2190)^{+}, \frac{1}{2}, 0} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1980)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{3}^{*}(1780)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K_{4}^{*}(2045)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1410)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(1680)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{+}, 1, \frac{1}{2}} &=& 1+0i \\ m_{0} &=& 4.101931740046389 \\ m_{1} &=& 1.115683 \\ m_{2} &=& 0.49367700000000003 \\ m_{3} &=& 0.1349768 \\ \end{array}\end{split}\]

Kinematic variable definitions

Hide code cell source

Math(aslatex(model.variables))

\[\begin{split}\displaystyle \begin{array}{rcl} \theta_{23} &=& \operatorname{acos}{\left(\frac{2 \sigma_{1} \left(- m_{1}^{2} - m_{2}^{2} + \sigma_{3}\right) - \left(m_{0}^{2} - m_{1}^{2} - \sigma_{1}\right) \left(m_{2}^{2} - m_{3}^{2} + \sigma_{1}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{1}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{1}, m_{2}^{2}, m_{3}^{2}\right)}} \right)} \\ \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)} \\ \theta_{12} &=& \operatorname{acos}{\left(\frac{2 \sigma_{3} \left(- m_{1}^{2} - m_{3}^{2} + \sigma_{2}\right) - \left(m_{0}^{2} - m_{3}^{2} - \sigma_{3}\right) \left(m_{1}^{2} - m_{2}^{2} + \sigma_{3}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{3}^{2}, \sigma_{3}\right)} \sqrt{\lambda\left(\sigma_{3}, m_{1}^{2}, m_{2}^{2}\right)}} \right)} \\ \zeta^0_{1(1)} &=& 0 \\ \zeta^1_{1(1)} &=& 0 \\ \zeta^0_{2(1)} &=& - \operatorname{acos}{\left(\frac{- 2 m_{0}^{2} \left(- m_{1}^{2} - m_{2}^{2} + \sigma_{3}\right) + \left(m_{0}^{2} + m_{1}^{2} - \sigma_{1}\right) \left(m_{0}^{2} + m_{2}^{2} - \sigma_{2}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{2}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(m_{0}^{2}, \sigma_{1}, m_{1}^{2}\right)}} \right)} \\ \zeta^1_{2(1)} &=& \operatorname{acos}{\left(\frac{2 m_{1}^{2} \left(- m_{0}^{2} - m_{3}^{2} + \sigma_{3}\right) + \left(m_{0}^{2} + m_{1}^{2} - \sigma_{1}\right) \left(- m_{1}^{2} - m_{3}^{2} + \sigma_{2}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{1}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{1}^{2}, m_{3}^{2}\right)}} \right)} \\ \zeta^0_{3(1)} &=& \operatorname{acos}{\left(\frac{- 2 m_{0}^{2} \left(- m_{1}^{2} - m_{3}^{2} + \sigma_{2}\right) + \left(m_{0}^{2} + m_{1}^{2} - \sigma_{1}\right) \left(m_{0}^{2} + m_{3}^{2} - \sigma_{3}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{1}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(m_{0}^{2}, \sigma_{3}, m_{3}^{2}\right)}} \right)} \\ \zeta^1_{3(1)} &=& - \operatorname{acos}{\left(\frac{2 m_{1}^{2} \left(- m_{0}^{2} - m_{2}^{2} + \sigma_{2}\right) + \left(m_{0}^{2} + m_{1}^{2} - \sigma_{1}\right) \left(- m_{1}^{2} - m_{2}^{2} + \sigma_{3}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{1}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{3}, m_{1}^{2}, m_{2}^{2}\right)}} \right)} \\ \end{array}\end{split}\]

Visualization

unfolded_expression = model.full_expression.doit()
intensity_func = create_parametrized_function(
    expression=unfolded_expression,
    parameters=model.parameter_defaults,
    backend="jax",
)

i, j = 3, 1
(k,) = {1, 2, 3} - {i, j}
σk, σk_expr = list(model.invariants.items())[k - 1]
Math(aslatex({σk: σk_expr}))

\[\begin{split}\displaystyle \begin{array}{rcl} \sigma_{2} &=& m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \sigma_{1} - \sigma_{3} \\ \end{array}\end{split}\]

Define meshgrid

Hide code cell source

resolution = 1_000
m = sorted(model.masses, key=str)
x_min = float(((m[j] + m[k]) ** 2).xreplace(model.masses))
x_max = float(((m[0] - m[i]) ** 2).xreplace(model.masses))
y_min = float(((m[i] + m[k]) ** 2).xreplace(model.masses))
y_max = float(((m[0] - m[j]) ** 2).xreplace(model.masses))
x_diff = x_max - x_min
y_diff = y_max - y_min
x_min -= 0.05 * x_diff
x_max += 0.05 * x_diff
y_min -= 0.05 * y_diff
y_max += 0.05 * y_diff
X, Y = jnp.meshgrid(
    jnp.linspace(x_min, x_max, num=resolution),
    jnp.linspace(y_min, y_max, num=resolution),
)

Define data transformer

Hide code cell source

definitions = dict(model.variables)
definitions[σk] = σk_expr
definitions = {
    symbol: expr.xreplace(definitions).xreplace(model.masses)
    for symbol, expr in definitions.items()
}

data_transformer = SympyDataTransformer.from_sympy(definitions, backend="jax")
dalitz_data = {
    f"sigma{i}": X,
    f"sigma{j}": Y,
}
dalitz_data.update(data_transformer(dalitz_data))

Collect resonance names per sub-system

Hide code cell source

resonances = defaultdict(set)
for transition in reaction.transitions:
    topology = transition.topology
    top_decay_products = topology.get_edge_ids_outgoing_from_node(0)
    (resonance_id, resonance), *_ = transition.intermediate_states.items()
    recoil_id, *_ = top_decay_products - {resonance_id}
    resonances[recoil_id].add(resonance.particle)
resonances = {k: sorted(v, key=lambda p: p.mass) for k, v in resonances.items()}
{k: [p.name for p in v] for k, v in resonances.items()}

Design slider UI

Hide code cell source

sliders = {}
categorized_sliders_m = defaultdict(list)
categorized_sliders_gamma = defaultdict(list)
categorized_cphi_pair = defaultdict(list)
couplings_name_root = R"\mathcal{H}^\mathrm{decay}"

for symbol, value in model.parameter_defaults.items():
    if symbol.name.startswith(R"\Gamma_{"):
        slider = w.FloatSlider(
            description=Rf"\({sp.latex(symbol)}\)",
            min=0.0,
            max=1.0,
            step=0.01,
            value=value,
            continuous_update=False,
        )
        sliders[symbol.name] = slider
        if symbol.name.startswith(R"\Gamma_{K"):
            categorized_sliders_gamma[1].append(slider)
        elif symbol.name.startswith(R"\Gamma_{\S"):
            categorized_sliders_gamma[2].append(slider)
        elif symbol.name.startswith(R"\Gamma_{N"):
            categorized_sliders_gamma[3].append(slider)

    if symbol.name.startswith("m_{"):
        slider = w.FloatSlider(
            description=Rf"\({sp.latex(symbol)}\)",
            min=0.63,
            max=4,
            step=0.01,
            value=value,
            continuous_update=False,
        )
        sliders[symbol.name] = slider
        if symbol.name.startswith("m_{K"):
            categorized_sliders_m[1].append(slider)
        elif symbol.name.startswith(R"m_{\S"):
            categorized_sliders_m[2].append(slider)
        elif symbol.name.startswith("m_{N"):
            categorized_sliders_m[3].append(slider)

    if symbol.name.startswith(couplings_name_root):
        c_latex = sp.latex(symbol)
        phi_latex = c_latex.replace(couplings_name_root, R"\phi", 1)

        slider_c = w.FloatSlider(
            description=Rf"\({c_latex}\)",
            min=0,
            max=10,
            step=0.01,
            value=abs(value),
            continuous_update=False,
        )
        slider_phi = w.FloatSlider(
            description=Rf"\({phi_latex}\)",
            min=-np.pi,
            max=+np.pi,
            step=0.01,
            value=np.angle(value),
            continuous_update=False,
        )

        sliders[symbol.name] = slider_c
        sliders[symbol.name.replace(couplings_name_root, "phi", 1)] = slider_phi

        cphi_hbox = w.HBox([slider_c, slider_phi])
        if "Sigma" in symbol.name:
            categorized_cphi_pair[2].append(cphi_hbox)
        elif "N" in symbol.name:
            categorized_cphi_pair[3].append(cphi_hbox)
        else:
            categorized_cphi_pair[1].append(cphi_hbox)


assert len(categorized_sliders_gamma) == 3
assert len(categorized_sliders_m) == 3
assert len(categorized_cphi_pair) == 3

subtabs = {}
for recoild_id, resonance_list in resonances.items():
    subtabs[recoild_id] = []
    for particle in resonance_list:
        m_sliders = [
            slider
            for slider in categorized_sliders_m[recoild_id]
            if particle.latex in slider.description
        ]
        gamma_sliders = [
            slider
            for slider in categorized_sliders_gamma[recoild_id]
            if particle.latex in slider.description
        ]
        cphi_pairs = [
            hbox
            for hbox in categorized_cphi_pair[recoild_id]
            if particle.latex in hbox.children[0].description
        ]
        pole_pair = w.HBox(m_sliders + gamma_sliders)
        resonance_tab = w.VBox([pole_pair, *cphi_pairs])
        subtabs[recoild_id].append(resonance_tab)
assert len(subtabs) == 3

main_tabs = []
for recoild_id, slider_boxes in subtabs.items():
    sub_tab_widget = w.Tab(children=slider_boxes)
    for i, particle in enumerate(resonances[recoild_id]):
        sub_tab_widget.set_title(i, particle.name)

    main_tabs.append(sub_tab_widget)

UI = w.Tab(children=main_tabs, titles=("K*", "Σ*", "N*"))

Define functions for inserting phase

Hide code cell source

def insert_phi(parameters: dict) -> dict:
    updated_parameters = {}
    for key, value in parameters.items():
        if key.startswith("phi"):
            continue
        if key.startswith(couplings_name_root):
            phi_key = key.replace(couplings_name_root, "phi")
            if phi_key in parameters:
                phi = parameters[phi_key]
                value *= np.exp(1j * phi)  # noqa:PLW2901
        updated_parameters[key] = value
    return updated_parameters

Show code cell source

Hide code cell source

%matplotlib widget
%config InlineBackend.figure_formats = ['png']
fig_2d, ax_2d = plt.subplots(dpi=200)
ax_2d.set_title("Model-weighted Phase space Dalitz Plot")
ax_2d.set_xlabel(R"$m^2(\Lambda K^+)\;\left[\mathrm{GeV}^2\right]$")
ax_2d.set_ylabel(R"$m^2(K^+ \pi^0)\;\left[\mathrm{GeV}^2\right]$")

mesh = None


def update_dalitz_plot(**parameters):
    global mesh
    parameters = insert_phi(parameters)
    intensity_func.update_parameters(parameters)
    intensities = intensity_func(dalitz_data)  # z
    intensities /= jnp.nansum(intensities)  # normalization

    if mesh is None:
        mesh = ax_2d.pcolormesh(X, Y, intensities, cmap="jet", vmax=3e-5)
    else:
        mesh.set_array(intensities)
    fig_2d.canvas.draw_idle()


interactive_plot = w.interactive_output(update_dalitz_plot, sliders)
fig_2d.tight_layout()
fig_2d.colorbar(mesh, ax=ax_2d)

if STATIC_PAGE:
    filename = "dalitz-plot-dpd.png"
    fig_2d.savefig(filename)
    plt.close(fig_2d)
    display(UI, Image(filename))
else:
    display(UI, interactive_plot)

Show code cell source

Hide code cell source

%matplotlib widget
%config InlineBackend.figure_formats = ['svg']

fig, axes = plt.subplots(figsize=(11, 3.5), ncols=2, sharey=True)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False
ax1, ax2 = axes

for ax in axes:
    recoil_id = 3 if ax is ax1 else 1
    decay_products = sorted(set(reaction.final_state) - {recoil_id})
    product_latex = " ".join([reaction.final_state[i].latex for i in decay_products])
    ax.set_xlabel(f"$m({product_latex})$ [GeV]")

LINES = defaultdict(lambda: None)
RESONANCE_LINE = defaultdict(lambda: None)


def update_plot(**parameters):
    parameters = insert_phi(parameters)
    intensity_func.update_parameters(parameters)
    intensities = intensity_func(dalitz_data)  # z
    intensities /= jnp.nansum(intensities)  # normalization

    max_value = 0
    color_id = 0
    for ax in axes:
        if ax is ax1:
            x = jnp.sqrt(X[0])
            y = jnp.nansum(intensities, axis=0)
        else:
            x = jnp.sqrt(Y[:, 0])
            y = jnp.nansum(intensities, axis=1)

        max_value = max(max_value, y.max())
        recoil_id = 3 if ax is ax1 else 1
        if LINES[recoil_id] is None:
            LINES[recoil_id] = ax.plot(x, y, alpha=0.5)[0]
        else:
            LINES[recoil_id].set_ydata(y)

        for resonance in resonances[recoil_id]:
            key = f"m_{{{resonance.latex}}}"
            val = parameters.get(key, resonance.mass)
            if RESONANCE_LINE[color_id] is None:
                RESONANCE_LINE[color_id] = ax.axvline(
                    val,
                    c=f"C{color_id}",
                    ls="dotted",
                    label=resonance.name,
                )
            else:
                RESONANCE_LINE[color_id].set_xdata([val, val])
            color_id += 1

    for ax in axes:
        ax.set_ylim(0, max_value * 1.1)


interactive_plot = w.interactive_output(update_plot, sliders)
for ax in axes:
    ax.legend(fontsize="small")

if STATIC_PAGE:
    filename = "histogram-dpd.svg"
    fig.savefig(filename)
    plt.close(fig)
    display(UI, SVG(filename))
else:
    display(UI, interactive_plot)