Formulate amplitude model¶
Generate transitions¶
In Generate transitions, we used generate_transitions()
to create a list of allowed StateTransitionGraph
instances for a specific decay channel:
import qrules as q
result = q.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 = q.io.asdot(result, collapse_graphs=True)
graphviz.Source(dot)
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 get_builder()
helps to get the correct amplitude builder class for this formalism_type
:
from ampform import get_builder
model_builder = get_builder(result)
type(model_builder)
ampform.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]
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. AmpForm 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 ampform.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
And we use the reconfigured HelicityAmplitudeBuilder
to generate another HelicityModel
, this time with relativistic Breit-Wigner functions and form factors:
model = model_builder.generate()
Export¶
There is no special export function to export an HelicityModel
. However, we can just use the built-in pickle
method to write the model to disk:
import pickle
with open("helicity_model.pickle", "wb") as stream:
pickle.dump(model, stream)
This model will be imported again in Inspect model interactively.
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()
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)
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
assert len(no_dynamics_substituted.free_symbols) == 1
import sympy as sp
theta = next(iter(no_dynamics_substituted.free_symbols))
sp.plot(
no_dynamics_substituted,
(theta, 0, sp.pi),
axis_center=(0, 0),
ylabel="$I$",
ylim=(0, 16),
);
For this decay channel, the amplitude model is built up of four components:
no_dynamics.subs(zip(no_dynamics.args, sp.symbols("I_{:4}")))
This can be nicely visualized as follows:
import sympy as sp
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(
sp.plot(
total,
(theta, 0, sp.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()
Plot the model¶
Tip
See Inspect model interactively for a much more didactic way to visualize 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 qrules import load_pdg
pi0 = load_pdg()["pi0"]
plotted_intensity = evaluated_angle_intensity.doit().subs(
{
sp.Symbol("m_1", real=True): pi0.mass,
sp.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
…and we can plot it with with sympy.plot
:
sp.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",
);
The expression itself looks like this (after some rounding of the float
values in this expression using tree traversal):
from sympy import preorder_traversal
def round_nested(expression: sp.Expr, n_decimals: int = 2) -> sp.Expr:
for node in preorder_traversal(expression):
if node.free_symbols:
continue
if isinstance(node, (float, sp.Float)):
expression = expression.subs(node, round(node, n_decimals))
if isinstance(node, sp.Pow) and node.args[1] == 1 / 2:
expression = expression.subs(node, round(node.n(), n_decimals))
return expression
rounded = round_nested(plotted_intensity)
round_nested(rounded)