Generate transitions

A Partial Wave Analysis starts by defining an amplitude model that describes the reaction process that is to be analyzed. Such a model is generally very complex and requires a fair amount of effort by the analyst (you). This gives a lot of room for mistakes.

QRules 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 produced by AmpForm can be exported into 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. For more info on that see Formulate amplitude model.

Note

Simple channels can be treated with the generate_transitions() façade function. This notebook shows how to treat more complicated cases with the StateTransitionManager.

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 qrules.

By default, the StateTransitionManager loads all particles from the PDG. The qrules 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.

How does the StateTransitionManager know what a "J/psi(1S)" is? Upon construction, the StateTransitionManager loads default definitions from the PDG. See Particle database for how to add custom particle definitions.

from qrules import InteractionType, StateTransitionManager

stm = StateTransitionManager(
    initial_state=["J/psi(1S)"],
    final_state=["gamma", "pi0", "pi0"],
    formalism="helicity",
    max_angular_momentum=2,
)

State Transition Manager

The StateTransitionManager (STM) is the main user interface class of qrules. 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 the StateTransitionManager instance. In total 3 steps are performed. The creation of reaction topologies. The creation of InitialFacts, 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, and WEAK) are used in the preparation stage. However, it is also possible to choose the allowed interaction types globally via set_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, qrules 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 instance allowed_intermediate_particles.

Tip

Custom topologies shows how to provide custom Topology instances to the STM, so that you generate more than just isobar decays.

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()
sorted(problem_sets, reverse=True)

[60.0, 1.0, 0.0001]

To get an idea of what these ProblemSets represent, you can use asdot() and Graphviz to visualize one of them (see Visualize solutions):

import graphviz

from qrules import io

some_problem_set = problem_sets[60.0][0]
dot = io.asdot(some_problem_set, render_node=True)
graphviz.Source(dot)

Each ProblemSet provides a mapping of initial_facts that represent the initial and final states with spin projections. The nodes and edges in between these initial_facts are still to be generated. This will be done from the provided solving_settings. There are two mechanisms there:

1. One the one hand, the EdgeSettings.qn_domains and NodeSettings.qn_domains contained in the solving_settings define the domain over which quantum number sets can be generated.

2. On the other, the EdgeSettings.rule_priorities and NodeSettings.rule_priorities in solving_settings define which conservation_rules are used to determine which of the sets of generated quantum numbers are valid.

Together, these two constraints allow the StateTransitionManager to generate a number of MutableTransitions that comply with the selected conservation_rules.

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.

reaction = stm.find_solutions(problem_sets)

Tip

See Quantum number solutions for a visualization of the intermediate steps.

The find_solutions() method returns a ReactionInfo object from which you can extract the transitions. Now, you can use get_intermediate_particles() to print the names of the intermediate states that the StateTransitionManager found:

print("found", len(reaction.transitions), "solutions!")
reaction.get_intermediate_particles().names

found 420 solutions!

['a(0)(980)0',
 'a(1)(1260)0',
 'a(2)(1320)0',
 'a(0)(1450)0',
 'a(1)(1640)0',
 'a(2)(1700)0',
 'b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(1)(1285)',
 'f(0)(1370)',
 'f(1)(1420)',
 "f(2)'(1525)",
 'f(0)(1500)',
 'f(0)(1710)',
 'f(2)(1950)',
 'f(0)(2020)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(782)',
 'omega(1420)',
 'omega(1650)',
 'phi(1020)',
 'phi(1680)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)0']

About the number of solutions

The “number of transitions” is the total number of allowed MutableTransition 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 ReactionInfo.get_intermediate_particles(), is just a set of the state names on the intermediate edges of the list of transitions, regardless of spin projection.

Now we have a lot of solutions that are actually heavily suppressed (they involve two weak decays).

Select interaction types

In general, you can modify the ProblemSets 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 InteractionType.WEAK interactions:

stm.set_allowed_interaction_types([InteractionType.STRONG])
problem_sets = stm.create_problem_sets()
reaction = stm.find_solutions(problem_sets)

print("found", len(reaction.transitions), "solutions!")
reaction.get_intermediate_particles().names

found 198 solutions!

['b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(0)(1370)',
 "f(2)'(1525)",
 'f(0)(1500)',
 'f(0)(1710)',
 'f(2)(1950)',
 'f(0)(2020)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)0']

Now note that, since a \(\gamma\) particle appears in one of the interaction nodes, qrules 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.

stm.set_allowed_interaction_types([InteractionType.STRONG, InteractionType.EM])
problem_sets = stm.create_problem_sets()
reaction = stm.find_solutions(problem_sets)

print("found", len(reaction.transitions), "solutions!")
reaction.get_intermediate_particles().names

found 324 solutions!

['a(0)(980)0',
 'a(2)(1320)0',
 'a(0)(1450)0',
 'a(2)(1700)0',
 'b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(0)(1370)',
 "f(2)'(1525)",
 'f(0)(1500)',
 'f(0)(1710)',
 'f(2)(1950)',
 'f(0)(2020)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(782)',
 'omega(1420)',
 'omega(1650)',
 'phi(1020)',
 'phi(1680)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)0']

This automatic selection of conservation rules can be switched of the StateTransitionManager.interaction_determinators. Here’s an example where we deselect the check that causes makes detects the existence of a photon in the decay chain. Note, however, that for \(J/\psi \to \gamma\pi^0\pi^0\), this results in non-executed node rules:

from qrules.system_control import GammaCheck

stm_no_check = StateTransitionManager(
    initial_state=["J/psi(1S)"],
    final_state=["gamma", "pi0", "pi0"],
)
stm_no_check.interaction_determinators = [
    check
    for check in stm_no_check.interaction_determinators
    if not isinstance(check, GammaCheck)
]
stm_no_check.set_allowed_interaction_types([InteractionType.STRONG])
problem_sets_no_check = stm_no_check.create_problem_sets()
try:
    reaction_no_check = stm_no_check.find_solutions(problem_sets_no_check)
except RuntimeError as e:
    msg, execution_info = e.args
execution_info

ExecutionInfo(
  not_executed_node_rules=defaultdict(set,
              {0: {'g_parity_conservation', 'isospin_conservation'},
               1: {'g_parity_conservation', 'isospin_conservation'}}),
  violated_node_rules=defaultdict(set, {0: set(), 1: set()}),
  not_executed_edge_rules=defaultdict(set, {0: {'isospin_validity'}}),
  violated_edge_rules=defaultdict(set, {}),
)

Be aware that after calling set_allowed_interaction_types(), the EM interaction is now selected for all nodes, for each node in the decay topology. Hence, there now might be solutions in which both nodes are electromagnetic. This is fine for the decay \(J/\psi \to \gamma \pi^0 \pi^0\), but for decays that require the WEAK interaction type, you want to set the interaction type per specific nodes. Take for example the decay \(\Lambda_c^+ \to p K^- \pi^+\), which has a production node that is mediated by the weak force and a decay node that goes via the strong force. In this case, only limit the decay node to the STRONG force:

lc2pkpi_stm = StateTransitionManager(
    initial_state=["Lambda(c)+"],
    final_state=["p", "K-", "pi+"],
    mass_conservation_factor=0.6,
)
lc2pkpi_stm.set_allowed_interaction_types([InteractionType.STRONG], node_id=1)
lc2pkpi_problem_sets = lc2pkpi_stm.create_problem_sets()
lc2pkpi_reaction = lc2pkpi_stm.find_solutions(lc2pkpi_problem_sets)

Show code cell source Hide code cell source

dot = io.asdot(lc2pkpi_reaction, collapse_graphs=True)
graphviz.Source(dot)

Select intermediate particles

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 qrules, but it is possible to hand qrules a list of allowed intermediate states,

stm.set_allowed_intermediate_particles(["f(0)", "f(2)"])
reaction = stm.find_solutions(problem_sets)
reaction.get_intermediate_particles().names

Show code cell output Hide code cell output

['f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(0)(1370)',
 "f(2)'(1525)",
 'f(0)(1500)',
 'f(0)(1710)',
 'f(2)(1950)',
 'f(0)(2020)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)']

or, using regular expressions,

stm.set_allowed_intermediate_particles(r"f\([02]\)", regex=True)
reaction = stm.find_solutions(problem_sets)
assert len(reaction.get_intermediate_particles().names) == 12

Now we have selected all amplitudes that involve f states:

Show code cell source Hide code cell source

dot = io.asdot(reaction, collapse_graphs=True, render_node=False)
graphviz.Source(dot)

See also

Visualize solutions

4. Export generated transitions

The ReactionInfo, MutableTransition, and Topology can be serialized to and from a dict with io.asdict() and io.fromdict():

io.asdict(reaction.transitions[0].topology)

{'nodes': frozenset({0, 1}),
 'edges': {-1: {'ending_node_id': 0},
  3: {'originating_node_id': 0, 'ending_node_id': 1},
  0: {'originating_node_id': 0},
  1: {'originating_node_id': 1},
  2: {'originating_node_id': 1}}}

YAML is more human-readable than JSON, but reading and writing JSON is faster.

This also means that the ReactionInfo can be written to JSON or YAML format with io.write() and loaded again with io.load():

io.write(reaction, "transitions.json")
imported_reaction = io.load("transitions.json")
assert imported_reaction == reaction

Handy if it takes a lot of computation time to re-generate the transitions!

Tip

The ampform package can formulate amplitude models based on the state transitions created by qrules. See Formulate amplitude model.