\(\Lambda_c^+ \to p K^- \pi^+\)

Model definition: lc2ppik-lhcb-2683025.json.

This notebooks illustrates the use of the ampform_dpd.io.serialization module for the decay \(\Lambda_c^+ \to p K^- \pi^+\). The corresponding model was optimized to a data sample of roughly half a million \(\Lambda_c^{\pm}\) decay candidates by the LHCb collaboration, INSPIRE-HEP 2683025.

Warning

The ampform_dpd.io.serialization module is a preview feature. This notebook illustrates the deserialization of the amplitude model JSON file to symbolic expressions. Keep an eye on ComPWA/ampform-dpd#133 for a list of tracked issues.

Import model

from __future__ import annotations
from pathlib import Path
import json
import logging
import os

import jax.numpy as jnp
import matplotlib.pyplot as plt
import pandas as pd
import sympy as sp
from ampform.dynamics import BlattWeisskopfSquared
from ampform.dynamics.form_factor import FormFactor
from ampform.dynamics.phasespace import BreakupMomentumSquared
from ampform.kinematics.phasespace import Kallen
from ampform.sympy import perform_cached_doit
from ampform_dpd import DefinedExpression
from ampform_dpd.decay import FinalStateID, State, ThreeBodyDecay
from ampform_dpd.dynamics import (
    BreitWigner,
    ChannelArguments,
    EnergyDependentWidth,
    MultichannelBreitWigner,
    P,
    SimpleBreitWigner,
)
from ampform_dpd.io import aslatex, perform_cached_lambdify, simplify_latex_rendering
from ampform_dpd.io.serialization.amplitude import (
    HelicityRecoupling,
    LSRecoupling,
    ParityRecoupling,
    formulate,
    formulate_aligned_amplitude,
    formulate_chain_amplitude,
    formulate_recoupling,
)
from ampform_dpd.io.serialization.decay import get_final_state, to_decay
from ampform_dpd.io.serialization.dynamics import (
    formulate_breit_wigner,
    formulate_dynamics,
    formulate_form_factor,
    formulate_multichannel_breit_wigner,
    to_mandelstam_symbol,
    to_mass_symbol,
)
from ampform_dpd.io.serialization.format import (
    ModelDefinition,
    Propagator,
    get_decay_chains,
    get_function_definition,
)
from IPython.display import JSON, Math
from tqdm.auto import tqdm

THIS_DIR = Path(".").absolute()
logging.getLogger("ampform.sympy").setLevel(logging.ERROR)
simplify_latex_rendering()

with open(THIS_DIR.parent.parent / "models" / "lc2ppik-lhcb-2683025.json") as stream:
    MODEL_DEFINITION = json.load(stream)

Construct ThreeBodyDecay

def to_latex(name: str) -> str:
    latex = {
        "Lc": R"\Lambda_c^+",
        "pi": R"\pi^+",
        "K": "K^-",
        "p": "p",
    }.get(name)
    if latex is not None:
        return latex
    mass_str = name[1:].strip("(").strip(")")
    subsystem_letter = name[0]
    subsystem = {"D": "D", "K": "K", "L": R"\Lambda"}.get(subsystem_letter)
    if subsystem is None:
        return name
    return f"{subsystem}({mass_str})"

DECAY = to_decay(MODEL_DEFINITION, to_latex=to_latex)
Math(aslatex(DECAY, with_jp=True))

\(\displaystyle \begin{array}{c} \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1232)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1600)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(D(1700)\left[J=\frac{3}{2}\right] \to p\left[J=\frac{1}{2}\right] \pi^+\left[J=0\right]\right) K^-\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(1430)\left[J=0\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(700)\left[J=0\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(K(892)\left[J=1\right] \to \pi^+\left[J=0\right] K^-\left[J=0\right]\right) p\left[J=\frac{1}{2}\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1405)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1520)\left[J=\frac{3}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1600)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1670)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(1690)\left[J=\frac{3}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \Lambda_c^+\left[J=\frac{1}{2}\right] \to \left(\Lambda(2000)\left[J=\frac{1}{2}\right] \to K^-\left[J=0\right] p\left[J=\frac{1}{2}\right]\right) \pi^+\left[J=0\right] \\ \end{array}\)

Dynamics

See also RUB-EP1/amplitude-serialization#22 about serialization of custom lineshapes.

CHAIN_DEFS = get_decay_chains(MODEL_DEFINITION)

Vertices

Blatt-Weisskopf form factor

z = sp.Symbol("z", nonnegative=True)
s, m1, m2, L, d = sp.symbols("s m1 m2 L R", nonnegative=True)
exprs = [
    FormFactor(s, m1, m2, L, d),
    BlattWeisskopfSquared(z, L),
    BreakupMomentumSquared(s, m1, m2),
]
Math(aslatex({e: e.doit(deep=False) for e in exprs}))

\(\displaystyle \begin{array}{rcl} \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right) &=& \sqrt{B_{L}^2\left(R^{2} q^2\left(s\right)\right)} \\ B_{L}^2\left(z\right) &=& \frac{\left|{h_{L}^{(1)}\left(1\right)}\right|^{2}}{z \left|{h_{L}^{(1)}\left(\sqrt{z}\right)}\right|^{2}} \\ q^2\left(s\right) &=& \frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{4 s} \\ \end{array}\)

ff_L1520 = formulate_form_factor(
    vertex=CHAIN_DEFS[2]["vertices"][0],
    model=MODEL_DEFINITION,
)
Math(aslatex(ff_L1520))

\(\displaystyle \begin{array}{rcl} \mathcal{F}_{1}\left(m_{0}, \sigma_{2}, m_{2}\right) &=& \mathcal{F}_{1}\left(m_{0}, \sigma_{2}, m_{2}\right) \\ R_{Lc} &=& 5.0 \\ \end{array}\)

Propagators

Breit-Wigner

x, y, z = sp.symbols("x:z")
s, m0, Γ0, m1, m2, L, d = sp.symbols("s m0 Gamma0 m1 m2 L R", nonnegative=True)
exprs = [
    BreitWigner(s, m0, Γ0, m1, m2, L, d),
    SimpleBreitWigner(s, m0, Γ0),
    EnergyDependentWidth(s, m0, Γ0, m1, m2, L, d),
    FormFactor(s, m1, m2, L, d),
    P(s, m1, m2),
    Kallen(x, y, z),
]
Math(aslatex({e: e.doit(deep=False) for e in exprs}))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_{L}\left(s; m_{0}, \Gamma_{0}\right) &=& \mathcal{R}^\mathrm{BW}\left(s; m_{0}, \Gamma_{L}\left(s\right)\right) \\ \mathcal{R}^\mathrm{BW}\left(s; m_{0}, \Gamma_{0}\right) &=& \frac{1}{- 1.0 i \Gamma_{0} m_{0} + m_{0}^{2} - s} \\ \Gamma_{L}\left(s\right) &=& \Gamma_{0} \frac{m_{0}}{\sqrt{s}} \frac{F_{L}\left(R p_{_{m_{1},m_{2}}}\left(s\right)\right)^{2}}{F_{L}\left(R p_{_{m_{1},m_{2}}}\left(m_{0}^{2}\right)\right)^{2}} \left(\frac{p_{_{m_{1},m_{2}}}\left(s\right)}{p_{_{m_{1},m_{2}}}\left(m_{0}^{2}\right)}\right)^{2 L + 1} \\ \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right) &=& \sqrt{B_{L}^2\left(R^{2} q^2\left(s\right)\right)} \\ p_{_{m_{1},m_{2}}}\left(s\right) &=& \frac{\sqrt{\lambda\left(s, m_{1}^{2}, m_{2}^{2}\right)}}{2 \sqrt{s}} \\ \lambda\left(x, y, z\right) &=& x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \\ \end{array}\)

K892_BW = formulate_breit_wigner(
    propagator=CHAIN_DEFS[20]["propagators"][0],
    resonance=to_latex(CHAIN_DEFS[20]["name"]),
    model=MODEL_DEFINITION,
)
Math(aslatex(K892_BW))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{K(892)}, \Gamma_{K(892)}\right) &=& \mathcal{R}^\mathrm{BW}\left(\sigma_{1}; m_{K(892)}, \Gamma_{1}\left(\sigma_{1}\right)\right) \\ m_{K(892)} &=& 0.8955 \\ \Gamma_{K(892)} &=& 0.047299999999999995 \\ m_{2} &=& 0.13957018 \\ m_{3} &=& 0.493677 \\ R_\mathrm{res} &=& 1.5 \\ \end{array}\)

Multi-channel Breit-Wigner

x, y, z = sp.symbols("x:z")
s, m0, Γ0, m1, m2, L, d = sp.symbols("s m0 Gamma0 m1 m2 L R", nonnegative=True)
channels = tuple(
    ChannelArguments(
        s,
        m0,
        width=sp.Symbol(f"Gamma{i}", nonnegative=True),
        m1=sp.Symbol(f"m_{{a,{i}}}", nonnegative=True),
        m2=sp.Symbol(f"m_{{b,{i}}}", nonnegative=True),
        angular_momentum=sp.Symbol(f"L{i}", integer=True, nonnegative=True),
        meson_radius=d,
    )
    for i in [1, 2]
)
exprs = [
    MultichannelBreitWigner(s, m0, channels),
    BreitWigner(s, m0, Γ0, m1, m2, L, d),
    BreitWigner(s, m0, Γ0),
    EnergyDependentWidth(s, m0, Γ0, m1, m2, L, d),
    FormFactor(s, m1, m2, L, d),
    P(s, m1, m2),
    Kallen(x, y, z),
]
Math(aslatex({e: e.doit(deep=False) for e in exprs}))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(s; \Gamma_{1}, \Gamma_{2}\right) &=& \mathcal{R}^\mathrm{BW}_{L=0}\left(s; m_{0}, \frac{\Gamma_{1} m_{0} \mathcal{F}_{L_{1}}\left(s, m_{a,1}, m_{b,1}\right)^{2}}{\sqrt{s}} + \frac{\Gamma_{2} m_{0} \mathcal{F}_{L_{2}}\left(s, m_{a,2}, m_{b,2}\right)^{2}}{\sqrt{s}}\right) \\ \mathcal{R}^\mathrm{BW}_{L}\left(s; m_{0}, \Gamma_{0}\right) &=& \mathcal{R}^\mathrm{BW}\left(s; m_{0}, \Gamma_{L}\left(s\right)\right) \\ \mathcal{R}^\mathrm{BW}_{L=0}\left(s; m_{0}, \Gamma_{0}\right) &=& \frac{1}{- 1.0 i \Gamma_{0} m_{0} + m_{0}^{2} - s} \\ \Gamma_{L}\left(s\right) &=& \Gamma_{0} \frac{m_{0}}{\sqrt{s}} \frac{F_{L}\left(R p_{_{m_{1},m_{2}}}\left(s\right)\right)^{2}}{F_{L}\left(R p_{_{m_{1},m_{2}}}\left(m_{0}^{2}\right)\right)^{2}} \left(\frac{p_{_{m_{1},m_{2}}}\left(s\right)}{p_{_{m_{1},m_{2}}}\left(m_{0}^{2}\right)}\right)^{2 L + 1} \\ \mathcal{F}_{L}\left(s, m_{1}, m_{2}\right) &=& \sqrt{B_{L}^2\left(R^{2} q^2\left(s\right)\right)} \\ p_{_{m_{1},m_{2}}}\left(s\right) &=& \frac{\sqrt{\lambda\left(s, m_{1}^{2}, m_{2}^{2}\right)}}{2 \sqrt{s}} \\ \lambda\left(x, y, z\right) &=& x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \\ \end{array}\)

L1405_Flatte = formulate_multichannel_breit_wigner(
    propagator=CHAIN_DEFS[0]["propagators"][0],
    resonance=to_latex(CHAIN_DEFS[0]["name"]),
    model=MODEL_DEFINITION,
)
Math(aslatex(L1405_Flatte))

\(\displaystyle \begin{array}{rcl} \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(\sigma_{2}; \Gamma_{\Lambda(1405)}, \Gamma_{\Lambda(1405)}^\text{ch. 2}\right) &=& \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{\Lambda(1405)}, \frac{\Gamma_{\Lambda(1405)} m_{\Lambda(1405)} \mathcal{F}_{0}\left(\sigma_{2}, m_{3}, m_{1}\right)^{2}}{\sqrt{\sigma_{2}}} + \frac{\Gamma_{\Lambda(1405)}^\text{ch. 2} m_{\Lambda(1405)} \mathcal{F}_{0}\left(\sigma_{2}, m_{a,2}, m_{b,2}\right)^{2}}{\sqrt{\sigma_{2}}}\right) \\ m_{\Lambda(1405)} &=& 1.4051 \\ \Gamma_{\Lambda(1405)} &=& 0.23395150538434703 \\ m_{3} &=& 0.938272046 \\ m_{1} &=& 0.493677 \\ R_{\Lambda(1405)} &=& 0 \\ m_{a,2} &=& 1.18937 \\ m_{b,2} &=& 0.13957018 \\ \Gamma_{\Lambda(1405)}^\text{ch. 2} &=& 0.23395150538434703 \\ \end{array}\)

Construct AmplitudeModel

Unpolarized intensity

λ0, λ1, λ2, λ3 = sp.symbols("lambda(:4)", rational=True)
amplitude_expr, _ = formulate_aligned_amplitude(MODEL_DEFINITION, λ0, λ1, λ2, λ3)
amplitude_expr.cleanup()

\(\displaystyle \sum_{\lambda_0^{\prime}=-1/2}^{1/2} \sum_{\lambda_1^{\prime}=-1/2}^{1/2}{A^{1}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{1(2)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{1(2)}\right) + A^{2}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{2(2)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{2(2)}\right) + A^{3}_{\lambda_0^{\prime}, \lambda_1^{\prime}, 0, 0} d^{\frac{1}{2}}_{\lambda_1^{\prime},\lambda_{1}}\left(\zeta^1_{3(2)}\right) d^{\frac{1}{2}}_{\lambda_{0},\lambda_0^{\prime}}\left(\zeta^0_{3(2)}\right)}\)

Amplitude for the decay chain

Helicity recouplings

λa = sp.Symbol(R"\lambda_a", rational=True)
λb = sp.Symbol(R"\lambda_b", rational=True)
λa0 = sp.Symbol(R"\lambda_a^0", rational=True)
λb0 = sp.Symbol(R"\lambda_b^0", rational=True)
f = sp.Symbol("f", integer=True)
l = sp.Symbol("l", integer=True, nonnegative=True)
s = sp.Symbol("s", nonnegative=True, rational=True)
ja = sp.Symbol("j_a", nonnegative=True, rational=True)
jb = sp.Symbol("j_b", nonnegative=True, rational=True)
j = sp.Symbol("j", nonnegative=True, rational=True)
exprs = [
    HelicityRecoupling(λa, λb, λa0, λb0),
    ParityRecoupling(λa, λb, λa0, λb0, f),
    LSRecoupling(λa, λb, l, s, ja, jb, j),
]
Math(aslatex({e: e.doit(deep=False) for e in exprs}))

\(\displaystyle \begin{array}{rcl} \mathcal{H}^\text{helicity}\left(\lambda_{a},\lambda_{b}|\lambda^{0}_{a},\lambda^{0}_{b}\right) &=& \delta_{\lambda_{a} \lambda^{0}_{a}} \delta_{\lambda_{b} \lambda^{0}_{b}} \\ \mathcal{H}^\text{parity}\left(\lambda_{a},\lambda_{b}|\lambda^{0}_{a},\lambda^{0}_{b},f\right) &=& f \delta_{\lambda_{a}, - \lambda^{0}_{a}} \delta_{\lambda_{b}, - \lambda^{0}_{b}} + \delta_{\lambda_{a} \lambda^{0}_{a}} \delta_{\lambda_{b} \lambda^{0}_{b}} \\ \mathcal{H}^\text{parity}\left(\lambda_{a},\lambda_{b}|l,s,j_{a},j_{b},j\right) &=& \frac{\sqrt{2 l + 1} C^{s,\lambda_{a} - \lambda_{b}}_{j_{a},\lambda_{a},j_{b},- \lambda_{b}} C^{j,\lambda_{a} - \lambda_{b}}_{l,0,s,\lambda_{a} - \lambda_{b}}}{\sqrt{2 j + 1}} \\ \end{array}\)

Recoupling deserialization

recouplings = [
    formulate_recoupling(MODEL_DEFINITION, chain_idx=0, vertex_idx=i) for i in range(2)
]
Math(aslatex({e: e.doit(deep=False) for e in recouplings}))

\(\displaystyle \begin{array}{rcl} \mathcal{H}^\text{helicity}\left(\lambda_{R},\lambda_{2}|\frac{1}{2},0\right) &=& \delta_{0 \lambda_{2}} \delta_{\frac{1}{2} \lambda_{R}} \\ \mathcal{H}^\text{parity}\left(\lambda_{3},\lambda_{1}|0,\frac{1}{2},1\right) &=& \delta_{- \frac{1}{2} \lambda_{1}} \delta_{0 \lambda_{3}} + \delta_{0 \lambda_{3}} \delta_{\frac{1}{2} \lambda_{1}} \\ \end{array}\)

Chain amplitudes

definitions = formulate_chain_amplitude(λ0, λ1, λ2, λ3, MODEL_DEFINITION, chain_idx=0)
Math(aslatex(definitions))

\(\displaystyle \begin{array}{rcl} A^{2}_{\lambda_{0}, \lambda_{1}, \lambda_{2}, \lambda_{3}} &=& c^{L1405[1/2]}_{\frac{1}{2}, 0, 0} \mathcal{H}^\text{helicity}\left(\frac{1}{2},\lambda_{2}|\frac{1}{2},0\right) \mathcal{R}^\mathrm{BW}_\mathrm{multi}\left(\sigma_{2}; \Gamma_{L1405}, \Gamma_{L1405}^\text{ch. 2}\right) \mathcal{H}^\text{parity}\left(\lambda_{3},\lambda_{1}|0,\frac{1}{2},1\right) d^{\frac{1}{2}}_{\frac{1}{2},- \frac{1}{2}}\left(\theta_{31}\right) \\ c^{L1405[1/2]}_{\frac{1}{2}, 0, 0} &=& 7.38649400481717+1.971018433257411i \\ m_{L1405} &=& 1.4051 \\ \Gamma_{L1405} &=& 0.23395150538434703 \\ m_{3} &=& 0.938272046 \\ m_{1} &=& 0.493677 \\ R_{L1405} &=& 0 \\ m_{a,2} &=& 1.18937 \\ m_{b,2} &=& 0.13957018 \\ \Gamma_{L1405}^\text{ch. 2} &=& 0.23395150538434703 \\ \theta_{31} &=& \operatorname{acos}{\left(\frac{2 \sigma_{2} \left(- m_{2}^{2} - m_{3}^{2} + \sigma_{1}\right) - \left(m_{0}^{2} - m_{2}^{2} - \sigma_{2}\right) \left(- m_{1}^{2} + m_{3}^{2} + \sigma_{2}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{2}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{3}^{2}, m_{1}^{2}\right)}} \right)} \\ \end{array}\)

Validation

Warning

The following serves as a numerical check on whether the amplitude model has been deserialized correctly. For now, this is not the case, see ComPWA/ampform-dpd#133 for updates.

checksums = {
    misc_key: {checksum["point"]: checksum["value"] for checksum in misc_value}
    for misc_key, misc_value in MODEL_DEFINITION["misc"].items()
    if "checksum" in misc_key
}
checksums

{'amplitude_model_checksums': {'validation_point': 9345.853380852355,
  'validation_point_m31sq': '-2.383989898030206 + 0.34006358781864954i',
  'validation_point_m12sq': '-0.32420527030886465 + 0.9152170090261661i',
  'validation_point_m23sq': '-1.6160420072107298 + 0.30076446839110044i'}}

checksum_points = {
    point["name"]: {par["name"]: par["value"] for par in point["parameters"]}
    for point in MODEL_DEFINITION["parameter_points"]
}
checksum_points

{'validation_point': {'cos_theta_31': -0.2309352648098208,
  'phi_31': 0.0,
  'm_31': 1.9101377207489973,
  'cos_theta_31_2': 0.0,
  'phi_31_2': 0.0,
  'm_31_2': 2.28646},
 'validation_point_m12sq': {'m_12_sq': 3.2},
 'validation_point_m23sq': {'m_23_sq': 1.4},
 'validation_point_m31sq': {'m_31_sq': 3.2}}