diff --git a/.cspell.json b/.cspell.json index a49073f36..03f4f2aaf 100644 --- a/.cspell.json +++ b/.cspell.json @@ -244,6 +244,7 @@ "hankel", "helicities", "helicity", + "Hippel", "htmlcov", "imap", "ipynb", diff --git a/docs/_extend_docstrings.py b/docs/_extend_docstrings.py index a8a89dcf6..b87a48746 100644 --- a/docs/_extend_docstrings.py +++ b/docs/_extend_docstrings.py @@ -63,12 +63,19 @@ def extend_docstrings() -> None: def extend_BlattWeisskopfSquared() -> None: - from ampform.dynamics import BlattWeisskopfSquared + from ampform.dynamics.form_factor import BlattWeisskopfSquared, SphericalHankel1 - z = sp.Symbol("z", real=True) - L = sp.Symbol("L", integer=True) + z = sp.Symbol("z", nonnegative=True, real=True) + L = sp.Symbol("L", integer=True, nonnegative=True) expr = BlattWeisskopfSquared(z, angular_momentum=L) - _append_latex_doit_definition(expr, deep=True, full_width=True) + h1lz = SphericalHankel1(L, z) + _append_latex_doit_definition(expr) + _append_to_docstring( + BlattWeisskopfSquared, + f"""\n + where :math:`{sp.latex(h1lz)}` is defined by :eq:`SphericalHankel1`. + """, + ) def extend_BoostMatrix() -> None: @@ -472,6 +479,15 @@ def extend_RotationZMatrix() -> None: ) +def extend_SphericalHankel1() -> None: + from ampform.dynamics.form_factor import SphericalHankel1 + + z = sp.Symbol("z", nonnegative=True, real=True) + ell = sp.Symbol(R"\ell", integer=True, nonnegative=True) + expr = SphericalHankel1(ell, z) + _append_latex_doit_definition(expr) + + def extend_Theta() -> None: from ampform.kinematics.angles import Theta diff --git a/docs/bibliography.bib b/docs/bibliography.bib index ffff591d3..401058569 100755 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -173,15 +173,6 @@ @incollection{ParticleDataGroup:2020ssz url = {https://pdg.lbl.gov/2010/reviews/rpp2010-rev-dalitz-analysis-formalism.pdf} } -@phdthesis{pychyGekoppeltePartialwellenanalyseAnnihilationen2016, - title = {{Gekoppelte Partialwellenanalyse von 𝑝̅𝑝-Annihilationen im Fluge in die Endzustände 𝐾⁺𝐾⁻𝜋⁰, 𝜋⁰𝜋⁰𝜂 und 𝜋⁰𝜂𝜂}}, - author = {Pychy, Julian}, - year = {2016}, - month = apr, - url = {https://hss-opus.ub.ruhr-uni-bochum.de/opus4/frontdoor/deliver/index/docId/4943/file/diss.pdf}, - school = {Ruhr-Universität Bochum} -} - @misc{Richman:1984gh, title = {An {{Experimenter}}'s {{Guide}} to the {{Helicity Formalism}}}, author = {Richman, Jeffrey D.}, @@ -189,3 +180,17 @@ @misc{Richman:1984gh month = jun, url = {https://inspirehep.net/literature/202987} } + +@article{VonHippel:1972fg, + title = {Centrifugal-{{Barrier Effects}} in {{Resonance Partial Decay Widths}}, {{Shapes}}, and {{Production Amplitudes}}}, + author = {{von Hippel}, Frank and Quigg, C.}, + year = {1972}, + month = feb, + journal = {Physical Review D}, + volume = {5}, + number = {3}, + pages = {624--638}, + issn = {0556-2821}, + doi = {10.1103/PhysRevD.5.624}, + url = {https://link.aps.org/doi/10.1103/PhysRevD.5.624} +} 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/__init__.py b/src/ampform/dynamics/__init__.py index 42b0be98d..fcc3d40ea 100644 --- a/src/ampform/dynamics/__init__.py +++ b/src/ampform/dynamics/__init__.py @@ -3,14 +3,14 @@ .. seealso:: :doc:`/usage/dynamics` and :doc:`/usage/dynamics/analytic-continuation` """ -# cspell:ignore asner mhash from __future__ import annotations -from typing import TYPE_CHECKING, Any, ClassVar +from typing import TYPE_CHECKING, Any import sympy as sp # pyright: reportUnusedImport=false +from ampform.dynamics.form_factor import BlattWeisskopfSquared from ampform.dynamics.phasespace import ( BreakupMomentumSquared, EqualMassPhaseSpaceFactor, # noqa: F401 @@ -27,114 +27,9 @@ from sympy.printing.latex import LatexPrinter -@unevaluated -class BlattWeisskopfSquared(sp.Expr): - # cspell:ignore pychyGekoppeltePartialwellenanalyseAnnihilationen - r"""Blatt-Weisskopf function :math:`B_L^2(z)`, up to :math:`L \leq 8`. - - Args: - z: Argument of the Blatt-Weisskopf function :math:`B_L^2(z)`. A usual - choice is :math:`z = (d q)^2` with :math:`d` the impact parameter and - :math:`q` the breakup-momentum (see `.BreakupMomentumSquared`). - - angular_momentum: Angular momentum :math:`L` of the decaying particle. - - Note that equal powers of :math:`z` appear in the nominator and the denominator, - while some sources have nominator :math:`1`, instead of :math:`z^L`. Compare for - instance Equation (50.27) in :pdg-review:`2021; Resonances; p.9`. - - Each of these cases for :math:`L` has been taken from - :cite:`pychyGekoppeltePartialwellenanalyseAnnihilationen2016`, p.59, - :cite:`Chung:1995dx`, p.415, and :cite:`Chung:1995dx`. For a good overview of where - to use these Blatt-Weisskopf functions, see :cite:`ParticleDataGroup:2020ssz`. - - See also :ref:`usage/dynamics:Form factor`. - """ - - z: Any - angular_momentum: Any - _latex_repr_ = R"B_{{{angular_momentum}}}^2\left({z}\right)" - - max_angular_momentum: ClassVar[int | None] = None - """Limit the maximum allowed angular momentum :math:`L`. - - This improves performance when :math:`L` is a `~sympy.core.symbol.Symbol` and you - are note interested in higher angular momenta. - """ - - def evaluate(self) -> sp.Expr: - z: sp.Expr = self.args[0] # type: ignore[assignment] - angular_momentum: sp.Expr = self.args[1] # type: ignore[assignment] - cases: dict[int, sp.Expr] = { - 0: sp.S.One, - 1: 2 * z / (z + 1), - 2: 13 * z**2 / ((z - 3) * (z - 3) + 9 * z), - 3: 277 * z**3 / (z * (z - 15) * (z - 15) + 9 * (2 * z - 5) * (2 * z - 5)), - 4: ( - 12746 - * z**4 - / ( - (z**2 - 45 * z + 105) * (z**2 - 45 * z + 105) - + 25 * z * (2 * z - 21) * (2 * z - 21) - ) - ), - 5: ( - 998881 - * z**5 - / (z**5 + 15 * z**4 + 315 * z**3 + 6300 * z**2 + 99225 * z + 893025) - ), - 6: ( - 118394977 - * z**6 - / ( - z**6 - + 21 * z**5 - + 630 * z**4 - + 18900 * z**3 - + 496125 * z**2 - + 9823275 * z - + 108056025 - ) - ), - 7: ( - 19727003738 - * z**7 - / ( - z**7 - + 28 * z**6 - + 1134 * z**5 - + 47250 * z**4 - + 1819125 * z**3 - + 58939650 * z**2 - + 1404728325 * z - + 18261468225 - ) - ), - 8: ( - 4392846440677 - * z**8 - / ( - z**8 - + 36 * z**7 - + 1890 * z**6 - + 103950 * z**5 - + 5457375 * z**4 - + 255405150 * z**3 - + 9833098275 * z**2 - + 273922023375 * z - + 4108830350625 - ) - ), - } - return sp.Piecewise(*[ - (expression, sp.Eq(angular_momentum, value)) - for value, expression in cases.items() - if self.max_angular_momentum is None or value <= self.max_angular_momentum - ]) - - @unevaluated class EnergyDependentWidth(sp.Expr): + # cspell:ignore asner r"""Mass-dependent width, coupled to the pole position of the resonance. See Equation (50.28) in :pdg-review:`2021; Resonances; p.9` and diff --git a/src/ampform/dynamics/form_factor.py b/src/ampform/dynamics/form_factor.py new file mode 100644 index 000000000..e54b0b6ce --- /dev/null +++ b/src/ampform/dynamics/form_factor.py @@ -0,0 +1,102 @@ +"""Implementations of the form factor, or barrier factor.""" + +from __future__ import annotations + +from functools import lru_cache +from typing import Any + +import sympy as sp + +from ampform.sympy import unevaluated + + +@unevaluated +class BlattWeisskopfSquared(sp.Expr): + r"""Normalized Blatt-Weisskopf function :math:`B_L^2(z)`, with :math:`B_L^2(1)=1`. + + Args: + z: Argument of the Blatt-Weisskopf function :math:`B_L^2(z)`. A usual + choice is :math:`z = (d q)^2` with :math:`d` the impact parameter and + :math:`q` the breakup-momentum (see `.BreakupMomentumSquared`). + + angular_momentum: Angular momentum :math:`L` of the decaying particle. + + Note that equal powers of :math:`z` appear in the nominator and the denominator, + while some sources define an *non-normalized* form factor :math:`F_L` with :math:`1` + in the nominator, instead of :math:`z^L`. See for instance Equation (50.27) in + :pdg-review:`2021; Resonances; p.9`. We normalize the form factor such that + :math:`B_L^2(1)=1` and that :math:`B_L^2` is unitless no matter what :math:`z` is. + + .. seealso:: :ref:`usage/dynamics:Form factor`, :doc:`TR-029`, + and :cite:`chungFormulasAngularMomentumBarrier2015`. + + With this, the implementation becomes + """ + + z: Any + angular_momentum: Any + _latex_repr_ = R"B_{{{angular_momentum}}}^2\left({z}\right)" + + def evaluate(self) -> sp.Expr: + 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 +class SphericalHankel1(sp.Expr): + r"""Spherical Hankel function of the first kind for real-valued :math:`z`. + + See :cite:`VonHippel:1972fg`, Equation (A12), and :doc:`TR-029` + for more info. `This page + `_ + explains the difference with the *general* Hankel function of the first kind, + :math:`H_\ell^{(1)}`. + + This expression class assumes that :math:`z` is real and evaluates to the following + series: + """ + + l: Any # noqa: E741 + z: Any + _latex_repr_ = R"h_{{{l}}}^{{(1)}}\left({z}\right)" + + 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) + * _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_deprecated.py b/tests/dynamics/test_deprecated.py index 93dc5a2ac..0315b30f6 100644 --- a/tests/dynamics/test_deprecated.py +++ b/tests/dynamics/test_deprecated.py @@ -2,11 +2,8 @@ import sympy as sp -from ampform.dynamics import ( - BlattWeisskopfSquared, - EnergyDependentWidth, - EqualMassPhaseSpaceFactor, -) +from ampform.dynamics import EnergyDependentWidth, EqualMassPhaseSpaceFactor +from ampform.dynamics.form_factor import BlattWeisskopfSquared from ampform.sympy import UnevaluatedExpression diff --git a/tests/dynamics/test_dynamics.py b/tests/dynamics/test_dynamics.py index fff816996..3473c6f33 100644 --- a/tests/dynamics/test_dynamics.py +++ b/tests/dynamics/test_dynamics.py @@ -6,13 +6,13 @@ import sympy as sp from ampform.dynamics import ( - BlattWeisskopfSquared, EnergyDependentWidth, EqualMassPhaseSpaceFactor, PhaseSpaceFactor, PhaseSpaceFactorSWave, relativistic_breit_wigner_with_ff, ) +from ampform.dynamics.form_factor import BlattWeisskopfSquared if TYPE_CHECKING: from qrules import ParticleCollection @@ -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: