diff --git a/.cspell.json b/.cspell.json index 75f27e90..e221f539 100644 --- a/.cspell.json +++ b/.cspell.json @@ -65,6 +65,7 @@ "breit", "compwa", "conda", + "Dalitz", "defaultdict", "eval", "flatté", @@ -73,6 +74,7 @@ "helicity", "itertools", "kwargs", + "Källén", "lambdification", "lambdified", "lambdifies", @@ -115,6 +117,7 @@ "astype", "autolaunch", "autonumbering", + "autoscale", "autoupdate", "axhline", "axvline", @@ -122,6 +125,7 @@ "bdist", "bgcolor", "boldsymbol", + "byckling", "cbff", "celltoolbar", "clim", @@ -155,6 +159,7 @@ "filterwarnings", "fontcolor", "fontsize", + "framealpha", "getitem", "getsource", "graphviz", @@ -162,6 +167,7 @@ "hasattr", "heatmap", "histtype", + "hspace", "imag", "infty", "iplt", @@ -179,6 +185,7 @@ "kutschke", "lambdifier", "lambdifygenerated", + "lightgray", "linecap", "linejoin", "linestyle", @@ -198,6 +205,7 @@ "maxsize", "meshgrid", "multiline", + "nansum", "nbconvert", "nbformat", "nbmake", @@ -233,6 +241,7 @@ "quadpy", "rankdir", "redeboer", + "relim", "repr", "richman", "rpartition", @@ -241,6 +250,8 @@ "rtfd", "rules's", "savefig", + "scalex", + "scaley", "scimath", "sdist", "seealso", @@ -265,6 +276,7 @@ "venv", "vmax", "vmin", + "wspace", "xdata", "xlabel", "xlim", diff --git a/docs/bibliography.bib b/docs/bibliography.bib index ccaacb47..b81a5385 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -38,6 +38,16 @@ @book{beckTestDrivenDevelopmentExample2003 lccn = {QA76.76.T48 B43 2003} } +@book{bycklingParticleKinematics1973, + title = {Particle {{Kinematics}}}, + author = {Byckling, Eero and Kajantie, Keijo}, + year = {1973}, + publisher = {{Wiley}}, + address = {{London, New York}}, + isbn = {978-0-471-12885-4}, + lccn = {QC794.6.K5 B95} +} + @article{chungPartialWaveAnalysis1995, title = {{Partial wave analysis in 𝐾-matrix formalism}}, author = {Chung, Suh-Urk and Brose, J. and Hackmann, R. and Klempt, E. and Spanier, S. and Strassburger, C.}, diff --git a/docs/conf.py b/docs/conf.py index f79ab980..17665ed7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -279,6 +279,7 @@ def get_minor_version(package_name: str) -> str: "report/014*", "report/015*", "report/016*", + "report/017*", ] nb_output_stderr = "remove" nb_render_priority = { diff --git a/docs/report/017.ipynb b/docs/report/017.ipynb new file mode 100644 index 00000000..468d5e9a --- /dev/null +++ b/docs/report/017.ipynb @@ -0,0 +1,752 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hideCode": true, + "hideOutput": true, + "hidePrompt": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "skip" + }, + "tags": [ + "remove-cell" + ] + }, + "outputs": [], + "source": [ + "%config InlineBackend.figure_formats = ['svg']\n", + "import os\n", + "\n", + "STATIC_WEB_PAGE = {\"EXECUTE_NB\", \"READTHEDOCS\"}.intersection(os.environ)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{autolink-concat}\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "````{margin}\n", + "```{spec} Symbolic phase space boundary for a three-body decay\n", + ":id: TR-017\n", + ":status: Finished\n", + ":tags: sympy;physics\n", + "\n", + "This reports shows how define the physical phase space region on a Dalitz plot using a Kibble function.\n", + "```\n", + "````\n", + "\n", + "# Phase space for a three-body decay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "remove-cell" + ] + }, + "outputs": [], + "source": [ + "%pip install -q ampform==0.14.0 sympy==1.10.1 tensorwaves==0.4.5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "from __future__ import annotations\n", + "\n", + "import warnings\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import sympy as sp\n", + "from ampform.sympy import (\n", + " UnevaluatedExpression,\n", + " create_expression,\n", + " implement_doit_method,\n", + " make_commutative,\n", + ")\n", + "from IPython.display import Math\n", + "from ipywidgets import FloatSlider, VBox, interactive_output\n", + "from matplotlib.axis import Axis\n", + "from matplotlib.contour import QuadContourSet\n", + "from matplotlib.lines import Line2D\n", + "from matplotlib.patches import Patch\n", + "from tensorwaves.function.sympy import create_parametrized_function\n", + "\n", + "warnings.filterwarnings(\"ignore\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Kinematics for a three-body decay $0 \\to 123$ can be fully described by two **Mandelstam variables** $\\sigma_1, \\sigma_2$, because the third variable $\\sigma_3$ can be expressed in terms $\\sigma_1, \\sigma_2$, the mass $m_0$ of the initial state, and the masses $m_1, m_2, m_3$ of the final state. As can be seen, the roles of $\\sigma_1, \\sigma_2, \\sigma_3$ are interchangeable.\n", + "\n", + "```{margin}\n", + "See Eq. (1.2) in {cite}`bycklingParticleKinematics1973`\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\sigma_{3} = m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \\sigma_{1} - \\sigma_{2}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3) -> sp.Expr:\n", + " return m0**2 + m1**2 + m2**2 + m3**2 - sigma1 - sigma2\n", + "\n", + "\n", + "m0, m1, m2, m3 = sp.symbols(\"m:4\")\n", + "s1, s2, s3 = sp.symbols(\"sigma1:4\")\n", + "computed_s3 = compute_third_mandelstam(s1, s2, m0, m1, m2, m3)\n", + "Math(Rf\"{sp.latex(s3)} = {sp.latex(computed_s3)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The phase space is defined by the closed area that satisfies the condition $\\phi(\\sigma_1,\\sigma_2) \\leq 0$, where $\\phi$ is a **Kibble function**:\n", + "\n", + "\n", + "```{margin}\n", + "See §V.2 in {cite}`bycklingParticleKinematics1973`\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "scroll-input", + "hide-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\phi\\left(\\sigma_{1}, \\sigma_{2}\\right) = \\lambda\\left(\\lambda\\left(\\sigma_{2}, m_{2}^{2}, m_{0}^{2}\\right), \\lambda\\left(\\sigma_{3}, m_{3}^{2}, m_{0}^{2}\\right), \\lambda\\left(\\sigma_{1}, m_{1}^{2}, m_{0}^{2}\\right)\\right)$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "@make_commutative\n", + "@implement_doit_method\n", + "class Kibble(UnevaluatedExpression):\n", + " def __new__(\n", + " cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints\n", + " ) -> Kibble:\n", + " return create_expression(\n", + " cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints\n", + " )\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " sigma1, sigma2, sigma3, m0, m1, m2, m3 = self.args\n", + " return Källén(\n", + " Källén(sigma2, m2**2, m0**2),\n", + " Källén(sigma3, m3**2, m0**2),\n", + " Källén(sigma1, m1**2, m0**2),\n", + " )\n", + "\n", + " def _latex(self, printer, *args):\n", + " sigma1, sigma2, *_ = map(printer._print, self.args)\n", + " return Rf\"\\phi\\left({sigma1}, {sigma2}\\right)\"\n", + "\n", + "\n", + "@make_commutative\n", + "@implement_doit_method\n", + "class Källén(UnevaluatedExpression):\n", + " def __new__(cls, x, y, z, **hints) -> Källén:\n", + " return create_expression(cls, x, y, z, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " x, y, z = self.args\n", + " return x**2 + y**2 + z**2 - 2 * x * y - 2 * y * z - 2 * z * x\n", + "\n", + " def _latex(self, printer, *args):\n", + " x, y, z = map(printer._print, self.args)\n", + " return Rf\"\\lambda\\left({x}, {y}, {z}\\right)\"\n", + "\n", + "\n", + "kibble = Kibble(s1, s2, s3, m0, m1, m2, m3)\n", + "Math(Rf\"{sp.latex(kibble)} = {sp.latex(kibble.doit(deep=False))}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and $\\lambda$ is the **Källén function**:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\lambda\\left(x, y, z\\right) = x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x, y, z = sp.symbols(\"x:z\")\n", + "expr = Källén(x, y, z)\n", + "Math(f\"{sp.latex(expr)} = {sp.latex(expr.doit())}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Any distribution over the phase space can now be defined using a two-dimensional grid over a Mandelstam pair $\\sigma_1,\\sigma_2$ of choice, with the condition $\\phi(\\sigma_1,\\sigma_2)<0$ selecting the values that are physically allowed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{cases} 1 & \\text{for}\\: \\phi\\left(\\sigma_{1}, \\sigma_{2}\\right) \\leq 0 \\\\\\text{NaN} & \\text{otherwise} \\end{cases}$" + ], + "text/plain": [ + "Piecewise((1, Kibble(sigma1, sigma2, m0**2 + m1**2 + m2**2 + m3**2 - sigma1 - sigma2, m0, m1, m2, m3) <= 0), (nan, True))" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def is_within_phasespace(\n", + " sigma1, sigma2, m0, m1, m2, m3, outside_value=sp.nan\n", + ") -> sp.Piecewise:\n", + " sigma3 = compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3)\n", + " kibble = Kibble(sigma1, sigma2, sigma3, m0, m1, m2, m3)\n", + " return sp.Piecewise(\n", + " (1, sp.LessThan(kibble, 0)),\n", + " (outside_value, True),\n", + " )\n", + "\n", + "\n", + "is_within_phasespace(s1, s2, m0, m1, m2, m3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "phsp_expr = is_within_phasespace(s1, s2, m0, m1, m2, m3, outside_value=0)\n", + "phsp_func = create_parametrized_function(\n", + " phsp_expr.doit(),\n", + " parameters={m0: 2.2, m1: 0.2, m2: 0.4, m3: 0.4},\n", + " backend=\"numpy\",\n", + ")" + ] + }, + { + "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", + "hide-input", + "remove-output" + ] + }, + "outputs": [], + "source": [ + "sliders = {\n", + " \"m0\": FloatSlider(description=\"m0\", max=3, value=2.1, step=0.01),\n", + " \"m1\": FloatSlider(description=\"m1\", max=2, value=0.2, step=0.01),\n", + " \"m2\": FloatSlider(description=\"m2\", max=2, value=0.4, step=0.01),\n", + " \"m3\": FloatSlider(description=\"m3\", max=2, value=0.4, step=0.01),\n", + "}\n", + "\n", + "resolution = 300\n", + "X, Y = np.meshgrid(\n", + " np.linspace(0, 4, num=resolution),\n", + " np.linspace(0, 4, num=resolution),\n", + ")\n", + "data = {\"sigma1\": X, \"sigma2\": Y}\n", + "\n", + "sidebar_ratio = 0.15\n", + "fig, ((ax1, _), (ax, ax2)) = plt.subplots(\n", + " figsize=(7, 7),\n", + " ncols=2,\n", + " nrows=2,\n", + " gridspec_kw=dict(\n", + " height_ratios=[sidebar_ratio, 1],\n", + " width_ratios=[1, sidebar_ratio],\n", + " ),\n", + ")\n", + "_.remove()\n", + "ax.set_xlim(0, 4)\n", + "ax.set_ylim(0, 4)\n", + "ax.set_xlabel(R\"$\\sigma_1$\")\n", + "ax.set_ylabel(R\"$\\sigma_2$\")\n", + "ax.set_xticks(range(5))\n", + "ax.set_yticks(range(5))\n", + "ax1.set_xlim(0, 4)\n", + "ax2.set_ylim(0, 4)\n", + "for a in [ax1, ax2]:\n", + " a.set_xticks([])\n", + " a.set_yticks([])\n", + " a.axis(\"off\")\n", + "fig.tight_layout()\n", + "fig.subplots_adjust(wspace=0, hspace=0)\n", + "\n", + "fig.canvas.toolbar_visible = False\n", + "fig.canvas.header_visible = False\n", + "fig.canvas.footer_visible = False\n", + "\n", + "MESH: QuadContourSet | None = None\n", + "PROJECTIONS: tuple[Line2D, Line2D] = None\n", + "BOUNDARIES: list[Line2D] | None = None\n", + "\n", + "\n", + "def plot(**parameters):\n", + " draw_boundaries(\n", + " parameters[\"m0\"],\n", + " parameters[\"m1\"],\n", + " parameters[\"m2\"],\n", + " parameters[\"m3\"],\n", + " )\n", + " global MESH, PROJECTIONS\n", + " if MESH is not None:\n", + " for coll in MESH.collections:\n", + " ax.collections.remove(coll)\n", + " phsp_func.update_parameters(parameters)\n", + " Z = phsp_func(data)\n", + " MESH = ax.contour(X, Y, Z, colors=\"black\")\n", + " contour = MESH.collections[0]\n", + " contour.set_facecolor(\"lightgray\")\n", + " x = X[0]\n", + " y = Y[:, 0]\n", + " Zx = np.nansum(Z, axis=0)\n", + " Zy = np.nansum(Z, axis=1)\n", + " if PROJECTIONS is None:\n", + " PROJECTIONS = (\n", + " ax1.plot(x, Zx, c=\"black\", lw=2)[0],\n", + " ax2.plot(Zy, y, c=\"black\", lw=2)[0],\n", + " )\n", + " else:\n", + " PROJECTIONS[0].set_data(x, Zx)\n", + " PROJECTIONS[1].set_data(Zy, y)\n", + " ax1.relim()\n", + " ax2.relim()\n", + " ax1.autoscale_view(scalex=False)\n", + " ax2.autoscale_view(scaley=False)\n", + " create_legend(ax)\n", + " fig.canvas.draw()\n", + "\n", + "\n", + "def draw_boundaries(m0, m1, m2, m3) -> None:\n", + " global BOUNDARIES\n", + " s1_min = (m2 + m3) ** 2\n", + " s1_max = (m0 - m1) ** 2\n", + " s2_min = (m1 + m3) ** 2\n", + " s2_max = (m0 - m2) ** 2\n", + " if BOUNDARIES is None:\n", + " BOUNDARIES = [\n", + " ax.axvline(s1_min, c=\"red\", ls=\"dotted\", label=\"$(m_2+m_3)^2$\"),\n", + " ax.axhline(s2_min, c=\"blue\", ls=\"dotted\", label=\"$(m_1+m_3)^2$\"),\n", + " ax.axvline(s1_max, c=\"red\", ls=\"dashed\", label=\"$(m_0-m_1)^2$\"),\n", + " ax.axhline(s2_max, c=\"blue\", ls=\"dashed\", label=\"$(m_0-m_2)^2$\"),\n", + " ]\n", + " else:\n", + " BOUNDARIES[0].set_data(get_line_data(s1_min))\n", + " BOUNDARIES[1].set_data(get_line_data(s2_min, horizontal=True))\n", + " BOUNDARIES[2].set_data(get_line_data(s1_max))\n", + " BOUNDARIES[3].set_data(get_line_data(s2_max, horizontal=True))\n", + "\n", + "\n", + "def create_legend(ax: Axis):\n", + " if ax.get_legend() is not None:\n", + " return\n", + " label = Rf\"${sp.latex(kibble)}\\leq0$\"\n", + " ax.legend(\n", + " handles=[\n", + " Patch(label=label, ec=\"black\", fc=\"lightgray\"),\n", + " *BOUNDARIES,\n", + " ],\n", + " loc=\"upper right\",\n", + " facecolor=\"white\",\n", + " framealpha=1,\n", + " )\n", + "\n", + "\n", + "def get_line_data(value, horizontal: bool = False) -> np.ndarray:\n", + " pair = (value, value)\n", + " if horizontal:\n", + " return np.array([(0, 1), pair])\n", + " return np.array([pair, (0, 1)])\n", + "\n", + "\n", + "output = interactive_output(plot, controls=sliders)\n", + "VBox([output, *sliders.values()])" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "{{ run_interactive }}\n", + "\n", + "```{image} https://user-images.githubusercontent.com/29308176/165257660-af66c9c1-2da9-475f-b4c4-4a52950acd00.gif\n", + ":alt: Cell output - interactive Dalitz plot\n", + ":align: center\n", + "```\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The phase space boundary can be described analytically in terms of $\\sigma_1$ or $\\sigma_2$, in which case there are two solutions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sol1, sol2 = sp.solve(kibble.doit().subs(s3, computed_s3), s2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "remove-input", + "keep_output", + "full-width" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{array}{c}\n", + " \\frac{- m_{0}^{2} m_{2}^{2} + m_{0}^{2} m_{3}^{2} + m_{0}^{2} \\sigma_{1} + m_{1}^{2} m_{2}^{2} - m_{1}^{2} m_{3}^{2} + m_{1}^{2} \\sigma_{1} + m_{2}^{2} \\sigma_{1} + m_{3}^{2} \\sigma_{1} - \\sigma_{1}^{2} - \\sqrt{\\left(m_{0}^{2} - 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{0}^{2} + 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} - 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} + 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right)}}{2 \\sigma_{1}} \\\\\n", + " \\frac{- m_{0}^{2} m_{2}^{2} + m_{0}^{2} m_{3}^{2} + m_{0}^{2} \\sigma_{1} + m_{1}^{2} m_{2}^{2} - m_{1}^{2} m_{3}^{2} + m_{1}^{2} \\sigma_{1} + m_{2}^{2} \\sigma_{1} + m_{3}^{2} \\sigma_{1} - \\sigma_{1}^{2} + \\sqrt{\\left(m_{0}^{2} - 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{0}^{2} + 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} - 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} + 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right)}}{2 \\sigma_{1}} \\\\\n", + "\\end{array}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "latex = R\"\\begin{array}{c}\" + \"\\n\"\n", + "for sol in [sol1, sol2]:\n", + " latex += Rf\" {sp.latex(sol)} \\\\\" + \"\\n\"\n", + "latex += R\"\\end{array}\"\n", + "Math(latex)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The boundary cannot be parametrized analytically in polar coordinates, but there is a numeric solution. The idea is to solve the condition $\\phi(\\sigma_1,\\sigma_2)=0$ after the following substitutions:\n", + "\n", + "```{margin}\n", + "See {cite}`bycklingParticleKinematics1973`, pp. 109–112\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "T0, T1, T2, T3 = sp.symbols(\"T0:4\")\n", + "r, theta = sp.symbols(\"r theta\", nonnegative=True)\n", + "substitutions = {\n", + " s1: (m2 + m3) ** 2 + T1,\n", + " s2: (m1 + m3) ** 2 + T2,\n", + " s3: (m1 + m2) ** 2 + T3,\n", + " T1: T0 / 3 - r * sp.cos(theta),\n", + " T2: T0 / 3 - r * sp.cos(theta + 2 * sp.pi / 3),\n", + " T3: T0 / 3 - r * sp.cos(theta + 4 * sp.pi / 3),\n", + " T0: m0**2\n", + " + m1**2\n", + " + m2**2\n", + " + m3**2\n", + " - (m2 + m3) ** 2\n", + " - (m1 + m3) ** 2\n", + " - (m1 + m2) ** 2,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "remove-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{array}{rcl}\n", + " \\sigma_{1} &=& T_{1} + \\left(m_{2} + m_{3}\\right)^{2} \\\\\n", + " \\sigma_{2} &=& T_{2} + \\left(m_{1} + m_{3}\\right)^{2} \\\\\n", + " \\sigma_{3} &=& T_{3} + \\left(m_{1} + m_{2}\\right)^{2} \\\\\n", + " T_{1} &=& \\frac{T_{0}}{3} - r \\cos{\\left(\\theta \\right)} \\\\\n", + " T_{2} &=& \\frac{T_{0}}{3} + r \\sin{\\left(\\theta + \\frac{\\pi}{6} \\right)} \\\\\n", + " T_{3} &=& \\frac{T_{0}}{3} + r \\cos{\\left(\\theta + \\frac{\\pi}{3} \\right)} \\\\\n", + " T_{0} &=& m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \\left(m_{1} + m_{2}\\right)^{2} - \\left(m_{1} + m_{3}\\right)^{2} - \\left(m_{2} + m_{3}\\right)^{2} \\\\\n", + "\\end{array}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "latex = R\"\\begin{array}{rcl}\" + \"\\n\"\n", + "for item in substitutions.items():\n", + " lhs, rhs = map(sp.latex, item)\n", + " latex += Rf\" {lhs} &=& {rhs} \\\\\" + \"\\n\"\n", + "latex += R\"\\end{array}\"\n", + "Math(latex)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For every value of $\\theta \\in [0, 2\\pi)$, the value of $r$ can now be found by solving the condition $\\phi(r, \\theta)=0$. Note that $\\phi(r, \\theta)$ is a **cubic polynomial of $r$**. For instance, if we take $m_0=5, m_1=2, m_{2,3}=1$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input", + "remove-output" + ] + }, + "outputs": [], + "source": [ + "phi_r = (\n", + " kibble.doit()\n", + " .subs(substitutions) # substitute sigmas\n", + " .subs(substitutions) # substitute T123\n", + " .subs(substitutions) # substitute T0\n", + " .subs({m0: 5, m1: 2, m2: 1, m3: 1})\n", + " .simplify()\n", + " .collect(r)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "remove-input", + "keep_output" + ] + }, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\begin{eqnarray}\n", + "\\phi(r, \\theta) & = & r^{3} \\cdot \\left(56 \\sqrt{3} \\sin{\\left(\\theta \\right)} + 400 \\cos^{3}{\\left(\\theta \\right)} - 356 \\cos{\\left(\\theta \\right)} + 112 \\cos{\\left(\\theta + \\frac{\\pi}{3} \\right)}\\right) \\nonumber\\\\\n", + "& & + r^{2} \\cdot \\left(2000 \\cos^{2}{\\left(\\theta \\right)} + 2100\\right) \\nonumber\\\\\n", + "& & - 4800 r \\cos{\\left(\\theta \\right)} \\nonumber\\\\\n", + "& & + -25200 \n", + "\\end{eqnarray}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lhs = sp.Symbol(R\"\\phi(r, \\theta)\")\n", + "latex = sp.multiline_latex(lhs, phi_r, environment=\"eqnarray\")\n", + "Math(latex)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The lowest value of $r$ that satisfies $\\phi(r,\\theta)=0$ defines the phase space boundary." + ] + } + ], + "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.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}