Inspect model interactively

In this notebook, we illustrate how to interactively inspect a HelicityModel using Sympy Plotting Backends and the ampform.sympy.slider module. The procedure should work for any sympy.Expr.

First, we create some HelicityModel. We could also have used pickle to load() the HelicityModel that we created in Formulate amplitude model, but the cell below allows running this notebook independently.

import qrules

from ampform import get_builder
from ampform.dynamics.builder import (
    create_non_dynamic_with_ff,
    create_relativistic_breit_wigner_with_ff,
)

reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["gamma", "pi0", "pi0"],
    allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"],
    allowed_interaction_types=["strong", "EM"],
    formalism="canonical-helicity",
)
builder = get_builder(reaction)
builder.config.stable_final_state_ids = {0, 1, 2}
builder.config.scalar_initial_state_mass = True
initial_state_particle = reaction.initial_state[-1]
builder.dynamics.assign(initial_state_particle, create_non_dynamic_with_ff)
for name in reaction.get_intermediate_particles().names:
    builder.dynamics.assign(name, create_relativistic_breit_wigner_with_ff)
model = builder.formulate()

In this case, as we saw, the overall model contains just one intensity term \(I = |\sum_i A_i|^2\), with \(\sum_i A_i\) some coherent sum of amplitudes. We can extract \(\sum_i A_i\) as follows:

import sympy as sp

amplitude = model.expression.args[0].args[0].args[0]
assert isinstance(amplitude, sp.Add)

Show code cell source

Hide code cell source

from IPython.display import Math

from ampform.io import aslatex

Math(aslatex(amplitude, terms_per_line=1))

\[\begin{split}\displaystyle \begin{aligned} & \frac{C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \Gamma_{f_{0}(1500)} m_{f_{0}(1500)} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,1}_{0,0,1,1} C^{1,1}_{1,1,0,0} \mathcal{F}_{0}\left(m_{012}^{2}, m_{0}, m_{12}\right) \mathcal{F}_{0}\left(m_{12}^{2}, m_{1}, m_{2}\right) D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{1,1}\left(- \phi_{0},\theta_{0},0\right)}{m_{12}^{2} - m_{f_{0}(1500)}^{2} + i m_{f_{0}(1500)} \Gamma\left(m_{12}^{2}\right)} \\ & \;+\; \frac{C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \Gamma_{f_{0}(980)} m_{f_{0}(980)} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,1}_{0,0,1,1} C^{1,1}_{1,1,0,0} \mathcal{F}_{0}\left(m_{012}^{2}, m_{0}, m_{12}\right) \mathcal{F}_{0}\left(m_{12}^{2}, m_{1}, m_{2}\right) D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{1,1}\left(- \phi_{0},\theta_{0},0\right)}{m_{12}^{2} - m_{f_{0}(980)}^{2} + i m_{f_{0}(980)} \Gamma\left(m_{12}^{2}\right)} \\ & \;+\; \frac{C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \Gamma_{f_{0}(1500)} m_{f_{0}(1500)} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,1}_{1,1,0,0} C^{1,1}_{2,0,1,1} \mathcal{F}_{2}\left(m_{012}^{2}, m_{0}, m_{12}\right) \mathcal{F}_{0}\left(m_{12}^{2}, m_{1}, m_{2}\right) D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{1,1}\left(- \phi_{0},\theta_{0},0\right)}{m_{12}^{2} - m_{f_{0}(1500)}^{2} + i m_{f_{0}(1500)} \Gamma\left(m_{12}^{2}\right)} \\ & \;+\; \frac{C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}} \Gamma_{f_{0}(980)} m_{f_{0}(980)} \left(C^{0,0}_{0,0,0,0}\right)^{2} C^{1,1}_{1,1,0,0} C^{1,1}_{2,0,1,1} \mathcal{F}_{2}\left(m_{012}^{2}, m_{0}, m_{12}\right) \mathcal{F}_{0}\left(m_{12}^{2}, m_{1}, m_{2}\right) D^{0}_{0,0}\left(- \phi^{12}_{1},\theta^{12}_{1},0\right) D^{1}_{1,1}\left(- \phi_{0},\theta_{0},0\right)}{m_{12}^{2} - m_{f_{0}(980)}^{2} + i m_{f_{0}(980)} \Gamma\left(m_{12}^{2}\right)} \\ \end{aligned}\end{split}\]

Substitute some of the boring parameters with the provided parameter_defaults:

Substitute certain parameters

Hide code cell source

amplitude = amplitude.doit().subs({
    symbol: value
    for symbol, value in model.parameter_defaults.items()
    if not str(symbol).startswith((R"\Gamma_", "m_{", "m_1", "m_2"))
})
amplitude = amplitude.subs({
    symbol: 0
    for symbol in amplitude.free_symbols
    if str(symbol).startswith(("phi", "theta"))
})
amplitude.free_symbols

{\Gamma_{f_{0}(1500)},
 \Gamma_{f_{0}(980)},
 m_1,
 m_12,
 m_2,
 m_{f_{0}(1500)},
 m_{f_{0}(980)}}

The ampform.sympy.slider module contains some handy functions for creating ipywidgets sliders from SymPy symbols. Here is an example that will be useful when using the spb module with interactive parameter sliders.

from ampform.sympy.slider import create_slider

sliders = {
    s: create_slider(s, value=v)
    for s, v in model.parameter_defaults.items()
    if s in amplitude.free_symbols
}
for symbol, slider in sliders.items():
    if str(symbol).startswith("m"):
        slider.max = 2.3
    if str(symbol).startswith(R"\Gamma"):
        slider.max = 1
    slider.step = 0.01
display(*sliders.values())

Finally, we can can use the spb module to create a few interactive plots!

%matplotlib widget

Create widget with Sympy Plotting Backends

Hide code cell source

import matplotlib.pyplot as plt
import spb

t = sp.Symbol("t")
m1, m2, m12 = sp.symbols("m_1 m_2 m_12", nonnegative=True)

line_kwargs = dict(
    expr_y=t,
    params=sliders,
    n=2,
    range_p=(t, 0, 10),
    rendering_kw=dict(ls="dotted"),
    use_cm=False,
    name="threshold",
)
xlim = (0.2, 2.5)
plt.rc("font", size=14)
fig, axes = plt.subplots(figsize=(12, 4), ncols=2, gridspec_kw={"width_ratios": [1, 2]})
ax1, ax2 = axes.ravel()
fig.subplots_adjust(bottom=0.15, left=0.1, right=0.98, top=0.95)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False
spb.graphics(
    spb.line_parametric_2d(
        expr_x=sp.re(amplitude),
        expr_y=sp.im(amplitude),
        params=sliders,
        range_p=(m12, 0.01, 3),
        use_cm=False,
    ),
    aspect="equal",
    ax=ax1,
    ncols=3,
    xlabel=R"$\text{Re}\,A$",
    ylabel=R"$\text{Im}\,A$",
)
exprs = [m1 + m2, *(s for s in sliders if str(s).startswith("m_{"))]
UI = spb.graphics(
    *(
        spb.line_parametric_2d(x, **line_kwargs, label=f"${sp.latex(x)}$")
        for x in exprs
    ),
    spb.line_abs_arg_colored(amplitude, (m12, *xlim), params=sliders, label="A"),
    ax=ax2,
    ncols=3,
    xlabel=R"$m_{\pi^0\pi^0}$",
    ylabel="$|A|$",
    xlim=xlim,
    ylim=(0, 2),
)
UI

Export widget as GIF

Hide code cell source

if STATIC_WEB_PAGE:
    from IPython.display import Image, display
    from matplotlib.animation import FuncAnimation, PillowWriter
    from tqdm.auto import tqdm

    def frame_generator():
        while True:
            if not forward and animated_slider.value <= start:
                return
            yield

    def animate(_):
        global forward
        if animated_slider.value >= end:
            forward = False
        if forward:
            animated_slider.value += step * animated_slider.step
        else:
            animated_slider.value -= step * animated_slider.step
        fig.canvas.draw_idle()
        pbar.set_postfix_str(f"value={animated_slider.value:g}")
        pbar.update()

    animated_slider, *_ = sliders.values()
    start = max(1.0, animated_slider.min)
    end = min(2.0, animated_slider.max)
    step = 3

    forward = True
    animated_slider.value = start
    pbar = tqdm(desc="Exporting animation", leave=False)
    output_path = "animation.gif"
    animation = FuncAnimation(
        fig,
        animate,
        cache_frame_data=False,
        frames=frame_generator(),
        repeat=False,
    )
    animation.save(output_path, writer=PillowWriter(fps=10))
    pbar.close()
    with open(output_path, "rb") as f:
        display(UI.children[0], Image(data=f.read(), format="png"))

This figure is an animation over just one of the parameters. Run the notebook itself to play around with all parameters!

Tip

See K-matrix for why \(\boldsymbol{K}\)-matrix dynamics are better than simple Breit–Wigners when resonances are close to each other.