docs/report/010.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{autolink-concat}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "physics",
     "sympy"
    ]
   },
   "source": [
    "::::{margin}\n",
    ":::{card} P-vector\n",
    "TR-010\n",
    "^^^\n",
    "This report is a sequel to [TR-005](005.ipynb) and [TR-009](009.ipynb).\n",
    "+++\n",
    "✅&nbsp;[ampform#131](https://github.com/ComPWA/ampform/issues/131)\n",
    ":::\n",
    "::::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# P-vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "%pip install -q ampform==0.10.5 sympy==1.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "import os\n",
    "import re\n",
    "import warnings\n",
    "from typing import TYPE_CHECKING\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import mpl_interactions.ipyplot as iplt\n",
    "import numpy as np\n",
    "import symplot\n",
    "import sympy as sp\n",
    "from ampform.dynamics import (\n",
    "    BlattWeisskopfSquared,\n",
    "    breakup_momentum_squared,\n",
    "    coupled_width,\n",
    "    phase_space_factor_complex,\n",
    ")\n",
    "from ampform.dynamics.decorator import (\n",
    "    UnevaluatedExpression,\n",
    "    create_expression,\n",
    "    implement_doit_method,\n",
    ")\n",
    "from IPython.display import display\n",
    "from ipywidgets import widgets as ipywidgets\n",
    "from matplotlib import cm\n",
    "from mpl_interactions.controller import Controls\n",
    "\n",
    "if TYPE_CHECKING:\n",
    "    from sympy.printing.latex import LatexPrinter\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "STATIC_WEB_PAGE = {\"EXECUTE_NB\", \"READTHEDOCS\"}.intersection(os.environ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "editable": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "if \"COLAB_RELEASE_TAG\" in os.environ:\n",
    "    import subprocess\n",
    "\n",
    "    from google.colab import output\n",
    "\n",
    "    output.enable_custom_widget_manager()\n",
    "    subprocess.run(\"pip install -q ipympl\".split(), check=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "scroll-input"
    ]
   },
   "source": [
    "As described in [TR-005](005.ipynb), the $\\boldsymbol{K}$-matrix describes **scattering processes** of the type $cd \\to ab$. The $P$-vector approach is one of two generalizations for **production processes** of the type $c \\to ab$. For more details on this approach, {cite}`chungPartialWaveAnalysis1995` refers to {cite}`aitchisonMatrixFormalismOverlapping1972`.\n",
    "\n",
    "If we take the production vector $P$ to be:\n",
    "\n",
    "$$\n",
    "P_i = \\sum_R \\frac{\\beta^0_R\\,g_{R,i}(m)}{m_R^2-m^2}\n",
    "$$ (P-vector)\n",
    "\n",
    "and, in its invariant form,\n",
    "\n",
    "$$\n",
    "\\hat{P}_i = \\sum_R \\frac{\\beta^0_R\\,g_{R,i}(m)}{\\left(m_R^2-m^2\\right)\\sqrt{\\rho_i}}\n",
    "$$ (invariant P-vector)\n",
    "\n",
    "with $g_{R,i}(m)$ given by Eq. {eq}`residue-function` (possibly with [`coupled_width()`](https://ampform.readthedocs.io/en/0.10.5/api/ampform.dynamics.html#ampform.dynamics.coupled_width)), then the vector $F$ describes the resulting amplitudes by\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "F & = & \\left(\\boldsymbol{I}-i\\boldsymbol{K}\\right)^{-1}P \\\\\n",
    "\\hat{F} & = & \\left(\\boldsymbol{I}-i\\boldsymbol{\\hat{K}\\boldsymbol{\\rho}}\\right)^{-1}\\hat{P}\n",
    "\\end{eqnarray}\n",
    "$$ (F-vector)\n",
    "\n",
    "with, from Eq. {eq}`covariant-K-matrix`:\n",
    "\n",
    "$$\n",
    "\\hat{\\boldsymbol{K}} = \\sqrt{\\left(\\boldsymbol{\\rho}^\\dagger\\right)^{-1}} \\boldsymbol{K} \\sqrt{\\boldsymbol{\\rho}^{-1}}\n",
    "$$ (K-hat in terms of K)\n",
    "\n",
    "\n",
    "Just like with the residue functions in [TR-005](005.ipynb) and [TR-009](009.ipynb), $\\beta$ is often expressed in terms of resonance mass and 'width':\n",
    "\n",
    "$$\n",
    "\\beta^0_R = \\beta_R\\sqrt{m_R\\Gamma^0_R}\n",
    "$$ (beta functions)\n",
    "\n",
    "When in addition, we use a [`coupled_width()`](https://ampform.readthedocs.io/en/0.10.5/api/ampform.dynamics.html#ampform.dynamics.coupled_width), the $\\hat{P}$-vector becomes:\n",
    "\n",
    "$$\n",
    "\\hat{P}_i = \\sum_R \\frac{\\beta_R\\gamma_{R,i}m_R\\Gamma^0_R B_{R,i}(m)}{m_R^2-m^2}\n",
    "$$\n",
    "\n",
    "with $B_{R,i}(m)$ the ratio of Blatt-Weisskopf barrier factors ({class}`~ampform.dynamics.BlattWeisskopfSquared`) for channel $i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "@implement_doit_method()\n",
    "class PhaseSpaceFactor(UnevaluatedExpression):\n",
    "    is_commutative = True\n",
    "\n",
    "    def __new__(\n",
    "        cls,\n",
    "        s: sp.Symbol,\n",
    "        m_a: sp.Symbol,\n",
    "        m_b: sp.Symbol,\n",
    "        i: int,\n",
    "        **hints,\n",
    "    ) -> PhaseSpaceFactor:\n",
    "        return create_expression(cls, s, m_a, m_b, i, **hints)\n",
    "\n",
    "    def evaluate(self) -> sp.Expr:\n",
    "        s, m_a, m_b, *_ = self.args\n",
    "        return phase_space_factor_complex(s, m_a, m_b)\n",
    "\n",
    "    def _latex(self, printer: LatexPrinter, *args) -> str:\n",
    "        s = printer._print(self.args[0])\n",
    "        i = self.args[-1]\n",
    "        return Rf\"\\rho_{{{i}}}({s})\"\n",
    "\n",
    "\n",
    "@implement_doit_method()\n",
    "class CoupledWidth(UnevaluatedExpression):\n",
    "    is_commutative = True\n",
    "\n",
    "    def __new__(\n",
    "        cls,\n",
    "        s: sp.Symbol,\n",
    "        mass0: sp.IndexedBase,\n",
    "        gamma0: sp.IndexedBase,\n",
    "        m_a: sp.IndexedBase,\n",
    "        m_b: sp.IndexedBase,\n",
    "        angular_momentum: int,\n",
    "        R: int | sp.Symbol,\n",
    "        i: int,\n",
    "        **hints,\n",
    "    ) -> CoupledWidth:\n",
    "        return create_expression(\n",
    "            cls, s, mass0, gamma0, m_a, m_b, angular_momentum, R, i, **hints\n",
    "        )\n",
    "\n",
    "    def evaluate(self) -> sp.Expr:\n",
    "        s, mass0, gamma0, m_a, m_b, angular_momentum, R, i = self.args\n",
    "\n",
    "        def phsp_factor(s, m_a, m_b):\n",
    "            return PhaseSpaceFactor(s, m_a, m_b, i)\n",
    "\n",
    "        return coupled_width(\n",
    "            s,\n",
    "            mass0[R],\n",
    "            gamma0[R, i],\n",
    "            m_a[i],\n",
    "            m_b[i],\n",
    "            angular_momentum=angular_momentum,\n",
    "            meson_radius=1,\n",
    "            phsp_factor=phsp_factor,\n",
    "        )\n",
    "\n",
    "    def _latex(self, printer: LatexPrinter, *args) -> str:\n",
    "        s = printer._print(self.args[0])\n",
    "        R = self.args[-2]\n",
    "        i = self.args[-1]\n",
    "        return Rf\"{{\\Gamma_{{{R},{i}}}}}({s})\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "scroll-input"
    ]
   },
   "outputs": [],
   "source": [
    "def Pi_relativistic(\n",
    "    m: sp.Symbol,\n",
    "    M: sp.IndexedBase,\n",
    "    Gamma: sp.IndexedBase,\n",
    "    gamma: sp.IndexedBase,\n",
    "    beta: sp.IndexedBase,\n",
    "    m_a: sp.IndexedBase,\n",
    "    m_b: sp.IndexedBase,\n",
    "    R: int | sp.Symbol,\n",
    "    i: int,\n",
    "    n_resonances: int | sp.Symbol,\n",
    "    angular_momentum: int | sp.Symbol = 0,\n",
    ") -> sp.Expr:\n",
    "    q_squared = breakup_momentum_squared(m**2, m_a[i], m_b[i])\n",
    "    ff2 = BlattWeisskopfSquared(z=q_squared, angular_momentum=angular_momentum)\n",
    "    parametrization = (beta[R] * gamma[R, i] * M[R] * Gamma[R, i] * sp.sqrt(ff2)) / (\n",
    "        M[R] ** 2 - m**2\n",
    "    )\n",
    "    return sp.Sum(parametrization, (R, 0, n_resonances - 1))\n",
    "\n",
    "\n",
    "def Kij_relativistic(\n",
    "    m: sp.Symbol,\n",
    "    M: sp.IndexedBase,\n",
    "    Gamma: sp.IndexedBase,\n",
    "    gamma: sp.IndexedBase,\n",
    "    m_a: sp.IndexedBase,\n",
    "    m_b: sp.IndexedBase,\n",
    "    R: sp.IndexedBase,\n",
    "    i: int,\n",
    "    j: int,\n",
    "    n_resonances: int | sp.Symbol,\n",
    "    angular_momentum: int | sp.Symbol = 0,\n",
    ") -> sp.Expr:\n",
    "    def residue_function(i):\n",
    "        return gamma[R, i] * sp.sqrt(\n",
    "            M[R] * CoupledWidth(m**2, M, Gamma, m_a, m_b, angular_momentum, R, i)\n",
    "        )\n",
    "\n",
    "    g_i = residue_function(i)\n",
    "    g_j = residue_function(j)\n",
    "    parametrization = (g_i * g_j) / (M[R] ** 2 - m**2)\n",
    "    return sp.Sum(parametrization, (R, 0, n_resonances - 1))\n",
    "\n",
    "\n",
    "def f_vector(\n",
    "    n_resonances: int,\n",
    "    n_channels: int,\n",
    "    angular_momentum: int | sp.Symbol = 0,\n",
    ") -> sp.Matrix:\n",
    "    # Define symbols\n",
    "    R = sp.Symbol(\"R\")\n",
    "    m = sp.Symbol(\"m\", real=True)\n",
    "    M = sp.IndexedBase(\"m\", shape=(n_resonances,))\n",
    "    Gamma = sp.IndexedBase(\"Gamma\", shape=(n_resonances, n_channels))\n",
    "    gamma = sp.IndexedBase(\"gamma\", shape=(n_resonances, n_channels))\n",
    "    beta = sp.IndexedBase(\"beta\", shape=(n_resonances,))\n",
    "    m_a = sp.IndexedBase(\"m_a\", shape=(n_channels,))\n",
    "    m_b = sp.IndexedBase(\"m_b\", shape=(n_channels,))\n",
    "    # Define phase space matrix\n",
    "    rho = sp.zeros(n_channels, n_channels)\n",
    "    for i in range(n_channels):\n",
    "        rho[i, i] = PhaseSpaceFactor(m**2, m_a[i], m_b[i], i)\n",
    "    # Define P-vector, K-matrix and T-matrix\n",
    "    P = create_symbol_matrix(\"P\", n_channels, 1)\n",
    "    K = create_symbol_matrix(\"K\", n_channels, n_channels)\n",
    "    F = (sp.eye(n_channels) - sp.I * K * rho).inv() * P\n",
    "    # Substitute elements\n",
    "    return F.subs({\n",
    "        K[i, j]: Kij_relativistic(\n",
    "            m=m,\n",
    "            M=M,\n",
    "            Gamma=Gamma,\n",
    "            gamma=gamma,\n",
    "            m_a=m_a,\n",
    "            m_b=m_b,\n",
    "            i=i,\n",
    "            j=j,\n",
    "            R=R,\n",
    "            n_resonances=n_resonances,\n",
    "            angular_momentum=angular_momentum,\n",
    "        )\n",
    "        for i in range(n_channels)\n",
    "        for j in range(n_channels)\n",
    "    }).subs({\n",
    "        P[i]: Pi_relativistic(\n",
    "            m=m,\n",
    "            M=M,\n",
    "            Gamma=Gamma,\n",
    "            gamma=gamma,\n",
    "            beta=beta,\n",
    "            i=i,\n",
    "            m_a=m_a,\n",
    "            m_b=m_b,\n",
    "            R=R,\n",
    "            n_resonances=n_resonances,\n",
    "            angular_momentum=angular_momentum,\n",
    "        )\n",
    "        for i in range(n_channels)\n",
    "    })\n",
    "\n",
    "\n",
    "def create_symbol_matrix(name: str, m: int, n: int) -> sp.Matrix:\n",
    "    symbol = sp.IndexedBase(name, shape=(m, n))\n",
    "    return sp.Matrix([[symbol[i, j] for j in range(n)] for i in range(m)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sqrt{B_{L}^2\\left(\\frac{\\left(m^{2} - \\left({m_{a}}_{0} - {m_{b}}_{0}\\right)^{2}\\right) \\left(m^{2} - \\left({m_{a}}_{0} + {m_{b}}_{0}\\right)^{2}\\right)}{4 m^{2}}\\right)} {\\Gamma}_{0,0} {\\beta}_{0} {\\gamma}_{0,0} {m}_{0}}{\\left(1 - \\frac{i {\\Gamma_{0,0}}(m^{2}) {\\gamma}_{0,0}^{2} {m}_{0} \\rho_{0}(m^{2})}{- m^{2} + {m}_{0}^{2}}\\right) \\left(- m^{2} + {m}_{0}^{2}\\right)}$"
      ],
      "text/plain": [
       "sqrt(BlattWeisskopfSquared(L, (m**2 - (m_a[0] - m_b[0])**2)*(m**2 - (m_a[0] + m_b[0])**2)/(4*m**2)))*Gamma[0, 0]*beta[0]*gamma[0, 0]*m[0]/((1 - I*CoupledWidth(m**2, m, Gamma, m_a, m_b, L, 0, 0)*gamma[0, 0]**2*m[0]*PhaseSpaceFactor(m**2, m_a[0], m_b[0], 0)/(-m**2 + m[0]**2))*(-m**2 + m[0]**2))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = sp.Symbol(\"L\", integer=True)\n",
    "symplot.partial_doit(\n",
    "    f_vector(n_resonances=1, n_channels=1, angular_momentum=L)[0, 0], sp.Sum\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "\n",
    "Compare {func}`~ampform.dynamics.relativistic_breit_wigner_with_ff`\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "scroll-input"
    ]
   },
   "outputs": [],
   "source": [
    "def plot_f_vector(\n",
    "    n_channels: int,\n",
    "    n_resonances: int,\n",
    "    angular_momentum: int | sp.Symbol = 0,\n",
    "    title: str = \"\",\n",
    ") -> None:\n",
    "    # Convert to Symbol: symplot cannot handle\n",
    "    m = sp.Symbol(\"m\", real=True)\n",
    "    epsilon = sp.Symbol(\"epsilon\", real=True)\n",
    "    i = sp.Symbol(\"i\", integer=True, negative=False)\n",
    "    expr = f_vector(n_resonances, n_channels, angular_momentum=angular_momentum).doit()[\n",
    "        i, 0\n",
    "    ]\n",
    "    expr = symplot.substitute_indexed_symbols(expr)\n",
    "    expr = expr.subs(m, m + epsilon * sp.I)\n",
    "    np_expr, sliders = symplot.prepare_sliders(expr, m)\n",
    "    symbol_to_arg = {symbol: arg for arg, symbol in sliders._arg_to_symbol.items()}\n",
    "\n",
    "    # Set plot domain\n",
    "    x_min, x_max = 1e-3, 3\n",
    "    y_min, y_max = -0.5, +0.5\n",
    "\n",
    "    plot_domain = np.linspace(x_min, x_max, num=500)\n",
    "    x_values = np.linspace(x_min, x_max, num=160)\n",
    "    y_values = np.linspace(y_min, y_max, num=80)\n",
    "    X, Y = np.meshgrid(x_values, y_values)\n",
    "    plot_domain_complex = X + Y * 1j\n",
    "\n",
    "    # Set slider values and ranges\n",
    "    m0_values = np.linspace(x_min, x_max, num=n_resonances + 2)\n",
    "    m0_values = m0_values[1:-1]\n",
    "\n",
    "    def set_default_values():\n",
    "        if \"L\" in sliders:\n",
    "            sliders.set_ranges(L=(0, 8))\n",
    "        sliders.set_ranges(\n",
    "            i=(0, n_channels - 1),\n",
    "            epsilon=(y_min * 0.2, y_max * 0.2, 0.01),\n",
    "        )\n",
    "        for R in range(n_resonances):\n",
    "            for i in range(n_channels):\n",
    "                sliders.set_ranges({\n",
    "                    f\"m{R}\": (0, 3, 100),\n",
    "                    f\"beta{R}\": (0, 5, 0.1),\n",
    "                    Rf\"\\Gamma_{{{R},{i}}}\": (-5, +5, 100),\n",
    "                    Rf\"\\gamma_{{{R},{i}}}\": (0, 20, 100),\n",
    "                    f\"m_a{i}\": (0, 1, 0.01),\n",
    "                    f\"m_b{i}\": (0, 1, 0.01),\n",
    "                })\n",
    "                sliders.set_values({\n",
    "                    f\"m{R}\": m0_values[R],\n",
    "                    f\"beta{R}\": 1,\n",
    "                    Rf\"\\Gamma_{{{R},{i}}}\": 3 * (0.4 + R * 0.2 - i * 0.3),\n",
    "                    Rf\"\\gamma_{{{R},{i}}}\": 0.2 * (10 - R + i),\n",
    "                    f\"m_a{i}\": (i + 1) * 0.25,\n",
    "                    f\"m_b{i}\": (i + 1) * 0.25,\n",
    "                })\n",
    "\n",
    "    set_default_values()\n",
    "\n",
    "    # Create interactive plots\n",
    "    controls = Controls(**sliders)\n",
    "    nrows = 2  # set to 3 for imag+real\n",
    "    fig, axes = plt.subplots(\n",
    "        nrows=nrows,\n",
    "        figsize=(8, nrows * 3.0),\n",
    "        sharex=True,\n",
    "        tight_layout=True,\n",
    "    )\n",
    "    fig.canvas.toolbar_visible = False\n",
    "    fig.canvas.header_visible = False\n",
    "    fig.canvas.footer_visible = False\n",
    "    for ax in axes:\n",
    "        ax.set_xlim(x_min, x_max)\n",
    "    if not title:\n",
    "        title = f\"{n_channels}-channel $F$-vector with {n_resonances} resonances\"\n",
    "    fig.suptitle(title)\n",
    "\n",
    "    # 2D plot\n",
    "    axes[0].set_ylabel(\"$|T|^{2}$\")\n",
    "    axes[0].set_yticks([])\n",
    "\n",
    "    def plot(channel: int):\n",
    "        def wrapped(*args, **kwargs) -> sp.Expr:\n",
    "            kwargs[\"i\"] = channel\n",
    "            return np.abs(np_expr(*args, **kwargs)) ** 2\n",
    "\n",
    "        return wrapped\n",
    "\n",
    "    for i in range(n_channels):\n",
    "        iplt.plot(\n",
    "            plot_domain,\n",
    "            plot(i),\n",
    "            ax=axes[0],\n",
    "            controls=controls,\n",
    "            ylim=\"auto\",\n",
    "            label=f\"channel {i}\",\n",
    "        )\n",
    "    if n_channels > 1:\n",
    "        axes[0].legend(loc=\"upper right\")\n",
    "    mass_line_style = {\n",
    "        \"c\": \"red\",\n",
    "        \"alpha\": 0.3,\n",
    "    }\n",
    "    for name in controls.params:\n",
    "        if not re.match(r\"^m[0-9]+$\", name):\n",
    "            continue\n",
    "        iplt.axvline(controls[name], ax=axes[0], **mass_line_style)\n",
    "\n",
    "    # 3D plot\n",
    "    def plot3(**kwargs):\n",
    "        z_cutoff = kwargs.pop(\"z_cutoff\")\n",
    "        epsilon = kwargs[\"epsilon\"]\n",
    "        kwargs[\"epsilon\"] = 0\n",
    "        imag_real = kwargs.pop(\"imag_real\")\n",
    "        Z = np_expr(plot_domain_complex, **kwargs)\n",
    "        if imag_real == \"imag\":\n",
    "            Z_values = Z.imag\n",
    "            ax_title = \"Re $T$\"\n",
    "        elif imag_real == \"real\":\n",
    "            Z_values = Z.real\n",
    "            ax_title = \"Im $T$\"\n",
    "        elif imag_real == \"abs\":\n",
    "            Z_values = np.abs(Z)\n",
    "            ax_title = \"$|T|$\"\n",
    "        else:\n",
    "            raise NotImplementedError\n",
    "        for ax in axes[1:]:\n",
    "            ax.clear()\n",
    "        axes[-1].pcolormesh(\n",
    "            X, Y, Z_values, cmap=cm.coolwarm, vmin=-z_cutoff, vmax=+z_cutoff\n",
    "        )\n",
    "        i = kwargs[\"i\"]\n",
    "        if n_channels == 1:\n",
    "            axes[-1].set_title(ax_title)\n",
    "        else:\n",
    "            axes[-1].set_title(f\"{ax_title}, channel {i}\")\n",
    "        for ax in axes[1:]:\n",
    "            ax.axhline(0, linewidth=0.5, c=\"black\", linestyle=\"dotted\")\n",
    "            if epsilon != 0.0:\n",
    "                ax.axhline(\n",
    "                    epsilon,\n",
    "                    linewidth=0.5,\n",
    "                    c=\"blue\",\n",
    "                    linestyle=\"dotted\",\n",
    "                    label=R\"$\\epsilon$\",\n",
    "                )\n",
    "                axes[-1].text(\n",
    "                    x=x_min + 0.008,\n",
    "                    y=epsilon + 0.01,\n",
    "                    s=R\"$\\epsilon$\",\n",
    "                    c=\"blue\",\n",
    "                )\n",
    "            for R in range(n_resonances):\n",
    "                mass = kwargs[f\"m{R}\"]\n",
    "                ax.axvline(mass, **mass_line_style)\n",
    "            if \"m_a0\" in kwargs:\n",
    "                colors = cm.plasma(np.linspace(0, 1, n_channels))\n",
    "                for i, color in enumerate(colors):\n",
    "                    m_a = kwargs[f\"m_a{i}\"]\n",
    "                    m_b = kwargs[f\"m_b{i}\"]\n",
    "                    s_thr = m_a + m_b\n",
    "                    ax.axvline(s_thr, c=color, linestyle=\"dotted\")\n",
    "                    ax.text(\n",
    "                        x=s_thr,\n",
    "                        y=0.95 * y_min,\n",
    "                        s=f\"$m_{{a{i}}}+m_{{b{i}}}$\",\n",
    "                        c=color,\n",
    "                        rotation=-90,\n",
    "                    )\n",
    "                    m_diff = m_a - m_b\n",
    "                    x_offset = (x_max - x_min) * 0.015\n",
    "                    if m_diff > x_offset + 0.01 and s_thr - abs(m_diff) > x_offset:\n",
    "                        ax.axvline(\n",
    "                            m_diff,\n",
    "                            c=color,\n",
    "                            linestyle=\"dashed\",\n",
    "                            alpha=0.5,\n",
    "                        )\n",
    "                        ax.text(\n",
    "                            x=m_diff - x_offset,\n",
    "                            y=0.95 * y_min,\n",
    "                            s=f\"$m_{{a{i}}}-m_{{b{i}}}$\",\n",
    "                            c=color,\n",
    "                            rotation=+90,\n",
    "                        )\n",
    "            ax.set_ylabel(\"Im $m$\")\n",
    "            ax.set_xticks([])\n",
    "            ax.set_yticks([])\n",
    "            ax.set_facecolor(\"white\")\n",
    "        for R in range(n_resonances):\n",
    "            mass = kwargs[f\"m{R}\"]\n",
    "            axes[-1].text(\n",
    "                x=mass + (x_max - x_min) * 0.008,\n",
    "                y=0.95 * y_min,\n",
    "                s=f\"$m_{R}$\",\n",
    "                c=\"red\",\n",
    "            )\n",
    "        axes[-1].set_xlabel(\"Re $m$\")\n",
    "        fig.canvas.draw_idle()\n",
    "\n",
    "    # Create switch for imag/real/abs\n",
    "    name = \"imag_real\"\n",
    "    sliders._sliders[name] = ipywidgets.RadioButtons(\n",
    "        options=[\"imag\", \"real\", \"abs\"],\n",
    "        description=R\"\\(s\\)-plane plot\",\n",
    "    )\n",
    "    sliders._arg_to_symbol[name] = name\n",
    "\n",
    "    # Create cut-off slider for z-direction\n",
    "    name = \"z_cutoff\"\n",
    "    sliders._sliders[name] = ipywidgets.IntSlider(\n",
    "        value=10,\n",
    "        min=+1,\n",
    "        max=+50,\n",
    "        description=R\"\\(z\\)-cutoff\",\n",
    "    )\n",
    "    sliders._arg_to_symbol[name] = name\n",
    "\n",
    "    # Create GUI\n",
    "    sliders_copy = dict(sliders)\n",
    "    h_boxes = []\n",
    "    for R in range(n_resonances):\n",
    "        buttons = [sliders_copy.pop(f\"m{R}\")]\n",
    "        if n_channels == 1:\n",
    "            buttons.append(sliders_copy.pop(symbol_to_arg[Rf\"\\Gamma_{{{R},0}}\"]))\n",
    "            buttons.append(sliders_copy.pop(symbol_to_arg[Rf\"\\gamma_{{{R},0}}\"]))\n",
    "        h_box = ipywidgets.HBox(buttons)\n",
    "        h_boxes.append(h_box)\n",
    "    remaining_sliders = sorted(\n",
    "        sliders_copy.values(), key=lambda s: (str(type(s)), s.description)\n",
    "    )\n",
    "    if n_channels == 1:\n",
    "        remaining_sliders.remove(sliders[\"i\"])\n",
    "    ui = ipywidgets.VBox(h_boxes + remaining_sliders)\n",
    "    output = ipywidgets.interactive_output(plot3, controls=sliders)\n",
    "    display(ui, output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "if STATIC_WEB_PAGE:\n",
    "    L = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false,
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_f_vector(\n",
    "    n_resonances=2,\n",
    "    n_channels=1,\n",
    "    angular_momentum=L,\n",
    "    title=\"Relativistic $F$-vector, single channel\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "{{ run_interactive }}\n",
    "\n",
    "![](https://user-images.githubusercontent.com/29308176/164993778-1f5987c2-4ff6-45e3-9ef4-cabbcb27b70a.gif)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "toc_visible": true
  },
  "keep_output": 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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}