Analytic continuation

Show code cell content

Hide code cell content

import logging
import warnings

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

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

logging.getLogger("tensorwaves.data").setLevel(logging.ERROR)  # hide progress bars
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:

Show code cell source

Hide code cell source

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,  # ty:ignore[invalid-argument-type]
)
model_builder.dynamics.assign("D0", create_non_dynamic_with_ff)
model_builder.dynamics.assign("a(0)(980)0", analytic_breit_wigner_builder)
model_builder.dynamics.assign("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 \(a_2(1320)^0\) resonance.

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 symbol in model.parameter_defaults:
    if not str(symbol).startswith("C") or "a_{0}(980)^{0}" not in str(symbol):
        continue
    model.parameter_defaults[symbol] = 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)

Show code cell source

Hide code cell source

%config InlineBackend.figure_formats = ['svg']

import numpy as np

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 = [plt.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, strict=True):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.set_yticks([])
ax.hist(phsp["m_01"], bins=100, alpha=0.5, weights=intensities)
ax.legend()
plt.show()