Create amplitude models¶
Your analysis starts by defining an amplitude model that describes the reaction process you want to study. Such a model is generally very complex and requires a fair amount of effort by the analyst (you). This also gives a lot of room for mistakes.
The ‘expert system’ is responsible to give you advice on the form of an amplitude model based on the problem set you define (initial state, final state, allowed interactions, intermediate states, etc.). Internally, the system propagates the quantum numbers through the reaction graph while satisfying the specified conservation rules. How to control this procedure is explained in more detail below.
Afterwards, the amplitude model of the expert system can be exported into ComPWA or Tensorwaves. The model can for instance be used to generate a data set (toy Monte Carlo) for this reaction and to optimize its parameters to resemble an actual data set as good as possible.
1.1. Define the problem set¶
We first define the boundary conditions of our physics problem, such as initial state, final state, formalism type, etc. and pass all of that information to the StateTransitionManager
. This is the main user interface class of the expertsystem
.
By default, the StateTransitionManager
loads all particles from the PDG. The expertsystem
would take a long time to check the quantum numbers of all these particles, so in this notebook, we use a smaller subset of relatively common particles.
from expertsystem.reaction import InteractionTypes, StateTransitionManager
stm = StateTransitionManager(
initial_state=["J/psi(1S)"],
final_state=["gamma", "pi0", "pi0"],
formalism_type="helicity",
)
The StateTransitionManager
(STM) is the main user interface class of the {mod}`expertsystem`. The boundary conditions of your physics problem, such as the initial state, final state, formalism type, etc., are defined through this interface.
create_problem_sets
of the STM creates all problem sets ― using the boundary conditions of theStateTransitionManager
instance. In total 3 steps are performed. The creation of reaction topologies. The creation ofInitialFacts
, based on a topology and the initial and final state information. And finally the solving settings such as the conservation laws and quantum number domains to use at which point of the topology.By default, all three interaction types (
EM
,Strong
, andWeak
) are used in the preparation stage. However, it is also possible to choose the allowed interaction types globally viaset_allowed_interaction_types
.
After the preparation step, you can modify the problem sets returned by create_problem_sets
to your liking. Since the output of this function contains quite a lot of information, the expertsystem
aids in the configuration (especially the STM).
A subset of particles that are allowed as intermediate states can also be specified: either through the
STM's constructor
or by setting the instanceallowed_intermediate_particles
.
1.2. Prepare Problem Sets¶
Create all ~.ProblemSet
’s using the boundary conditions of the ~.StateTransitionManager
instance. By default it uses the isobar model (tree of two-body
decays) to build ~.Topology
’s. Various ~.InitialFacts
are created for each topology based on the initial and final state. Lastly some reasonable default settings for the solving process are chosen. Remember that each interaction node defines its own set of conservation laws.
The StateTransitionManager
(STM) defines three interaction types:
Interaction |
Strength |
---|---|
strong |
\(60\) |
electromagnetic (EM) |
\(1\) |
weak |
\(10^{-4}\) |
By default, all three are used in the preparation stage. The create_problem_sets()
method of the STM generates graphs with all possible combinations of interaction nodes. An overall interaction strength is assigned to each graph and they are grouped according to this strength.
problem_sets = stm.create_problem_sets()
1.3. Find solutions¶
If you are happy with the default settings generated by the StateTransitionManager
, just start with solving directly!
This step takes about 23 sec on an Intel(R) Core(TM) i7-6820HQ CPU of 2.70GHz running, multi-threaded.
result = stm.find_solutions(problem_sets)
The find_solutions()
method returns a Result
object from which you can extract the solutions
and any violated_edge_rules
in case solutions were found. Now, you can use get_intermediate_particles()
to print the names of the intermediate states that the StateTransitionManager
found:
print("found", len(result.solutions), "solutions!")
result.get_intermediate_particles().names
found 414 solutions!
{'a(0)(1450)0',
'a(0)(980)0',
'a(1)(1260)0',
'a(1)(1640)0',
'a(2)(1320)0',
'a(2)(1700)0',
'b(1)(1235)0',
'f(0)(1370)',
'f(0)(1500)',
'f(0)(1710)',
'f(0)(500)',
'f(0)(980)',
'f(1)(1285)',
'f(1)(1420)',
"f(2)'(1525)",
'f(2)(1270)',
'f(2)(1950)',
'f(2)(2010)',
'f(2)(2300)',
'f(2)(2340)',
'h(1)(1170)',
'omega(1420)',
'omega(1650)',
'omega(782)',
'phi(1020)',
'phi(1680)',
'rho(1450)0',
'rho(1700)0',
'rho(770)0'}
About the number of solutions
The “number of solutions
” is the total number of allowed StateTransitionGraph
instances that the StateTransitionManager
has found. This also includes all allowed spin projection combinations. In this channel, we for example consider a \(J/\psi\) with spin projection \(\pm1\) that decays into a \(\gamma\) with spin projection \(\pm1\), which already gives us four possibilities.
On the other hand, the intermediate state names that was extracted with Result.get_intermediate_particles()
, is just a set
of the state names on the intermediate edges of the list of solutions
, regardless of spin projection.
Tip
Have a look at Result.get_particle_graphs()
, if you want to extract a list
of StateTransitionGraph
s of which the spin projections have been stripped. This is especially useful when visualizing the decay topologies.
Now we have a lot of solutions that are actually heavily suppressed (they involve two weak decays).
In general, you can modify the ~.ProblemSet
s returned by create_problem_sets()
directly, but the STM also comes with functionality to globally choose the allowed interaction types. So, go ahead and disable the EM
and Weak
interactions:
stm.set_allowed_interaction_types([InteractionTypes.Strong])
problem_sets = stm.create_problem_sets()
result = stm.find_solutions(problem_sets)
print("found", len(result.solutions), "solutions!")
result.get_intermediate_particles().names
found 192 solutions!
{'b(1)(1235)0',
'f(0)(1370)',
'f(0)(1500)',
'f(0)(1710)',
'f(0)(500)',
'f(0)(980)',
"f(2)'(1525)",
'f(2)(1270)',
'f(2)(1950)',
'f(2)(2010)',
'f(2)(2300)',
'f(2)(2340)',
'rho(1450)0',
'rho(1700)0',
'rho(770)0'}
Now note that, since a \(\gamma\) particle appears in one of the interaction nodes, the expert system knows that this node must involve EM interactions! Because the node can be an effective interaction, the weak interaction cannot be excluded, as it contains only a subset of conservation laws.
Since only the strong interaction was supposed to be used, this results in a warning and the STM automatically corrects the mistake.
Once the EM interaction is included, this warning disappears. Be aware, however, that the EM interaction is now available globally. Hence, there now might be solutions in which both nodes are electromagnetic.
stm.set_allowed_interaction_types([InteractionTypes.Strong, InteractionTypes.EM])
problem_sets = stm.create_problem_sets()
result = stm.find_solutions(problem_sets)
print("found", len(result.solutions), "solutions!")
result.get_intermediate_particles().names
found 318 solutions!
{'a(0)(1450)0',
'a(0)(980)0',
'a(2)(1320)0',
'a(2)(1700)0',
'b(1)(1235)0',
'f(0)(1370)',
'f(0)(1500)',
'f(0)(1710)',
'f(0)(500)',
'f(0)(980)',
"f(2)'(1525)",
'f(2)(1270)',
'f(2)(1950)',
'f(2)(2010)',
'f(2)(2300)',
'f(2)(2340)',
'h(1)(1170)',
'omega(1420)',
'omega(1650)',
'omega(782)',
'phi(1020)',
'phi(1680)',
'rho(1450)0',
'rho(1700)0',
'rho(770)0'}
Great! Now we selected only the strongest contributions. Be aware, though, that there are more effects that can suppress certain decays, like small branching ratios. In this example, the initial state \(J/\Psi\) can decay into \(\pi^0 + \rho^0\) or \(\pi^0 + \omega\).
decay |
branching ratio |
---|---|
\(\omega\to\gamma+\pi^0\) |
0.0828 |
\(\rho^0\to\gamma+\pi^0\) |
0.0006 |
Unfortunately, the \(\rho^0\) mainly decays into \(\pi^0+\pi^0\), not \(\gamma+\pi^0\) and is therefore suppressed. This information is currently not known to the expert system, but it is possible to hand the expert system a list of allowed intermediate states.
# particles are found by name comparison,
# i.e. f2 will find all f2's and f all f's independent of their spin
stm.set_allowed_intermediate_particles(["f(0)", "f(2)"])
result = stm.find_solutions(problem_sets)
print("found", len(result.solutions), "solutions!")
result.get_intermediate_particles().names
found 138 solutions!
{'f(0)(1370)',
'f(0)(1500)',
'f(0)(1710)',
'f(0)(500)',
'f(0)(980)',
"f(2)'(1525)",
'f(2)(1270)',
'f(2)(1950)',
'f(2)(2010)',
'f(2)(2300)',
'f(2)(2340)'}
Now we have selected all amplitudes that involve f states. The warnings appear only to notify the user that the list of solutions is not exhaustive: for certain edges in the graph, no suitable particle was found (since only f states were allowed).
1.4. Generate an amplitude model¶
Now that we are satisfied with the intermediate resonances, we are all set convert the solutions that the STM found to an amplitude model! This can be done with generate
. This function deduces the amplitude formalism from the formalism_type
that you specified when constructing the STM. The result is an AmplitudeModel
from which you can extract all information about Particle
instances, FitParameters
, Kinematics
etc.
import expertsystem as es
amplitude_model = es.amplitude.generate(result)
list(amplitude_model.parameters)
['MesonRadius_J/psi(1S)',
'Magnitude_J/psi(1S)_to_f(2)(2300)_2+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2300)_2+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2010)_2+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2010)_2+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2340)_2+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2340)_2+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(1270)_2+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1270)_2+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
"Magnitude_J/psi(1S)_to_f(2)'(1525)_2+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
"Phase_J/psi(1S)_to_f(2)'(1525)_2+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
'Magnitude_J/psi(1S)_to_f(2)(1950)_2+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1950)_2+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2300)_1+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2300)_1+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2010)_1+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2010)_1+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2340)_1+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2340)_1+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(1270)_1+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1270)_1+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
"Magnitude_J/psi(1S)_to_f(2)'(1525)_1+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
"Phase_J/psi(1S)_to_f(2)'(1525)_1+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
'Magnitude_J/psi(1S)_to_f(2)(1950)_1+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1950)_1+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2300)_0+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2300)_0+gamma_1;f(2)(2300)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2010)_0+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2010)_0+gamma_1;f(2)(2010)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(2340)_0+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(2340)_0+gamma_1;f(2)(2340)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(2)(1270)_0+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1270)_0+gamma_1;f(2)(1270)_to_pi0_0+pi0_0;',
"Magnitude_J/psi(1S)_to_f(2)'(1525)_0+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
"Phase_J/psi(1S)_to_f(2)'(1525)_0+gamma_1;f(2)'(1525)_to_pi0_0+pi0_0;",
'Magnitude_J/psi(1S)_to_f(2)(1950)_0+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(2)(1950)_0+gamma_1;f(2)(1950)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(0)(500)_0+gamma_1;f(0)(500)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(0)(500)_0+gamma_1;f(0)(500)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(0)(1370)_0+gamma_1;f(0)(1370)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(0)(1370)_0+gamma_1;f(0)(1370)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(0)(1710)_0+gamma_1;f(0)(1710)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(0)(1710)_0+gamma_1;f(0)(1710)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(0)(980)_0+gamma_1;f(0)(980)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(0)(980)_0+gamma_1;f(0)(980)_to_pi0_0+pi0_0;',
'Magnitude_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;',
'Phase_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;',
'strength_incoherent']
The AmplitudeModel
also contains a dynamics
section. By default, it only contains the dynamics for the initial state ― NonDynamic
:
list(amplitude_model.dynamics)
['J/psi(1S)']
type(amplitude_model.dynamics["J/psi(1S)"])
expertsystem.amplitude.model.NonDynamic
To set dynamics for specific resonances, use e.g. ParticleDynamics.set_breit_wigner()
:
for name in result.get_intermediate_particles().names:
amplitude_model.dynamics.set_breit_wigner(name)
{name: type(dyn) for name, dyn in amplitude_model.dynamics.items()}
{'J/psi(1S)': expertsystem.amplitude.model.NonDynamic,
'f(0)(1500)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(0)(500)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(0)(1710)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(0)(980)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(2)(1950)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(2)(2300)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(2)(2010)': expertsystem.amplitude.model.RelativisticBreitWigner,
"f(2)'(1525)": expertsystem.amplitude.model.RelativisticBreitWigner,
'f(0)(1370)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(2)(2340)': expertsystem.amplitude.model.RelativisticBreitWigner,
'f(2)(1270)': expertsystem.amplitude.model.RelativisticBreitWigner}
Finally, you can use the expertsystem.io
module to write()
the AmplitudeModel
to a file (either XML or YAML):
es.io.write(amplitude_model, filename="model.xml")
es.io.write(amplitude_model, filename="model.yml")
This allows you to cache file AmplitudeModel
to disk and then load it back using load_amplitude_model()
:
imported_model_xml = es.io.load_amplitude_model(filename="model.xml")
imported_model_yml = es.io.load_amplitude_model(filename="model.yml")
assert imported_model_xml == amplitude_model
assert imported_model_yml == amplitude_model
Have a look through the sections of the resulting XML or YAML recipe file to see what you recognize from the problem set defined above. There may also be some things you want to change in there manually, so make sure you store this recipe file carefully (e.g. track it with Git) as to avoid overwriting it your changes after rerunning the expert system.
Now you can use this recipe file as an amplitude model in an external PWA fitter package!