{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "%config Completer.use_jedi = False\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "\n",
    "# Install on Google Colab\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "from IPython import get_ipython\n",
    "\n",
    "install_packages = \"google.colab\" in str(get_ipython())\n",
    "if install_packages:\n",
    "    for package in [\"expertsystem\", \"graphviz\"]:\n",
    "        subprocess.check_call(\n",
    "            [sys.executable, \"-m\", \"pip\", \"install\", package]\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard lineshapes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{warning}\n",
    "The {doc}`PWA Expert System <index>` has been split up into {doc}`QRules <qrules:index>` and {doc}`AmpForm <ampform:index>`. Please use these packages instead!\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "The formulas and figures on this page have been generated with the lineshapes in the {mod}`.lineshape` module, so as to [glue](https://myst-nb.readthedocs.io/en/latest/use/glue.html) them back in to the API of that module.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import myst_nb\n",
    "import sympy as sp\n",
    "\n",
    "from expertsystem.amplitude.dynamics.lineshape import (\n",
    "    BlattWeisskopf,\n",
    "    relativistic_breit_wigner,\n",
    "    relativistic_breit_wigner_with_ff,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.figure import Figure\n",
    "\n",
    "from expertsystem.reaction.combinatorics import arange\n",
    "\n",
    "\n",
    "def plot_real_imag(\n",
    "    expression: sp.Expr,\n",
    "    variable: sp.Symbol,\n",
    "    x_min: float,\n",
    "    x_max: float,\n",
    "    resolution: int = 100,\n",
    ") -> Figure:\n",
    "    delta = (x_max - x_min) / resolution\n",
    "    x = list(arange(x_min, x_max, delta))\n",
    "    y_real = list(map(lambda x: sp.Abs(expression).subs(variable, x), x))\n",
    "    y_imag = list(map(lambda x: sp.arg(expression).subs(variable, x), x))\n",
    "    fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(8, 8))\n",
    "    for a in ax:\n",
    "        a.xaxis.set_ticks([])\n",
    "        a.yaxis.set_ticks([])\n",
    "    ax_real, ax_imag = ax\n",
    "    ax_imag.set(xlabel=f\"${variable.name}$\")\n",
    "    ax_imag.set(ylabel=f\"imag $f({variable.name})$\")\n",
    "    ax_real.set(ylabel=f\"real $f({variable.name})$\")\n",
    "    ax_real.yaxis.set_ticks([])\n",
    "    ax_imag.yaxis.set_ticks([0, float(sp.pi)])\n",
    "    ax_imag.yaxis.set_ticklabels([0, R\"$\\pi$\"])\n",
    "    ax_imag.set_ylim([0, float(sp.pi)])\n",
    "    ax_real.plot(x, y_real)\n",
    "    ax_imag.plot(x, y_imag)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Form factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The {mod}`expertsystem` uses {class}`.BlattWeisskopf` functions $B_L$ as _barrier factors_ (also called _form factors_):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_L\\left(q\\right)$"
      ],
      "text/plain": [
       "BlattWeisskopf(q, d, L)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{cases} 1 & \\text{for}\\: L = 0 \\\\\\frac{2 d^{2} q^{2}}{d^{2} q^{2} + 1} & \\text{for}\\: L = 1 \\\\\\frac{13 d^{4} q^{4}}{9 d^{2} q^{2} + \\left(d^{2} q^{2} - 3\\right)^{2}} & \\text{for}\\: L = 2 \\\\\\frac{277 d^{6} q^{6}}{d^{2} q^{2} \\left(d^{2} q^{2} - 15\\right)^{2} + \\left(2 d^{2} q^{2} - 5\\right) \\left(18 d^{2} q^{2} - 45\\right)} & \\text{for}\\: L = 3 \\\\\\frac{12746 d^{8} q^{8}}{25 d^{2} q^{2} \\left(2 d^{2} q^{2} - 21\\right)^{2} + \\left(d^{4} q^{4} - 45 d^{2} q^{2} + 105\\right)^{2}} & \\text{for}\\: L = 4 \\end{cases}$"
      ],
      "text/plain": [
       "Piecewise((1, Eq(L, 0)), (2*d**2*q**2/(d**2*q**2 + 1), Eq(L, 1)), (13*d**4*q**4/(9*d**2*q**2 + (d**2*q**2 - 3)**2), Eq(L, 2)), (277*d**6*q**6/(d**2*q**2*(d**2*q**2 - 15)**2 + (2*d**2*q**2 - 5)*(18*d**2*q**2 - 45)), Eq(L, 3)), (12746*d**8*q**8/(25*d**2*q**2*(2*d**2*q**2 - 21)**2 + (d**4*q**4 - 45*d**2*q**2 + 105)**2), Eq(L, 4)))"
      ]
     },
     "metadata": {
      "scrapbook": {
       "mime_prefix": "",
       "name": "BlattWeisskopf"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "q, d, L = sp.symbols(\"q, d, L\", real=True)\n",
    "ff = BlattWeisskopf(q, d, L)\n",
    "display(ff)\n",
    "myst_nb.glue(\"BlattWeisskopf\", ff.doit())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relativistic Breit-Wigner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Without_ form factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See {func}`.relativistic_breit_wigner`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\Gamma m_{0}}{- i \\Gamma m_{0} - m^{2} + m_{0}^{2}}$"
      ],
      "text/plain": [
       "Gamma*m0/(-I*Gamma*m0 - m**2 + m0**2)"
      ]
     },
     "metadata": {
      "scrapbook": {
       "mime_prefix": "",
       "name": "relativistic_breit_wigner"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m, m0, w0 = sp.symbols(\"m m0 Gamma\", real=True)\n",
    "myst_nb.glue(\n",
    "    \"relativistic_breit_wigner\",\n",
    "    relativistic_breit_wigner(m, m0, w0),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plot_real_imag(relativistic_breit_wigner(m, 1.0, 0.3), m, x_min=0, x_max=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _With_ form factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See {func}`.relativistic_breit_wigner_with_ff`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\Gamma m_{0}^{2} \\sqrt{\\frac{\\left(m^{2} - \\left(m_{a} - m_{b}\\right)^{2}\\right) \\left(m^{2} - \\left(m_{a} + m_{b}\\right)^{2}\\right)}{m^{2}}}}{m \\sqrt{\\frac{\\left(m_{0}^{2} - \\left(m_{a} - m_{b}\\right)^{2}\\right) \\left(m_{0}^{2} - \\left(m_{a} + m_{b}\\right)^{2}\\right)}{m_{0}^{2}}} \\left(- \\frac{i \\Gamma m_{0}^{2} \\sqrt{\\frac{\\left(m^{2} - \\left(m_{a} - m_{b}\\right)^{2}\\right) \\left(m^{2} - \\left(m_{a} + m_{b}\\right)^{2}\\right)}{m^{2}}}}{m \\sqrt{\\frac{\\left(m_{0}^{2} - \\left(m_{a} - m_{b}\\right)^{2}\\right) \\left(m_{0}^{2} - \\left(m_{a} + m_{b}\\right)^{2}\\right)}{m_{0}^{2}}}} - m^{2} + m_{0}^{2}\\right)}$"
      ],
      "text/plain": [
       "Gamma*m0**2*sqrt((m**2 - (m_a - m_b)**2)*(m**2 - (m_a + m_b)**2)/m**2)/(m*sqrt((m0**2 - (m_a - m_b)**2)*(m0**2 - (m_a + m_b)**2)/m0**2)*(-I*Gamma*m0**2*sqrt((m**2 - (m_a - m_b)**2)*(m**2 - (m_a + m_b)**2)/m**2)/(m*sqrt((m0**2 - (m_a - m_b)**2)*(m0**2 - (m_a + m_b)**2)/m0**2)) - m**2 + m0**2))"
      ]
     },
     "metadata": {
      "scrapbook": {
       "mime_prefix": "",
       "name": "relativistic_breit_wigner_with_ff"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 0\n",
    "m, m0, w0, ma, mb, meson_radius = sp.symbols(\n",
    "    \"m m0 Gamma m_a m_b q_r\", real=True\n",
    ")\n",
    "myst_nb.glue(\n",
    "    \"relativistic_breit_wigner_with_ff\",\n",
    "    relativistic_breit_wigner_with_ff(\n",
    "        mass=m,\n",
    "        mass0=m0,\n",
    "        gamma0=w0,\n",
    "        m_a=ma,\n",
    "        m_b=mb,\n",
    "        angular_momentum=L,\n",
    "        meson_radius=meson_radius,\n",
    "    ).doit(),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "ma = 0.2\n",
    "mb = 0.3\n",
    "complex_bw_ff = relativistic_breit_wigner_with_ff(\n",
    "    mass=m,\n",
    "    mass0=1.0,\n",
    "    gamma0=0.3,\n",
    "    m_a=ma,\n",
    "    m_b=mb,\n",
    "    angular_momentum=0,\n",
    "    meson_radius=1,\n",
    ")\n",
    "plot_real_imag(complex_bw_ff.doit(), m, x_min=ma + mb, x_max=2);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}