Skip to content

Commit

Permalink
ENH: formulate Blatt–Weisskopf with spherical Hankel (#418)
Browse files Browse the repository at this point in the history
* FEAT: simplify Hankel polynomial internally
* MAINT: create separate `form_factor` submodule
  • Loading branch information
redeboer authored May 16, 2024
1 parent df3710a commit c1480fe
Show file tree
Hide file tree
Showing 8 changed files with 158 additions and 154 deletions.
1 change: 1 addition & 0 deletions .cspell.json
Original file line number Diff line number Diff line change
Expand Up @@ -244,6 +244,7 @@
"hankel",
"helicities",
"helicity",
"Hippel",
"htmlcov",
"imap",
"ipynb",
Expand Down
24 changes: 20 additions & 4 deletions docs/_extend_docstrings.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand Down
23 changes: 14 additions & 9 deletions docs/bibliography.bib
Original file line number Diff line number Diff line change
Expand Up @@ -173,19 +173,24 @@ @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.},
year = {1984},
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}
}
22 changes: 14 additions & 8 deletions docs/usage/dynamics.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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)**):"
]
},
{
Expand All @@ -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)"
]
},
Expand All @@ -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}))"
]
},
{
Expand All @@ -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)"
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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"
]
Expand Down Expand Up @@ -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",
Expand Down
111 changes: 3 additions & 108 deletions src/ampform/dynamics/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
102 changes: 102 additions & 0 deletions src/ampform/dynamics/form_factor.py
Original file line number Diff line number Diff line change
@@ -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<compwa:report/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<compwa:report/029>`
for more info. `This page
<https://mathworld.wolfram.com/SphericalHankelFunctionoftheFirstKind.html>`_
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)}
Loading

0 comments on commit c1480fe

Please sign in to comment.