Analytic continuation

Analytic continuation#

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)
../../_images/ddfcf6a9039ce28ebb5ad570ebc3e683c34bd30d9f64584f351ca41cc74d3563.svg

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:

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,
)
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()

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)
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):
    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()
../../_images/7a7c2b7bceaa026fe92d026a40f611133a1f6d432e4ef70d8687ecca291d87ae.svg