diff --git a/docs/usage/dynamics.ipynb b/docs/usage/dynamics.ipynb index 2f3a957e7..629aaa20d 100644 --- a/docs/usage/dynamics.ipynb +++ b/docs/usage/dynamics.ipynb @@ -148,7 +148,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "AmpForm uses Blatt-Weisskopf functions $B_L$ as _barrier factors_ (also called _form factors_, see {class}`.BlattWeisskopfSquared`):" + "AmpForm uses Blatt-Weisskopf functions $B_L$ as _barrier factors_ (also called _form factors_, see {class}`.BlattWeisskopfSquared` and **[TR-029](https://compwa.github.io/report/029)**):" ] }, { @@ -161,8 +161,8 @@ "source": [ "from ampform.dynamics import BlattWeisskopfSquared\n", "\n", - "L = sp.Symbol(\"L\", integer=True)\n", - "z = sp.Symbol(\"z\", real=True)\n", + "L = sp.Symbol(\"L\", integer=True, nonnegative=True)\n", + "z = sp.Symbol(\"z\", nonnegative=True, real=True)\n", "ff2 = BlattWeisskopfSquared(z, L)" ] }, @@ -183,9 +183,15 @@ }, "outputs": [], "source": [ + "from ampform.dynamics.form_factor import SphericalHankel1\n", "from ampform.io import aslatex\n", "\n", - "Math(aslatex({ff2: ff2.doit()}))" + "ell = sp.Symbol(R\"\\ell\", integer=True, nonnegative=True)\n", + "exprs = [\n", + " BlattWeisskopfSquared(z, L),\n", + " SphericalHankel1(ell, z),\n", + "]\n", + "Math(aslatex({e: e.doit(deep=False) for e in exprs}))" ] }, { @@ -205,7 +211,7 @@ "source": [ "from ampform.dynamics import BreakupMomentumSquared\n", "\n", - "m, m_a, m_b, d = sp.symbols(\"m, m_a, m_b, d\")\n", + "m, m_a, m_b, d = sp.symbols(\"m, m_a, m_b, d\", nonnegative=True)\n", "s = m**2\n", "q_squared = BreakupMomentumSquared(s, m_a, m_b)\n", "ff2 = BlattWeisskopfSquared(q_squared * d**2, angular_momentum=L)" @@ -243,7 +249,7 @@ "source": [ "plot_domain = np.linspace(0.01, 4, 500)\n", "sliders.set_ranges(\n", - " L=(0, 8),\n", + " L=(0, 10),\n", " m_a=(0, 2, 200),\n", " m_b=(0, 2, 200),\n", " d=(0, 5),\n", @@ -359,7 +365,7 @@ "source": [ "from ampform.dynamics import relativistic_breit_wigner\n", "\n", - "m, m0, w0 = sp.symbols(\"m, m0, Gamma0\")\n", + "m, m0, w0 = sp.symbols(\"m, m0, Gamma0\", nonnegative=True)\n", "rel_bw = relativistic_breit_wigner(s=m**2, mass0=m0, gamma0=w0)\n", "rel_bw" ] @@ -537,7 +543,7 @@ "sliders.set_ranges(\n", " m0=(0, 4, 200),\n", " Gamma0=(0, 1, 100),\n", - " L=(0, 8),\n", + " L=(0, 10),\n", " m_a=(0, 2, 200),\n", " m_b=(0, 2, 200),\n", " d=(0, 5),\n", diff --git a/src/ampform/dynamics/form_factor.py b/src/ampform/dynamics/form_factor.py index d17af8a8d..e54b0b6ce 100644 --- a/src/ampform/dynamics/form_factor.py +++ b/src/ampform/dynamics/form_factor.py @@ -2,6 +2,7 @@ from __future__ import annotations +from functools import lru_cache from typing import Any import sympy as sp @@ -37,15 +38,19 @@ class BlattWeisskopfSquared(sp.Expr): _latex_repr_ = R"B_{{{angular_momentum}}}^2\left({z}\right)" def evaluate(self) -> sp.Expr: - z, angular_momentum = self.args - return ( - sp.Abs(SphericalHankel1(angular_momentum, 1)) ** 2 - / sp.Abs(SphericalHankel1(angular_momentum, sp.sqrt(z))) ** 2 + ell = self.angular_momentum + z = sp.Dummy("z", nonnegative=True, real=True) + expr = ( + sp.Abs(SphericalHankel1(ell, 1)) ** 2 + / sp.Abs(SphericalHankel1(ell, sp.sqrt(z))) ** 2 / z ) + if not ell.free_symbols: + expr = expr.doit().simplify() + return expr.xreplace({z: self.z}) -@unevaluated(implement_doit=False) +@unevaluated class SphericalHankel1(sp.Expr): r"""Spherical Hankel function of the first kind for real-valued :math:`z`. @@ -63,22 +68,35 @@ class SphericalHankel1(sp.Expr): z: Any _latex_repr_ = R"h_{{{l}}}^{{(1)}}\left({z}\right)" - def doit(self, deep: bool = True, **kwargs): - expr = self.evaluate() - if deep and isinstance(self.l, sp.Integer): - return expr.doit() - return expr - def evaluate(self) -> sp.Expr: l, z = self.args # noqa: E741 k = sp.Dummy("k", integer=True, nonnegative=True) return ( (-sp.I) ** (1 + l) # type:ignore[operator] * (sp.exp(z * sp.I) / z) - * sp.Sum( + * _SymbolicSum( sp.factorial(l + k) / (sp.factorial(l - k) * sp.factorial(k)) * (sp.I / (2 * z)) ** k, # type:ignore[operator] (k, 0, l), ) ) + + +class _SymbolicSum(sp.Sum): + """See [TR-029](https://compwa.github.io/report/029.html) for why this class is needed.""" + + def doit(self, deep: bool = True, **kwargs) -> sp.Expr: + if _get_indices(self): + expression = self.args[0] + indices = self.args[1:] + return _SymbolicSum(expression.doit(deep=deep, **kwargs), *indices) + return super().doit(deep=deep, **kwargs) + + +@lru_cache(maxsize=None) +def _get_indices(expr: sp.Sum) -> set[sp.Basic]: + free_symbols = set() + for index in expr.args[1:]: + free_symbols.update(index.free_symbols) + return {s for s in free_symbols if not isinstance(s, sp.Dummy)} diff --git a/tests/dynamics/test_dynamics.py b/tests/dynamics/test_dynamics.py index 3ad8f561e..3473c6f33 100644 --- a/tests/dynamics/test_dynamics.py +++ b/tests/dynamics/test_dynamics.py @@ -21,7 +21,7 @@ class TestBlattWeisskopfSquared: - def test_max_angular_momentum(self): + def test_factorials(self): z = sp.Symbol("z") angular_momentum = sp.Symbol("L", integer=True) form_factor = BlattWeisskopfSquared(z, angular_momentum) @@ -29,24 +29,6 @@ def test_max_angular_momentum(self): factor, z_power, _ = form_factor_9.args assert factor == 4392846440677 assert z_power == z**8 - assert BlattWeisskopfSquared.max_angular_momentum is None - BlattWeisskopfSquared.max_angular_momentum = 1 - assert form_factor.evaluate() == sp.Piecewise( - (1, sp.Eq(angular_momentum, 0)), - (2 * z / (z + 1), sp.Eq(angular_momentum, 1)), - ) - BlattWeisskopfSquared.max_angular_momentum = None - - def test_unevaluated_expression(self): - z = sp.Symbol("z") - ff1 = BlattWeisskopfSquared(z, angular_momentum=1) - ff2 = BlattWeisskopfSquared(z, angular_momentum=2) - assert ff1.max_angular_momentum is None - assert ff2.max_angular_momentum is None - BlattWeisskopfSquared.max_angular_momentum = 3 - assert ff1.max_angular_momentum is 3 # noqa: F632 - assert ff2.max_angular_momentum is 3 # noqa: F632 - BlattWeisskopfSquared.max_angular_momentum = None class TestEnergyDependentWidth: