5. Uncertainties#

Hide code cell content
from __future__ import annotations

import itertools
import json
import logging
import os
import tarfile
from collections import defaultdict
from textwrap import wrap

import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
import sympy as sp
import yaml
from ampform.kinematics.phasespace import is_within_phasespace
from ampform.sympy import PoolSum
from IPython.display import Latex, Markdown
from matplotlib import cm
from matplotlib.ticker import FuncFormatter, MultipleLocator
from plotly.subplots import make_subplots
from scipy.interpolate import griddata
from tensorwaves.function.sympy import create_function
from tensorwaves.interface import DataSample, ParametrizedFunction
from tqdm.auto import tqdm

from polarimetry import formulate_polarimetry
from polarimetry.amplitude import AmplitudeModel
from polarimetry.data import (
    create_data_transformer,
    generate_meshgrid_sample,
    generate_phasespace_sample,
)
from polarimetry.function import integrate_intensity, sub_intensity
from polarimetry.io import (
    export_polarimetry_field,
    mute_jax_warnings,
    perform_cached_doit,
    perform_cached_lambdify,
)
from polarimetry.lhcb import (
    ParameterBootstrap,
    load_model_builder,
    load_model_parameters,
)
from polarimetry.lhcb.particle import load_particles
from polarimetry.plot import add_watermark, use_mpl_latex_fonts

logging.getLogger("polarimetry.io").setLevel(logging.INFO)
mute_jax_warnings()
FUNCTION_CACHE: dict[sp.Expr, ParametrizedFunction] = {}
UNFOLDED_POOLSUM_CACHE: dict[PoolSum, sp.Expr] = {}

NO_TQDM = "EXECUTE_NB" in os.environ
if NO_TQDM:
    logging.getLogger().setLevel(logging.ERROR)
    logging.getLogger("polarimetry.io").setLevel(logging.ERROR)

5.1. Model loading#

Hide code cell source
model_file = "../data/model-definitions.yaml"
particles = load_particles("../data/particle-definitions.yaml")
reference_subsystem = 1
with open(model_file) as f:
    model_titles = list(yaml.safe_load(f))

models = {}
for title in tqdm(model_titles, desc="Formulating models", disable=NO_TQDM):
    amplitude_builder = load_model_builder(model_file, particles, model_id=title)
    model = amplitude_builder.formulate(reference_subsystem)
    imported_parameter_values = load_model_parameters(
        model_file, model.decay, title, particles
    )
    model.parameter_defaults.update(imported_parameter_values)
    models[title] = model
Hide code cell source
def unfold_poolsum(expr: PoolSum) -> sp.Expr:
    unfolded_expr = UNFOLDED_POOLSUM_CACHE.get(expr)
    if unfolded_expr is None:
        unfolded_expr = expr.doit()
        UNFOLDED_POOLSUM_CACHE[expr] = unfolded_expr
        return unfolded_expr
    return unfolded_expr


nominal_model_title = "Default amplitude model"
nominal_model = models[nominal_model_title]
unfolded_exprs = {}
for title, model in tqdm(
    models.items(), desc="Unfolding expressions", disable=NO_TQDM
):
    amplitude_builder = load_model_builder(model_file, particles, model_id=title)
    polarimetry_exprs = formulate_polarimetry(amplitude_builder, reference_subsystem)
    folded_exprs = {f"alpha_{i}": expr for i, expr in zip("xyz", polarimetry_exprs)}
    folded_exprs["intensity"] = model.intensity
    unfolded_exprs[title] = {
        key: perform_cached_doit(unfold_poolsum(expr).xreplace(model.amplitudes))
        for key, expr in tqdm(
            folded_exprs.items(), disable=NO_TQDM, leave=False, postfix=title
        )
    }
Hide code cell source
expression_hashes = {
    title: hash(exprs["intensity"]) for title, exprs in unfolded_exprs.items()
}
table = R"""
$$
\begin{array}{rllr}
  && \textbf{Model description} & n\textbf{ ops.} \\
  \hline
"""
for i, (title, expressions) in enumerate(unfolded_exprs.items()):
    expr = expressions["intensity"]
    h = hash(expr)
    for j, v in enumerate(expression_hashes.values()):
        if h == v:
            break
    same_as = ""
    if i != j:
        same_as = f"= {j}"
    ops = sp.count_ops(expr)
    title = "".join(Rf"\text{{{t}}}\\ " for t in wrap(title, width=56))
    title = Rf"\begin{{array}}{{l}}{title}\end{{array}}"
    title = title.replace("^", R"\^{}")
    table += Rf"  \mathbf{{{i}}} & {same_as} & {title} & {ops:,} \\ \hline" "\n"
table += R"""\end{array}
$$
"""

n_models = len(models)
n_unique_hashes = len(set(expression_hashes.values()))
src = Rf"""
Of the {n_models} models, there are {n_unique_hashes} with a unique expression tree.
"""
if NO_TQDM:
    src += "\n:::{dropdown} Show number of mathematical operations per model\n"
src += table
if NO_TQDM:
    src += "\n:::"
Markdown(src.strip())

Of the 18 models, there are 9 with a unique expression tree.

Show number of mathematical operations per model
\[\begin{split} \begin{array}{rllr} && \textbf{Model description} & n\textbf{ ops.} \\ \hline \mathbf{0} & & \begin{array}{l}\text{Default amplitude model}\\ \end{array} & 43,198 \\ \hline \mathbf{1} & = 0 & \begin{array}{l}\text{Alternative amplitude model with K(892) with free mass}\\ \text{and width}\\ \end{array} & 43,198 \\ \hline \mathbf{2} & = 0 & \begin{array}{l}\text{Alternative amplitude model with L(1670) with free mass}\\ \text{and width}\\ \end{array} & 43,198 \\ \hline \mathbf{3} & = 0 & \begin{array}{l}\text{Alternative amplitude model with L(1690) with free mass}\\ \text{and width}\\ \end{array} & 43,198 \\ \hline \mathbf{4} & = 0 & \begin{array}{l}\text{Alternative amplitude model with D(1232) with free mass}\\ \text{and width}\\ \end{array} & 43,198 \\ \hline \mathbf{5} & = 0 & \begin{array}{l}\text{Alternative amplitude model with L(1600), D(1600),}\\ \text{D(1700) with free mass and width}\\ \end{array} & 43,198 \\ \hline \mathbf{6} & = 0 & \begin{array}{l}\text{Alternative amplitude model with free L(1405) Flatt'e}\\ \text{widths, indicated as G1 (pK channel) and G2 (Sigmapi)}\\ \end{array} & 43,198 \\ \hline \mathbf{7} & & \begin{array}{l}\text{Alternative amplitude model with L(1800) contribution}\\ \text{added with free mass and width}\\ \end{array} & 44,222 \\ \hline \mathbf{8} & & \begin{array}{l}\text{Alternative amplitude model with L(1810) contribution}\\ \text{added with free mass and width}\\ \end{array} & 46,782 \\ \hline \mathbf{9} & & \begin{array}{l}\text{Alternative amplitude model with D(1620) contribution}\\ \text{added with free mass and width}\\ \end{array} & 44,222 \\ \hline \mathbf{10} & & \begin{array}{l}\text{Alternative amplitude model in which a Relativistic}\\ \text{Breit-Wigner is used for the K(700) contribution}\\ \end{array} & 43,470 \\ \hline \mathbf{11} & = 0 & \begin{array}{l}\text{Alternative amplitude model with K(700) with free mass}\\ \text{and width}\\ \end{array} & 43,198 \\ \hline \mathbf{12} & & \begin{array}{l}\text{Alternative amplitude model with K(1410) contribution}\\ \text{added with mass and width from PDG2020}\\ \end{array} & 46,780 \\ \hline \mathbf{13} & & \begin{array}{l}\text{Alternative amplitude model in which a Relativistic}\\ \text{Breit-Wigner is used for the K(1430) contribution}\\ \end{array} & 43,470 \\ \hline \mathbf{14} & = 0 & \begin{array}{l}\text{Alternative amplitude model with K(1430) with free width}\\ \end{array} & 43,198 \\ \hline \mathbf{15} & & \begin{array}{l}\text{Alternative amplitude model with an additional overall}\\ \text{exponential form factor exp(-alpha q\^{}2) multiplying Bugg}\\ \text{lineshapes. The exponential parameter is indicated as}\\ \text{``alpha''}\\ \end{array} & 43,582 \\ \hline \mathbf{16} & = 0 & \begin{array}{l}\text{Alternative amplitude model with free radial parameter d}\\ \text{for the Lc resonance, indicated as dLc}\\ \end{array} & 43,198 \\ \hline \mathbf{17} & & \begin{array}{l}\text{Alternative amplitude model obtained using LS couplings}\\ \end{array} & 110,839 \\ \hline \end{array} \end{split}\]
Hide code cell source
def cached_lambdify(expr: sp.Expr, model: AmplitudeModel) -> ParametrizedFunction:
    func = FUNCTION_CACHE.get(expr)
    if func is None:
        func = perform_cached_lambdify(
            expr,
            parameters=model.parameter_defaults,
            backend="jax",
        )
        FUNCTION_CACHE[expr] = func
    str_parameters = {str(k): v for k, v in model.parameter_defaults.items()}
    func.update_parameters(str_parameters)
    return func


jax_functions = {}
original_parameters: dict[str, dict[str, complex | float | int]] = {}
progress_bar = tqdm(desc="Lambdifying to JAX", disable=NO_TQDM, total=len(models))
for title, model in models.items():
    progress_bar.set_postfix_str(title)
    jax_functions[title] = {
        key: cached_lambdify(expr, model)
        for key, expr in unfolded_exprs[title].items()
    }
    original_parameters[title] = dict(jax_functions[title]["intensity"].parameters)
    progress_bar.update()
progress_bar.set_postfix_str("")
progress_bar.close()

5.2. Statistical uncertainties#

5.2.1. Parameter bootstrapping#

n_bootstraps = 100
nominal_functions = jax_functions[nominal_model_title]
bootstrap = ParameterBootstrap(model_file, nominal_model.decay, nominal_model_title)
bootstrap_parameters = bootstrap.create_distribution(n_bootstraps, seed=0)
n_events = 100_000
transformer = create_data_transformer(nominal_model)
phsp_sample = generate_phasespace_sample(nominal_model.decay, n_events, seed=0)
phsp_sample = transformer(phsp_sample)
Hide code cell source
resonances = [chain.resonance for chain in nominal_model.decay.chains]
nominal_parameters = dict(original_parameters[nominal_model_title])
stat_grids = defaultdict(list)
stat_decay_rates = defaultdict(list)
for i in tqdm(
    range(n_bootstraps),
    desc="Computing polarimetry and intensities for parameter combinations",
    disable=NO_TQDM,
):
    new_parameters = {k: v[i] for k, v in bootstrap_parameters.items()}
    for key, func in nominal_functions.items():
        func.update_parameters(nominal_parameters)
        func.update_parameters(new_parameters)
        stat_grids[key].append(func(phsp_sample).real)
    I_tot = integrate_intensity(stat_grids["intensity"][-1])
    for resonance in resonances:
        res_filter = resonance.name.replace("(", R"\(").replace(")", R"\)")
        I_sub = sub_intensity(
            nominal_functions["intensity"], phsp_sample, [res_filter]
        )
        stat_decay_rates[resonance.name].append(I_sub / I_tot)

stat_intensities = jnp.array(stat_grids["intensity"])
stat_polarimetry = jnp.array([stat_grids[f"alpha_{i}"] for i in "xyz"])
stat_polarimetry_norm = jnp.sqrt(jnp.sum(stat_polarimetry**2, axis=0))
stat_decay_rates = {k: jnp.array(v) for k, v in stat_decay_rates.items()}

5.2.2. Mean and standard deviations#

assert stat_intensities.shape == (n_bootstraps, n_events)
assert stat_polarimetry.shape == (3, n_bootstraps, n_events)
assert stat_polarimetry_norm.shape == (n_bootstraps, n_events)
n_bootstraps, n_events
(100, 100000)
Hide code cell source
stat_alpha_mean = [
    jnp.mean(stat_polarimetry_norm, axis=0),
    *jnp.mean(stat_polarimetry, axis=1),
]
stat_alpha_std = [
    jnp.std(stat_polarimetry_norm, axis=0),
    *jnp.std(stat_polarimetry, axis=1),
]

stat_alpha_times_I_mean = [
    jnp.mean(stat_polarimetry_norm * stat_intensities, axis=0),
    *jnp.mean(stat_polarimetry * stat_intensities, axis=1),
]
stat_alpha_times_I_std = [
    jnp.std(stat_polarimetry_norm * stat_intensities, axis=0),
    *jnp.std(stat_polarimetry * stat_intensities, axis=1),
]
stat_alpha_times_I_mean = jnp.array(stat_alpha_times_I_mean)
stat_alpha_times_I_std = jnp.array(stat_alpha_times_I_std)

stat_intensity_mean = jnp.mean(stat_intensities, axis=0)
stat_intensity_std = jnp.std(stat_intensities, axis=0)

5.2.3. Distributions#

Hide code cell source
def interpolate_to_grid(values: np.ndarray, method: str = "linear"):
    return griddata(POINTS, values, (X, Y))


resolution = 200
POINTS = np.transpose(
    [
        phsp_sample["sigma1"],
        phsp_sample["sigma2"],
    ]
)
grid_sample = generate_meshgrid_sample(nominal_model.decay, resolution)
X = np.array(grid_sample["sigma1"])
Y = np.array(grid_sample["sigma2"])
Hide code cell content
def create_contour(phsp: DataSample) -> jnp.ndarray:
    # See also https://compwa-org.rtfd.io/report/017.html
    m0, m1, m2, m3, s1, s2 = sp.symbols("m(:4) sigma(1:3)", nonnegative=True)
    filter_expr = is_within_phasespace(
        s1,
        s2,
        nominal_model.parameter_defaults[m0],
        nominal_model.parameter_defaults[m1],
        nominal_model.parameter_defaults[m2],
        nominal_model.parameter_defaults[m3],
        outside_value=0,
    ).doit()
    filter_func = create_function(filter_expr, backend="jax")
    return filter_func(contour_sample)  # noqa: F821


def draw_dalitz_contour(
    ax, color: str = "black", width: float = 0.1, **kwargs
) -> None:
    ax.contour(
        X_CONTOUR,
        Y_CONTOUR,
        Z_CONTOUR,
        colors=color,
        linewidths=width,
        **kwargs,
    )


contour_sample = generate_meshgrid_sample(nominal_model.decay, resolution=500)
X_CONTOUR = np.array(contour_sample["sigma1"])
Y_CONTOUR = np.array(contour_sample["sigma2"])
Z_CONTOUR = create_contour(contour_sample)
del contour_sample, create_contour
Hide code cell source
%config InlineBackend.figure_formats = ['png']
plt.rcdefaults()
use_mpl_latex_fonts()
plt.rc("font", size=18)
fig, axes = plt.subplots(
    dpi=200,
    figsize=(16.7, 8),
    gridspec_kw={"width_ratios": [1, 1, 1, 1.18]},
    ncols=4,
    nrows=2,
    sharex=True,
    sharey=True,
)
plt.subplots_adjust(hspace=0.02, wspace=0.02)
fig.suptitle(R"Polarimetry field $\vec\alpha$ (statistical \& systematic)")
s1_label = R"$m^2\left(K^-\pi^+\right)$ [GeV$^2$]"
s2_label = R"$m^2\left(pK^-\right)$ [GeV$^2$]"
s3_label = R"$m^2\left(p\pi^+\right)$ [GeV$^2$]"
axes[0, 0].set_ylabel(s2_label)
axes[1, 0].set_ylabel(s2_label)

global_max_std = max(map(jnp.nanmax, stat_alpha_std))
for i in range(4):
    if i != 0:
        title = Rf"$\alpha_{'xyz'[i-1]}$"
    else:
        title = R"$\left|\vec\alpha\right|$"
    axes[0, i].set_title(title)

    draw_dalitz_contour(axes[0, i], zorder=10)
    Z = interpolate_to_grid(stat_alpha_mean[i])
    mesh = axes[0, i].pcolormesh(X, Y, Z, cmap=cm.RdYlGn_r)
    mesh.set_clim(vmin=-1, vmax=+1)
    if axes[0, i] is axes[0, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[0, i], pad=0.01)
        c_bar.ax.set_ylabel(Rf"$\alpha$ averaged with {n_bootstraps} bootstraps")
        c_bar.ax.set_yticks([-1, 0, +1])
        c_bar.ax.set_yticklabels(["-1", "0", "+1"])

    draw_dalitz_contour(axes[1, i], zorder=10)
    Z = interpolate_to_grid(stat_alpha_std[i])
    mesh = axes[1, i].pcolormesh(X, Y, Z)
    mesh.set_clim(vmin=0, vmax=global_max_std)
    axes[1, i].set_xlabel(s1_label)
    if axes[1, i] is axes[1, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[1, i], pad=0.01)
        c_bar.ax.set_ylabel("standard deviation")
plt.show()
_images/8dbb18a62825321312188fa9f23b899dfb6a07e0018aeef9ef4db22ce83db787.png
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, axes = plt.subplots(
    dpi=200,
    figsize=(16.7, 8),
    gridspec_kw={"width_ratios": [1, 1, 1, 1.18]},
    ncols=4,
    nrows=2,
    sharex=True,
    sharey=True,
)
plt.subplots_adjust(hspace=0.02, wspace=0.02)
fig.suptitle(R"$\vec\alpha \cdot I$ distributions (statistical \& systematic)")
axes[0, 0].set_ylabel(s2_label)
axes[1, 0].set_ylabel(s2_label)

global_max_mean = jnp.nanmax(jnp.abs(stat_alpha_times_I_mean))
global_max_std = jnp.nanmax(stat_alpha_times_I_std / stat_intensity_mean)
for i in range(4):
    if i != 0:
        title = Rf"$\alpha_{'xyz'[i-1]}$"
    else:
        title = R"$\left|\vec\alpha\right|$"
    axes[0, i].set_title(title)
    axes[1, i].set_xlabel(s1_label)

    draw_dalitz_contour(axes[0, i], zorder=10)
    Z = interpolate_to_grid(stat_alpha_times_I_mean[i])
    mesh = axes[0, i].pcolormesh(X, Y, Z, cmap=cm.RdYlGn_r)
    mesh.set_clim(vmin=-global_max_mean, vmax=+global_max_mean)
    if axes[0, i] is axes[0, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[0, i], pad=0.01)
        c_bar.ax.set_ylabel(
            Rf"$\alpha \cdot I$ averaged with {n_bootstraps} bootstraps"
        )

    draw_dalitz_contour(axes[1, i], zorder=10)
    Z = interpolate_to_grid(stat_alpha_times_I_std[i] / stat_intensity_mean)
    mesh = axes[1, i].pcolormesh(X, Y, Z)
    mesh.set_clim(vmin=0, vmax=global_max_std)
    if axes[1, i] is axes[1, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[1, i], pad=0.01)
        c_bar.ax.set_ylabel("standard deviation / intensity")
plt.show()
_images/9b7c51bd7e09763b97a7d67afbca16d7cab4bf3a92404d8c5cc06515241708ea.png
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, (ax1, ax2) = plt.subplots(
    dpi=200,
    figsize=(12, 6.2),
    ncols=2,
    sharey=True,
)
fig.suptitle(R"Intensity distribution (statistical \& systematics)", y=0.95)
ax1.set_xlabel(s1_label)
ax2.set_xlabel(s1_label)
ax1.set_ylabel(s2_label)

Z = interpolate_to_grid(stat_intensity_mean)
mesh = ax1.pcolormesh(X, Y, Z, cmap=cm.Reds)
fig.colorbar(mesh, ax=ax1, pad=0.01)
draw_dalitz_contour(ax1, width=0.2)
ax1.set_title(f"average of {n_bootstraps} bootstraps")

Z = interpolate_to_grid(stat_intensity_std / stat_intensity_mean)
mesh = ax2.pcolormesh(X, Y, Z)
fig.colorbar(mesh, ax=ax2, pad=0.01)
draw_dalitz_contour(ax2, width=0.2)
ax2.set_title("standard deviation / intensity")
fig.tight_layout()
plt.show()
_images/37cfbcf39e7401e0f7befed30d58226e747e34418a97dc087d730922404551b2.png

5.2.4. Comparison with nominal values#

Hide code cell source
for func in nominal_functions.values():
    func.update_parameters(nominal_parameters)
nominal_intensity = nominal_functions["intensity"](phsp_sample)
nominal_polarimetry = jnp.array(
    [nominal_functions[f"alpha_{i}"](phsp_sample).real for i in "xyz"]
)
nominal_polarimetry_norm = jnp.sqrt(jnp.sum(nominal_polarimetry**2, axis=0))
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, axes = plt.subplots(
    dpi=200,
    figsize=(17.3, 4),
    gridspec_kw={"width_ratios": [1, 1, 1, 1.2]},
    ncols=4,
    sharey=True,
)
plt.subplots_adjust(hspace=0.2, wspace=0.05)
fig.suptitle("Comparison with nominal values", y=1.04)
axes[0].set_ylabel(s2_label)
for ax in axes:
    ax.set_xlabel(s1_label)

vmax = 5.0  # %
for i in range(4):
    if i != 0:
        title = Rf"$\alpha_{'xyz'[i-1]}$"
        z_values = jnp.abs(
            (stat_alpha_mean[i] - nominal_polarimetry[i - 1])
            / nominal_polarimetry[i - 1]
        )
    else:
        title = "$I$"
        z_values = 100 * jnp.abs(
            (stat_intensity_mean - nominal_intensity) / nominal_intensity
        )
    axes[i].set_title(title)

    draw_dalitz_contour(axes[i], zorder=10)
    Z = interpolate_to_grid(z_values)
    mesh = axes[i].pcolormesh(X, Y, Z, cmap=cm.Reds)
    mesh.set_clim(vmin=0, vmax=vmax)
    if axes[i] is axes[-1]:
        c_bar = fig.colorbar(mesh, ax=axes[i], pad=0.02)
        c_bar.ax.set_ylabel(R"difference with nominal value (\%)")

plt.show()
_images/41ca9dde584ff16080c095ad554cb2a26e65674aae7fb7580afc15c56cf82436.png

5.3. Systematic uncertainties#

Hide code cell content
syst_grids = defaultdict(list)
syst_decay_rates = defaultdict(list)
progress_bar = tqdm(desc="Computing systematics", disable=NO_TQDM, total=len(models))
for title, model in models.items():
    progress_bar.set_postfix_str(title)
    for key, func in jax_functions[title].items():
        func.update_parameters(original_parameters[title])
        syst_grids[key].append(func(phsp_sample).real)
    I_tot = integrate_intensity(syst_grids["intensity"][-1])
    intensity_func = jax_functions[title]["intensity"]
    for resonance in resonances:
        res_filter = resonance.name.replace("(", R"\(").replace(")", R"\)")
        I_sub = sub_intensity(intensity_func, phsp_sample, [res_filter])
        syst_decay_rates[resonance.name].append(I_sub / I_tot)
    progress_bar.update()
progress_bar.set_postfix_str("")
progress_bar.close()

syst_intensities = jnp.array(syst_grids["intensity"])
syst_polarimetry = jnp.array([syst_grids[f"alpha_{i}"] for i in "xyz"])
syst_polarimetry_norm = jnp.sqrt(jnp.sum(syst_polarimetry**2, axis=0))
syst_decay_rates = {k: jnp.array(v) for k, v in syst_decay_rates.items()}

5.3.1. Mean and standard deviations#

n_models = len(models)
assert syst_intensities.shape == (n_models, n_events)
assert syst_polarimetry.shape == (3, n_models, n_events)
assert syst_polarimetry_norm.shape == (n_models, n_events)
n_models, n_events
(18, 100000)
Hide code cell source
syst_alpha_mean = [
    jnp.mean(syst_polarimetry_norm, axis=0),
    *jnp.mean(syst_polarimetry, axis=1),
]
alpha_diff_with_model_0 = [
    syst_polarimetry_norm - syst_polarimetry_norm[0],
    *(syst_polarimetry - syst_polarimetry[:, None, 0]),
]

syst_alpha_mean = jnp.array(syst_alpha_mean)
alpha_diff_with_model_0 = jnp.array(alpha_diff_with_model_0)

assert alpha_diff_with_model_0.shape == (4, n_models, n_events)
assert jnp.nanmax(alpha_diff_with_model_0[:, 0]) == 0.0
alpha_syst_extrema = jnp.abs(alpha_diff_with_model_0).max(axis=1)
syst_polarimetry_times_I = [
    syst_polarimetry_norm * syst_intensities,
    *(syst_polarimetry * syst_intensities),
]
syst_polarimetry_times_I = jnp.array(syst_polarimetry_times_I)


syst_alpha_times_I_mean = syst_polarimetry_times_I.mean(axis=1)
syst_alpha_times_I_diff = (
    syst_polarimetry_times_I - syst_polarimetry_times_I[:, None, 0]
)
assert syst_alpha_times_I_diff.shape == (4, n_models, n_events)
assert jnp.nanmax(syst_alpha_times_I_diff[:, 0]) == 0.0
syst_alpha_times_I_extrema = jnp.abs(syst_alpha_times_I_diff).max(axis=1)
intensity_diff_with_model_0 = syst_intensities - syst_intensities[0]
intensity_extrema = jnp.nanmax(intensity_diff_with_model_0, axis=0)

5.3.2. Distributions#

Hide code cell content
def plot_intensity_distributions(sigma: int) -> None:
    original_font_size = plt.rcParams["font.size"]
    use_mpl_latex_fonts()
    plt.rc("font", size=10)
    n_subplots = n_models - 1
    nrows = int(np.floor(np.sqrt(n_subplots)))
    ncols = int(np.ceil(n_subplots / nrows))
    fig, axes = plt.subplots(
        figsize=2.0 * np.array([ncols, nrows]),
        dpi=200,
        ncols=ncols,
        nrows=nrows,
        sharex=True,
        sharey=True,
    )
    fig.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.06)
    x_label = {1: s1_label, 2: s2_label, 3: s3_label}[sigma]
    fig.text(0.5, 0.04, x_label, ha="center")
    fig.text(0.04, 0.5, "$I$ (normalized)", va="center", rotation="vertical")
    for ax in axes.flatten()[n_subplots:]:
        fig.delaxes(ax)
    n_bins = 80
    nominal_bin_values, bin_edges = jnp.histogram(
        phsp_sample[f"sigma{sigma}"],
        weights=syst_intensities[0],
        bins=n_bins,
        density=True,
    )
    for i, (ax, intensities) in enumerate(
        zip(axes.flatten(), syst_intensities[1:]), 1
    ):
        ax.set_title(f"Model {i}", y=0.01)
        ax.set_yticks([])
        bin_values, _ = jnp.histogram(
            phsp_sample[f"sigma{sigma}"],
            weights=intensities,
            bins=n_bins,
            density=True,
        )
        ax.fill_between(
            bin_edges[:-1],
            bin_values,
            alpha=0.5,
            step="pre",
        )
        ax.step(
            x=bin_edges[:-1],
            y=nominal_bin_values,
            linewidth=0.3,
            color="red",
        )
    fig.savefig(f"_images/intensity-distributions-sigma{sigma}.svg")
    plt.show()
    use_mpl_latex_fonts()
    plt.rc("font", size=original_font_size)


%config InlineBackend.figure_formats = ['svg']
plot_intensity_distributions(1)
plot_intensity_distributions(2)
plot_intensity_distributions(3)
_images/79839bfe8eb5da921c6f2d0d21970d5603aad570cc11639d684c778fae95edae.svg_images/df5bffaef62d07c2a8b5508331fa4b94c13dc40823e0119566957b46ca1117b1.svg_images/8fddb1ac9389ddd34531c7619326a1e5d0c43e2462eed5a6ea49d146a48ee225.svg
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, axes = plt.subplots(
    dpi=200,
    figsize=(16.7, 8),
    gridspec_kw={"width_ratios": [1, 1, 1, 1.18]},
    ncols=4,
    nrows=2,
    sharex=True,
    sharey=True,
)
plt.subplots_adjust(hspace=0.02, wspace=0.02)
fig.suptitle(R"Polarimetry field $\vec\alpha$ (model)")
axes[0, 0].set_ylabel(s2_label)
axes[1, 0].set_ylabel(s2_label)

global_max_std = jnp.nanmax(alpha_syst_extrema)
for i in range(4):
    if i != 0:
        title = Rf"$\alpha_{'xyz'[i-1]}$"
    else:
        title = R"$\left|\vec\alpha\right|$"
    axes[0, i].set_title(title)

    draw_dalitz_contour(axes[0, i], zorder=10)
    Z = interpolate_to_grid(syst_alpha_mean[i])
    mesh = axes[0, i].pcolormesh(X, Y, Z, cmap=cm.RdYlGn_r)
    mesh.set_clim(vmin=-1, vmax=+1)
    if axes[0, i] is axes[0, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[0, i], pad=0.01)
        c_bar.ax.set_ylabel(Rf"$\alpha$ averaged with {n_models} models")
        c_bar.ax.set_yticks([-1, 0, +1])
        c_bar.ax.set_yticklabels(["-1", "0", "+1"])

    draw_dalitz_contour(axes[1, i], zorder=10)
    Z = interpolate_to_grid(alpha_syst_extrema[i])
    mesh = axes[1, i].pcolormesh(X, Y, Z)
    mesh.set_clim(vmin=0, vmax=global_max_std)
    axes[1, i].set_xlabel(s1_label)
    if axes[1, i] is axes[1, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[1, i], pad=0.01)
        c_bar.ax.set_ylabel("maximum absolute difference\nwith the default model")
plt.show()
_images/bb1937fec40033d1c56e7e1fbcb0f5197607118646a5fdb51e449ef622fb32b5.png
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, axes = plt.subplots(
    dpi=200,
    figsize=(16.7, 8),
    gridspec_kw={"width_ratios": [1, 1, 1, 1.18]},
    ncols=4,
    nrows=2,
    sharex=True,
    sharey=True,
)
plt.subplots_adjust(hspace=0.02, wspace=0.02)
axes[0, 0].set_ylabel(s2_label)
axes[1, 0].set_ylabel(s2_label)

syst_intensity_mean = jnp.mean(syst_intensities, axis=0)
global_max_mean = jnp.nanmax(jnp.abs(syst_alpha_times_I_mean))
global_max_diff = jnp.nanmax(syst_alpha_times_I_extrema / syst_intensity_mean)
for i in range(4):
    if i != 0:
        title = Rf"$\alpha_{'xyz'[i-1]}$"
    else:
        title = R"$\left|\vec\alpha\right|$"
    axes[0, i].set_title(title)
    axes[1, i].set_xlabel(s1_label)

    draw_dalitz_contour(axes[0, i], zorder=10)
    Z = interpolate_to_grid(syst_alpha_times_I_mean[i])
    mesh = axes[0, i].pcolormesh(X, Y, Z, cmap=cm.RdYlGn_r)
    mesh.set_clim(vmin=-global_max_mean, vmax=+global_max_mean)
    if axes[0, i] is axes[0, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[0, i], pad=0.01)
        c_bar.ax.set_ylabel(Rf"$\alpha \cdot I$ averaged with {n_models} models")

    draw_dalitz_contour(axes[1, i], zorder=10)
    Z = interpolate_to_grid(syst_alpha_times_I_extrema[i] / syst_intensity_mean)
    mesh = axes[1, i].pcolormesh(X, Y, Z)
    mesh.set_clim(vmin=0, vmax=global_max_diff)
    if axes[1, i] is axes[1, -1]:
        c_bar = fig.colorbar(mesh, ax=axes[1, i], pad=0.01)
        c_bar.ax.set_ylabel(
            "maximum absolute difference\n"
            "with the default model\n"
            "divided by nominal intensity"
        )
plt.show()
_images/6353937fb5c8f111ed66c910d65a92157636e62f47b44b11062e175714146f1c.png
Hide code cell source
%config InlineBackend.figure_formats = ['png']
fig, (ax1, ax2) = plt.subplots(
    dpi=200,
    figsize=(12, 6.9),
    ncols=2,
    sharey=True,
)
fig.suptitle("Intensity distribution (model)", y=0.95)
ax1.set_xlabel(s1_label)
ax2.set_xlabel(s1_label)
ax1.set_ylabel(s2_label)

Z = interpolate_to_grid(syst_intensity_mean)
mesh = ax1.pcolormesh(X, Y, Z, cmap=cm.Reds)
fig.colorbar(mesh, ax=ax1, pad=0.01)
draw_dalitz_contour(ax1, width=0.2)
ax1.set_title(f"average of {n_models} models")

Z = interpolate_to_grid(intensity_extrema / syst_intensity_mean)
mesh = ax2.pcolormesh(X, Y, Z)
fig.colorbar(mesh, ax=ax2, pad=0.01)
draw_dalitz_contour(ax2, width=0.2)
ax2.set_title("maximum absolute difference\n" R"with the default model (\%)")
fig.tight_layout()
plt.show()
_images/4e84166c26e485124fa255466d10ef07bcee1d6d1ae327cf9636aa119477f72c.png

5.4. Uncertainty on polarimetry#

For each bootstrap or alternative model \(i\), we compute the angle between each aligned polarimeter vector \(\vec\alpha_i\) and the one from the nominal model, \(\vec\alpha_0\):

\[ \cos\theta_i = \frac{\vec \alpha_i \cdot \vec \alpha_0}{|\alpha_i||\alpha_0|}. \]

The solid angle can then be computed as:

\[ \delta\Omega = \int_0^{2\pi}\int_0^{\theta} \mathrm{d}\phi \, \mathrm{d}\cos\theta = 2\pi\left(1-\cos\theta\right). \]

The statistical uncertainty is given by taking the standard deviation on the \(\delta\Omega\) distribution and the systematic uncertainty is given by taking finding \(\theta_\mathrm{max} = \max\theta_i\) and computing \(\delta\Omega_\mathrm{max}\) from that.

Hide code cell source
def dot(array1: jnp.ndarray, array2: jnp.ndarray, axis=0) -> jnp.ndarray:
    return jnp.sum(array1 * array2, axis=axis)


def norm(array: jnp.ndarray, axis=0) -> jnp.ndarray:
    return jnp.sqrt(dot(array, array, axis=axis))


def compute_theta(polarimetry: jnp.ndarray) -> jnp.ndarray:
    return jnp.arccos(
        dot(polarimetry, nominal_polarimetry[:, None])
        / (norm(polarimetry) * norm(nominal_polarimetry))
    )


def compute_solid_angle(theta: jnp.ndarray) -> jnp.ndarray:
    return 2 * jnp.pi * (1 - jnp.cos(theta))


stat_theta = compute_theta(stat_polarimetry)
syst_theta = compute_theta(syst_polarimetry)
assert stat_theta.shape == (n_bootstraps, n_events)
assert syst_theta.shape == (n_models, n_events)

stat_solid_angle = jnp.nanstd(compute_solid_angle(stat_theta), axis=0)
syst_solid_angle = jnp.nanmax(compute_solid_angle(syst_theta), axis=0)


def plot_angle_uncertainties(watermark: bool, titles: bool) -> None:
    plt.ioff()
    fig, axes = plt.subplots(
        dpi=200,
        figsize=(11, 4),
        ncols=2,
    )
    fig.subplots_adjust(wspace=0.3)
    ax1, ax2 = axes
    if titles:
        fig.suptitle(R"Uncertainty over $\vec\alpha$ polar angle", y=1.04)
        ax1.set_title(R"Statistical \& systematic")
        ax2.set_title("Model")
    for ax in axes:
        ax.set_box_aspect(1)
        ax.set_xlabel(s1_label)
        ax.set_ylabel(s2_label)
        draw_dalitz_contour(ax, width=0.2, zorder=10)

    Z = interpolate_to_grid(stat_solid_angle)
    mesh = ax1.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    mesh.set_clim(0, 4 * np.pi)
    c_bar = fig.colorbar(mesh, ax=ax1, pad=0.02)
    c_bar.ax.set_ylabel(R"Std. of $\delta\Omega$")
    c_bar.ax.set_yticks([0, 2 * np.pi, 4 * np.pi])
    c_bar.ax.set_yticklabels(["$0$", R"$2\pi$", R"$4\pi$"])

    Z = interpolate_to_grid(syst_solid_angle)
    mesh = ax2.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    mesh.set_clim(0, 4 * np.pi)
    c_bar = fig.colorbar(mesh, ax=ax2, pad=0.02)
    c_bar.ax.set_ylabel(R"Max. of $\delta\Omega$")
    c_bar.ax.set_yticks([0, 2 * np.pi, 4 * np.pi])
    c_bar.ax.set_yticklabels(["$0$", R"$2\pi$", R"$4\pi$"])
    output_file = "polarimetry-field-angle-uncertainties"
    if watermark:
        output_file += "-watermark"
        add_watermark_to_uncertainties(ax1)
        add_watermark_to_uncertainties(ax2)
    if titles:
        output_file += "-with-titles"
    fig.savefig(f"_static/images/{output_file}.png", bbox_inches="tight")
    if watermark and titles:
        plt.show()
    plt.close(fig)
    plt.ion()


def add_watermark_to_uncertainties(ax) -> None:
    add_watermark(ax, 0.70, 0.82, fontsize=16)


%config InlineBackend.figure_formats = ['png']
plt.rcdefaults()
use_mpl_latex_fonts()
plt.rc("font", size=16)

boolean_combinations = list(itertools.product(*2 * [[False, True]]))
for watermark, titles in tqdm(boolean_combinations, disable=NO_TQDM):
    plot_angle_uncertainties(watermark, titles)
_images/1e490232e3d285faa979e3bd57a98eca60b8a788c2a6e4fc0f6bb1f7ad7d12bf.png
Hide code cell source
stat_alpha_difference = norm(stat_polarimetry - nominal_polarimetry[:, None])
syst_alpha_difference = norm(syst_polarimetry - nominal_polarimetry[:, None])
nominal_polarimetry_norm = norm(nominal_polarimetry[:, None])
stat_alpha_rel_difference = 100 * stat_alpha_difference / nominal_polarimetry_norm
syst_alpha_rel_difference = 100 * syst_alpha_difference / nominal_polarimetry_norm
assert stat_alpha_rel_difference.shape == (n_bootstraps, n_events)
assert syst_alpha_rel_difference.shape == (n_models, n_events)


def create_figure(titles: bool):
    fig, axes = plt.subplots(
        dpi=200,
        figsize=(11.5, 4),
        ncols=2,
    )
    fig.subplots_adjust(wspace=0.4)
    ax1, ax2 = axes
    if titles:
        fig.suptitle(R"Uncertainty over $\left|\vec\alpha\right|$", y=1.04)
        ax1.set_title(R"Statistical \& systematic")
        ax2.set_title("Model")
    for ax in axes:
        ax.set_box_aspect(1)
        ax.set_xlabel(s1_label)
        ax.set_ylabel(s2_label)
        draw_dalitz_contour(ax, width=0.2)
    return fig, (ax1, ax2)


def plot_norm_rel_uncertainties(watermark: bool, titles: bool) -> None:
    plt.ioff()
    fig, (ax1, ax2) = create_figure(titles)
    Z = interpolate_to_grid(jnp.nanstd(stat_alpha_rel_difference, axis=0))
    mesh = ax1.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    mesh.set_clim(0, 100)
    c_bar = fig.colorbar(mesh, ax=ax1, pad=0.02)
    c_bar.ax.set_ylabel(
        R"Std. of"
        R" $\frac{|\alpha^{(i)}-\alpha^\mathrm{default}|}{|\alpha^\mathrm{default}|}$"
    )
    c_bar.ax.set_yticks([0, 50, 100])
    c_bar.ax.set_yticklabels([R"0\%", R"50\%", R"100\%"])
    Z = interpolate_to_grid(jnp.nanmax(syst_alpha_rel_difference, axis=0))
    mesh = ax2.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    mesh.set_clim(0, 100)
    c_bar = fig.colorbar(mesh, ax=ax2, pad=0.02)
    c_bar.ax.set_ylabel(
        R"Max. of"
        R" $\frac{|\alpha^{(i)}-\alpha^\mathrm{default}|}{|\alpha^\mathrm{default}|}$"
    )
    c_bar.ax.set_yticks([0, 50, 100])
    c_bar.ax.set_yticklabels([R"0\%", R"50\%", R"100\%"])
    output_file = "polarimetry-field-norm-uncertainties-relative"
    if watermark:
        output_file += "-watermark"
        add_watermark_to_uncertainties(ax1)
        add_watermark_to_uncertainties(ax2)
    if titles:
        output_file += "-with-titles"
    fig.savefig(f"_static/images/{output_file}.png", bbox_inches="tight")
    if watermark and titles:
        plt.show()
    plt.close(fig)
    plt.ion()


def plot_norm_uncertainties(watermark: bool, titles: bool) -> None:
    plt.ioff()
    vmax = 0.5
    fig, (ax1, ax2) = create_figure(titles)
    Z = interpolate_to_grid(jnp.nanstd(stat_alpha_difference, axis=0))
    mesh = ax1.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    c_bar = fig.colorbar(mesh, ax=ax1, pad=0.02)
    mesh.set_clim(0, vmax)
    c_bar.ax.set_ylabel(R"Std. of $|\alpha^{(i)}-\alpha^\mathrm{default}|$")
    Z = interpolate_to_grid(jnp.nanmax(syst_alpha_difference, axis=0))
    mesh = ax2.pcolormesh(X, Y, Z, cmap=plt.cm.YlOrRd)
    mesh.set_clim(0, vmax)
    c_bar = fig.colorbar(mesh, ax=ax2, pad=0.02)
    c_bar.ax.set_ylabel(R"Max. of $|\alpha^{(i)}-\alpha^\mathrm{default}|$")
    output_file = "polarimetry-field-norm-uncertainties"
    if watermark:
        output_file += "-watermark"
        add_watermark_to_uncertainties(ax1)
        add_watermark_to_uncertainties(ax2)
    if titles:
        output_file += "-with-titles"
    fig.savefig(f"_static/images/{output_file}.png", bbox_inches="tight")
    if watermark and titles:
        plt.show()
    plt.close(fig)
    plt.ion()


%config InlineBackend.figure_formats = ['png']
plt.rcdefaults()
use_mpl_latex_fonts()
plt.rc("font", size=18)

for watermark, titles in tqdm(boolean_combinations, disable=NO_TQDM):
    plot_norm_rel_uncertainties(watermark, titles)
    plot_norm_uncertainties(watermark, titles)
_images/08c71828411ef998285df5358428d33571fa46a6446611f31e84f522abdc6266.png _images/3ad48ef0c2787dbec1f7eb7810205dc66a233e2f2ccab136624e80ea0ade6d31.png

5.5. Decay rates#

Hide code cell source
with open("../data/observable-references.yaml") as f:
    data = yaml.load(f, Loader=yaml.SafeLoader)
lhcb_values = {
    k: tuple(float(v) for v in row.split(" ± "))
    for k, row in data[nominal_model_title]["rate"].items()
}

src = R"""
\begin{array}{l|c|c}
\textbf{Resonance} & \textbf{Decay rate} & \textbf{LHCb} \\
\hline
"""
for resonance in resonances:
    ff_statistical = 100 * stat_decay_rates[resonance.name]
    ff_systematics = 100 * syst_decay_rates[resonance.name]
    ff_nominal = f"{ff_systematics[0]:.2f}"
    ff_stat = f"{ff_statistical.std():.2f}"
    ff_syst_min = f"{(ff_systematics[1:]-ff_systematics[0]).min():+.2f}"
    ff_syst_max = f"{(ff_systematics[1:]-ff_systematics[0]).max():+.2f}"
    src += f"{resonance.latex} & "
    src += Rf"{ff_nominal} \pm {ff_stat}_{{{ff_syst_min}}}^{{{ff_syst_max}}} & "
    lhcb_value, lhcb_stat, lhcb_syst, _ = lhcb_values[resonance.name]
    src += Rf"{lhcb_value} \pm {lhcb_stat} \pm {lhcb_syst} \\" "\n"
src += R"\end{array}"
Latex(src)
\[\begin{split}\begin{array}{l|c|c} \textbf{Resonance} & \textbf{Decay rate} & \textbf{LHCb} \\ \hline \Lambda(1405) & 7.78 \pm 0.43_{-2.53}^{+3.01} & 7.7 \pm 0.2 \pm 3.0 \\ \Lambda(1520) & 1.91 \pm 0.10_{-0.24}^{+0.04} & 1.86 \pm 0.09 \pm 0.23 \\ \Lambda(1600) & 5.16 \pm 0.28_{-1.93}^{+0.50} & 5.2 \pm 0.2 \pm 1.9 \\ \Lambda(1670) & 1.15 \pm 0.04_{-0.29}^{+0.06} & 1.18 \pm 0.06 \pm 0.32 \\ \Lambda(1690) & 1.16 \pm 0.01_{-0.33}^{+0.06} & 1.19 \pm 0.09 \pm 0.34 \\ \Lambda(2000) & 9.55 \pm 0.67_{-2.26}^{+0.83} & 9.58 \pm 0.27 \pm 0.93 \\ \Delta(1232) & 28.73 \pm 1.34_{-0.79}^{+1.76} & 28.6 \pm 0.29 \pm 0.76 \\ \Delta(1600) & 4.50 \pm 0.51_{-1.40}^{+0.93} & 4.5 \pm 0.3 \pm 1.5 \\ \Delta(1700) & 3.89 \pm 0.07_{-0.48}^{+0.94} & 3.9 \pm 0.2 \pm 0.94 \\ K(700) & 2.99 \pm 0.20_{-0.59}^{+0.91} & 3.02 \pm 0.16 \pm 0.92 \\ K(892) & 21.95 \pm 1.24_{-0.70}^{+0.59} & 22.14 \pm 0.23 \pm 0.64 \\ K(1430) & 14.70 \pm 0.80_{-2.67}^{+2.78} & 14.7 \pm 0.6 \pm 2.7 \\ \end{array}\end{split}\]
Hide code cell source
alignment = "c".join("" for _ in range(n_models) if i)
names = " & ".join(Rf"\textbf{{{i}}}" for i in range(n_models) if i)
src = Rf"""
$$
\begin{{array}}{{l|{alignment}}}
  \textbf{{Resonance}} & {names} \\
  \hline
"""
for resonance in resonances:
    ff_systematics = 100 * syst_decay_rates[resonance.name]
    src += f"  {resonance.latex:13s}"
    for ff_model in ff_systematics[1:]:
        diff = f"{ff_model-ff_systematics[0]:+.2f}"
        if ff_model == ff_systematics[1:].min():
            src += Rf" & \color{{blue}}{{{diff}}}"
        elif ff_model == ff_systematics[1:].max():
            src += Rf" & \color{{red}}{{{diff}}}"
        else:
            src += f" & {diff}"
    src += R" \\" "\n"
src += """\
\\end{array}
$$
"""
for i, title in enumerate(models):
    src += f"\n- **{i}**: {title}"
Markdown(src)
\[\begin{split} \begin{array}{l|ccccccccccccccccc} \textbf{Resonance} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} & \textbf{5} & \textbf{6} & \textbf{7} & \textbf{8} & \textbf{9} & \textbf{10} & \textbf{11} & \textbf{12} & \textbf{13} & \textbf{14} & \textbf{15} & \textbf{16} & \textbf{17} \\ \hline \Lambda(1405) & +0.11 & -0.14 & -0.01 & -0.33 & -0.99 & \color{red}{+3.01} & \color{blue}{-2.53} & -0.66 & -1.58 & -0.43 & -0.01 & -1.97 & -0.11 & -0.04 & -0.18 & +0.06 & -0.75 \\ \Lambda(1520) & +0.03 & +0.00 & +0.01 & +0.01 & \color{blue}{-0.24} & -0.01 & -0.04 & -0.08 & -0.06 & -0.06 & \color{red}{+0.04} & -0.15 & -0.00 & +0.00 & -0.09 & -0.00 & +0.03 \\ \Lambda(1600) & -0.02 & -0.09 & +0.13 & +0.22 & \color{red}{+0.50} & -0.09 & -0.30 & +0.23 & \color{blue}{-1.93} & -0.46 & +0.12 & -1.85 & -0.12 & -0.06 & +0.00 & -0.03 & +0.05 \\ \Lambda(1670) & -0.01 & \color{red}{+0.06} & +0.03 & +0.01 & -0.01 & -0.12 & \color{blue}{-0.29} & -0.03 & -0.11 & +0.05 & -0.01 & -0.03 & +0.01 & +0.00 & -0.03 & -0.01 & +0.02 \\ \Lambda(1690) & +0.00 & -0.00 & +0.04 & -0.13 & +0.01 & -0.06 & -0.04 & -0.26 & \color{blue}{-0.33} & -0.08 & -0.04 & \color{red}{+0.06} & +0.01 & +0.00 & -0.10 & -0.01 & -0.08 \\ \Lambda(2000) & +0.05 & +0.10 & -0.08 & -0.09 & +0.08 & -0.85 & \color{blue}{-2.26} & \color{red}{+0.83} & -0.93 & +0.31 & -0.23 & -0.86 & +0.35 & +0.19 & -1.40 & -0.02 & +0.30 \\ \Delta(1232) & -0.27 & +0.02 & +0.31 & \color{red}{+1.76} & -0.44 & -0.14 & +0.49 & -0.63 & -0.77 & +0.53 & -0.31 & +0.65 & +0.10 & +0.06 & \color{blue}{-0.79} & -0.25 & +0.24 \\ \Delta(1600) & +0.33 & -0.10 & -0.15 & -0.28 & +0.59 & -0.38 & \color{blue}{-1.40} & -0.29 & \color{red}{+0.93} & +0.03 & +0.05 & -0.58 & +0.07 & +0.04 & -0.10 & +0.01 & -0.26 \\ \Delta(1700) & -0.01 & +0.03 & -0.13 & +0.07 & +0.39 & \color{blue}{-0.48} & -0.15 & +0.82 & \color{red}{+0.94} & +0.18 & +0.05 & +0.75 & +0.03 & +0.01 & +0.23 & +0.01 & +0.10 \\ K(700) & +0.17 & -0.02 & +0.04 & +0.75 & +0.62 & \color{blue}{-0.59} & -0.31 & -0.06 & \color{red}{+0.91} & +0.56 & +0.25 & +0.42 & +0.10 & +0.07 & -0.06 & +0.01 & +0.26 \\ K(892) & -0.53 & -0.02 & -0.12 & -0.46 & -0.58 & +0.55 & \color{red}{+0.59} & -0.06 & +0.18 & +0.28 & -0.16 & +0.25 & -0.01 & +0.00 & -0.52 & +0.02 & \color{blue}{-0.70} \\ K(1430) & -0.29 & +0.07 & -0.50 & -0.03 & +0.18 & -0.76 & \color{red}{+2.78} & +2.40 & \color{blue}{-2.67} & +1.29 & +0.23 & -2.27 & +0.91 & +0.44 & +2.07 & -0.00 & +0.71 \\ \end{array} \end{split}\]
  • 0: Default amplitude model

  • 1: Alternative amplitude model with K(892) with free mass and width

  • 2: Alternative amplitude model with L(1670) with free mass and width

  • 3: Alternative amplitude model with L(1690) with free mass and width

  • 4: Alternative amplitude model with D(1232) with free mass and width

  • 5: Alternative amplitude model with L(1600), D(1600), D(1700) with free mass and width

  • 6: Alternative amplitude model with free L(1405) Flatt’e widths, indicated as G1 (pK channel) and G2 (Sigmapi)

  • 7: Alternative amplitude model with L(1800) contribution added with free mass and width

  • 8: Alternative amplitude model with L(1810) contribution added with free mass and width

  • 9: Alternative amplitude model with D(1620) contribution added with free mass and width

  • 10: Alternative amplitude model in which a Relativistic Breit-Wigner is used for the K(700) contribution

  • 11: Alternative amplitude model with K(700) with free mass and width

  • 12: Alternative amplitude model with K(1410) contribution added with mass and width from PDG2020

  • 13: Alternative amplitude model in which a Relativistic Breit-Wigner is used for the K(1430) contribution

  • 14: Alternative amplitude model with K(1430) with free width

  • 15: Alternative amplitude model with an additional overall exponential form factor exp(-alpha q^2) multiplying Bugg lineshapes. The exponential parameter is indicated as ``alpha’’

  • 16: Alternative amplitude model with free radial parameter d for the Lc resonance, indicated as dLc

  • 17: Alternative amplitude model obtained using LS couplings

5.6. Average polarimetry values#

The components of the averaged polarimeter vector \(\overline{\alpha}\) are defined as:

\[ \overline{\alpha}_j = \int I_0\left(\tau\right) \alpha_j\left(\tau\right) \mathrm{d}^n \tau \; \big / \int I_0\left(\tau\right)\,\mathrm{d}^n \tau \]

The averages of the norm of \(\vec\alpha\) are computed as follows:

  • \(\left|\overline{\alpha}\right| = \sqrt{\overline{\alpha_x}^2+\overline{\alpha_y}^2+\overline{\alpha_z}^2}\), with the statistical uncertainties added in quadrature and the systematic uncertainties by taking the same formula on the extrema values of each \(\overline{\alpha_j}\)

  • \(\overline{\left|\alpha\right|} = \sqrt{\int I_0\left(\tau\right) \left|\vec\alpha\left(\tau\right)\right|^2 \mathrm{d}^n \tau \; \big / \int I_0\left(\tau\right)\,\mathrm{d}^n \tau}\)

Hide code cell source
def compute_weighted_average(v: jnp.ndarray, weights: jnp.ndarray) -> jnp.ndarray:
    return jnp.nansum(v * weights, axis=-1) / jnp.nansum(weights, axis=-1)


def compute_polar_coordinates(cartesian_vector: jnp.ndarray) -> jnp.ndarray:
    x, y, z = cartesian_vector
    norm = jnp.sqrt(jnp.sum(cartesian_vector**2, axis=0))
    theta = np.arccos(z / norm)
    phi = np.pi - np.arctan2(y, -x)
    return jnp.array([norm, theta, phi])


stat_weighted_alpha_norm = jnp.sqrt(
    compute_weighted_average(stat_polarimetry_norm**2, weights=stat_intensities)
)
stat_weighted_alpha = compute_weighted_average(
    stat_polarimetry,
    weights=stat_intensities,
)
stat_weighted_alpha_polar = compute_polar_coordinates(stat_weighted_alpha)
assert stat_weighted_alpha_norm.shape == (n_bootstraps,)
assert stat_weighted_alpha.shape == (3, n_bootstraps)
assert stat_weighted_alpha_polar.shape == (3, n_bootstraps)

syst_weighted_alpha_norm = jnp.sqrt(
    compute_weighted_average(syst_polarimetry_norm**2, weights=syst_intensities)
)
syst_weighted_alpha = compute_weighted_average(
    syst_polarimetry,
    weights=syst_intensities,
)
syst_weighted_alpha_polar = compute_polar_coordinates(syst_weighted_alpha)
assert syst_weighted_alpha_norm.shape == (n_models,)
assert syst_weighted_alpha.shape == (3, n_models)
assert syst_weighted_alpha_polar.shape == (3, n_models)

nominal_weighted_alpha_norm = syst_weighted_alpha_norm[0]
nominal_weighted_alpha = syst_weighted_alpha[:, 0]
nominal_weighted_alpha_polar = syst_weighted_alpha_polar[:, 0]

syst_weighted_alpha_norm_diff = (
    syst_weighted_alpha_norm - nominal_weighted_alpha_norm
)
syst_weighted_alpha_diff = (syst_weighted_alpha.T - nominal_weighted_alpha).T
syst_weighted_alpha_diff_polar = (
    syst_weighted_alpha_polar.T - nominal_weighted_alpha_polar
).T

stat_weighted_alpha_std = stat_weighted_alpha.std(axis=1)
syst_weighted_alpha_min = syst_weighted_alpha_diff.min(axis=1)
syst_weighted_alpha_max = syst_weighted_alpha_diff.max(axis=1)
stat_weighted_alpha_polar_std = stat_weighted_alpha_polar.std(axis=1)
syst_weighted_alpha_polar_min = syst_weighted_alpha_diff_polar.min(axis=1)
syst_weighted_alpha_polar_max = syst_weighted_alpha_diff_polar.max(axis=1)

Cartesian coordinates:

Hide code cell source
def render_cartesian_uncertainties(
    value, stat, syst_min, syst_max, plus: bool = True
) -> str:
    value *= 1e3
    stat *= 1e3
    syst_min *= 1e3
    syst_max *= 1e3
    if plus:
        val = f"{value:+.1f}"
    else:
        val = f"{value:.1f}"
    stat = f"{stat:.1f}"
    syst_min = f"-{abs(syst_min):.1f}"
    syst_max = f"+{abs(syst_max):.1f}"
    return (
        Rf"\left({val} \pm {stat}_{{{syst_min}}}^{{{syst_max}}} \right) \times"
        R" 10^{-3}"
    )


src = R"""
\begin{array}{ccr}
"""
for i in range(3):
    value = render_cartesian_uncertainties(
        nominal_weighted_alpha[i],
        stat_weighted_alpha_std[i],
        syst_weighted_alpha_min[i],
        syst_weighted_alpha_max[i],
    )
    src += Rf"  \overline{{\alpha_{'xyz'[i]}}} & = & {value} \\" "\n"

value = render_cartesian_uncertainties(
    nominal_weighted_alpha_norm,
    stat_weighted_alpha_norm.std(),
    syst_weighted_alpha_norm_diff.min(),
    syst_weighted_alpha_norm_diff.max(),
    plus=False,
)
src += Rf"  \overline{{\left|\alpha\right|}} & = & {value} \\"
src += R"\end{array}"
Latex(src)
\[\begin{split}\begin{array}{ccr} \overline{\alpha_x} & = & \left(-62.6 \pm 4.5_{-14.8}^{+8.4} \right) \times 10^{-3} \\ \overline{\alpha_y} & = & \left(+8.9 \pm 8.9_{-12.7}^{+9.1} \right) \times 10^{-3} \\ \overline{\alpha_z} & = & \left(-278.0 \pm 23.7_{-40.4}^{+12.6} \right) \times 10^{-3} \\ \overline{\left|\alpha\right|} & = & \left(669.4 \pm 9.3_{-10.4}^{+15.3} \right) \times 10^{-3} \\\end{array}\end{split}\]

Polar coordinates:

\[\begin{split} \begin{array}{lcl} \theta\left(\vec\alpha\right) &=& \arccos\left(\alpha_z \, \big/ \left|\alpha\right|\right) \\ \phi\left(\vec\alpha\right) &=& \pi - \mathrm{atan2}\left(\alpha_y, -\alpha_x\right) \end{array} \end{split}\]
Hide code cell source
def render_radian_angle_uncertainties(value, stat, syst_min, syst_max) -> str:
    val = f"{value:+.2f}"
    stat = f"{stat:.2f}"
    syst_min = f"-{abs(syst_min):.2f}"
    syst_max = f"+{abs(syst_max):.2f}"
    return Rf"{val} \pm {stat}_{{{syst_min}}}^{{{syst_max}}}\;\mathrm{{rad}}"


def render_angle_uncertainties(value, stat, syst_min, syst_max) -> str:
    value /= np.pi
    stat /= np.pi
    syst_min /= np.pi
    syst_max /= np.pi
    val = f"{value:+.3f}"
    stat = f"{stat:.3f}"
    syst_min = f"-{abs(syst_min):.3f}"
    syst_max = f"+{abs(syst_max):.3f}"
    return (
        Rf"\left({val} \pm {stat}_{{{syst_min}}}^{{{syst_max}}} \right) \times \pi"
    )


src = R"""
\begin{array}{ccl}
"""
labels = [
    R"\left|\overline{\alpha}\right|",
    R"\theta\left(\overline{\alpha}\right)",
    R"\phi\left(\overline{\alpha}\right)",
]
for i, label in enumerate(labels):
    renderer = (
        render_cartesian_uncertainties
        if i == 0
        else render_radian_angle_uncertainties
    )
    args = (
        nominal_weighted_alpha_polar[i],
        stat_weighted_alpha_polar_std[i],
        syst_weighted_alpha_polar_min[i],
        syst_weighted_alpha_polar_max[i],
    )
    src += Rf"  {label} & = & {renderer(*args)} \\" "\n"
    if i > 0:
        src += Rf"  & = & {render_angle_uncertainties(*args)} \\" "\n"

src += R"\end{array}"
Latex(src)
\[\begin{split}\begin{array}{ccl} \left|\overline{\alpha}\right| & = & \left(+285.1 \pm 24.0_{-13.8}^{+37.9} \right) \times 10^{-3} \\ \theta\left(\overline{\alpha}\right) & = & +2.92 \pm 0.01_{-0.04}^{+0.05}\;\mathrm{rad} \\ & = & \left(+0.929 \pm 0.002_{-0.011}^{+0.017} \right) \times \pi \\ \phi\left(\overline{\alpha}\right) & = & +3.00 \pm 0.14_{-0.09}^{+0.21}\;\mathrm{rad} \\ & = & \left(+0.955 \pm 0.045_{-0.028}^{+0.067} \right) \times \pi \\ \end{array}\end{split}\]

Averaged polarimeter values for each model (and the difference with the nominal model):

Hide code cell source
src = R"""
\begin{array}{r|rrrr|rrrr}
  \textbf{Model}
  & \overline{\alpha}_x & \overline{\alpha}_y & \overline{\alpha}_z & \overline{\left|\alpha\right|}
  & \Delta\overline{\alpha}_x & \Delta\overline{\alpha}_y & \Delta\overline{\alpha}_z
  & \Delta\overline{\left|\alpha\right|} \\
  \hline
"""
for i, title in enumerate(models):
    α = 1e3 * syst_weighted_alpha[:, i]
    abs_α = 1e3 * syst_weighted_alpha_norm[i]
    Δα = 1e3 * syst_weighted_alpha_diff[:, i]
    abs_Δα = 1e3 * syst_weighted_alpha_norm_diff[i]
    src += Rf"  \textbf{{{i}}}"
    src += Rf"  & {α[0]:+.1f} & {α[1]:+.1f} & {α[2]:+.1f} & {abs_α:.1f}"
    if i != 0:
        src += Rf"  & {Δα[0]:+.1f} & {Δα[1]:+.1f} & {Δα[2]:+.1f} & {abs_Δα:+.1f}"
    src += R"  \\" "\n"
    del i, title, α, abs_α, Δα, abs_Δα
src += R"\end{array}"
Latex(src)
\[\begin{split}\begin{array}{r|rrrr|rrrr} \textbf{Model} & \overline{\alpha}_x & \overline{\alpha}_y & \overline{\alpha}_z & \overline{\left|\alpha\right|} & \Delta\overline{\alpha}_x & \Delta\overline{\alpha}_y & \Delta\overline{\alpha}_z & \Delta\overline{\left|\alpha\right|} \\ \hline \textbf{0} & -62.6 & +8.9 & -278.0 & 669.4 \\ \textbf{1} & -61.6 & +8.5 & -279.4 & 670.7 & +1.0 & -0.4 & -1.4 & +1.3 \\ \textbf{2} & -62.9 & +9.1 & -278.4 & 669.8 & -0.3 & +0.2 & -0.5 & +0.4 \\ \textbf{3} & -58.4 & +7.4 & -276.2 & 667.7 & +4.2 & -1.5 & +1.8 & -1.6 \\ \textbf{4} & -69.3 & +9.5 & -277.2 & 666.9 & -6.6 & +0.6 & +0.8 & -2.5 \\ \textbf{5} & -70.7 & +9.6 & -277.4 & 668.7 & -8.0 & +0.8 & +0.6 & -0.6 \\ \textbf{6} & -69.7 & +9.1 & -281.7 & 673.0 & -7.1 & +0.2 & -3.8 & +3.7 \\ \textbf{7} & -77.4 & +18.0 & -305.4 & 671.4 & -14.8 & +9.1 & -27.5 & +2.1 \\ \textbf{8} & -55.8 & +10.9 & -284.6 & 675.5 & +6.8 & +2.0 & -6.7 & +6.1 \\ \textbf{9} & -66.9 & +4.4 & -290.4 & 672.8 & -4.3 & -4.5 & -12.4 & +3.5 \\ \textbf{10} & -56.4 & +2.4 & -265.4 & 659.0 & +6.2 & -6.5 & +12.6 & -10.4 \\ \textbf{11} & -64.7 & +9.3 & -278.6 & 670.4 & -2.1 & +0.4 & -0.6 & +1.0 \\ \textbf{12} & -75.1 & +1.8 & -283.4 & 663.5 & -12.5 & -7.1 & -5.4 & -5.8 \\ \textbf{13} & -61.8 & +8.1 & -277.3 & 668.8 & +0.9 & -0.8 & +0.7 & -0.6 \\ \textbf{14} & -62.2 & +8.7 & -277.6 & 669.2 & +0.5 & -0.2 & +0.4 & -0.2 \\ \textbf{15} & -54.2 & -3.8 & -318.4 & 684.6 & +8.4 & -12.7 & -40.4 & +15.3 \\ \textbf{16} & -62.1 & +8.2 & -278.1 & 669.5 & +0.5 & -0.7 & -0.1 & +0.2 \\ \textbf{17} & -58.1 & +12.1 & -278.6 & 666.5 & +4.5 & +3.2 & -0.6 & -2.9 \\ \end{array}\end{split}\]
Hide code cell source
src = R"""
\begin{array}{r|rrr|rrr}
  \textbf{Model}
  & 10^3 \cdot \left|\overline{\alpha}\right|
  & \theta\left(\overline{\alpha}\right) / \pi
  & \phi\left(\overline{\alpha}\right) / \pi
  & 10^3 \cdot \Delta\left|\overline{\alpha}\right|
  & \Delta\theta\left(\overline{\alpha}\right) / \pi
  & \Delta\phi\left(\overline{\alpha}\right) / \pi \\
  \hline
"""
for i, title in enumerate(models):
    α, θ, φ = syst_weighted_alpha_polar[:, i]
    Δα, Δθ, Δφ = syst_weighted_alpha_diff_polar[:, i]
    src += Rf"  \textbf{{{i}}}"
    src += Rf"  & {1e3*α:+.1f} & {θ/np.pi:+.3f} & {φ/np.pi:+.3f}"
    if i != 0:
        src += Rf"  & {1e3*Δα:+.1f} & {Δθ/np.pi:+.3f} & {Δφ/np.pi:+.3f}"
    src += R"  \\" "\n"
    del i, title, α, θ, φ, Δα, Δθ, Δφ
src += R"\end{array}"
Latex(src)
\[\begin{split}\begin{array}{r|rrr|rrr} \textbf{Model} & 10^3 \cdot \left|\overline{\alpha}\right| & \theta\left(\overline{\alpha}\right) / \pi & \phi\left(\overline{\alpha}\right) / \pi & 10^3 \cdot \Delta\left|\overline{\alpha}\right| & \Delta\theta\left(\overline{\alpha}\right) / \pi & \Delta\phi\left(\overline{\alpha}\right) / \pi \\ \hline \textbf{0} & +285.1 & +0.929 & +0.955 \\ \textbf{1} & +286.2 & +0.930 & +0.956 & +1.1 & +0.001 & +0.001 \\ \textbf{2} & +285.6 & +0.929 & +0.954 & +0.5 & -0.000 & -0.001 \\ \textbf{3} & +282.4 & +0.933 & +0.960 & -2.7 & +0.004 & +0.005 \\ \textbf{4} & +285.8 & +0.921 & +0.956 & +0.8 & -0.007 & +0.001 \\ \textbf{5} & +286.4 & +0.920 & +0.957 & +1.4 & -0.009 & +0.002 \\ \textbf{6} & +290.4 & +0.922 & +0.959 & +5.3 & -0.007 & +0.004 \\ \textbf{7} & +315.6 & +0.919 & +0.927 & +30.5 & -0.010 & -0.028 \\ \textbf{8} & +290.3 & +0.937 & +0.939 & +5.2 & +0.008 & -0.017 \\ \textbf{9} & +298.0 & +0.928 & +0.979 & +12.9 & -0.001 & +0.024 \\ \textbf{10} & +271.3 & +0.933 & +0.987 & -13.8 & +0.004 & +0.031 \\ \textbf{11} & +286.2 & +0.927 & +0.955 & +1.1 & -0.002 & -0.000 \\ \textbf{12} & +293.2 & +0.918 & +0.992 & +8.1 & -0.011 & +0.037 \\ \textbf{13} & +284.2 & +0.930 & +0.958 & -0.9 & +0.001 & +0.003 \\ \textbf{14} & +284.6 & +0.929 & +0.956 & -0.5 & +0.000 & +0.001 \\ \textbf{15} & +323.0 & +0.946 & +1.022 & +37.9 & +0.017 & +0.067 \\ \textbf{16} & +285.1 & +0.929 & +0.958 & -0.0 & +0.001 & +0.003 \\ \textbf{17} & +284.8 & +0.933 & +0.935 & -0.2 & +0.004 & -0.021 \\ \end{array}\end{split}\]

Tip

These values can be downloaded in serialized JSON format under Exported distributions.

Hide code cell source
alpha_x_reversed = syst_weighted_alpha[0, ::-1]
alpha_y_reversed = syst_weighted_alpha[1, ::-1]
alpha_z_reversed = syst_weighted_alpha[2, ::-1]
alpha_norm_reversed = syst_weighted_alpha_polar[0, ::-1]
alpha_theta_reversed = syst_weighted_alpha_polar[1, ::-1]
alpha_phi_reversed = syst_weighted_alpha_polar[2, ::-1]
marker_options = dict(
    color=(n_models - 1) * ["blue"] + ["red"],
    size=(n_models - 1) * [7] + [10],
    opacity=0.7,
)
model_titles_reversed = [
    "<br>".join(wrap(title, width=60)) for title in reversed(models)
]

fig = make_subplots(cols=2, shared_yaxes=True)
fig.add_trace(
    go.Scatter(
        x=alpha_phi_reversed,
        y=alpha_theta_reversed,
        hovertemplate="<b>%{text}</b><br>"
        + "ϕ(ɑ̅) = %{x:.3f}, "
        + "θ(ɑ̅) = %{y:.3f}"
        + "<extra></extra>",  # hide trace box
        mode="markers",
        marker=marker_options,
        text=model_titles_reversed,
    ),
    col=1,
    row=1,
)
fig.add_trace(
    go.Scatter(
        x=alpha_norm_reversed,
        y=alpha_theta_reversed,
        hovertemplate="<b>%{text}</b><br>"
        + "|ɑ̅| = %{x:.3f}, "
        + "θ(ɑ̅) = %{y:.3f}"
        + "<extra></extra>",  # hide trace box
        mode="markers",
        marker=marker_options,
        text=model_titles_reversed,
    ),
    col=2,
    row=1,
)
fig.update_layout(
    height=600,
    showlegend=False,
    title_text="Averaged ɑ̅ <b>systematics</b> distribution (polar)",
    xaxis=dict(
        title="ϕ(ɑ̅)",
        tickmode="array",
        tickvals=np.array([0.9, 0.95, 1, 1.05]) * np.pi,
        ticktext=["0.9 π", "0.95 π", "π", "1.05 π"],
    ),
    yaxis=dict(
        title="θ(ɑ̅)",
        tickmode="array",
        tickvals=np.array([0.90, 0.91, 0.92, 0.93, 0.94, 0.95]) * np.pi,
        ticktext=["0.90 π", "0.91 π", "0.92 π", "0.93 π", "0.94 π", "0.95 π"],
    ),
    xaxis2=dict(title="|ɑ̅|"),
)
fig.show()

fig = go.Figure(
    data=go.Scatter3d(
        x=alpha_x_reversed,
        y=alpha_y_reversed,
        z=alpha_z_reversed,
        hovertemplate="<b>%{text}</b><br>"
        + "ɑ̅<sub>x</sub> = %{x:.3f}, "
        + "ɑ̅<sub>y</sub> = %{y:.3f}, "
        + "ɑ̅<sub>z</sub> = %{z:.3f}"
        + "<extra></extra>",  # hide trace box
        mode="markers",
        marker=marker_options,
        text=model_titles_reversed,
    ),
)
fig.update_layout(
    width=700,
    height=700,
    scene=dict(
        aspectmode="cube",
        xaxis_title="ɑ̅<sub>x</sub>",
        yaxis_title="ɑ̅<sub>y</sub>",
        zaxis_title="ɑ̅<sub>z</sub>",
    ),
    showlegend=False,
    title_text=R"Average ɑ̅ <b>systematics</b> distribution (cartesian)",
)
fig.show()