Custom dynamics¶

import inspect
from typing import Dict, Tuple

import graphviz
import qrules as q
import sympy as sp
from ampform.dynamics.builder import (
    ResonanceDynamicsBuilder,
    TwoBodyKinematicVariableSet,
    create_relativistic_breit_wigner,
)
from qrules.particle import Particle

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-2-3575a74e0649> in <module>
      2 from typing import Dict, Tuple
      3 
----> 4 import graphviz
      5 import qrules as q
      6 import sympy as sp

ModuleNotFoundError: No module named 'graphviz'

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:

print(inspect.getsource(ResonanceDynamicsBuilder))
print(inspect.getsource(create_relativistic_breit_wigner))

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 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:

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

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

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));

AmpForm

Custom dynamics¶