Spin alignment with data

Spin alignment with data#

Hide code cell content 
%config InlineBackend.figure_formats = ['svg']

import logging
import warnings
from pathlib import Path

from IPython.display import display

LOGGER = logging.getLogger()
LOGGER.setLevel(logging.ERROR)
warnings.filterwarnings("ignore")

Phase space sample#

Hide code cell content 
import qrules

PDG = qrules.load_pdg()

from tensorwaves.data import TFPhaseSpaceGenerator, TFUniformRealNumberGenerator

phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=PDG["Lambda(c)+"].mass,
    final_state_masses={
        0: PDG["p"].mass,
        1: PDG["K-"].mass,
        2: PDG["pi+"].mass,
    },
)
rng = TFUniformRealNumberGenerator(seed=0)
phsp_momenta = phsp_generator.generate(1_000_000, rng)

Generate transitions#

Hide code cell content 
from qrules.particle import ParticleCollection, create_particle

particle_db = ParticleCollection()
particle_db.add(PDG["Lambda(c)+"])
particle_db.add(PDG["p"])
particle_db.add(PDG["K-"])
particle_db.add(PDG["pi+"])

particle_db.add(
    create_particle(
        PDG["K*(892)0"],
        name="K*",
        latex="K^*",
    )
)
particle_db.add(
    create_particle(
        PDG["Lambda(1405)"],
        name="Lambda*",
        latex=R"\Lambda^*",
    )
)
particle_db.add(
    create_particle(
        PDG["Delta(1232)++"],
        name="Delta*++",
        latex=R"\Delta^*",
    )
)

reaction = qrules.generate_transitions(
    initial_state=("Lambda(c)+", [-0.5, +0.5]),
    final_state=["p", "K-", "pi+"],
    formalism="helicity",
    particle_db=particle_db,
)
Hide code cell source 
import graphviz

n = len(reaction.transitions)
for i, t in enumerate(reaction.transitions[:: n // 3]):
    dot = qrules.io.asdot([t], collapse_graphs=True, size=3.5)
    graph = graphviz.Source(dot)
    output_file = Path(f"graph{i}")
    graph.render(output_file, format="svg")
    output_file.unlink()
    display(graph)

013_12_0 013_12_1 013_12_2

Distribution without alignment#

Amplitude model formulated following Appendix C:

import ampform
from ampform.dynamics.builder import RelativisticBreitWignerBuilder

builder = ampform.get_builder(reaction)
builder.align_spin = False
builder.stable_final_state_ids = list(reaction.final_state)
builder.scalar_initial_state_mass = True
bw_builder = RelativisticBreitWignerBuilder()
for name in reaction.get_intermediate_particles().names:
    builder.set_dynamics(name, bw_builder)
standard_model = builder.formulate()
standard_model.intensity
\[\displaystyle \sum_{m_{A}=-1/2}^{1/2} \sum_{m_{0}=-1/2}^{1/2} \sum_{m_{1}=0} \sum_{m_{2}=0}{\left|{{A^{01}}_{m_{A},m_{0},m_{1},m_{2}} + {A^{02}}_{m_{A},m_{0},m_{1},m_{2}} + {A^{12}}_{m_{A},m_{0},m_{1},m_{2}}}\right|^{2}}\]
Hide code cell source 
import sympy as sp
from IPython.display import Math, display

for i, (symbol, expr) in enumerate(standard_model.amplitudes.items()):
    if i == 3:
        display(Math(R"\dots"))
        break
    latex = sp.multiline_latex(symbol, expr, environment="eqnarray")
    display(Math(latex))

Importing the parameter values given by Table 1:

Hide code cell content 
from ampform.helicity import HelicityModel

# fmt: off
parameter_table = {
    # K*
    R"C_{\Lambda_{c}^{+} \to K^*_{0} p_{+1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}}": 1,
    R"C_{\Lambda_{c}^{+} \to K^*_{+1} p_{+1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}}": 0.5 + 0.5j,
    R"C_{\Lambda_{c}^{+} \to K^*_{-1} p_{-1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}}": 1j,
    R"C_{\Lambda_{c}^{+} \to K^*_{0} p_{-1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}}": -0.5 - 0.5j,
    "m_{K^*}": 0.9,  # GeV
    R"\Gamma_{K^*}": 0.2,  # GeV
    # Λ*
    R"C_{\Lambda_{c}^{+} \to \Lambda^*_{-1/2} \pi^{+}_{0}; \Lambda^* \to K^{-}_{0} p_{+1/2}}": 1j,
    R"C_{\Lambda_{c}^{+} \to \Lambda^*_{+1/2} \pi^{+}_{0}; \Lambda^* \to K^{-}_{0} p_{+1/2}}": 0.8 - 0.4j,
    R"m_{\Lambda^*}": 1.6,  # GeV
    R"\Gamma_{\Lambda^*}": 0.2,  # GeV
    # Δ*
    R"C_{\Lambda_{c}^{+} \to \Delta^*_{+1/2} K^{-}_{0}; \Delta^* \to p_{+1/2} \pi^{+}_{0}}": 0.6 - 0.4j,
    R"C_{\Lambda_{c}^{+} \to \Delta^*_{-1/2} K^{-}_{0}; \Delta^* \to p_{+1/2} \pi^{+}_{0}}": 0.1j,
    R"m_{\Delta^*}": 1.4,  # GeV
    R"\Gamma_{\Delta^*}": 0.2,  # GeV
}
# fmt: on


def set_coefficients(model: HelicityModel) -> None:
    for name, value in parameter_table.items():
        model.parameter_defaults[name] = value
Hide code cell source 
set_coefficients(standard_model)

latex = R"\begin{array}{lc}" + "\n"
for par_name, value in parameter_table.items():
    value = str(value).lstrip("(").rstrip(")").replace("j", "i")
    symbol = sp.Symbol(par_name)
    latex += Rf"  {sp.latex(symbol)} & {value} \\" + "\n"
latex += R"\end{array}"
Math(latex)
\[\begin{split}\displaystyle \begin{array}{lc} C_{\Lambda_{c}^{+} \to K^*_{0} p_{+1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}} & 1 \\ C_{\Lambda_{c}^{+} \to K^*_{+1} p_{+1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}} & 0.5+0.5i \\ C_{\Lambda_{c}^{+} \to K^*_{-1} p_{-1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}} & 1i \\ C_{\Lambda_{c}^{+} \to K^*_{0} p_{-1/2}; K^* \to K^{-}_{0} \pi^{+}_{0}} & -0.5-0.5i \\ m_{K^*} & 0.9 \\ \Gamma_{K^*} & 0.2 \\ C_{\Lambda_{c}^{+} \to \Lambda^*_{-1/2} \pi^{+}_{0}; \Lambda^* \to K^{-}_{0} p_{+1/2}} & 1i \\ C_{\Lambda_{c}^{+} \to \Lambda^*_{+1/2} \pi^{+}_{0}; \Lambda^* \to K^{-}_{0} p_{+1/2}} & 0.8-0.4i \\ m_{\Lambda^*} & 1.6 \\ \Gamma_{\Lambda^*} & 0.2 \\ C_{\Lambda_{c}^{+} \to \Delta^*_{+1/2} K^{-}_{0}; \Delta^* \to p_{+1/2} \pi^{+}_{0}} & 0.6-0.4i \\ C_{\Lambda_{c}^{+} \to \Delta^*_{-1/2} K^{-}_{0}; \Delta^* \to p_{+1/2} \pi^{+}_{0}} & 0.1i \\ m_{\Delta^*} & 1.4 \\ \Gamma_{\Delta^*} & 0.2 \\ \end{array}\end{split}\]

Generate data#

Hide code cell content 
import matplotlib.pyplot as plt
import numpy as np
from tensorwaves.data import SympyDataTransformer
from tensorwaves.function.sympy import create_function


def compute_sub_intensities(
    model: HelicityModel, resonance_name: str, phsp, full_expression
) -> np.ndarray:
    parameter_values = {}
    for symbol, value in model.parameter_defaults.items():
        if resonance_name not in symbol.name and symbol.name.startswith("C"):
            parameter_values[symbol] = 0
        else:
            parameter_values[symbol] = value
    sub_expression = full_expression.subs(parameter_values)
    sub_intensity = create_function(sub_expression, backend="jax")
    return np.array(sub_intensity(phsp).real)


def plot_distributions(model: HelicityModel) -> None:
    helicity_transformer = SympyDataTransformer.from_sympy(
        model.kinematic_variables, backend="jax"
    )
    phsp = helicity_transformer(phsp_momenta)
    phsp = {k: v.real for k, v in phsp.items()}

    full_expression = model.expression.doit()
    substituted_expression = full_expression.xreplace(model.parameter_defaults)
    intensity_func = create_function(substituted_expression, backend="jax")
    intensities_all = np.array(intensity_func(phsp).real)
    intensities_k = compute_sub_intensities(model, "K^*", phsp, full_expression)
    intensities_delta = compute_sub_intensities(model, "Delta^*", phsp, full_expression)
    intensities_lambda = compute_sub_intensities(
        model, "Lambda^*", phsp, full_expression
    )

    fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
    hist_kwargs = {
        "bins": 80,
        "histtype": "step",
    }

    for x in ax.flatten():
        x.set_yticks([])

    ax[0, 0].set_xlabel("$m^2(pK^-)$ [GeV$^2/c^4$]")
    ax[0, 1].set_xlabel(R"$m^2(K^-\pi^+)$ [GeV$^2/c^4$]")
    ax[0, 2].set_xlabel(R"$m^2(p\pi^+)$ [GeV$^2/c^4$]")
    ax[1, 0].set_xlabel(R"$\cos\theta(p)$")
    ax[1, 1].set_xlabel(R"$\phi(p)$")
    ax[1, 2].set_xlabel(R"$\chi$")

    for x, xticks in {
        ax[0, 0]: [2, 2.5, 3, 3.5, 4, 4.5],
        ax[0, 1]: [0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2],
        ax[0, 2]: [1, 1.5, 2, 2.5, 3],
        ax[1, 0]: [-1, -0.5, 0, 0.5, 1],
        ax[1, 1]: [-3, -2, -1, 0, 1, 2, 3],
    }.items():
        x.set_xticks(xticks)
        x.set_xticklabels(xticks)

    for weights, color, label in [
        (intensities_all, "red", "Model"),
        (intensities_k, "orange", R"$K^*\to\,K^{^-}\pi^+$"),
        (intensities_delta, "brown", R"$\Delta^{*^{++}} \to\,p\pi^+$"),
        (intensities_lambda, "purple", R"$\Lambda^* \to\,p K^{^-}$"),
    ]:
        kwargs = dict(weights=weights, color=color, **hist_kwargs)
        ax[0, 0].hist(np.array(phsp["m_01"] ** 2), **kwargs)
        ax[0, 1].hist(np.array(phsp["m_12"] ** 2), **kwargs)
        ax[0, 2].hist(np.array(phsp["m_02"] ** 2), **kwargs)
        ax[1, 0].hist(np.array(np.cos(phsp["theta_01"])), **kwargs)
        ax[1, 1].hist(np.array(phsp["phi_01"]), **kwargs, label=label)

    ax[1, 2].remove()
    handles, labels = ax[1, 1].get_legend_handles_labels()
    fig.legend(handles, labels, loc="lower right")

    ax[0, 2].set_xlim(1, 3.4)
    ax[1, 0].set_xlim(-1, +1)
    ax[1, 1].set_xlim(-np.pi, +np.pi)

    fig.tight_layout()

    plt.show()

Warning

It takes several minutes to lambdify the full expression and expressions for the Wigner rotation angles.

Hide code cell source 
plot_distributions(standard_model)

013_23_0

Spin alignment sum#

Now, with the spin alignment sum from ampform#245 inserted:

builder.align_spin = True
aligned_model = builder.formulate()
set_coefficients(aligned_model)
aligned_model.intensity
\[\displaystyle \sum_{m_{A}=-1/2}^{1/2} \sum_{m_{0}=-1/2}^{1/2} \sum_{m_{1}=0} \sum_{m_{2}=0}{\left|{\sum_{\lambda^{01}_{0}=-1/2}^{1/2} \sum_{\mu^{01}_{0}=-1/2}^{1/2} \sum_{\nu^{01}_{0}=-1/2}^{1/2} \sum_{\lambda^{01}_{1}=0} \sum_{\mu^{01}_{1}=0} \sum_{\nu^{01}_{1}=0} \sum_{\lambda^{01}_{2}=0}{{A^{01}}_{m_{A},\lambda^{01}_{0},- \lambda^{01}_{1},- \lambda^{01}_{2}} D^{0}_{m_{1},\nu^{01}_{1}}\left(\alpha^{01}_{1},\beta^{01}_{1},\gamma^{01}_{1}\right) D^{0}_{m_{2},\lambda^{01}_{2}}\left(\phi_{01},\theta_{01},0\right) D^{0}_{\mu^{01}_{1},\lambda^{01}_{1}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{0}_{\nu^{01}_{1},\mu^{01}_{1}}\left(\phi_{01},\theta_{01},0\right) D^{\frac{1}{2}}_{m_{0},\nu^{01}_{0}}\left(\alpha^{01}_{0},\beta^{01}_{0},\gamma^{01}_{0}\right) D^{\frac{1}{2}}_{\mu^{01}_{0},\lambda^{01}_{0}}\left(\phi^{01}_{0},\theta^{01}_{0},0\right) D^{\frac{1}{2}}_{\nu^{01}_{0},\mu^{01}_{0}}\left(\phi_{01},\theta_{01},0\right)} + \sum_{\lambda^{02}_{0}=-1/2}^{1/2} \sum_{\mu^{02}_{0}=-1/2}^{1/2} \sum_{\nu^{02}_{0}=-1/2}^{1/2} \sum_{\lambda^{02}_{1}=0} \sum_{\lambda^{02}_{2}=0} \sum_{\mu^{02}_{2}=0} \sum_{\nu^{02}_{2}=0}{{A^{02}}_{m_{A},\lambda^{02}_{0},- \lambda^{02}_{1},- \lambda^{02}_{2}} D^{0}_{m_{1},\lambda^{02}_{1}}\left(\phi_{02},\theta_{02},0\right) D^{0}_{m_{2},\nu^{02}_{2}}\left(\alpha^{02}_{2},\beta^{02}_{2},\gamma^{02}_{2}\right) D^{0}_{\mu^{02}_{2},\lambda^{02}_{2}}\left(\phi^{02}_{0},\theta^{02}_{0},0\right) D^{0}_{\nu^{02}_{2},\mu^{02}_{2}}\left(\phi_{02},\theta_{02},0\right) D^{\frac{1}{2}}_{m_{0},\nu^{02}_{0}}\left(\alpha^{02}_{0},\beta^{02}_{0},\gamma^{02}_{0}\right) D^{\frac{1}{2}}_{\mu^{02}_{0},\lambda^{02}_{0}}\left(\phi^{02}_{0},\theta^{02}_{0},0\right) D^{\frac{1}{2}}_{\nu^{02}_{0},\mu^{02}_{0}}\left(\phi_{02},\theta_{02},0\right)} + \sum_{\lambda^{12}_{0}=-1/2}^{1/2} \sum_{\lambda^{12}_{1}=0} \sum_{\mu^{12}_{1}=0} \sum_{\nu^{12}_{1}=0} \sum_{\lambda^{12}_{2}=0} \sum_{\mu^{12}_{2}=0} \sum_{\nu^{12}_{2}=0}{{A^{12}}_{m_{A},\lambda^{12}_{0},\lambda^{12}_{1},- \lambda^{12}_{2}} D^{0}_{m_{1},\nu^{12}_{1}}\left(\alpha^{12}_{1},\beta^{12}_{1},\gamma^{12}_{1}\right) D^{0}_{m_{2},\nu^{12}_{2}}\left(\alpha^{12}_{2},\beta^{12}_{2},\gamma^{12}_{2}\right) D^{0}_{\mu^{12}_{1},\lambda^{12}_{1}}\left(\phi^{12}_{1},\theta^{12}_{1},0\right) D^{0}_{\mu^{12}_{2},\lambda^{12}_{2}}\left(\phi^{12}_{1},\theta^{12}_{1},0\right) D^{0}_{\nu^{12}_{1},\mu^{12}_{1}}\left(\phi_{0},\theta_{0},0\right) D^{0}_{\nu^{12}_{2},\mu^{12}_{2}}\left(\phi_{0},\theta_{0},0\right) D^{\frac{1}{2}}_{m_{0},\lambda^{12}_{0}}\left(\phi_{0},\theta_{0},0\right)}}\right|^{2}}\]
Hide code cell source 
plot_distributions(aligned_model)

013_28_0

Compare with Figure 2. Note that the distributions differ close to threshold, because the distributions in the paper are produced with form factors and an energy-dependent width.