Formulate amplitude model

Warning

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

Generate transitions

In Generate transitions, we use generate_transitions() to create a list of allowed StateTransitionGraph instances for a specific decay channel:

import expertsystem as es

result = es.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["gamma", "pi0", "pi0"],
    allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"],
    allowed_interaction_types=["strong", "EM"],
    formalism_type="canonical-helicity",
)
import graphviz

dot = es.io.asdot(result, collapse_graphs=True)
graphviz.Source(dot)
../_images/amplitude_6_0.svg

Build SymPy expression

We can now use the Result to formulate an amplitude model. The type of this amplitude model is dependent on the formalism_type. The function amplitude.get_builder() helps to get the correct amplitude builder class for this formalism_type:

model_builder = es.amplitude.get_builder(result)
type(model_builder)
expertsystem.amplitude.helicity.CanonicalAmplitudeBuilder

If we now use the HelicityAmplitudeBuilder.generate() method of this builder, we get an amplitude model without any dynamics:

model_no_dynamics = model_builder.generate()

The HelicityModel.expression is just a sympy.Expr, which we can pull apart by using its args (see Advanced Expression Manipulation):

intensities = model_no_dynamics.expression.args
intensity_1 = intensities[0]
base, power = intensity_1.args
abs_arg = base.args[0]
amplitudes = abs_arg.args
amplitudes[0]
\[\displaystyle C[J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1};f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}] D^{0}_{0,0}\left(- \phi_{1,1+2},\theta_{1,1+2},0\right) D^{1}_{1,1}\left(- \phi_{1+2},\theta_{1+2},0\right) {C^('0', '0')_('0', '0', '0', '0')}^{2} {C^('1', '-1')_('0', '0', '1', '-1')} {C^('1', '-1')_('1', '-1', '0', '0')}\]

To set dynamics for specific resonances, use set_dynamics() on the same HelicityAmplitudeBuilder instance. You can set the dynamics to be any kind of Expr, as long as you keep track of which Symbol names you use. The expertsystem does provide a few common lineshape functions however, which can be constructed as Expr with the correct Symbol names using set_dynamics(). This function takes specific builder functions, such as create_relativistic_breit_wigner(), which would create a relativistic_breit_wigner() for a specific resonance. Here’s an example for a relativistic Breit-Wigner with form factor for the intermediate resonances and use a Blatt-Weisskopf barrier factor for the production decay:

from expertsystem.amplitude.dynamics.builder import (
    create_non_dynamic_with_ff,
    create_relativistic_breit_wigner_with_ff,
)

initial_state_particle = result.get_initial_state()[0]
model_builder.set_dynamics(initial_state_particle, create_non_dynamic_with_ff)
for name in result.get_intermediate_particles().names:
    model_builder.set_dynamics(name, create_relativistic_breit_wigner_with_ff)

See also

Custom dynamics

And we use the reconfigured HelicityAmplitudeBuilder to generate another HelicityModel, this time with relativistic Breit-Wigner functions and form factors:

model = model_builder.generate()

Visualize

Mathematical formula

It’s possible to view the complete amplitude model as an expression. This would, however, clog the screen here, so we instead just view the formula of one of its components:

some_amplitude = list(model.components.values())[0]
some_amplitude.doit()
\[\displaystyle \frac{\sqrt{10} C[J/\psi(1S) \to f_{0}(980)_{0} \gamma_{+1};f_{0}(980) \to \pi^{0}_{0} \pi^{0}_{0}] \Gamma_{f(0)(980)} m_{f(0)(980)} \left(\frac{1}{2} - \frac{\cos{\left(\theta_{1+2} \right)}}{2}\right) e^{- i \phi_{1+2}}}{10 \left(- \frac{i \Gamma_{f(0)(980)} m_{f(0)(980)}^{2} \sqrt{\frac{\left(m_{12}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{12}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{12}^{2}}}}{m_{12} \sqrt{\frac{\left(m_{f(0)(980)}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{f(0)(980)}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{f(0)(980)}^{2}}}} - m_{12}^{2} + m_{f(0)(980)}^{2}\right)}\]

Note

We use doit() to evaluate the Wigner-\(D\) (Rotation.D) and Clebsch-Gordan (CG) functions in the full expression.

The HelicityModel.parameter_defaults attribute can be used to substitute all parameters with suggested values:

some_amplitude.doit().subs(model.parameter_defaults)
\[\displaystyle \frac{0.00594 \sqrt{10} \left(\frac{1}{2} - \frac{\cos{\left(\theta_{1+2} \right)}}{2}\right) e^{- i \phi_{1+2}}}{- m_{12}^{2} + 0.9801 - \frac{0.05821794 i \sqrt{\frac{\left(m_{12}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{12}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{12}^{2}}}}{m_{12} \sqrt{\left(0.9801 - \left(m_{1} - m_{2}\right)^{2}\right) \left(0.9801 - \left(m_{1} + m_{2}\right)^{2}\right)}}}\]

Plotting

In this case (\(J/\psi \to \gamma f_0, f_0 \to \pi^0\pi^0\)) without dynamics, the total intensity is only dependent on the \(\theta\) angle of \(\gamma\) in the center of mass frame (see get_helicity_angle_label()):

no_dynamics = model_no_dynamics.expression.doit()
no_dynamics_substituted = no_dynamics.subs(model.parameter_defaults)
no_dynamics_substituted
\[\displaystyle 0.8 \sqrt{10} \cos^{2}{\left(\theta_{1+2} \right)} + 4.4 \cos^{2}{\left(\theta_{1+2} \right)} + 0.8 \sqrt{10} + 4.4\]
assert len(no_dynamics_substituted.free_symbols) == 1
import sympy as sy

theta = next(iter(no_dynamics_substituted.free_symbols))
sy.plot(
    no_dynamics_substituted,
    (theta, 0, sy.pi),
    axis_center=(0, 0),
    ylabel="$I$",
    ylim=(0, 16),
);
../_images/amplitude_30_0.svg

For this decay channel, the amplitude model is built up of four components:

no_dynamics.subs(zip(no_dynamics.args, sy.symbols("I_{:4}")))
\[\displaystyle I_{0} + I_{1} + I_{2} + I_{3}\]

This can be nicely visualized as follows:

import sympy as sy

plots = []
colors = ["red", "blue", "green", "purple"]

total = 0
for i, intensity in enumerate(no_dynamics.args):
    total += intensity.subs(model.parameter_defaults).doit()
    plots.append(
        sy.plot(
            total,
            (theta, 0, sy.pi),
            axis_center=(0, 0),
            ylabel="$I$",
            ylim=(0, 16),
            line_color=colors[i],
            show=False,
            label=f"$I_{i}$",
            legend=True,
        )
    )
for i in range(1, 4):
    plots[0].extend(plots[i])
plots[0].show()
../_images/amplitude_34_0.svg

Plot the model

In the model with dynamics, we have several free symbols, such as the mass and width of the resonances. For the fitting package these will be considered parameters that are to be optimized and (kinematic) variables that represent the data set. Examples of parameters are mass (\(m_\text{particle}\)) and width (\(\Gamma_\text{particle}\)) of the resonances and certain amplitude coefficients (\(C\)). Examples of kinematic variables are the helicity angles \(\theta\) and \(\phi\) and the invariant masses (\(m_{ij...}\)).

sorted(model.expression.free_symbols, key=lambda s: s.name)
[C[J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1};f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}],
 C[J/\psi(1S) \to f_{0}(980)_{0} \gamma_{+1};f_{0}(980) \to \pi^{0}_{0} \pi^{0}_{0}],
 Gamma_f(0)(1500),
 Gamma_f(0)(980),
 d_f(0)(1500),
 d_f(0)(980),
 m_1,
 m_12,
 m_2,
 m_f(0)(1500),
 m_f(0)(980),
 phi_1+2,
 phi_1,1+2,
 theta_1+2,
 theta_1,1+2]

Let’s say we want to plot the amplitude model with respect to \(m_{3+4}\). We would have to substitute all other free symbols with some value.

First, we can use HelicityModel.parameter_defaults to substitute the parameters with suggested values:

suggested_expression = model.expression.subs(model.parameter_defaults)
suggested_expression.free_symbols
{m_1, m_12, m_2, phi_1+2, phi_1,1+2, theta_1+2, theta_1,1+2}

Ideally, we would now ‘integrate out’ the helicity angles. Here, we however just set these angles to \(0\), as computing the integral would take quite some time:

angle = 0
angle_substitutions = {
    s: angle
    for s in suggested_expression.free_symbols
    if s.name.startswith("phi") or s.name.startswith("theta")
}
evaluated_angle_intensity = suggested_expression.subs(angle_substitutions)
evaluated_angle_intensity.free_symbols
{m_1, m_12, m_2}

By now we are only left with the masses of the final state particles (\(m_1\) and \(m_2\)), since they appear as symbols in the relativistic_breit_wigner_with_ff(). Final state particles 3 and 4 are the \(\pi^0\)’s, so we can just substitute them with their masses.

from expertsystem.reaction.particle import load_pdg

pi0 = load_pdg()["pi0"]
plotted_intensity = evaluated_angle_intensity.doit().subs(
    {
        sy.Symbol("m_1", real=True): pi0.mass,
        sy.Symbol("m_2", real=True): pi0.mass,
    },
    simultaneous=True,
)

Tip

Use subs() with simultaneous=True, as that avoids a bug later on when plotting with matplotlib.pyplot.

That’s it! Now we are only left with the invariant mass \(m_{3+4}\) of the two pions:

assert len(plotted_intensity.free_symbols) == 1
m = next(iter(plotted_intensity.free_symbols))
m
\[\displaystyle m_{12}\]

…and we can plot it with with sympy.plot:

sy.plot(
    plotted_intensity,
    (m, 2 * pi0.mass, 2.5),
    axis_center=(2 * pi0.mass, 0),
    xlabel=fR"$m(\pi^{0}\pi^{0})$",
    ylabel="$I$",
    backend="matplotlib",
);
../_images/amplitude_49_0.svg

The expression itself looks like this (after some rounding of the float values in this expression:

from sympy import preorder_traversal

rounded_intensity = plotted_intensity
rounded_intensity = rounded_intensity.subs({sy.sqrt(10): sy.sqrt(10).n()})
for a in preorder_traversal(rounded_intensity):
    if isinstance(a, sy.Float):
        rounded_intensity = rounded_intensity.subs(a, round(a, 2))
rounded_intensity
\[\displaystyle 2 \left|{\frac{0.22}{- m_{12}^{2} + 2.27 - \frac{0.17 i \sqrt{m_{12}^{2} - 0.07}}{m_{12}}} + \frac{0.08}{- m_{12}^{2} + 0.98 - \frac{0.06 i \sqrt{m_{12}^{2} - 0.07}}{m_{12}}}}\right|^{2}\]