Analytic continuation#
Show 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
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}}"
),
)
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)
Show 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()