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)
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]
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
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()
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 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),
);
For this decay channel, the amplitude model is built up of four components:
no_dynamics.subs(zip(no_dynamics.args, sy.symbols("I_{:4}")))
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()
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
…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",
);
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