BREAK: replace EnergyDependentWidth with AmpForm's implementation (#190)

* FEAT: allow selecting phase space factor protocol

Author: redeboer

docs/conf.py

@@ -61,6 +61,7 @@
     "sp.Basic": "sympy.core.basic.Basic",
     "sp.Expr": "sympy.core.expr.Expr",
     "sp.Indexed": "sympy.tensor.indexed.Indexed",
+    "PhaseSpaceFactorProtocol": "ampform.dynamics.phasespace.PhaseSpaceFactorProtocol",
     "sp.Rational": "sympy.core.numbers.Rational",
     "sp.Symbol": "sympy.core.symbol.Symbol",
     "StateID": ("obj", "ampform_dpd.decay.StateID"),

docs/index.ipynb

@@ -472,12 +472,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from ampform.dynamics.phasespace import PhaseSpaceFactor\n",
+    "\n",
     "from ampform_dpd.dynamics.builder import BreitWignerBuilder\n",
     "\n",
     "dynamics_builder = BreitWignerBuilder(\n",
     "    energy_dependent_width=True,\n",
     "    production_form_factor=True,\n",
     "    decay_form_factor=True,\n",
+    "    phsp_factor=PhaseSpaceFactor,\n",
     ")\n",
     "for chain in model_builder.decay.chains:\n",
     "    model_builder.dynamics_choices.register_builder(chain, dynamics_builder)\n",

docs/serialization.ipynb

@@ -46,6 +46,7 @@
     "import matplotlib.pyplot as plt\n",
     "import pandas as pd\n",
     "import sympy as sp\n",
+    "from ampform.dynamics import EnergyDependentWidth\n",
     "from ampform.dynamics.form_factor import BlattWeisskopfSquared, FormFactor\n",
     "from ampform.dynamics.phasespace import BreakupMomentum, BreakupMomentumSquared\n",
     "from IPython.display import JSON, Math\n",
@@ -57,7 +58,6 @@
     "    BreitWigner,\n",
     "    BuggBreitWigner,\n",
     "    ChannelArguments,\n",
-    "    EnergyDependentWidth,\n",
     "    MultichannelBreitWigner,\n",
     "    SimpleBreitWigner,\n",
     ")\n",

src/ampform_dpd/dynamics/__init__.py

@@ -6,9 +6,14 @@
 from typing import TYPE_CHECKING, Any
 
 import sympy as sp
+from ampform.dynamics import EnergyDependentWidth
 from ampform.dynamics.form_factor import FormFactor
-from ampform.dynamics.phasespace import BreakupMomentum
-from ampform.sympy import unevaluated
+from ampform.dynamics.phasespace import (
+    BreakupMomentum,
+    PhaseSpaceFactor,
+    PhaseSpaceFactorProtocol,
+)
+from ampform.sympy import argument, unevaluated
 
 if TYPE_CHECKING:
     from sympy.printing.latex import LatexPrinter
@@ -23,13 +28,14 @@ class RelativisticBreitWigner(sp.Expr):
     m2: Any
     angular_momentum: Any
     meson_radius: Any
+    phsp_factor: PhaseSpaceFactorProtocol = argument(  # type:ignore[assignment]
+        default=PhaseSpaceFactor, sympify=False
+    )
     _latex_repr_ = (
         R"\mathcal{{R}}_{{{angular_momentum}}}\left({s}, {mass0}, {gamma0}\right)"
     )
 
     def evaluate(self):
-        from ampform.dynamics import EnergyDependentWidth  # noqa: PLC0415
-
         s, m0, w0, m1, m2, angular_momentum, meson_radius = self.args
         width = EnergyDependentWidth(
             s=s,
@@ -39,6 +45,7 @@ def evaluate(self):
             m_b=m2,
             angular_momentum=angular_momentum,
             meson_radius=meson_radius,
+            phsp_factor=self.phsp_factor,
             name=Rf"\Gamma_{{{sp.latex(angular_momentum)}}}",
         )
         return (m0 * w0) / (m0**2 - s - width * m0 * sp.I)
@@ -57,6 +64,9 @@ class BreitWignerMinL(sp.Expr):
     l_prod: Any
     R_dec: Any
     R_prod: Any
+    phsp_factor: PhaseSpaceFactorProtocol = argument(  # type:ignore[assignment]
+        default=PhaseSpaceFactor, sympify=False
+    )
     _latex_repr_ = R"\mathcal{{R}}^\mathrm{{BW}}_{{{l_dec},{l_prod}}}\left({s}\right)"
 
     def evaluate(self):  # noqa: PLR0914
@@ -65,7 +75,7 @@ def evaluate(self):  # noqa: PLR0914
         ff0_prod = FormFactor(m_top**2, m0, m_spec, l_prod, R_prod)
         ff_dec = FormFactor(s, m1, m2, l_dec, R_dec)
         ff0_dec = FormFactor(m0**2, m1, m2, l_dec, R_dec)
-        width = EnergyDependentWidth(s, m0, Γ0, m1, m2, l_dec, R_dec)
+        width = EnergyDependentWidth(s, m0, Γ0, m1, m2, l_dec, R_dec, self.phsp_factor)
         return sp.Mul(
             ff_prod / ff0_prod,
             1 / (m0**2 - s - sp.I * m0 * width),
@@ -117,37 +127,6 @@ def evaluate(self):
         return 1 / (m0**2 - s - sp.I * m0 * Γ)
 
 
-@unevaluated
-class EnergyDependentWidth(sp.Expr):
-    s: Any
-    m0: Any
-    Γ0: Any
-    m1: Any
-    m2: Any
-    L: Any
-    R: Any
-
-    def evaluate(self):
-        s, m0, Γ0, m1, m2, L, R = self.args
-        ff = FormFactor(s, m1, m2, L, R) ** 2
-        ff0 = FormFactor(m0**2, m1, m2, L, R) ** 2
-        return sp.Mul(
-            Γ0,
-            m0 / sp.sqrt(s),
-            ff / ff0,
-            evaluate=False,
-        )
-
-    def _latex_repr_(self, printer: LatexPrinter) -> str:
-        s = printer._print(self.s)  # pyright:ignore[reportPrivateUsage]
-        if self.L == 0:
-            if self.m1 == 0 and self.m2 == 0:
-                return printer._print(self.Γ0)  # pyright:ignore[reportPrivateUsage]
-            return Rf"\Gamma\left({s}\right)"
-        L = printer._print(self.L)  # pyright:ignore[reportPrivateUsage]
-        return Rf"\Gamma_{{{L}}}\left({s}\right)"
-
-
 @unevaluated
 class MultichannelBreitWigner(sp.Expr):
     s: Any
@@ -194,6 +173,9 @@ class BreitWigner(sp.Expr):
     m2: Any = 0
     angular_momentum: Any = 0
     meson_radius: Any = 1
+    phsp_factor: PhaseSpaceFactorProtocol = argument(  # type:ignore[assignment]
+        default=PhaseSpaceFactor, sympify=False
+    )
 
     def evaluate(self):
         width = self.energy_dependent_width()
@@ -206,7 +188,7 @@ def energy_dependent_width(self) -> sp.Expr:
         s, m0, Γ0, m1, m2, L, d = self.args
         if L == 0 and m1 == 0 and m2 == 0:
             return Γ0  # type:ignore[return-value]
-        return EnergyDependentWidth(s, m0, Γ0, m1, m2, L, d)
+        return EnergyDependentWidth(s, m0, Γ0, m1, m2, L, d, self.phsp_factor)
 
     def _latex_repr_(self, printer: LatexPrinter, *args) -> str:
         s = printer._print(self.s)

src/ampform_dpd/dynamics/builder.py

@@ -10,6 +10,7 @@
 
 import sympy as sp
 from ampform.dynamics.form_factor import FormFactor
+from ampform.dynamics.phasespace import PhaseSpaceFactor, PhaseSpaceFactorProtocol
 from attrs import define
 
 from ampform_dpd import DefinedExpression, create_mass_symbol, to_particle
@@ -22,13 +23,14 @@ class BreitWignerBuilder:
     energy_dependent_width: bool = True
     decay_form_factor: bool = True
     production_form_factor: bool = True
+    phsp_factor: PhaseSpaceFactorProtocol = PhaseSpaceFactor  # type:ignore[assignment]
 
     def __call__(self, decay_chain: ThreeBodyDecayChain) -> DefinedExpression:
         """Formulate a (relativistic) Breit-Wigner for this resonance."""
         decay_node = decay_chain.decay_node
         s = get_mandelstam_s(decay_node)
         if self.energy_dependent_width:
-            expression = _create_breit_wigner(s, decay_node)
+            expression = _create_breit_wigner(s, decay_node, self.phsp_factor)
         else:
             expression = _create_simple_breit_wigner(s, decay_node)
         if self.decay_form_factor:
@@ -66,7 +68,9 @@ def _create_form_factor(s: sp.Symbol, isobar: IsobarNode) -> DefinedExpression:
     return DefinedExpression(form_factor, parameter_defaults)
 
 
-def _create_breit_wigner(s: sp.Symbol, isobar: DecayNode) -> DefinedExpression:
+def _create_breit_wigner(
+    s: sp.Symbol, isobar: DecayNode, phsp_factor: PhaseSpaceFactorProtocol
+) -> DefinedExpression:
     outgoing_state_mass1 = create_mass_symbol(isobar.child1)
     outgoing_state_mass2 = create_mass_symbol(isobar.child2)
     angular_momentum = _get_angular_momentum(isobar)
@@ -81,6 +85,7 @@ def _create_breit_wigner(s: sp.Symbol, isobar: DecayNode) -> DefinedExpression:
         m2=outgoing_state_mass2,
         angular_momentum=angular_momentum,
         meson_radius=meson_radius,
+        phsp_factor=phsp_factor,
     )
     parameter_defaults: dict[sp.Symbol, complex | float] = {
         res_mass: isobar.parent.mass,