Analytic continuation

import logging
import warnings

import ampform
import graphviz
import matplotlib.pyplot as plt
import qrules
from IPython.display import Math

from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

logger = logging.getLogger()
logger.setLevel(logging.ERROR)
warnings.filterwarnings("ignore")

Sometimes qrules finds resonances that lie outside phase space, because resonances can ‘leak’ into phase space. The example below is simplified and the two selected resonances are a bit unusual, but it serves to illustrate how to handle these sub-threshold resonances.

reaction = qrules.generate_transitions(
    initial_state="D0",
    final_state=["K-", "K+", "K0"],
    allowed_intermediate_particles=["a(0)(980)0", "a(2)(1320)0"],
    formalism="canonical-helicity",
)
dot = qrules.io.asdot(
    reaction,
    collapse_graphs=True,
    render_node=False,
)
graphviz.Source(dot)

Of the two resonances, \(a_0(980)\) lies just below threshold ― it’s mass is smaller than the the masses of the two decay products combined:

pdg = qrules.load_pdg()
a_meson = pdg["a(0)(980)0"]
k_plus = pdg["K+"]
k_minus = pdg["K+"]
two_k_mass = round(k_minus.mass + k_plus.mass, 3)
display(
    Math(
        Rf"m_{{{a_meson.latex}}} = {a_meson.mass} \pm"
        Rf" {a_meson.width}\;\mathrm{{GeV}}"
    ),
    Math(
        Rf"m_{{{k_plus.latex}}} + m_{{{k_minus.latex}}} ="
        Rf" {two_k_mass}\;\mathrm{{GeV}}"
    ),
)

\[\displaystyle m_{a_{0}(980)^{0}} = 0.98 \pm 0.075\;\mathrm{GeV}\]

\[\displaystyle m_{K^{+}} + m_{K^{+}} = 0.987\;\mathrm{GeV}\]

To correctly describe the dynamics for this resonance, we should use make use of analytic continuation. As opposed to Step 1: Formulate model, we now construct a RelativisticBreitWignerBuilder where we set its phase space factor to PhaseSpaceFactorSWave:

from ampform.dynamics import PhaseSpaceFactorSWave
from ampform.dynamics.builder import (
    RelativisticBreitWignerBuilder,
    create_non_dynamic_with_ff,
    create_relativistic_breit_wigner_with_ff,
)

model_builder = ampform.get_builder(reaction)
analytic_breit_wigner_builder = RelativisticBreitWignerBuilder(
    form_factor=True,
    energy_dependent_width=True,
    phsp_factor=PhaseSpaceFactorSWave,
)
model_builder.set_dynamics(
    "J/psi(1S)",
    create_non_dynamic_with_ff,
)
model_builder.set_dynamics(
    "a(0)(980)0",
    analytic_breit_wigner_builder,
)
model_builder.set_dynamics(
    "a(2)(1320)0",
    create_relativistic_breit_wigner_with_ff,
)
model = model_builder.formulate()

Note that we have reduced the coupling for the sub-threshold resonance \(a_0(980)^0\), so that it doesn’t dominate over $$

The effect can be seen once we generate data. Despite the fact that the resonance lies outside phase space, there is still a contribution to the intensity:

for par in model.parameter_defaults:
    if not par.name.startswith("C") or "a_{0}(980)^{0}" not in par.name:
        continue
    model.parameter_defaults[par] = 0.1

intensity = create_parametrized_function(
    expression=model.expression.doit(),
    parameters=model.parameter_defaults,
    backend="jax",
)
helicity_transformer = SympyDataTransformer.from_sympy(
    model.kinematic_variables, backend="numpy"
)
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
phsp_momenta = phsp_generator.generate(100_000, rng)

import numpy as np
from matplotlib import cm

phsp = helicity_transformer(phsp_momenta)
intensities = np.array(intensity(phsp))

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]

fig, ax = plt.subplots()
ax.set_xlabel("$m_{02}$ [GeV]")
for p, color in zip(resonances, colors):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.set_yticks([])
ax.hist(phsp["m_02"], bins=100, alpha=0.5, weights=intensities)
ax.legend()
plt.show()