7.1. Dynamics lineshapes

Import Python libraries

Hide code cell content

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from ampform.dynamics.form_factor import (
    BlattWeisskopfSquared,
    FormFactor,
    SphericalHankel1,
)
from ampform.dynamics.phasespace import BreakupMomentum
from ampform_dpd import DynamicsBuilder, create_mass_symbol_mapping
from ampform_dpd.decay import ThreeBodyDecay
from ampform_dpd.dynamics import (
    BreitWignerMinL,
    BuggBreitWigner,
    EnergyDependentWidth,
    FlattéSWave,
)
from ampform_dpd.dynamics.builder import get_mandelstam_s
from ampform_dpd.io import cached
from matplotlib_inline.backend_inline import set_matplotlib_formats

from polarimetry.io import (
    display_doit,
    display_latex,
    mute_ampform_warnings,
    mute_jax_warnings,
)
from polarimetry.lhcb import load_model_parameters, load_three_body_decay
from polarimetry.lhcb.dynamics import (
    formulate_breit_wigner,
    formulate_bugg_breit_wigner,
    formulate_exponential_bugg_breit_wigner,
    formulate_flatte_1405,
)
from polarimetry.lhcb.particle import load_particles
from polarimetry.lhcb.symbol import create_meson_radius_symbol
from polarimetry.plot import use_mpl_latex_fonts


def load_parameters(
    decay: ThreeBodyDecay, model_id: int | str = 0
) -> dict[sp.Symbol, complex | float]:
    return load_model_parameters(
        decay=decay,
        filename="../../data/model-definitions.yaml",
        model_id=model_id,
        particle_definitions=particles,
    )


particles = load_particles("../../data/particle-definitions.yaml")
Sigma = particles["Sigma-"]
K = particles["K-"]
Lambda_c = particles["Lambda_c+"]
p = particles["p"]
pi = particles["pi+"]
mute_jax_warnings()
mute_ampform_warnings()
set_matplotlib_formats("svg")
use_mpl_latex_fonts()

ls_model_id = "Alternative amplitude model obtained using LS couplings"
exp_model_id = "Alternative amplitude model with an additional overall exponential form factor exp(-alpha q^2) multiplying Bugg lineshapes. The exponential parameter is indicated as alpha"

In the following, we consider a decay \(0 \xrightarrow{L} (r \xrightarrow{\ell} ij)k\), where the resonance \(r\) is produced at the production vertex with angular momentum \(L\) and decays at the decay vertex with angular momentum \(\ell\). The meson radii for both vertices are denoted \(R_\mathrm{prod}\) and \(R_\mathrm{dec}\), respectively.

7.1.1. Relativistic Breit–Wigner

Most resonances were parametrized with a standard relativistic Breit–Wigner with two form factors for the vertices.

Show code cell source

Hide code cell source

s, mr, Γr, mi, mj = sp.symbols("s m_r Gamma_R m_i m_j", nonnegative=True)
m0, mk = sp.symbols("m0 m_k", nonnegative=True)
R_prod = create_meson_radius_symbol("prod")
R_dec = create_meson_radius_symbol("dec")
L, ell = sp.symbols("L ell", integer=True, nonnegative=True)
display_doit(BreitWignerMinL(s, m0, mk, mr, Γr, mi, mj, ell, L, R_dec, R_prod))

\[\begin{split}\displaystyle \begin{aligned} \mathcal{R}^\mathrm{BW}_{\ell,L}\left(s\right) \;&=\; \frac{\frac{\mathcal{F}_{L}\left(m_{0}^{2}, \sqrt{s}, m_{k}\right)}{\mathcal{F}_{L}\left(m_{0}^{2}, m_{r}, m_{k}\right)} \frac{\mathcal{F}_{\ell}\left(s, m_{i}, m_{j}\right)}{\mathcal{F}_{\ell}\left(m_{r}^{2}, m_{i}, m_{j}\right)}}{m_{r}^{2} - i m_{r} \Gamma\left(s\right) - s} \\ \end{aligned}\end{split}\]

7.1.2. Bugg Breit–Wigner

The \(K(700)\) and \(K(1430)\) were modeled with a Bugg Breit–Wigner function. Here, \(s_A\) is the Adler zero.

Show code cell source

Hide code cell source

gamma, sA = sp.symbols("gamma s_A")
bugg = BuggBreitWigner(s, mr, Γr, mi, mj, gamma)
sA_expr = mi**2 - mj**2 / 2
display_latex({
    bugg: bugg.evaluate().subs(sA_expr, sA),
    sA: sA_expr,
})

\[\begin{split}\displaystyle \begin{aligned} \mathcal{R}^\mathrm{Bugg}\left(s\right) \;&=\; \frac{1}{- \frac{i \Gamma_{R} m_{r} \left(s - s_{A}\right) e^{- \gamma s}}{m_{r}^{2} - s_{A}} + m_{r}^{2} - s} \\ s_{A} \;&=\; m_{i}^{2} - \frac{m_{j}^{2}}{2} \\ \end{aligned}\end{split}\]

One of the models uses a Bugg Breit–Wigner with an exponential factor (see this paper, Eq. (4)):

Show code cell source

Hide code cell source

q = BreakupMomentum(m0**2, sp.sqrt(s), mk)
alpha = sp.Symbol("alpha")
bugg * sp.exp(-alpha * q**2)

\[\displaystyle e^{- \alpha q_{0}\left(m_{0}^{2}\right)^{2}} \mathcal{R}^\mathrm{Bugg}\left(s\right)\]

7.1.3. Flatté for S-waves

In the original amplitude analysis, the \(\varLambda(1405)\) baryon was modeled with a Flatté for \(S\)-waves (\(\ell=0\)), coupled to the channel \(\varLambda(1405) \to \varSigma^-\pi^+\). Note that both the \(\varSigma^-\pi^+\) threshold and the mass of the \(\varLambda(1405)\) lies below the \(pK^-\) threshold.

Show code cell source

Hide code cell source

Γ1, Γ2, mπ, mΣ = sp.symbols("Gamma1 Gamma2 m_pi m_Sigma")
display_doit(FlattéSWave(s, mr, (Γ1, Γ2), (mi, mj), (mπ, mΣ)))

\[\begin{split}\displaystyle \begin{aligned} \mathcal{R}^\mathrm{Flatté}\left(s\right) \;&=\; \frac{1}{m_{r}^{2} - i m_{r} \left(\frac{\Gamma_{1} m_{r} q\left(s\right)}{\sqrt{s} q\left(m_{r}^{2}\right)} + \frac{\Gamma_{2} m_{r} q\left(s\right)}{\sqrt{s} q\left(m_{r}^{2}\right)}\right) - s} \\ \end{aligned}\end{split}\]

7.1.4. Other definitions

Show code cell source

Hide code cell source

x = sp.Symbol("x", nonnegative=True)
exprs = [
    EnergyDependentWidth(s, mr, Γr, mi, mj, ell, R_dec),
    FormFactor(m0**2, sp.sqrt(s), mk, L, R_prod),
    FormFactor(s, mi, mj, ell, R_dec),
    BreakupMomentum(m0**2, sp.sqrt(s), mk),
    BreakupMomentum(s, mi, mj),
    BlattWeisskopfSquared(x**2, ell),
    SphericalHankel1(ell, x),
]
display_latex({x: x.doit(deep=False) for x in exprs})

\[\begin{split}\displaystyle \begin{aligned} \Gamma\left(s\right) \;&=\; \frac{\Gamma_{R} \mathcal{F}_{\ell}\left(s, m_{i}, m_{j}\right)^{2} \rho\left(s\right)}{\mathcal{F}_{\ell}\left(m_{r}^{2}, m_{i}, m_{j}\right)^{2} \rho\left(m_{r}^{2}\right)} \\ \mathcal{F}_{L}\left(m_{0}^{2}, \sqrt{s}, m_{k}\right) \;&=\; \sqrt{B_{L}^2\left(R_\mathrm{prod}^{2} q^2_{0}\left(m_{0}^{2}\right)\right)} \\ \mathcal{F}_{\ell}\left(s, m_{i}, m_{j}\right) \;&=\; \sqrt{B_{\ell}^2\left(R_\mathrm{dec}^{2} q^2\left(s\right)\right)} \\ q_{0}\left(m_{0}^{2}\right) \;&=\; \frac{\sqrt{\left(m_{0}^{2} - \left(- m_{k} + \sqrt{s}\right)^{2}\right) \left(m_{0}^{2} - \left(m_{k} + \sqrt{s}\right)^{2}\right)}}{2 m_{0}} \\ q\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{i} - m_{j}\right)^{2}\right) \left(s - \left(m_{i} + m_{j}\right)^{2}\right)}}{2 \sqrt{s}} \\ B_{\ell}^2\left(x^{2}\right) \;&=\; \frac{\left|{h_{\ell}^{(1)}\left(1\right)}\right|^{2}}{x^{2} \left|{h_{\ell}^{(1)}\left(x\right)}\right|^{2}} \\ h_{\ell}^{(1)}\left(x\right) \;&=\; \frac{\left(- i\right)^{\ell + 1} e^{i x} \sum_{k=0}^{\ell} \frac{\left(\frac{i}{2 x}\right)^{k} \left(k + \ell\right)!}{k! \left(- k + \ell\right)!}}{x} \\ \end{aligned}\end{split}\]

7.1.5. Lineshape visualizations

The lineshapes are visualized below for some of the relevant resonances. The section of lineshapes for each resonance can be found model-definitions.yaml

Function for plotting a dynamics builder

Hide code cell source

def plot_lineshape(
    dynamics_builder: DynamicsBuilder,
    title: str,
    resonances: str,
    x_range: tuple[float, float],
    y_title: float,
    legend_position: tuple[int, tuple[float, float]],
    model_id: int | str = 0,
) -> None:
    min_ls = len(resonances) != 1
    if not min_ls:
        resonance = particles[resonances[0]]
        title = f"{title} for ${resonance.latex}$"
    decay = load_three_body_decay(resonances, particles, min_ls)
    parameter_defaults = load_parameters(decay, model_id)
    parameter_defaults.update(create_mass_symbol_mapping(decay))
    display_latex({
        chain: dynamics_builder(chain).expression.doit(deep=False)
        for chain in decay.chains
    })
    fig, axes = plt.subplots(
        figsize=(6, 2 * len(decay.chains)),
        nrows=len(decay.chains),
        sharex=True,
    )
    fig.patch.set_facecolor("none")
    for ax in axes:
        ax.patch.set_facecolor("none")
        ax.spines["bottom"].set_position("zero")
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)
    fig.suptitle(title, y=y_title)
    x = np.linspace(*x_range, num=300)
    for ax, chain in zip(axes, decay.chains, strict=True):
        dynamics = dynamics_builder(chain)
        func = cached.lambdify(
            dynamics.expression.doit(),
            parameters=dynamics.parameters | parameter_defaults,
            backend="jax",
        )
        mass_name = f"m_{{{chain.resonance.latex}}}"
        resonance_mass = func.parameters[mass_name]
        child1, child2 = chain.decay_products
        phsp_min = child1.mass + child2.mass
        phsp_max = chain.parent.mass - chain.spectator.mass
        ax.axvspan(x.min(), phsp_min, alpha=0.15, color="gray", label="non-physical")
        ax.axvspan(phsp_max, x.max(), alpha=0.15, color="gray")
        ax.axvline(resonance_mass, c="red", ls="dotted", label=R"$m_\mathrm{res}$")
        sigma = get_mandelstam_s(chain.decay_node)
        z = func({sigma.name: x**2})
        ax.plot(x, np.abs(z), label="abs")
        ax.plot(x, z.imag, label="imag")
        ax.plot(x, z.real, label="real")
        if min_ls:
            axis_title = f"${chain.resonance.latex}$"
        else:
            axis_title = Rf"$L={chain.incoming_ls.L}, \ell={chain.outgoing_ls.L}$"
        ax.set_title(axis_title, y=0.86)
        ax.set_xlim(*x_range)
    legend_axis, anchor = legend_position
    axes[legend_axis].legend(
        bbox_to_anchor=anchor,
        framealpha=1,
        handlelength=1,
        handletextpad=0.5,
        loc="upper right",
    )
    fig.tight_layout()
    plt.show()

7.1.5.1. Breit–Wigner with form factors

Show code cell source

Hide code cell source

plot_lineshape(
    formulate_breit_wigner,
    title="Relativistic Breit–Wigner with form factors",
    resonances=["K(892)"],
    x_range=(0.5, 1.5),
    y_title=0.97,
    model_id=ls_model_id,
    legend_position=(0, (1.03, 1.03)),
)

\[\begin{split}\displaystyle \begin{aligned} \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(K(892) \xrightarrow[S=0]{L=1} \pi^+ K^-\right) p \;&=\; \frac{\frac{\mathcal{F}_{0}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right)}{\mathcal{F}_{0}\left(m_{0}^{2}, m_{K(892)}, m_{1}\right)} \frac{\mathcal{F}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right)}{\mathcal{F}_{1}\left(m_{K(892)}^{2}, m_{2}, m_{3}\right)}}{m_{K(892)}^{2} - i m_{K(892)} \Gamma_{892}\left(\sigma_{1}\right) - \sigma_{1}} \\ \Lambda_c^+ \xrightarrow[S=1/2]{L=1} \left(K(892) \xrightarrow[S=0]{L=1} \pi^+ K^-\right) p \;&=\; \frac{\frac{\mathcal{F}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right)}{\mathcal{F}_{1}\left(m_{0}^{2}, m_{K(892)}, m_{1}\right)} \frac{\mathcal{F}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right)}{\mathcal{F}_{1}\left(m_{K(892)}^{2}, m_{2}, m_{3}\right)}}{m_{K(892)}^{2} - i m_{K(892)} \Gamma_{892}\left(\sigma_{1}\right) - \sigma_{1}} \\ \Lambda_c^+ \xrightarrow[S=3/2]{L=1} \left(K(892) \xrightarrow[S=0]{L=1} \pi^+ K^-\right) p \;&=\; \frac{\frac{\mathcal{F}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right)}{\mathcal{F}_{1}\left(m_{0}^{2}, m_{K(892)}, m_{1}\right)} \frac{\mathcal{F}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right)}{\mathcal{F}_{1}\left(m_{K(892)}^{2}, m_{2}, m_{3}\right)}}{m_{K(892)}^{2} - i m_{K(892)} \Gamma_{892}\left(\sigma_{1}\right) - \sigma_{1}} \\ \Lambda_c^+ \xrightarrow[S=3/2]{L=2} \left(K(892) \xrightarrow[S=0]{L=1} \pi^+ K^-\right) p \;&=\; \frac{\frac{\mathcal{F}_{2}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right)}{\mathcal{F}_{2}\left(m_{0}^{2}, m_{K(892)}, m_{1}\right)} \frac{\mathcal{F}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right)}{\mathcal{F}_{1}\left(m_{K(892)}^{2}, m_{2}, m_{3}\right)}}{m_{K(892)}^{2} - i m_{K(892)} \Gamma_{892}\left(\sigma_{1}\right) - \sigma_{1}} \\ \end{aligned}\end{split}\]

7.1.5.2. Bugg Breit–Wigner

Show code cell source

Hide code cell source

plot_lineshape(
    formulate_bugg_breit_wigner,
    title="Bugg Breit–Wigner",
    resonances=["K(700)", "K(1430)"],
    x_range=(0.5, 1.5),
    y_title=0.92,
    legend_position=(0, (1.03, 1.15)),
)

\[\begin{split}\displaystyle \begin{aligned} \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(K(700) \xrightarrow[S=0]{L=0} \pi^+ K^-\right) p \;&=\; \frac{1}{- \frac{i \Gamma_{K(700)} m_{K(700)} \left(\frac{m_{2}^{2}}{2} - m_{3}^{2} + \sigma_{1}\right) e^{- \gamma_{K(700)} \sigma_{1}}}{\frac{m_{2}^{2}}{2} - m_{3}^{2} + m_{K(700)}^{2}} + m_{K(700)}^{2} - \sigma_{1}} \\ \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(K(1430) \xrightarrow[S=0]{L=0} \pi^+ K^-\right) p \;&=\; \frac{1}{- \frac{i \Gamma_{K(1430)} m_{K(1430)} \left(\frac{m_{2}^{2}}{2} - m_{3}^{2} + \sigma_{1}\right) e^{- \gamma_{K(1430)} \sigma_{1}}}{\frac{m_{2}^{2}}{2} - m_{3}^{2} + m_{K(1430)}^{2}} + m_{K(1430)}^{2} - \sigma_{1}} \\ \end{aligned}\end{split}\]

Show code cell source

Hide code cell source

plot_lineshape(
    formulate_exponential_bugg_breit_wigner,
    title="Exponential Bugg Breit–Wigner",
    resonances=["K(700)", "K(1430)"],
    x_range=(0.5, 1.5),
    y_title=0.92,
    legend_position=(0, (1.03, 1.15)),
    model_id=exp_model_id,
)

\[\begin{split}\displaystyle \begin{aligned} \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(K(700) \xrightarrow[S=0]{L=0} \pi^+ K^-\right) p \;&=\; e^{- \alpha_{K(700)} q_{0}\left(m_{0}^{2}\right)^{2}} \mathcal{R}^\mathrm{Bugg}\left(\sigma_{1}\right) \\ \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(K(1430) \xrightarrow[S=0]{L=0} \pi^+ K^-\right) p \;&=\; e^{- \alpha_{K(1430)} q_{0}\left(m_{0}^{2}\right)^{2}} \mathcal{R}^\mathrm{Bugg}\left(\sigma_{1}\right) \\ \end{aligned}\end{split}\]

7.1.5.3. Flatté for \(\varLambda(1405)\)

Show code cell source

Hide code cell source

plot_lineshape(
    formulate_flatte_1405,
    title="S-wave Flatté",
    resonances=["L(1405)"],
    x_range=(1.3, 2.3),
    y_title=0.94,
    legend_position=(1, (1.08, 1.02)),
    model_id=ls_model_id,
)

\[\begin{split}\displaystyle \begin{aligned} \Lambda_c^+ \xrightarrow[S=1/2]{L=0} \left(\Lambda(1405) \xrightarrow[S=1/2]{L=0} K^- p\right) \pi^+ \;&=\; \frac{\mathcal{R}^\mathrm{Flatté}\left(\sigma_{2}\right) \mathcal{F}_{0}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right)}{\mathcal{F}_{0}\left(m_{0}^{2}, m_{\Lambda(1405)}, m_{2}\right)} \\ \Lambda_c^+ \xrightarrow[S=1/2]{L=1} \left(\Lambda(1405) \xrightarrow[S=1/2]{L=0} K^- p\right) \pi^+ \;&=\; \frac{\mathcal{R}^\mathrm{Flatté}\left(\sigma_{2}\right) \mathcal{F}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right)}{\mathcal{F}_{1}\left(m_{0}^{2}, m_{\Lambda(1405)}, m_{2}\right)} \\ \end{aligned}\end{split}\]