Custom dynamics


The PWA Expert System has been split up into QRules and AmpForm. Please use these packages instead!

import inspect
from typing import Dict, Tuple

import graphviz
import sympy as sp

import expertsystem as es
from expertsystem.amplitude.dynamics.builder import (
from expertsystem.reaction.particle import Particle

We start by generating allowed transitions for a simple decay channel, just like in Formulate amplitude model:

result = es.generate_transitions(
    initial_state=("J/psi(1S)", [+1]),
    final_state=[("gamma", [+1]), "pi0", "pi0"],
    allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"],
    allowed_interaction_types=["strong", "EM"],
dot =, collapse_graphs=True)

Next, create a HelicityAmplitudeBuilder using get_builder():

model_builder = es.amplitude.get_builder(result)

In Formulate amplitude model, we used set_dynamics() with some standard lineshape builders from the builder module. These builders have a signature that follows the ResonanceDynamicsBuilder Protocol:

class ResonanceDynamicsBuilder(Protocol):
    """Protocol that is used by `.set_dynamics`.

    Follow this `~typing.Protocol` when defining a builder function that is to
    be used by `.set_dynamics`. For an example, see the source code
    `.create_relativistic_breit_wigner`, which creates a

    .. seealso:: :doc:`/usage/dynamics/custom`

    def __call__(
        self, particle: Particle, variable_pool: TwoBodyKinematicVariableSet
    ) -> Tuple[sp.Expr, Dict[sp.Symbol, float]]:

def create_relativistic_breit_wigner(
    particle: Particle, variable_pool: TwoBodyKinematicVariableSet
) -> Tuple[sp.Expr, Dict[sp.Symbol, float]]:
    inv_mass = variable_pool.in_edge_inv_mass
    res_mass = sp.Symbol(f"m_{}")
    res_width = sp.Symbol(f"Gamma_{}")
    return (
        {res_mass: particle.mass, res_width: particle.width},

A function that behaves like a ResonanceDynamicsBuilder should return a tuple of some Expr (which formulates your lineshape) and a dict of Symbols to some suggested initial values. This signature is required so that set_dynamics() knows how to extract the correct symbol names and their suggested initial values from a StateTransitionGraph.

The Expr you use for the lineshape can be anything. Here, we use a Gaussian function and wrap it in a function. As you can see, this function stands on its own, independent of the expertsystem:

def my_dynamics(x: sp.Symbol, mu: sp.Symbol, sigma: sp.Symbol) -> sp.Expr:
    return sp.exp(-((x - mu) ** 2) / sigma ** 2 / 2) / (
        sigma * sp.sqrt(2 * sp.pi)
x, mu, sigma = sp.symbols("x mu sigma")
sp.plot(my_dynamics(x, 0, 1), (x, -3, 3), axis_center=(0, 0))
my_dynamics(x, mu, sigma)
\[\displaystyle \frac{\sqrt{2} e^{- \frac{\left(- \mu + x\right)^{2}}{2 \sigma^{2}}}}{2 \sqrt{\pi} \sigma}\]

We can now follow the example of the create_relativistic_breit_wigner() to create a builder for this custom lineshape:

def create_my_dynamics(
    particle: Particle, variable_pool: TwoBodyKinematicVariableSet
) -> Tuple[sp.Expr, Dict[sp.Symbol, float]]:
    res_mass = sp.Symbol(f"m_{}")
    res_width = sp.Symbol(f"sigma_{}")
    return (
        {res_mass: particle.mass, res_width: particle.width},

Now, just like in Build SymPy expression, it’s simply a matter of plugging this builder into set_dynamics() and we can generate() a model with this custom lineshape:

for name in result.get_intermediate_particles().names:
    model_builder.set_dynamics(name, create_my_dynamics)
model = model_builder.generate()

As can be seen, the HelicityModel.parameter_defaults section has been updated with the some additional parameters for the custom parameter and there corresponding suggested initial values:

{m_f(0)(980): 0.99,
 sigma_f(0)(980): 0.06,
 C[J/\psi(1S) \to f_{0}(980)_{0} \gamma_{+1};f_{0}(980) \to \pi^{0}_{0} \pi^{0}_{0}]: (1+0j),
 m_f(0)(1500): 1.506,
 sigma_f(0)(1500): 0.112,
 C[J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1};f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}]: (1+0j)}

Let’s quickly have a look what this lineshape looks like. First, check which Symbols remain once we replace the parameters with their suggested initial values. These are the kinematic variables of the model:

expr = model.expression.doit().subs(model.parameter_defaults)
free_symbols = tuple(sorted(expr.free_symbols, key=lambda s:
(m_12, theta_1+2)

To create an invariant mass distribution, we should integrate out the \(\theta\) angle. This can be done with integrate():

m, theta = free_symbols
integrated_expr = sp.integrate(
    (theta, 0, sp.pi),
\[\displaystyle 9.0 \cdot 10^{1} e^{- 277.777777777778 \left(m_{12} - 0.99\right)^{2}} + 3.0 \cdot 10^{1} e^{- 180.806441326531 \left(0.664010624169987 m_{12} - 1\right)^{2}} + 1.0 \cdot 10^{2} e^{- 90.4032206632653 \left(0.664010624169987 m_{12} - 1\right)^{2}} e^{- 138.888888888889 \left(m_{12} - 0.99\right)^{2}}\]

Finally, here is the resulting expression as a function of the invariant mass, with custom dynamics!

x1, x2 = 0.6, 1.9
sp.plot(integrated_expr, (m, x1, x2), axis_center=(x1, 0));