From 40f69e5aabb7da677ba642691c872dc27aa3dca2 Mon Sep 17 00:00:00 2001 From: SVJ_Vitor Date: Fri, 12 Jul 2024 12:15:53 +0200 Subject: [PATCH] DOC: formulate pure SymPy amplitude model (#13) --- .github/workflows/ci-docs.yml | 4 + docs/index.md | 1 + docs/symbolics.ipynb | 409 ++++++++++++++++++++++++++++++++++ 3 files changed, 414 insertions(+) create mode 100644 docs/symbolics.ipynb diff --git a/.github/workflows/ci-docs.yml b/.github/workflows/ci-docs.yml index eeef184..13fadc0 100644 --- a/.github/workflows/ci-docs.yml +++ b/.github/workflows/ci-docs.yml @@ -1,4 +1,8 @@ name: Documentation +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: |- + ${{ github.ref != format('refs/heads/{0}', github.event.repository.default_branch) }} on: pull_request: branches: [main] diff --git a/docs/index.md b/docs/index.md index dcab849..c3dad8b 100644 --- a/docs/index.md +++ b/docs/index.md @@ -16,6 +16,7 @@ address to the issue ✅ [ComPWA/gluex-nstar#1](https://github.com/ComPWA/g ```{toctree} reaction +symbolics pgamma-state branching-fraction ampform/etapi0p diff --git a/docs/symbolics.ipynb b/docs/symbolics.ipynb new file mode 100644 index 0000000..bbb5f4f --- /dev/null +++ b/docs/symbolics.ipynb @@ -0,0 +1,409 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# pγ→ηπ⁰p with SymPy only" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This section is a follow-up to formulate the amplitude model for the $\\gamma p \\to \\eta\\pi^0 p$ channel example symbolically. See **[TR‑033](https://compwa.github.io/report/033)** for a purely numerical tutorial.\n", + "\n", + "The model we want to implement is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{array}{rcl}\n", + "I &=& \\left|A^{12} + A^{23} + A^{31}\\right|^2 \\\\\n", + "A^{12} &=& \\frac{\\sum a_m Y_2^m (\\Omega_1)}{s_{12}-m^2_{a_2}+im_{a_2} \\Gamma_{a_2}} \\\\\n", + "A^{23} &=& \\frac{\\sum b_m Y_1^m (\\Omega_2)}{s_{23}-m^2_{\\Delta}+im_{\\Delta} \\Gamma_{\\Delta}} \\\\\n", + "A^{31} &=& \\frac{c_0}{s_{31}-m^2_{N^*}+im_{N^*} \\Gamma_{N^*}} \\,,\n", + "\\end{array}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $1=\\eta$, $2=\\pi^0$, and $3=p$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Import Python libraries" + }, + "tags": [ + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import sympy as sp" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### $A^{12}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s12, m_a2, Gamma_a2 = sp.symbols(r\"s_{12} m_{a_2} \\Gamma_{a_2}\")\n", + "theta1, phi1 = sp.symbols(\"theta_1 phi_1\")\n", + "a = sp.IndexedBase(\"a\")\n", + "m = sp.symbols(\"m\", cls=sp.Idx)\n", + "A12 = sp.Sum(a[m] * sp.Ynm(2, m, theta1, phi1), (m, -2, 2)) / (\n", + " s12 - m_a2**2 + sp.I * m_a2 * Gamma_a2\n", + ")\n", + "A12" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sp.Ynm(2, 1, 1, 1).expand(func=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sp.Ynm(2, 1, 1, 1).expand(func=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A12.free_symbols" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A12_func = sp.lambdify(\n", + " [s12, a[-2], a[-1], a[0], a[1], a[2], m_a2, Gamma_a2, theta1, phi1],\n", + " A12.doit().expand(func=True),\n", + ")\n", + "A12_func" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### $A^{23}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s23, m_delta, gamma_delta = sp.symbols(\"s_{23} m_Delta Gamma_Delta\")\n", + "b = sp.IndexedBase(\"b\")\n", + "m = sp.symbols(\"m\", cls=sp.Idx)\n", + "\n", + "theta2, phi2 = sp.symbols(\"theta_2 phi_2\")\n", + "A23 = sp.Sum(b[m] * sp.Ynm(1, m, theta2, phi2), (m, -1, 1)) / (\n", + " s23 - m_delta**2 + sp.I * m_delta * gamma_delta\n", + ")\n", + "A23" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A23_func = sp.lambdify(\n", + " [s23, b[-1], b[0], b[1], m_delta, gamma_delta, theta2, phi2],\n", + " A23.doit().expand(func=True),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### $A^{31}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s31, c0, m_nstar, gamma_nstar = sp.symbols(r\"s_{31} c_0 m_{N^*} \\Gamma_{N^*}\")\n", + "A31 = c0 / (s31 - m_nstar**2 + sp.I * m_nstar * gamma_nstar)\n", + "A31" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A31_func = sp.lambdify([s31, c0, m_nstar, gamma_nstar], A31)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### $A = A^{12}+A^{23}+A^{31}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "intensity_expr = sp.Abs(A12 + A23 + A31) ** 2\n", + "intensity_expr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "scroll-input" + ] + }, + "outputs": [], + "source": [ + "intensity_func = sp.lambdify(\n", + " [\n", + " s12,\n", + " a[-2],\n", + " a[-1],\n", + " a[0],\n", + " a[1],\n", + " a[2],\n", + " m_a2,\n", + " Gamma_a2,\n", + " theta1,\n", + " phi1,\n", + " s23,\n", + " b[-1],\n", + " b[0],\n", + " b[1],\n", + " m_delta,\n", + " gamma_delta,\n", + " theta2,\n", + " phi2,\n", + " s31,\n", + " c0,\n", + " m_nstar,\n", + " gamma_nstar,\n", + " ],\n", + " intensity_expr.doit().expand(func=True),\n", + ")\n", + "intensity_func" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a_vals = [0, 0.5, 3.5, 4, 2.5]\n", + "b_vals = [-1.5, 4, 0.5]\n", + "c0_val = 2.5\n", + "m_a2_val = 1.32\n", + "m_delta_val = 1.54\n", + "m_nstar_val = 1.87\n", + "Gamma_a2_val = 0.1\n", + "gamma_delta_val = 0.1\n", + "gamma_nstar_val = 0.1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "phi = np.pi / 4\n", + "theta = np.pi / 4\n", + "s_values = np.linspace(0, 5, num=500)\n", + "A12_values = A12_func(s_values, *a_vals, m_a2_val, Gamma_a2_val, theta, phi)\n", + "A23_values = A23_func(s_values, *b_vals, m_delta_val, gamma_delta_val, theta, phi)\n", + "A31_values = A31_func(s_values, c0_val, m_nstar_val, gamma_nstar_val)" + ] + }, + { + "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", + "full-width" + ] + }, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(figsize=(10, 9), nrows=3)\n", + "angle_projection_text = R\"$\\phi=\\frac{\\pi}{4},\\theta=\\frac{\\pi}{4}$\"\n", + "for i, A_values in enumerate([A12_values, A23_values, A31_values]):\n", + " ax = axes[i]\n", + " recoil_id = i + 1\n", + " decay_products = sorted({1, 2, 3} - {recoil_id})\n", + " subsystem = \"\".join(map(str, decay_products))\n", + " ax.plot(s_values, A_values.imag, label=\"Imaginary Part\", c=\"blue\", ls=\"-\")\n", + " ax.plot(s_values, A_values.real, label=\"Real Part\", c=\"red\", ls=\"--\")\n", + " ax.plot(\n", + " s_values,\n", + " np.abs(A_values) ** 2,\n", + " label=f\"Intensity $I^{{{subsystem}}}=|A^{{{subsystem}}}|^2$\",\n", + " c=\"g\",\n", + " ls=\"-.\",\n", + " )\n", + " ax.plot(\n", + " s_values,\n", + " np.angle(A_values),\n", + " label=\"Phase (angle)\",\n", + " c=\"m\",\n", + " ls=\":\",\n", + " )\n", + " ax.set_title(Rf\"Components of $A^{{{subsystem}}}$ vs s ({angle_projection_text})\")\n", + " ax.set_xlabel(f\"$s_{{{subsystem}}}$\")\n", + " ax.set_ylabel(Rf\"$A^{{{subsystem}}}$ components\")\n", + " ax.legend()\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unitarity is preserved in each of the subsystems (assuming fixed $\\phi,\\theta$), because we assume there is only one resonance in the subsystem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "full-width" + ] + }, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(figsize=(12, 4), ncols=3)\n", + "for i, A_values in enumerate([A12_values, A23_values, A31_values]):\n", + " ax = axes[i]\n", + " recoil_id = i + 1\n", + " decay_products = sorted({1, 2, 3} - {recoil_id})\n", + " subsystem = \"\".join(map(str, decay_products))\n", + " ax.set_xlabel(Rf\"$\\mathrm{{Re}}\\,A^{{{subsystem}}}$\")\n", + " ax.set_ylabel(Rf\"$\\mathrm{{Im}}\\,A^{{{subsystem}}}$\")\n", + " ax.set_title(Rf\"$A^{{{subsystem}}}$ ({angle_projection_text})\")\n", + " ax.plot(A_values.real, A_values.imag)\n", + "fig.suptitle(\"Argand diagrams\")\n", + "fig.tight_layout()\n", + "plt.show(fig)" + ] + } + ], + "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 +}