implement LS-coupling generator functions

docs/lambda-decay.ipynb

@@ -60,6 +60,11 @@
     "    RelativisticBreitWigner,\n",
     ")\n",
     "from polarization.io import as_latex\n",
+    "from polarization.spin import (\n",
+    "    create_spin_range,\n",
+    "    filter_parity_violating_ls,\n",
+    "    generate_ls_couplings,\n",
+    ")\n",
     "\n",
     "LOGGER = logging.getLogger()\n",
     "LOGGER.setLevel(logging.ERROR)\n",
@@ -173,44 +178,17 @@
     "\n",
     "    @staticmethod\n",
     "    def generate_ls(particle: Particle, chain_id: int) -> Resonance:\n",
-    "        LS_prod = generate_ls(Λc, particle, siblings[chain_id], strong=False)\n",
+    "        LS_prod = generate_ls_couplings(Λc.spin, particle.spin, siblings[chain_id].spin)\n",
     "        LS_prod = [L for L, _ in LS_prod]\n",
-    "        LS_dec = generate_ls(particle, *decay_products[chain_id])\n",
+    "        child1, child2 = decay_products[chain_id]\n",
+    "        LS_dec = generate_ls_couplings(particle.spin, child1.spin, child2.spin)\n",
+    "        LS_dec = filter_parity_violating_ls(\n",
+    "            LS_dec, particle.parity, child1.parity, child2.parity\n",
+    "        )\n",
     "        LS_dec = [L for L, _ in LS_dec]\n",
     "        return Resonance(particle, l_R=min(LS_dec), l_Λc=min(LS_prod))\n",
     "\n",
     "\n",
-    "def generate_ls(\n",
-    "    parent: Particle,\n",
-    "    child1: Particle,\n",
-    "    child2: Particle,\n",
-    "    strong: bool = True,\n",
-    "    max_L: int = 3,\n",
-    "):\n",
-    "    s1 = child1.spin\n",
-    "    s2 = child2.spin\n",
-    "    s_values = arange(abs(s1 - s2), s1 + s2)\n",
-    "    LS_values = set()\n",
-    "    for S in s_values:\n",
-    "        for L in arange(0, max_L):\n",
-    "            if not abs(L - S) <= parent.spin <= L + S:\n",
-    "                continue\n",
-    "            η0, η1, η2 = [\n",
-    "                int(parent.parity),\n",
-    "                int(child1.parity),\n",
-    "                int(child2.parity),\n",
-    "            ]\n",
-    "            if strong and η0 != η1 * η2 * (-1) ** L:\n",
-    "                continue\n",
-    "            LS_values.add((L, S))\n",
-    "    return sorted(LS_values)\n",
-    "\n",
-    "\n",
-    "def arange(x1, x2):\n",
-    "    spin_range = np.arange(float(x1), +float(x2) + 0.5)\n",
-    "    return list(map(sp.Rational, spin_range))\n",
-    "\n",
-    "\n",
     "resonance_particles = {\n",
     "    chain_id: [PDG[name] for name in names]\n",
     "    for chain_id, names in resonance_names.items()\n",
@@ -506,10 +484,6 @@
     "    return Str(particle.latex)\n",
     "\n",
     "\n",
-    "def create_spin_range(j):\n",
-    "    return arange(-j, +j)\n",
-    "\n",
-    "\n",
     "display(\n",
     "    formulate_chain_amplitude(1, ν, λ),\n",
     "    formulate_chain_amplitude(2, ν, λ),\n",

setup.cfg

@@ -13,6 +13,7 @@ setup_requires =
 install_requires =
     attrs >=20.1.0  # on_setattr and https://www.attrs.org/en/stable/api.html#next-gen
     ampform >=0.14.0
+    numpy
     sympy >=1.10  # module sympy.printing.numpy and array expressions with shape kwarg
 packages = find:
 package_dir =

src/polarization/spin.py

@@ -0,0 +1,69 @@
+from __future__ import annotations
+
+from typing import SupportsFloat, SupportsInt
+
+import numpy as np
+import sympy as sp
+
+
+def generate_ls_couplings(
+    parent_spin: SupportsFloat,
+    child1_spin: SupportsFloat,
+    child2_spin: SupportsFloat,
+    max_L: int = 3,
+) -> list[tuple[int, sp.Rational]]:
+    r"""
+    >>> generate_ls_couplings(1.5, 0.5, 0)
+    [(1, 1/2), (2, 1/2)]
+    """
+    s1 = float(child1_spin)
+    s2 = float(child2_spin)
+    angular_momenta = create_rational_range(0, max_L)
+    coupled_spins = create_rational_range(abs(s1 - s2), s1 + s2)
+    ls_couplings = {
+        (int(L), S)
+        for L in angular_momenta
+        for S in coupled_spins
+        if abs(L - S) <= parent_spin <= L + S
+    }
+    return sorted(ls_couplings)
+
+
+def filter_parity_violating_ls(
+    ls_couplings: list[tuple[int, sp.Rational]],
+    parent_parity: SupportsInt,
+    child1_parity: SupportsInt,
+    child2_parity: SupportsInt,
+) -> list[tuple[int, sp.Rational]]:
+    r"""
+    >>> LS = generate_ls_couplings(0.5, 1.5, 0)  # Λc → Λ(1520)π
+    >>> LS
+    [(1, 3/2), (2, 3/2)]
+    >>> filter_parity_violating_ls(LS, +1, -1, -1)
+    [(2, 3/2)]
+    """
+    η0, η1, η2 = (
+        int(parent_parity),
+        int(child1_parity),
+        int(child2_parity),
+    )
+    return [(L, S) for L, S in ls_couplings if η0 == η1 * η2 * (-1) ** L]
+
+
+def create_spin_range(spin: SupportsFloat) -> list[sp.Rational]:
+    """
+    >>> create_spin_range(1.5)
+    [-3/2, -1/2, 1/2, 3/2]
+    """
+    return create_rational_range(-spin, spin)
+
+
+def create_rational_range(
+    __from: SupportsFloat, __to: SupportsFloat
+) -> list[sp.Rational]:
+    """
+    >>> create_rational_range(-0.5, +1.5)
+    [-1/2, 1/2, 3/2]
+    """
+    spin_range = np.arange(float(__from), +float(__to) + 0.5)
+    return list(map(sp.Rational, spin_range))