diff --git a/docs/lambda-k-pi/ampform-dpd.ipynb b/docs/lambda-k-pi/ampform-dpd.ipynb
new file mode 100644
index 0000000..adbd634
--- /dev/null
+++ b/docs/lambda-k-pi/ampform-dpd.ipynb
@@ -0,0 +1,759 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Amplitude model with `ampform-dpd` (Draft)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "PWA study on $p \\gamma \\to \\Lambda K^+ \\pi^0$.\n",
+    "We formulate the helicity amplitude model symbolically using AmpForm-DPD here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Import Python libraries"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from __future__ import annotations\n",
+    "\n",
+    "import logging\n",
+    "import os\n",
+    "import warnings\n",
+    "from collections import defaultdict\n",
+    "from fractions import Fraction\n",
+    "from textwrap import dedent\n",
+    "\n",
+    "import ampform\n",
+    "import graphviz\n",
+    "import ipywidgets as w\n",
+    "import jax\n",
+    "import jax.numpy as jnp\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import qrules\n",
+    "import sympy as sp\n",
+    "from ampform.dynamics.builder import RelativisticBreitWignerBuilder\n",
+    "from ampform.io import aslatex, improve_latex_rendering\n",
+    "from IPython.display import Markdown, Math, display\n",
+    "from qrules.particle import Particle, Spin, create_particle, load_pdg\n",
+    "from tensorwaves.data import (\n",
+    "    SympyDataTransformer,\n",
+    "    TFPhaseSpaceGenerator,\n",
+    "    TFUniformRealNumberGenerator,\n",
+    ")\n",
+    "from tensorwaves.function.sympy import create_parametrized_function\n",
+    "\n",
+    "os.environ[\"TF_CPP_MIN_LOG_LEVEL\"] = \"3\"\n",
+    "logging.disable(logging.WARNING)\n",
+    "warnings.filterwarnings(\"ignore\")\n",
+    "\n",
+    "improve_latex_rendering()\n",
+    "particle_db = load_pdg()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Decay definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Particle definitions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "scroll-output",
+     "full-width",
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "def generate_markdown_table(particles: list[str]):\n",
+    "    src = dedent(r\"\"\"\n",
+    "    | Particle | Name | PID | $J^{PC} (I^G)$ | $I_3$ | $M$ | $\\Gamma$ | $Q$ | $S$ | $B$ |\n",
+    "    | :------- |------|-----|----------------|-------|-----|----------|-----|-----|-----|\n",
+    "    \"\"\")\n",
+    "    for name in particles:\n",
+    "        p = particle_db[name]\n",
+    "        src += f\"| ${p.latex}$ | `{p.name}` | {p.pid} | {jpc_ig(p)} | {i_3(p)} |  {p.mass:.3g}| {p.width:g} | {p.charge} |{p.strangeness} | {p.baryon_number}|\\n\"\n",
+    "    return src\n",
+    "\n",
+    "\n",
+    "def jpc_ig(particle: Particle) -> str:\n",
+    "    j = format_fraction(particle.spin)\n",
+    "    p = format_parity(particle.parity)\n",
+    "    c = format_parity(particle.c_parity)\n",
+    "    if particle.isospin is None:\n",
+    "        return f\"${j}^{{{p}{c}}}$\"\n",
+    "    i = format_fraction(particle.isospin.magnitude)\n",
+    "    g = format_parity(particle.g_parity)\n",
+    "    return rf\"${j}^{{{p}{c}}} \\; ({i}^{{{g}}})$\"\n",
+    "\n",
+    "\n",
+    "def i_3(particle: Particle) -> str:\n",
+    "    if particle.isospin is None:\n",
+    "        return \"N/A\"\n",
+    "    return f\"${format_fraction(particle.isospin.projection)}$\"\n",
+    "\n",
+    "\n",
+    "def format_fraction(value: float) -> str:\n",
+    "    value = Fraction(value)\n",
+    "    if value.denominator == 1:\n",
+    "        return str(value.numerator)\n",
+    "    return rf\"\\frac{{{value.numerator}}}{{{value.denominator}}}\"\n",
+    "\n",
+    "\n",
+    "def format_parity(parity: int | None) -> str:\n",
+    "    if parity is None:\n",
+    "        return \" \"\n",
+    "    if parity == -1:\n",
+    "        return \"-\"\n",
+    "    if parity == 1:\n",
+    "        return \"+\"\n",
+    "    raise NotImplementedError\n",
+    "\n",
+    "\n",
+    "particles = [\"Lambda\", \"K+\", \"pi0\", \"gamma\", \"p\"]\n",
+    "src = generate_markdown_table(particles)\n",
+    "Markdown(src)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the table above, PID is the PDG ID from PDG particle numbering scheme,\n",
+    "$J$ is the spin, $P$ is the parity, $C$ is the C parity, $I$ is the isospin (magnitude), $G$ is the G parity.\n",
+    "$I_3$ is the isospin projection (or the 3rd component),\n",
+    "$M$ is the mass,\n",
+    "$\\Gamma$ is the width,\n",
+    "$Q$ is the charge,\n",
+    "$S$ is the strangeness number,\n",
+    "and $B$ is the baryon number."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initial state definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Mass for $p \\gamma$ system"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "E_lab_gamma = 8.5\n",
+    "m_proton = 0.938\n",
+    "m_0 = np.sqrt(2 * E_lab_gamma * m_proton + m_proton**2)\n",
+    "m_eta = 0.548\n",
+    "m_pi = 0.135\n",
+    "m_0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Add custom particle $p \\gamma$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pgamma1 = Particle(\n",
+    "    name=\"pgamma1\",\n",
+    "    latex=r\"p\\gamma (s1/2)\",\n",
+    "    spin=0.5,\n",
+    "    mass=m_0,\n",
+    "    charge=1,\n",
+    "    isospin=Spin(1 / 2, +1 / 2),\n",
+    "    baryon_number=1,\n",
+    "    parity=-1,\n",
+    "    pid=99990,\n",
+    ")\n",
+    "pgamma2 = create_particle(\n",
+    "    template_particle=pgamma1,\n",
+    "    name=\"pgamma2\",\n",
+    "    latex=R\"p\\gamma (s3/2)\",\n",
+    "    spin=1.5,\n",
+    "    pid=pgamma1.pid + 1,\n",
+    ")\n",
+    "particle_db.update([pgamma1, pgamma2])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generate transitions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For simplicity, we use the initial state $p \\gamma$ (with spin-$\\frac{1}{2}$), \n",
+    "and set the allowed interaction type to be strong only,\n",
+    "the formalism is selected to be helicity formalism instead of canonical."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{seealso}\n",
+    "[Helicity versus canonical](https://ampform.readthedocs.io/stable/usage/helicity/formalism.html)\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "reaction = qrules.generate_transitions(\n",
+    "    initial_state=(\"pgamma1\"),\n",
+    "    final_state=[\"Lambda\", \"K+\", \"pi0\"],\n",
+    "    allowed_interaction_types=[\"strong\"],\n",
+    "    formalism=\"helicity\",\n",
+    "    particle_db=particle_db,\n",
+    "    max_angular_momentum=4,\n",
+    "    max_spin_magnitude=4,\n",
+    "    mass_conservation_factor=0,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "dot = qrules.io.asdot(reaction, collapse_graphs=True)\n",
+    "graphviz.Source(dot)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Formulate amplitude model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_builder = ampform.get_builder(reaction)\n",
+    "model_builder.config.scalar_initial_state_mass = True\n",
+    "model_builder.config.stable_final_state_ids = 0, 1, 2\n",
+    "bw_builder = RelativisticBreitWignerBuilder(\n",
+    "    energy_dependent_width=False,\n",
+    "    form_factor=False,\n",
+    ")\n",
+    "for name in reaction.get_intermediate_particles().names:\n",
+    "    model_builder.dynamics.assign(name, bw_builder)\n",
+    "model = model_builder.formulate()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.intensity"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The first term in the amplitude model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "scroll-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "(symbol, expr), *_ = model.amplitudes.items()\n",
+    "Math(aslatex({symbol: expr}, terms_per_line=1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "Model parameters"
+    },
+    "tags": [
+     "scroll-output",
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Math(aslatex(model.parameter_defaults))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "Kinematic variable definitions"
+    },
+    "tags": [
+     "hide-input",
+     "scroll-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Math(aslatex(model.kinematic_variables))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "unfolded_expression = model.expression.doit()\n",
+    "intensity_func = create_parametrized_function(\n",
+    "    expression=unfolded_expression,\n",
+    "    parameters=model.parameter_defaults,\n",
+    "    backend=\"jax\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "phsp_event = 500_000\n",
+    "rng = TFUniformRealNumberGenerator(seed=0)\n",
+    "phsp_generator = TFPhaseSpaceGenerator(\n",
+    "    initial_state_mass=reaction.initial_state[-1].mass,\n",
+    "    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},\n",
+    ")\n",
+    "phsp_momenta = phsp_generator.generate(phsp_event, rng)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "helicity_transformer = SympyDataTransformer.from_sympy(\n",
+    "    model.kinematic_variables,\n",
+    "    backend=\"jax\",\n",
+    ")\n",
+    "phsp = helicity_transformer(phsp_momenta)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-output",
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "resonances = defaultdict(set)\n",
+    "for transition in reaction.transitions:\n",
+    "    topology = transition.topology\n",
+    "    top_decay_products = topology.get_edge_ids_outgoing_from_node(0)\n",
+    "    (resonance_id, resonance), *_ = transition.intermediate_states.items()\n",
+    "    recoil_id, *_ = top_decay_products - {resonance_id}\n",
+    "    resonances[recoil_id].add(resonance.particle)\n",
+    "resonances = {k: sorted(v, key=lambda p: p.mass) for k, v in resonances.items()}\n",
+    "{k: [p.name for p in v] for k, v in resonances.items()}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Design slider UI"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "sliders = {}\n",
+    "categorized_sliders_m = defaultdict(list)\n",
+    "categorized_sliders_gamma = defaultdict(list)\n",
+    "categorized_cphi_pair = defaultdict(list)\n",
+    "\n",
+    "for symbol, value in model.parameter_defaults.items():\n",
+    "    if symbol.name.startswith(R\"\\Gamma_{\"):\n",
+    "        slider = w.FloatSlider(\n",
+    "            description=Rf\"\\({sp.latex(symbol)}\\)\",\n",
+    "            min=0.0,\n",
+    "            max=1.0,\n",
+    "            value=value,\n",
+    "            continuous_update=False,\n",
+    "        )\n",
+    "        sliders[symbol.name] = slider\n",
+    "        if symbol.name.startswith(R\"\\Gamma_{K\"):\n",
+    "            categorized_sliders_gamma[0].append(slider)\n",
+    "        elif symbol.name.startswith(R\"\\Gamma_{\\S\"):\n",
+    "            categorized_sliders_gamma[1].append(slider)\n",
+    "        elif symbol.name.startswith(R\"\\Gamma_{N\"):\n",
+    "            categorized_sliders_gamma[2].append(slider)\n",
+    "\n",
+    "    if symbol.name.startswith(\"m_{\"):\n",
+    "        slider = w.FloatSlider(\n",
+    "            description=Rf\"\\({sp.latex(symbol)}\\)\",\n",
+    "            min=0.63,\n",
+    "            max=4,\n",
+    "            value=value,\n",
+    "            continuous_update=False,\n",
+    "        )\n",
+    "        sliders[symbol.name] = slider\n",
+    "        if symbol.name.startswith(\"m_{K\"):\n",
+    "            categorized_sliders_m[0].append(slider)\n",
+    "        elif symbol.name.startswith(R\"m_{\\S\"):\n",
+    "            categorized_sliders_m[1].append(slider)\n",
+    "        elif symbol.name.startswith(\"m_{N\"):\n",
+    "            categorized_sliders_m[2].append(slider)\n",
+    "\n",
+    "    if symbol.name.startswith(\"C_{\"):\n",
+    "        c_latex = sp.latex(symbol)\n",
+    "        phi_latex = c_latex.replace(\"C\", R\"\\phi\", 1)\n",
+    "\n",
+    "        slider_c = w.FloatSlider(\n",
+    "            description=Rf\"\\({c_latex}\\)\",\n",
+    "            min=0,\n",
+    "            max=10,\n",
+    "            value=abs(value),\n",
+    "            continuous_update=False,\n",
+    "        )\n",
+    "        slider_phi = w.FloatSlider(\n",
+    "            description=Rf\"\\({phi_latex}\\)\",\n",
+    "            min=-np.pi,\n",
+    "            max=+np.pi,\n",
+    "            value=np.angle(value),\n",
+    "            continuous_update=False,\n",
+    "        )\n",
+    "\n",
+    "        sliders[symbol.name] = slider_c\n",
+    "        sliders[symbol.name.replace(\"C\", \"phi\", 1)] = slider_phi\n",
+    "\n",
+    "        cphi_hbox = w.HBox([slider_c, slider_phi])\n",
+    "        if \"Sigma\" in symbol.name:\n",
+    "            categorized_cphi_pair[1].append(cphi_hbox)\n",
+    "        elif R\"\\to N\" in symbol.name:\n",
+    "            categorized_cphi_pair[2].append(cphi_hbox)\n",
+    "        else:\n",
+    "            categorized_cphi_pair[0].append(cphi_hbox)\n",
+    "\n",
+    "\n",
+    "assert len(categorized_sliders_gamma) == 3\n",
+    "assert len(categorized_sliders_m) == 3\n",
+    "assert len(categorized_cphi_pair) == 3\n",
+    "\n",
+    "subtabs = {}\n",
+    "for category, resonance_list in resonances.items():\n",
+    "    subtabs[category] = []\n",
+    "    for particle in resonance_list:\n",
+    "        m_sliders = [\n",
+    "            slider\n",
+    "            for slider in categorized_sliders_m[category]\n",
+    "            if particle.latex in slider.description\n",
+    "        ]\n",
+    "        gamma_sliders = [\n",
+    "            slider\n",
+    "            for slider in categorized_sliders_gamma[category]\n",
+    "            if particle.latex in slider.description\n",
+    "        ]\n",
+    "        cphi_pairs = [\n",
+    "            hbox\n",
+    "            for hbox in categorized_cphi_pair[category]\n",
+    "            if particle.latex in hbox.children[0].description\n",
+    "        ]\n",
+    "        pole_pair = w.HBox(m_sliders + gamma_sliders)\n",
+    "        resonance_tab = w.VBox([pole_pair, *cphi_pairs])\n",
+    "        subtabs[category].append(resonance_tab)\n",
+    "assert len(subtabs) == 3\n",
+    "\n",
+    "main_tabs = []\n",
+    "for category, slider_boxes in subtabs.items():\n",
+    "    sub_tab_widget = w.Tab(children=slider_boxes)\n",
+    "    for i, particle in enumerate(resonances[category]):\n",
+    "        sub_tab_widget.set_title(i, particle.name)\n",
+    "\n",
+    "    main_tabs.append(sub_tab_widget)\n",
+    "\n",
+    "UI = w.Tab(children=main_tabs, titles=(\"K*\", \"Σ*\", \"N*\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def insert_phi(parameters: dict) -> dict:\n",
+    "    updated_parameters = {}\n",
+    "    for key, value in parameters.items():\n",
+    "        if key.startswith(\"phi_\"):\n",
+    "            continue\n",
+    "        if key.startswith(\"C_\"):\n",
+    "            phi_key = key.replace(\"C_\", \"phi_\")\n",
+    "            if phi_key in parameters:\n",
+    "                phi = parameters[phi_key]\n",
+    "                value *= np.exp(1j * phi)  # noqa:PLW2901\n",
+    "        updated_parameters[key] = value\n",
+    "    return updated_parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib widget\n",
+    "%config InlineBackend.figure_formats = ['png']\n",
+    "fig_2d, ax_2d = plt.subplots(dpi=200)\n",
+    "ax_2d.set_title(\"Model-weighted Phase space Dalitz Plot\")\n",
+    "ax_2d.set_xlabel(R\"$m^2(\\Lambda K^+)\\;\\left[\\mathrm{GeV}^2\\right]$\")\n",
+    "ax_2d.set_ylabel(R\"$m^2(K^+ \\pi^0)\\;\\left[\\mathrm{GeV}^2\\right]$\")\n",
+    "\n",
+    "mesh = None\n",
+    "\n",
+    "\n",
+    "def update_histogram(**parameters):\n",
+    "    global mesh\n",
+    "    parameters = insert_phi(parameters)\n",
+    "    intensity_func.update_parameters(parameters)\n",
+    "    intensity_weights = intensity_func(phsp)\n",
+    "    bin_values, xedges, yedges = jnp.histogram2d(\n",
+    "        phsp[\"m_01\"].real ** 2,\n",
+    "        phsp[\"m_12\"].real ** 2,\n",
+    "        bins=200,\n",
+    "        weights=intensity_weights,\n",
+    "        density=True,\n",
+    "    )\n",
+    "    bin_values = jnp.where(bin_values < 1e-6, jnp.nan, bin_values)\n",
+    "    x, y = jnp.meshgrid(xedges[:-1], yedges[:-1])\n",
+    "    if mesh is None:\n",
+    "        mesh = ax_2d.pcolormesh(x, y, bin_values.T, cmap=\"jet\", vmax=0.15)\n",
+    "    else:\n",
+    "        mesh.set_array(bin_values.T)\n",
+    "    fig_2d.canvas.draw_idle()\n",
+    "\n",
+    "\n",
+    "interactive_plot = w.interactive_output(update_histogram, sliders)\n",
+    "fig_2d.tight_layout()\n",
+    "fig_2d.colorbar(mesh, ax=ax_2d)\n",
+    "display(UI, interactive_plot)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "full-width",
+     "scroll-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib widget\n",
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "\n",
+    "fig, axes = plt.subplots(figsize=(11, 3.5), ncols=3, sharey=True)\n",
+    "fig.canvas.toolbar_visible = False\n",
+    "fig.canvas.header_visible = False\n",
+    "fig.canvas.footer_visible = False\n",
+    "ax1, ax2, ax3 = axes\n",
+    "\n",
+    "for recoil_id, ax in enumerate(axes):\n",
+    "    decay_products = sorted({0, 1, 2} - {recoil_id})\n",
+    "    product_latex = \" \".join([reaction.final_state[i].latex for i in decay_products])\n",
+    "    ax.set_xlabel(f\"$m({product_latex})$ [GeV]\")\n",
+    "\n",
+    "LINES = 3 * [None]\n",
+    "RESONANCE_LINE = [None] * len(reaction.get_intermediate_particles())\n",
+    "\n",
+    "\n",
+    "def update_plot(**parameters):\n",
+    "    parameters = insert_phi(parameters)\n",
+    "    intensity_func.update_parameters(parameters)\n",
+    "    intensities = intensity_func(phsp)\n",
+    "    max_value = 0\n",
+    "    color_id = 0\n",
+    "    for recoil_id, ax in enumerate(axes):\n",
+    "        decay_products = sorted({0, 1, 2} - {recoil_id})\n",
+    "        key = f\"m_{''.join(str(i) for i in decay_products)}\"\n",
+    "        bin_values, bin_edges = jax.numpy.histogram(\n",
+    "            phsp[key].real,\n",
+    "            bins=120,\n",
+    "            density=True,\n",
+    "            weights=intensities,\n",
+    "        )\n",
+    "        max_value = max(max_value, bin_values.max())\n",
+    "\n",
+    "        if LINES[recoil_id] is None:\n",
+    "            LINES[recoil_id] = ax.step(bin_edges[:-1], bin_values, alpha=0.5)[0]\n",
+    "        else:\n",
+    "            LINES[recoil_id].set_ydata(bin_values)\n",
+    "\n",
+    "        for resonance in resonances[recoil_id]:\n",
+    "            key = f\"m_{{{resonance.latex}}}\"\n",
+    "            val = parameters.get(key, resonance.mass)\n",
+    "            if RESONANCE_LINE[color_id] is None:\n",
+    "                RESONANCE_LINE[color_id] = ax.axvline(\n",
+    "                    val,\n",
+    "                    c=f\"C{color_id}\",\n",
+    "                    ls=\"dotted\",\n",
+    "                    label=resonance.name,\n",
+    "                )\n",
+    "            else:\n",
+    "                RESONANCE_LINE[color_id].set_xdata([val, val])\n",
+    "            color_id += 1\n",
+    "\n",
+    "    for ax in axes:\n",
+    "        ax.set_ylim(0, max_value * 1.1)\n",
+    "\n",
+    "\n",
+    "interactive_plot = w.interactive_output(update_plot, sliders)\n",
+    "for ax in axes:\n",
+    "    ax.legend(fontsize=\"small\")\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "display(UI, interactive_plot)"
+   ]
+  }
+ ],
+ "metadata": {
+  "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.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/lambda-k-pi/automated.ipynb b/docs/lambda-k-pi/ampform.ipynb
similarity index 100%
rename from docs/lambda-k-pi/automated.ipynb
rename to docs/lambda-k-pi/ampform.ipynb