Analytic continuation¶
import logging
import warnings
import ampform
import graphviz
import matplotlib.pyplot as plt
import pandas as pd
import qrules
from IPython.display import Math
from tensorwaves.data import generate_data, generate_phsp
from tensorwaves.data.transform import HelicityTransformer
from tensorwaves.model import LambdifiedFunction, SympyModel
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)+"],
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(
fR"m_{{{a_meson.latex}}} = {a_meson.mass} \pm"
fR" {a_meson.width}\;\mathrm{{GeV}}"
),
Math(
fR"m_{{{k_plus.latex}}} + m_{{{k_minus.latex}}} ="
fR" {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: Create amplitude model, we now use create_analytic_breit_wigner() (a relativistic Breit-Wigner with analytic continuation) instead of create_relativistic_breit_wigner_with_ff():
from ampform.dynamics.builder import (
create_analytic_breit_wigner,
create_non_dynamic_with_ff,
create_relativistic_breit_wigner_with_ff,
)
model_builder = ampform.get_builder(reaction)
model_builder.set_dynamics(
"J/psi(1S)",
create_non_dynamic_with_ff,
)
model_builder.set_dynamics(
"a(0)(980)0",
create_analytic_breit_wigner,
)
model_builder.set_dynamics(
"a(2)(1320)+",
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:
sympy_model = SympyModel(
expression=model.expression.doit(),
parameters=model.parameter_defaults,
)
intensity = LambdifiedFunction(sympy_model, backend="jax")
data_converter = HelicityTransformer(model.adapter)
reaction_info = model.adapter.reaction_info
initial_state_mass = reaction_info.initial_state[-1].mass
final_state_masses = {i: p.mass for i, p in reaction_info.final_state.items()}
phsp_sample = generate_phsp(100_000, initial_state_mass, final_state_masses)
data_sample = generate_data(
2_000,
initial_state_mass,
final_state_masses,
data_converter,
intensity,
)
import numpy as np
from matplotlib import cm
reaction_info = model.adapter.reaction_info
intermediate_states = sorted(
(
p
for p in model.particles
if p not in reaction_info.final_state.values()
and p not in reaction_info.initial_state.values()
),
key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(intermediate_states))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]
def indicate_masses():
plt.xlabel("$m_{12}$ [GeV]")
for i, p in enumerate(intermediate_states):
plt.gca().axvline(
x=p.mass, linestyle="dotted", label=p.name, color=colors[i]
)
phsp_set = data_converter.transform(phsp_sample)
data_set = data_converter.transform(data_sample)
data_frame = pd.DataFrame(data_set)
data_frame["m_12"].hist(bins=50, alpha=0.5, density=True)
indicate_masses()
plt.legend();
