diff --git a/docs/index.md b/docs/index.md
index 8b1fb26c..d46e898b 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -170,6 +170,7 @@ maxdepth: 1
 ---
 lc2pkpi
 jpsi2ksp
+xib2pkk
 ```
 
 ```{toctree}
diff --git a/docs/xib2pkk.ipynb b/docs/xib2pkk.ipynb
new file mode 100644
index 00000000..27e197ca
--- /dev/null
+++ b/docs/xib2pkk.ipynb
@@ -0,0 +1,598 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "# Ξ<sub>b</sub>⁻ → p K⁻ K⁻\n",
+    "\n",
+    "```{autolink-concat}\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-cell",
+     "scroll-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from __future__ import annotations\n",
+    "\n",
+    "import logging\n",
+    "import os\n",
+    "import warnings\n",
+    "from typing import TYPE_CHECKING\n",
+    "\n",
+    "import graphviz\n",
+    "import jax.numpy as jnp\n",
+    "import matplotlib.pyplot as plt\n",
+    "import qrules\n",
+    "import sympy as sp\n",
+    "from ampform.dynamics import BlattWeisskopfSquared, EnergyDependentWidth\n",
+    "from ampform.sympy import perform_cached_doit\n",
+    "from IPython.display import Latex, Markdown\n",
+    "from tensorwaves.data.transform import SympyDataTransformer\n",
+    "from tqdm.auto import tqdm\n",
+    "\n",
+    "from ampform_dpd import DalitzPlotDecompositionBuilder\n",
+    "from ampform_dpd.adapter.qrules import (\n",
+    "    load_particles,\n",
+    "    normalize_state_ids,\n",
+    "    permute_equal_final_states,\n",
+    "    to_three_body_decay,\n",
+    ")\n",
+    "from ampform_dpd.decay import State\n",
+    "from ampform_dpd.dynamics import FormFactor, RelativisticBreitWigner\n",
+    "from ampform_dpd.dynamics.builder import formulate_breit_wigner_with_form_factor\n",
+    "from ampform_dpd.io import (\n",
+    "    as_markdown_table,\n",
+    "    aslatex,\n",
+    "    perform_cached_lambdify,\n",
+    "    simplify_latex_rendering,\n",
+    ")\n",
+    "\n",
+    "simplify_latex_rendering()\n",
+    "logging.getLogger(\"absl\").setLevel(logging.ERROR)  # mute JAX\n",
+    "warnings.simplefilter(\"ignore\")\n",
+    "\n",
+    "NO_TQDM = \"EXECUTE_NB\" in os.environ\n",
+    "if NO_TQDM:\n",
+    "    logging.getLogger(\"ampform_dpd.io\").setLevel(logging.ERROR)\n",
+    "if TYPE_CHECKING:\n",
+    "    from tensorwaves.interface import DataSample, ParametrizedFunction"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
+   "source": [
+    "## Decay definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "See [DOI:10.1103/PhysRevD.104.052010](https://doi.org/10.1103/PhysRevD.104.052010) [[pdf](https://journals.aps.org/prd/pdf/10.1103/PhysRevD.104.052010)], _Search for CP violation in $\\Xi_b^- \\to p K^- K^-$ decays_ by LHCb. It found six asymmetry parameters, for $\\Lambda(1405)$, $\\Lambda(1520)$, $\\Lambda(1670)$, $\\Sigma(1385)$, $\\Sigma(1775)$, and $\\Sigma(1915)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "PARTICLES = load_particles()\n",
+    "excluded_particles = [\n",
+    "    \"Lambda(1600)\",\n",
+    "    \"Lambda(1690)\",\n",
+    "    \"Lambda(1800)\",\n",
+    "    \"Lambda(1810)\",\n",
+    "    \"Lambda(1890)\",\n",
+    "    \"Lambda(2000)\",\n",
+    "    \"Sigma(1660)0\",\n",
+    "    \"Sigma(1670)0\",\n",
+    "    \"Sigma(1750)0\",\n",
+    "    \"Sigma(1910)0\",\n",
+    "    \"Sigma(c)(2455)0\",\n",
+    "    \"Sigma(c)(2520)0\",\n",
+    "]\n",
+    "for name in excluded_particles:\n",
+    "    PARTICLES.remove(PARTICLES[name])\n",
+    "\n",
+    "REACTION = qrules.generate_transitions(\n",
+    "    initial_state=\"Xi(b)-\",\n",
+    "    final_state=[\"K-\", \"K-\", \"p\"],\n",
+    "    formalism=\"canonical-helicity\",\n",
+    "    allowed_intermediate_particles=[\"Lambda\", \"Sigma\"],\n",
+    "    max_angular_momentum=2,\n",
+    "    particle_db=PARTICLES,\n",
+    ")\n",
+    "REACTION = normalize_state_ids(REACTION)\n",
+    "REACTION = permute_equal_final_states(REACTION)\n",
+    "dot = qrules.io.asdot(REACTION, collapse_graphs=True)\n",
+    "graphviz.Source(dot)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "DECAY = to_three_body_decay(REACTION.transitions, min_ls=True)\n",
+    "Markdown(as_markdown_table([DECAY.initial_state, *DECAY.final_state.values()]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "resonances = sorted(\n",
+    "    {t.resonance for t in DECAY.chains},\n",
+    "    key=lambda p: (p.name[0], p.mass),\n",
+    ")\n",
+    "resonance_names = [p.name for p in resonances]\n",
+    "Markdown(as_markdown_table(resonances))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "def get_decay_identifier(decay):\n",
+    "    return (decay.resonance, *(particle.name for particle in decay.decay_products))\n",
+    "\n",
+    "\n",
+    "chains = {get_decay_identifier(c): c for c in DECAY.chains}\n",
+    "Latex(aslatex(chains.values(), with_jp=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true
+   },
+   "source": [
+    "## Model formulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "model_builder = DalitzPlotDecompositionBuilder(DECAY, min_ls=False)\n",
+    "for chain in model_builder.decay.chains:\n",
+    "    model_builder.dynamics_choices.register_builder(\n",
+    "        chain, formulate_breit_wigner_with_form_factor\n",
+    "    )\n",
+    "model = model_builder.formulate(reference_subsystem=1)\n",
+    "model.intensity"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Latex(aslatex(model.variables))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each **unaligned** amplitude is defined as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide-input",
+     "full-width",
+     "hide-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Latex(aslatex({k: v for k, v in model.amplitudes.items() if v}))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "s, m0, w0, m1, m2, L, R, z = sp.symbols(\"s m0 Gamma0 m1 m2 L R z\")\n",
+    "exprs = [\n",
+    "    RelativisticBreitWigner(s, m0, w0, m1, m2, L, R),\n",
+    "    EnergyDependentWidth(s, m0, w0, m1, m2, L, R),\n",
+    "    FormFactor(s, m1, m2, L, R),\n",
+    "    BlattWeisskopfSquared(z, L),\n",
+    "]\n",
+    "Latex(aslatex({e: e.doit(deep=False) for e in exprs}))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true
+   },
+   "source": [
+    "## Preparing for input data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "i, j = (1, 2)\n",
+    "k, *_ = {1, 2, 3} - {i, j}\n",
+    "σk, σk_expr = list(model.invariants.items())[k - 1]\n",
+    "Latex(aslatex({σk: σk_expr}))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Define meshgrid for Dalitz plot"
+    },
+    "tags": [
+     "hide-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "resolution = 1_000\n",
+    "m = sorted(model.masses, key=str)\n",
+    "x_min = float(((m[j] + m[k]) ** 2).xreplace(model.masses))\n",
+    "x_max = float(((m[0] - m[i]) ** 2).xreplace(model.masses))\n",
+    "y_min = float(((m[i] + m[k]) ** 2).xreplace(model.masses))\n",
+    "y_max = float(((m[0] - m[j]) ** 2).xreplace(model.masses))\n",
+    "x_diff = x_max - x_min\n",
+    "y_diff = y_max - y_min\n",
+    "x_min -= 0.05 * x_diff\n",
+    "x_max += 0.05 * x_diff\n",
+    "y_min -= 0.05 * y_diff\n",
+    "y_max += 0.05 * y_diff\n",
+    "X, Y = jnp.meshgrid(\n",
+    "    jnp.linspace(x_min, x_max, num=resolution),\n",
+    "    jnp.linspace(y_min, y_max, num=resolution),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Create data converter for Dalitz coordinates"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "definitions = dict(model.variables)\n",
+    "definitions[σk] = σk_expr\n",
+    "definitions = {\n",
+    "    symbol: expr.xreplace(definitions).xreplace(model.masses)\n",
+    "    for symbol, expr in definitions.items()\n",
+    "}\n",
+    "data_transformer = SympyDataTransformer.from_sympy(definitions, backend=\"jax\")\n",
+    "dalitz_data = {\n",
+    "    f\"sigma{i}\": X,\n",
+    "    f\"sigma{j}\": Y,\n",
+    "}\n",
+    "dalitz_data.update(data_transformer(dalitz_data))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "remove-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "for key, array in dalitz_data.items():\n",
+    "    assert not jnp.all(jnp.isnan(array)), f\"All values for {key} are NaN\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Prepare parametrized numerical function"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "unfolded_intensity_expr = perform_cached_doit(model.intensity)\n",
+    "full_intensity_expr = perform_cached_doit(\n",
+    "    unfolded_intensity_expr.xreplace(model.amplitudes)\n",
+    ")\n",
+    "free_parameters = {\n",
+    "    k: v\n",
+    "    for k, v in model.parameter_defaults.items()\n",
+    "    if isinstance(k, sp.Indexed)\n",
+    "    if \"production\" in str(k) or \"decay\" in str(k)\n",
+    "}\n",
+    "fixed_parameters = {\n",
+    "    k: v for k, v in model.parameter_defaults.items() if k not in free_parameters\n",
+    "}\n",
+    "intensity_func = perform_cached_lambdify(\n",
+    "    full_intensity_expr.xreplace(fixed_parameters),\n",
+    "    parameters=free_parameters,\n",
+    "    backend=\"jax\",\n",
+    ")\n",
+    "intensities = intensity_func(dalitz_data)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "remove-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "assert not jnp.all(jnp.isnan(intensities)), \"All intensities are NaN\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "## Dalitz plot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['png']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "def get_decay_products(subsystem_id: int) -> tuple[State, State]:\n",
+    "    return tuple(s for s in DECAY.final_state.values() if s.index != subsystem_id)\n",
+    "\n",
+    "\n",
+    "plt.rc(\"font\", size=18)\n",
+    "I_tot = jnp.nansum(intensities)\n",
+    "normalized_intensities = intensities / I_tot\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(14, 10))\n",
+    "mesh = ax.pcolormesh(X, Y, normalized_intensities)\n",
+    "ax.set_aspect(\"equal\")\n",
+    "c_bar = plt.colorbar(mesh, ax=ax, pad=0.01)\n",
+    "c_bar.ax.set_ylabel(\"Normalized intensity (a.u.)\")\n",
+    "sigma_labels = {\n",
+    "    i: Rf\"$\\sigma_{i} = M^2\\left({' '.join(p.latex for p in get_decay_products(i))}\\right)$\"\n",
+    "    for i in (1, 2, 3)\n",
+    "}\n",
+    "ax.set_xlabel(sigma_labels[i])\n",
+    "ax.set_ylabel(sigma_labels[j])\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "scroll-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "def compute_sub_intensity(\n",
+    "    func: ParametrizedFunction, phsp: DataSample, resonance_latex: str\n",
+    ") -> jnp.ndarray:\n",
+    "    original_parameters = dict(func.parameters)\n",
+    "    zero_parameters = {\n",
+    "        k: 0\n",
+    "        for k, v in func.parameters.items()\n",
+    "        if R\"\\mathcal{H}\" in k\n",
+    "        if resonance_latex not in k\n",
+    "    }\n",
+    "    func.update_parameters(zero_parameters)\n",
+    "    intensities = func(phsp)\n",
+    "    func.update_parameters(original_parameters)\n",
+    "    return intensities\n",
+    "\n",
+    "\n",
+    "plt.rc(\"font\", size=16)\n",
+    "fig, ax = plt.subplots(figsize=(10, 6), sharey=True)\n",
+    "fig.subplots_adjust(wspace=0.02)\n",
+    "x = jnp.sqrt(X[0])\n",
+    "y = jnp.sqrt(Y[:, 0])\n",
+    "ax.fill_between(x, jnp.nansum(normalized_intensities, axis=0), alpha=0.5)\n",
+    "_, y_max = ax.get_ylim()\n",
+    "ax.set_ylim(0, y_max)\n",
+    "ax.autoscale(enable=False, axis=\"x\")\n",
+    "ax.set_ylabel(\"Normalized intensity (a.u.)\")\n",
+    "ax.set_xlabel(sigma_labels[i])\n",
+    "resonance_counter = 0\n",
+    "for chain in tqdm(model.decay.chains, disable=NO_TQDM):\n",
+    "    if {p.index for p in chain.decay_products} != {1, 3}:\n",
+    "        continue\n",
+    "    resonance = chain.resonance\n",
+    "    sub_intensities = compute_sub_intensity(\n",
+    "        intensity_func, dalitz_data, resonance.latex\n",
+    "    )\n",
+    "    color = f\"C{resonance_counter}\"\n",
+    "    ax.plot(x, jnp.nansum(sub_intensities / I_tot, axis=0), c=color)\n",
+    "    ax.axvline(resonance.mass, label=f\"${resonance.latex}$\", c=color, ls=\"dashed\")\n",
+    "    resonance_counter += 1\n",
+    "ax.legend(fontsize=12)\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "colab": {
+   "toc_visible": true
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.19"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}