Ξb⁻ → p K⁻ K⁻

Ξb⁻ → p K⁻ K⁻#

Hide code cell content 
from __future__ import annotations

import logging
import os
import warnings
from typing import TYPE_CHECKING

import graphviz
import jax.numpy as jnp
import matplotlib.pyplot as plt
import qrules
import sympy as sp
from ampform.dynamics import EnergyDependentWidth
from ampform.dynamics.form_factor import BlattWeisskopfSquared, FormFactor
from ampform.sympy import perform_cached_doit
from IPython.display import Latex, Markdown
from tensorwaves.data.transform import SympyDataTransformer
from tqdm.auto import tqdm

from ampform_dpd import DalitzPlotDecompositionBuilder
from ampform_dpd.adapter.qrules import (
    load_particles,
    normalize_state_ids,
    permute_equal_final_states,
    to_three_body_decay,
)
from ampform_dpd.decay import State
from ampform_dpd.dynamics import RelativisticBreitWigner
from ampform_dpd.dynamics.builder import formulate_breit_wigner_with_form_factor
from ampform_dpd.io import (
    as_markdown_table,
    aslatex,
    perform_cached_lambdify,
    simplify_latex_rendering,
)

simplify_latex_rendering()
logging.getLogger("absl").setLevel(logging.ERROR)  # mute JAX
warnings.simplefilter("ignore", category=RuntimeWarning)

NO_TQDM = "EXECUTE_NB" in os.environ
if NO_TQDM:
    logging.getLogger("ampform_dpd.io").setLevel(logging.ERROR)
if TYPE_CHECKING:
    from tensorwaves.interface import DataSample, ParametrizedFunction

Decay definition#

See DOI:10.1103/PhysRevD.104.052010 [pdf], Search for CP violation in \(\Xi_b^- \to p K^- K^-\) decays by LHCb. It found six asymmetry parameters, for \(\Lambda(1405)\), \(\Lambda(1520)\), \(\Lambda(1670)\), \(\Sigma(1385)\), \(\Sigma(1775)\), and \(\Sigma(1915)\).

Hide code cell source 
PARTICLES = load_particles()
excluded_particles = [
    "Lambda(1600)",
    "Lambda(1690)",
    "Lambda(1800)",
    "Lambda(1810)",
    "Lambda(1890)",
    "Lambda(2000)",
    "Sigma(1660)0",
    "Sigma(1670)0",
    "Sigma(1750)0",
    "Sigma(1910)0",
    "Sigma(c)(2455)0",
    "Sigma(c)(2520)0",
]
for name in excluded_particles:
    PARTICLES.remove(PARTICLES[name])

REACTION = qrules.generate_transitions(
    initial_state="Xi(b)-",
    final_state=["K-", "K-", "p"],
    formalism="canonical-helicity",
    allowed_intermediate_particles=["Lambda", "Sigma"],
    max_angular_momentum=2,
    particle_db=PARTICLES,
)
REACTION = normalize_state_ids(REACTION)
REACTION = permute_equal_final_states(REACTION)
dot = qrules.io.asdot(REACTION, collapse_graphs=True)
graphviz.Source(dot)
Hide code cell source 
DECAY = to_three_body_decay(REACTION.transitions, min_ls=True)
Markdown(as_markdown_table([DECAY.initial_state, *DECAY.final_state.values()]))
Hide code cell source 
resonances = sorted(
    {t.resonance for t in DECAY.chains},
    key=lambda p: (p.name[0], p.mass),
)
resonance_names = [p.name for p in resonances]
Markdown(as_markdown_table(resonances))
Hide code cell source 
def get_decay_identifier(decay):
    return (decay.resonance, *(particle.name for particle in decay.decay_products))


chains = {get_decay_identifier(c): c for c in DECAY.chains}
Latex(aslatex(chains.values(), with_jp=True))

Model formulation#

model_builder = DalitzPlotDecompositionBuilder(DECAY, min_ls=False)
for chain in model_builder.decay.chains:
    model_builder.dynamics_choices.register_builder(
        chain, formulate_breit_wigner_with_form_factor
    )
model = model_builder.formulate(reference_subsystem=1)
model.intensity
Hide code cell source 
Latex(aslatex(model.variables))

Each unaligned amplitude is defined as follows:

Hide code cell source 
Latex(aslatex(model.amplitudes))
Hide code cell source 
s, m0, w0, m1, m2, L, R, z = sp.symbols("s m0 Gamma0 m1 m2 L R z", nonnegative=True)
exprs = [
    RelativisticBreitWigner(s, m0, w0, m1, m2, L, R),
    EnergyDependentWidth(s, m0, w0, m1, m2, L, R),
    FormFactor(s, m1, m2, L, R),
    BlattWeisskopfSquared(z, L),
]
Latex(aslatex({e: e.doit(deep=False) for e in exprs}))

Preparing for input data#

Hide code cell source 
i, j = (1, 2)
k, *_ = {1, 2, 3} - {i, j}
σk, σk_expr = list(model.invariants.items())[k - 1]
Latex(aslatex({σk: σk_expr}))
Hide code cell source 
resolution = 1_000
m = sorted(model.masses, key=str)
x_min = float(((m[j] + m[k]) ** 2).xreplace(model.masses))
x_max = float(((m[0] - m[i]) ** 2).xreplace(model.masses))
y_min = float(((m[i] + m[k]) ** 2).xreplace(model.masses))
y_max = float(((m[0] - m[j]) ** 2).xreplace(model.masses))
x_diff = x_max - x_min
y_diff = y_max - y_min
x_min -= 0.05 * x_diff
x_max += 0.05 * x_diff
y_min -= 0.05 * y_diff
y_max += 0.05 * y_diff
X, Y = jnp.meshgrid(
    jnp.linspace(x_min, x_max, num=resolution),
    jnp.linspace(y_min, y_max, num=resolution),
)
Hide code cell source 
definitions = dict(model.variables)
definitions[σk] = σk_expr
definitions = {
    symbol: expr.xreplace(definitions).xreplace(model.masses)
    for symbol, expr in definitions.items()
}
data_transformer = SympyDataTransformer.from_sympy(definitions, backend="jax")
dalitz_data = {
    f"sigma{i}": X,
    f"sigma{j}": Y,
}
dalitz_data.update(data_transformer(dalitz_data))
Hide code cell source 
unfolded_intensity_expr = perform_cached_doit(model.intensity)
full_intensity_expr = perform_cached_doit(
    unfolded_intensity_expr.xreplace(model.amplitudes)
)
free_parameters = {
    k: v
    for k, v in model.parameter_defaults.items()
    if isinstance(k, sp.Indexed)
    if "production" in str(k) or "decay" in str(k)
}
fixed_parameters = {
    k: v for k, v in model.parameter_defaults.items() if k not in free_parameters
}
intensity_func = perform_cached_lambdify(
    full_intensity_expr.xreplace(fixed_parameters),
    parameters=free_parameters,
    backend="jax",
)
intensities = intensity_func(dalitz_data)

Dalitz plot#

Hide code cell source 
def get_decay_products(subsystem_id: int) -> tuple[State, State]:
    return tuple(s for s in DECAY.final_state.values() if s.index != subsystem_id)


plt.rc("font", size=18)
I_tot = jnp.nansum(intensities)
normalized_intensities = intensities / I_tot

fig, ax = plt.subplots(figsize=(14, 10))
mesh = ax.pcolormesh(X, Y, normalized_intensities)
ax.set_aspect("equal")
c_bar = plt.colorbar(mesh, ax=ax, pad=0.01)
c_bar.ax.set_ylabel("Normalized intensity (a.u.)")
sigma_labels = {
    i: Rf"$\sigma_{i} = M^2\left({' '.join(p.latex for p in get_decay_products(i))}\right)$"
    for i in (1, 2, 3)
}
ax.set_xlabel(sigma_labels[i])
ax.set_ylabel(sigma_labels[j])
plt.show()
Hide code cell source 
def compute_sub_intensity(
    func: ParametrizedFunction, phsp: DataSample, resonance_latex: str
) -> jnp.ndarray:
    original_parameters = dict(func.parameters)
    zero_parameters = {
        k: 0
        for k, v in func.parameters.items()
        if R"\mathcal{H}" in k
        if resonance_latex not in k
    }
    func.update_parameters(zero_parameters)
    intensities = func(phsp)
    func.update_parameters(original_parameters)
    return intensities


plt.rc("font", size=16)
fig, ax = plt.subplots(figsize=(10, 6), sharey=True)
fig.subplots_adjust(wspace=0.02)
x = jnp.sqrt(X[0])
y = jnp.sqrt(Y[:, 0])
ax.fill_between(x, jnp.nansum(normalized_intensities, axis=0), alpha=0.5)
_, y_max = ax.get_ylim()
ax.set_ylim(0, y_max)
ax.autoscale(enable=False, axis="x")
ax.set_ylabel("Normalized intensity (a.u.)")
ax.set_xlabel(sigma_labels[i])
resonance_counter = 0
for chain in tqdm(model.decay.chains, disable=NO_TQDM):
    if {p.index for p in chain.decay_products} != {1, 3}:
        continue
    resonance = chain.resonance
    sub_intensities = compute_sub_intensity(
        intensity_func, dalitz_data, resonance.latex
    )
    color = f"C{resonance_counter}"
    ax.plot(x, jnp.nansum(sub_intensities / I_tot, axis=0), c=color)
    ax.axvline(resonance.mass, label=f"${resonance.latex}$", c=color, ls="dashed")
    resonance_counter += 1
ax.legend(fontsize=12)
plt.show()