Custom dynamics¶
Warning
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 (
ResonanceDynamicsBuilder,
TwoBodyKinematicVariableSet,
create_relativistic_breit_wigner,
)
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"],
formalism_type="canonical-helicity",
)
dot = es.io.asdot(result, collapse_graphs=True)
graphviz.Source(dot)
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:
print(inspect.getsource(ResonanceDynamicsBuilder))
print(inspect.getsource(create_relativistic_breit_wigner))
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
`.relativistic_breit_wigner`.
.. 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_{particle.name}")
res_width = sp.Symbol(f"Gamma_{particle.name}")
return (
relativistic_breit_wigner(
inv_mass,
res_mass,
res_width,
),
{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)
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_{particle.name}")
res_width = sp.Symbol(f"sigma_{particle.name}")
return (
my_dynamics(
x=variable_pool.in_edge_inv_mass,
mu=res_mass,
sigma=res_width,
),
{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:
model.parameter_defaults
{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: s.name))
free_symbols
(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(
expr,
(theta, 0, sp.pi),
meijerg=True,
conds="piecewise",
risch=None,
heurisch=None,
manual=None,
)
integrated_expr.n(1)
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));
