Inspect model interactively#

In this notebook, we illustrate how to interactively inspect a HelicityModel. The procedure should in fact work for any sympy.Expr.

Create amplitude model#

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",
)
model_builder = get_builder(reaction)
initial_state_particle = reaction.initial_state[-1]
model_builder.dynamics.assign(initial_state_particle, create_non_dynamic_with_ff)
for name in reaction.get_intermediate_particles().names:
    model_builder.dynamics.assign(name, create_relativistic_breit_wigner_with_ff)
model = 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)

Substitute some of the boring parameters with the provided parameter_defaults:

substitutions = {
    symbol: value
    for symbol, value in model.parameter_defaults.items()
    if not symbol.name.startswith(R"\Gamma_") and not symbol.name.startswith("m_")
}
amplitude = amplitude.doit().subs(substitutions)

Lambdify#

We now need to identify the Symbol over which the amplitude is to be plotted. The remaining symbols will be turned into slider parameters.

plot_symbol = sp.Symbol("m_12", nonnegative=True)
slider_symbols = sorted(amplitude.free_symbols, key=lambda s: s.name)
slider_symbols.remove(plot_symbol)
slider_symbols
[\Gamma_{f_{0}(1500)},
 \Gamma_{f_{0}(980)},
 m_0,
 m_012,
 m_1,
 m_2,
 m_{f_{0}(1500)},
 m_{f_{0}(980)},
 phi_0,
 theta_0]

Next, lambdify() the expression:

np_amplitude = sp.lambdify(
    (plot_symbol, *slider_symbols),
    amplitude,
    "numpy",
)

We also have to define some functions that formulate what we want to plot. A pure amplitude won’t do, because we can only plot real values:

import numpy as np


def intensity(plot_variable, **kwargs):
    values = np_amplitude(plot_variable, **kwargs)
    return np.abs(values) ** 2


def argand(**kwargs):
    values = np_amplitude(plot_domain, **kwargs)
    argand = np.array([values.real, values.imag])
    return argand.T

Prepare sliders#

Tip

This procedure has been extracted to the façade function symplot.prepare_sliders().

We also need to define the domain over which to plot, as well sliders for each of the remaining symbols. The function create_slider() helps creating an ipywidgets slider for each Symbol:

from symplot import create_slider

plot_domain = np.linspace(0.2, 2.5, num=400)
sliders_mapping = {symbol.name: create_slider(symbol) for symbol in slider_symbols}

If the name of a Symbol is not a valid name for a Python variable (see str.isidentifier()), lambdify() ‘dummifies’ it, so it can be used as argument for the lambdified function. Since interactive_plot() works with these (dummified) arguments as identifiers for the sliders, we need some mapping between the Symbol and their dummified name. This can be constructed with the help of inspect.signature():

import inspect

symbols_names = [s.name for s in (plot_symbol, *slider_symbols)]
arg_names = list(inspect.signature(np_amplitude).parameters)
arg_to_symbol = dict(zip(arg_names, symbols_names))

The class SliderKwargs comes in as a handy manager for this dict of sliders:

from symplot import SliderKwargs

sliders = SliderKwargs(sliders_mapping, arg_to_symbol)
sliders
SliderKwargs(
  sliders={'\\Gamma_{f_{0}(1500)}': FloatSlider(value=0.0, description='\\(\\Gamma_{f_{0}(1500)}\\)'),
   '\\Gamma_{f_{0}(980)}': FloatSlider(value=0.0, description='\\(\\Gamma_{f_{0}(980)}\\)'),
   'm_0': FloatSlider(value=0.0, description='\\(m_{0}\\)'),
   'm_012': FloatSlider(value=0.0, description='\\(m_{012}\\)'),
   'm_1': FloatSlider(value=0.0, description='\\(m_{1}\\)'),
   'm_2': FloatSlider(value=0.0, description='\\(m_{2}\\)'),
   'm_{f_{0}(1500)}': FloatSlider(value=0.0, description='\\(m_{f_{0}(1500)}\\)'),
   'm_{f_{0}(980)}': FloatSlider(value=0.0, description='\\(m_{f_{0}(980)}\\)'),
   'phi_0': FloatSlider(value=0.0, description='\\(\\phi_{0}\\)'),
   'theta_0': FloatSlider(value=0.0, description='\\(\\theta_{0}\\)')},
  arg_to_symbol={'Dummy_138': '\\Gamma_{f_{0}(1500)}',
   'Dummy_137': '\\Gamma_{f_{0}(980)}',
   'm_0': 'm_0',
   'm_012': 'm_012',
   'm_1': 'm_1',
   'm_2': 'm_2',
   'Dummy_136': 'm_{f_{0}(1500)}',
   'Dummy_135': 'm_{f_{0}(980)}',
   'phi_0': 'phi_0',
   'theta_0': 'theta_0'},
)

SliderKwargs also provides convenient methods for setting ranges and initial values for the sliders.

n_steps = 100
sliders.set_ranges({
    "m_{f_{0}(980)}": (0.3, 1.8, n_steps),
    "m_{f_{0}(1500)}": (0.3, 1.8, n_steps),
    R"\Gamma_{f_{0}(980)}": (0.01, 1, n_steps),
    R"\Gamma_{f_{0}(1500)}": (0.01, 1, n_steps),
    "m_1": (0.01, 1, n_steps),
    "m_2": (0.01, 1, n_steps),
    "phi_0": (0, 2 * np.pi, 40),
    "theta_0": (-np.pi, np.pi, 40),
})

The method set_values() is designed in particular for parameter_defaults, but also works well with both argument names (that may have been dummified) and the original symbol names:

import qrules

pdg = qrules.load_pdg()
sliders.set_values(model.parameter_defaults)
sliders.set_values(
    {  # symbol names
        "phi_0": 0,
        "theta_0": np.pi / 2,
    },
    # argument names
    m_1=pdg["pi0"].mass,
    m_2=pdg["pi0"].mass,
)

Interactive Argand plot#

Finally, we can use mpl-interactions to plot the functions defined above with regard to the parameter values:

Hide code cell source
%config InlineBackend.figure_formats = ['svg']

import matplotlib.pyplot as plt
import mpl_interactions.ipyplot as iplt

# Create figure
fig, axes = plt.subplots(1, 2, figsize=0.9 * np.array((8, 3.8)), tight_layout=True)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False
fig.suptitle(R"$J/\psi \to \gamma f_0, f_0 \to \pi^0\pi^0$")
ax_intensity, ax_argand = axes
m_label = R"$m_{\pi^0\pi^0}$ (GeV)"
ax_intensity.set_xlabel(m_label)
ax_intensity.set_ylabel("$I = |A|^2$")
ax_argand.set_xlabel("Re($A$)")
ax_argand.set_ylabel("Im($A$)")

# Fill plots
controls = iplt.plot(
    plot_domain,
    intensity,
    **sliders,
    slider_formats=dict.fromkeys(arg_names, "{:.3f}"),
    ylim="auto",
    ax=ax_intensity,
)
iplt.scatter(
    argand,
    controls=controls,
    xlim="auto",
    ylim="auto",
    parametric=True,
    c=plot_domain,
    s=1,
    ax=ax_argand,
)
plt.colorbar(label=m_label, ax=ax_argand, aspect=30, pad=0.01)
plt.winter()
../../_images/d28e7668d367c00093b8fd7c4c42a7cdbcef40e829d5b1bc9195be2862f92676.png

Tip

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