Custom dynamics¶
import graphviz
import qrules as q
import sympy as sp
We start by generating allowed transitions for a simple decay channel, just like in Formulate amplitude model:
result = q.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 = q.io.asdot(result, collapse_graphs=True)
graphviz.Source(dot)
Next, create a HelicityAmplitudeBuilder
using get_builder()
:
from ampform import get_builder
model_builder = 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
:
import inspect
from ampform.dynamics.builder import (
ResonanceDynamicsBuilder,
create_relativistic_breit_wigner,
)
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, resonance: Particle, variable_pool: TwoBodyKinematicVariableSet
) -> BuilderReturnType:
...
def create_relativistic_breit_wigner(
resonance: Particle, variable_pool: TwoBodyKinematicVariableSet
) -> BuilderReturnType:
"""Create a `.relativistic_breit_wigner` for a two-body decay."""
inv_mass = variable_pool.in_edge_inv_mass
res_mass = sp.Symbol(f"m_{resonance.name}")
res_width = sp.Symbol(f"Gamma_{resonance.name}")
expression = relativistic_breit_wigner(
s=inv_mass ** 2,
mass0=res_mass,
gamma0=res_width,
)
parameter_defaults = {
res_mass: resonance.mass,
res_width: resonance.width,
}
return expression, parameter_defaults
A function that behaves like a ResonanceDynamicsBuilder
should return a tuple
of some Expr
(which formulates your lineshape) and a dict
of Symbol
s 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 ampform
:
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:
from typing import Dict, Tuple
from qrules.particle import Particle
from ampform.dynamics.builder import TwoBodyKinematicVariableSet
def create_my_dynamics(
resonance: Particle, variable_pool: TwoBodyKinematicVariableSet
) -> Tuple[sp.Expr, Dict[sp.Symbol, float]]:
res_mass = sp.Symbol(f"m_{resonance.name}")
res_width = sp.Symbol(f"sigma_{resonance.name}")
expression = my_dynamics(
x=variable_pool.in_edge_inv_mass,
mu=res_mass,
sigma=res_width,
)
parameter_defaults = {
res_mass: resonance.mass,
res_width: resonance.width,
}
return expression, parameter_defaults
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 Symbol
s 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));