Formulate amplitude model#
Show code cell content
import sympy as sp
from IPython.display import Math
from ampform.io import aslatex, improve_latex_rendering
improve_latex_rendering()
Generate transitions#
In Generate transitions, we used generate_transitions()
to create a list of allowed StateTransition
s for a specific decay channel:
import qrules
reaction = qrules.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="canonical-helicity",
)
Show code cell source
import graphviz
dot = qrules.io.asdot(reaction, collapse_graphs=True)
graphviz.Source(dot)
Build model#
We can now use the ReactionInfo
to formulate an amplitude model. The type of this amplitude model is dependent on the formalism
. The function get_builder()
helps to get the correct amplitude builder class for this formalism
:
from ampform import get_builder
model_builder = get_builder(reaction)
type(model_builder)
ampform.helicity.CanonicalAmplitudeBuilder
If we now use the HelicityAmplitudeBuilder.formulate()
method of this builder, we get a HelicityModel
without any dynamics:
See also
Main expressions#
A HelicityModel
has a few attributes. The expression for the total intensity is given by intensity
:
This shows that the main intensity is an incoherent sum of the amplitude for each spin projection combination of the initial and final states. The expressions for each of these amplitudes are provided with the amplitudes
attribute. This is an OrderedDict
, so we can inspect the first of these amplitudes as follows:
symbol, expression = list(model_no_dynamics.amplitudes.items())[0]
The intensity and its amplitudes are recombined through the expression
attribute. This is just a sympy.Expr
, which we can pull apart by using its args
(see Advanced Expression Manipulation). Here’s an example:
intensities = model_no_dynamics.expression.args
intensity_1 = intensities[0]
base, power = intensity_1.args
abs_arg = base.args[0]
amplitude_sum = abs_arg.args
some_amplitude = amplitude_sum[0]
some_amplitude
some_amplitude.doit()
Parameters and kinematic variables#
As can be seen, the expression contains several Symbol
s. Some of these represent (kinematic) variables, such as the helicity angles \(\phi_0\) and \(\theta_0\) (see get_helicity_angle_symbols()
for the meaning of their subscripts). Others will later on be interpreted parameters when fitting the model to data.
The HelicityModel
comes with expressions for these kinematic_variables
, so that it’s possible to compute them from 4-momentum data.
kinematic_variables = set(model_no_dynamics.kinematic_variables)
kinematic_variables
{m_0, m_012, m_1, m_12, m_2, phi_0, phi_1^12, theta_0, theta_1^12}
theta = sp.Symbol("theta_0", real=True)
expr = model_no_dynamics.kinematic_variables[theta]
expr.doit()
The remaining symbols in the HelicityModel
are parameters. Each of these parameters comes with suggested parameter values (parameter_defaults
), that have been extracted from the ReactionInfo
object where possible:
Helicity couplings#
If you prefer to characterize the strength of each partial two-body decay, set HelicityAmplitudeBuilder.use_helicity_couplings
to True
and formulate a model:
model_builder.use_helicity_couplings = True
model_with_couplings = model_builder.formulate()
Math(aslatex(model_with_couplings.parameter_defaults))
Scalar masses#
By default, the HelicityAmplitudeBuilder
creates sympy.Expr
s for each kinematic variable, including all ‘invariant’ final state masses (\(m_0, m_1, \dots\)). However, it often happens that certain particles in a final state are stable. In that case, you may want to substitute these symbols with scalar values. This can be achieved by specifying which final state IDs are to be considered stable. Their corresponding mass symbols will then be considered parameters and a scalar suggested parameter value will be provided.
{m_012, m_12, phi_0, phi_1^12, theta_0, theta_1^12}
We can reset this option as follows:
Similarly, it can happen that the initial state mass is fixed (for instance when analyzing generated data or when applying a kinematic fit while reconstructing measured data). In that case, set HelicityAmplitudeBuilder.scalar_initial_state_mass
to True
.
Extend kinematic variables#
The HelicityAmplitudeBuilder
by default only generates kinematic_variables
(helicity angles and invariant masses) for the topologies that are available in the ReactionInfo
object that it was created with. If you want to calculate more kinematic variables, you can use the method register_topology()
of its helicity HelicityAmplitudeBuilder.adapter
to register more topologies and generate more kinematic variables. This is especially useful when generating data later on with TensorWaves.
To make this a bit easier, there is permutate_registered_topologies()
, which is a small convenience function makes it possible to generate permutations of all registered_topologies
and register them as well. Note that initially, only one Topology
is registered in the HelicityAmplitudeBuilder.adapter
, namely the one for the decay \(J/\psi \to \gamma f_0, f_0 \to \pi^0\pi^0\):
dot = qrules.io.asdot(model_builder.adapter.registered_topologies)
graphviz.Source(dot)
We now permutate_registered_topologies()
to register permutations of this Topology
:
model_builder.adapter.permutate_registered_topologies()
There are now three registered_topologies
―one for each permutation:
len(model_builder.adapter.registered_topologies)
3
assert len(model_builder.adapter.registered_topologies) == 3
dot = qrules.io.asdot(model_builder.adapter.registered_topologies)
graphviz.Source(dot)
And if we formulate()
a new HelicityModel
, we see that it has many more kinematic_variables
:
{m_0,
m_01,
m_02,
m_1,
m_12,
m_2,
phi_0,
phi_01,
phi_02,
phi_0^01,
phi_0^02,
phi_1^12,
theta_0,
theta_01,
theta_02,
theta_0^01,
theta_0^02,
theta_1^12}
Tip
To register even more topologies, use e.g. create_isobar_topologies()
to generate other, non-isomorphic topologies that cannot be created with permutations. This is relevant for more than three final states.
Set dynamics#
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 (see Custom dynamics).
AmpForm does provide a few common dynamics
functions, which can be constructed as Expr
with the correct Symbol
names using set_dynamics()
. This function takes specific dynamics.builder
functions and classes, such as RelativisticBreitWignerBuilder
, which can create relativistic_breit_wigner()
functions for specific resonances. 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 import PhaseSpaceFactor # optional
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
bw_builder = RelativisticBreitWignerBuilder(
energy_dependent_width=True,
form_factor=True,
phsp_factor=PhaseSpaceFactor, # optional
)
for name in reaction.get_intermediate_particles().names:
model_builder.set_dynamics(name, bw_builder)
Note that this RelativisticBreitWignerBuilder
can also be initialized with a different PhaseSpaceFactorProtocol
. This allows us to insert different phase space factors (see Analytic continuation and create_analytic_breit_wigner()
).
See also
And we use the reconfigured HelicityAmplitudeBuilder
to generate another HelicityModel
, this time with relativistic Breit-Wigner functions and form factors:
Export#
There is no special export function to export an HelicityModel
. However, we can just use the built-in pickle
module to write the model to disk and load it back:
import pickle
with open("helicity_model.pickle", "wb") as stream:
pickle.dump(model, stream)
with open("helicity_model.pickle", "rb") as stream:
model = pickle.load(stream)
Cached expression ‘unfolding’#
Amplitude model expressions can be extremely large. AmpForm can formulate such expressions relatively fast, but sympy
has to ‘unfold’ these expressions with doit()
, which can take a long time. AmpForm provides a function that can cache the ‘unfolded’ expression to disk, so that the expression unfolding runs faster upon the next run.
from ampform.sympy import perform_cached_doit
full_expression = perform_cached_doit(model.expression)
sp.count_ops(full_expression)
815
See perform_cached_doit()
for some tips on how to improve performance.
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)
Tip
To view the full expression for the amplitude model without crashing Jupyter Lab, install jupyterlab-katex
.
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_symbols()
):
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
Show code cell source
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:
Show code cell source
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()
Tip
Use HelicityModel.sum_components()
for adding up separate components of the model.
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) \xrightarrow[S=1]{L=0} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(1500) \gamma; f_{0}(1500) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}},
C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma; f_{0}(980) \xrightarrow[S=0]{L=0} \pi^{0} \pi^{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_0,
phi_1^12,
theta_0,
theta_1^12]
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_0, phi_1^12, theta_0, theta_1^12}
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. (Alternatively, see Scalar masses.)
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:
…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=Rf"$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):
Show code cell source
def round_nested(expression: sp.Expr, n_decimals: int = 2) -> sp.Expr:
for node in sp.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
Show code cell source
rounded = round_nested(plotted_intensity)
round_nested(rounded)