Welcome to QRules!¶

QRules is a Python package for validating and generating particle reactions using quantum number conservation rules. The user only has to provide a certain set of boundary conditions (initial and final state, allowed interaction types, expected decay topologies, etc.). QRules will then span the space of allowed quantum numbers over all allowed decay topologies and particle instances that correspond with the sets of allowed quantum numbers it has found.

The resulting state transition objects are particularly useful for amplitude analysis / Partial Wave Analysis as they contain all information (such as expected masses, widths, and spin projections) that is needed to formulate an amplitude model.

The Usage pages illustrate several features of qrules. You can run each of them as Jupyter notebooks with the launch button in the top-right corner. Enjoy!

Internal design

QRules consists of three major components:

State transition graphs

A StateTransitionGraph is a directed graph that consists of nodes and edges. In a directed graph, each edge must be connected to at least one node (in correspondence to Feynman graphs). This way, a graph describes the transition from one state to another.
- Edges correspond to states (particles with spin). In other words, edges are a collection of properties such as the quantum numbers that characterize a state that the particle is in.
- Nodes represents interactions and contain all information for the transition of this specific step. Most importantly, a node contains a collection of conservation rules that have to be satisfied. An interaction node has \(M\) ingoing lines and \(N\) outgoing lines, where \(M,N \in \mathbb{Z}\), \(M > 0, N > 0\).
Conservation rules

The central component are the conservation_rules. They belong to individual nodes and receive properties about the node itself, as well as properties of the ingoing and outgoing edges of that node. Based on those properties the conservation rules determine whether edges pass or not.
Solvers

The determination of the correct state properties in the graph is done by solvers. New properties are set for intermediate edges and interaction nodes and their validity is checked with the conservation rules.

QRules workflow

Preparation

1.1. Build all possible topologies. A topology is represented by a StateTransitionGraph, in which the edges and nodes are empty (no particle information).

1.2. Fill the topology graphs with the user provided information. Typically these are the graph’s ingoing edges (initial state) and outgoing edges (final state).
Solving

2.1. Propagate quantum number information through the complete graph while respecting the specified conservation laws. Information like mass is not used in this first solving step.

2.2. Clone graphs while inserting concrete matching particles for the intermediate edges (mainly adds the mass variable).

2.3. Validate the complete graphs, so run all conservation law check that were postponed from the first step.

Table of Contents

Installation¶

The fastest way of installing this package is through PyPI or Conda:

PyPI

python3 -m pip install qrules

Conda

conda install -c conda-forge qrules

This installs the latest, stable release that you can find on the stable branch.

The latest version on the main branch can be installed as follows:

python3 -m pip install git+https://github.com/ComPWA/qrules@main

In that case, however, we highly recommend using the more dynamic ‘editable installation’ instead. This goes as follows:

Get the source code:

git clone https://github.com/ComPWA/qrules.git
cd qrules

[Recommended] Create a virtual environment (see here).
Install the project as an ‘editable installation’ and install additional packages for the developer:
```
python3 -m pip install -e .[dev]
```
Pinning dependency versions
In order to install the exact same versions of the dependencies with which the framework has been tested, use the provided constraints files for the specific Python version 3.x you are using:
python3 -m pip install -c .constraints/py3.x.txt -e .[dev]
See also

Pinning dependency versions

That’s all! Have a look at the Usage page to try out the package. You can also have a look at the Help developing page for tips on how to work with this ‘editable’ developer setup!

Usage¶

Main interface¶

Here are some quick examples of how to use qrules. For more fine-grained control, have a look at Advanced.

Investigate intermediate resonances¶

import qrules

reaction = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", "Sigma+", "p~"],
    allowed_interaction_types="strong",
    formalism="canonical-helicity",
)

import graphviz

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

Next, you use the ampform package to convert these transitions into a mathematical description that you can use to fit your data and perform Partial Wave Analysis!

Quantum number search¶

The load_pdg() function creates a ParticleCollection containing the latest PDG info. Its find() and filter() methods allows you to quickly look up the quantum numbers of a particle and, vice versa, look up particle candidates based on a set of quantum numbers.

import qrules

pdg = qrules.load_pdg()
pdg.find(22)  # by pid
pdg.find("Delta(1920)++")

Particle(
  name='Delta(1920)++',
  pid=22224,
  latex='\\Delta(1920)^{++}',
  spin=1.5,
  mass=1.92,
  width=0.3,
  charge=2,
  isospin=Spin(3/2, +3/2),
  baryon_number=1,
  parity=+1,
)

subset = pdg.filter(
    lambda p: p.spin in [2.5, 3.5, 4.5] and p.name.startswith("N")
)
subset.names

['N(1675)~-',
 'N(1675)~0',
 'N(1675)0',
 'N(1675)+',
 'N(1680)~-',
 'N(1680)~0',
 'N(1680)0',
 'N(1680)+',
 'N(2190)~-',
 'N(2190)0',
 'N(2190)~0',
 'N(2190)+']

Tip

See Particle database

Check allowed reactions¶

qrules can be used to check whether a transition between an initial and final state is violated by any conservation rules:

qrules.check_reaction_violations(
    initial_state="pi0",
    final_state=["gamma", "gamma", "gamma"],
)

{frozenset({'c_parity_conservation'})}

Advanced¶

Each of the qrules’s sub-modules offer functionality to handle more advanced reaction types. The following notebooks illustrate how use them.

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.

from qrules import InteractionType, StateTransitionManager

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

StateTransitionManager

The StateTransitionManager (STM) is the main user interface class of {mod}`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, {mod}`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 decay topologies):

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 (GraphSettings). There are two mechanisms there:

One the one hand, the EdgeSettings.qn_domains and NodeSettings.qn_domains contained in the GraphSettings define the domain over which quantum number sets can be generated.
On the other, the EdgeSettings.rule_priorities and NodeSettings.rule_priorities in GraphSettings 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 StateTransitionGraphs 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)

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 270 solutions!

['a(0)(980)0',
 'a(1)(1260)0',
 'a(0)(1450)0',
 'a(1)(1640)0',
 'b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(1)(1285)',
 'f(0)(1370)',
 'f(1)(1420)',
 'f(0)(1500)',
 'f(0)(1710)',
 '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 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 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).

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 84 solutions!

['b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 '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. 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([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 174 solutions!

['a(0)(980)0',
 'a(0)(1450)0',
 'b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'h(1)(1170)',
 'omega(782)',
 'omega(1420)',
 'omega(1650)',
 'phi(1020)',
 'phi(1680)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)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 qrules, but it is possible to hand qrules 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)"])

reaction = stm.find_solutions(problem_sets)

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

found 30 solutions!

['f(0)(500)', 'f(0)(980)', 'f(0)(1370)', 'f(0)(1500)', 'f(0)(1710)']

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

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

See also

Visualize decay topologies

4. Export generated transitions¶

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

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

{'nodes': [0, 1],
 'edges': {-1: {'ending_node_id': 0},
  0: {'originating_node_id': 0},
  1: {'originating_node_id': 1},
  2: {'originating_node_id': 1},
  3: {'originating_node_id': 0, 'ending_node_id': 1}}}

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.

Particle database¶

In PWA, you usually want to search for special resonances, possibly even some not listed in the PDG. In this notebook, we go through a few ways to add or overwrite Particle instances in the database with your own particle definitions.

Loading the default database¶

In Generate transitions, we made use of the StateTransitionManager. By default, if you do not specify the particle_db argument, the StateTransitionManager calls the function load_default_particles(). This functions returns a ParticleCollection instance with Particle definitions from the PDG, along with additional definitions that are provided in the file additional_definitions.yml.

Here, we call this method directly to illustrate what happens (we use load_pdg(), which loads a subset):

from qrules.particle import load_pdg

particle_db = load_pdg()
print("Number of loaded particles:", len(particle_db))

Number of loaded particles: 531

In the following, we illustrate how to use the methods of the ParticleCollection class to find and ‘modify’ Particles and add() them back to the ParticleCollection.

Finding particles¶

The ParticleCollection class offers some methods to search for particles by name or by PID (see find()):

particle_db.find(333)

Particle(
  name='phi(1020)',
  pid=333,
  latex='\\phi(1020)',
  spin=1.0,
  mass=1.019461,
  width=0.004248999999999999,
  isospin=Spin(0, 0),
  parity=-1,
  c_parity=-1,
  g_parity=-1,
)

With filter(), you can perform more sophisticated searches. This is done by either passing a function or lambda.

subset = particle_db.filter(lambda p: p.name.startswith("f(2)"))
subset.names

['f(2)(1270)',
 "f(2)'(1525)",
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)']

subset = particle_db.filter(
    lambda p: p.strangeness == 1
    and p.spin >= 1
    and p.mass > 1.8
    and p.mass < 1.9
)
subset.names

['K(2)(1820)0',
 'K(2)(1820)+',
 'Lambda(1820)~',
 'Lambda(1830)~',
 'Lambda(1890)~']

subset = particle_db.filter(lambda p: p.is_lepton())
subset.names

['e-',
 'e+',
 'mu-',
 'mu+',
 'nu(e)',
 'nu(tau)',
 'nu(mu)',
 'nu(tau)~',
 'nu(mu)~',
 'nu(e)~',
 'tau-',
 'tau+']

Note that in each of these examples, we call the names property. This is just to only display the names, sorted alphabetically, otherwise the output becomes a bit of a mess:

particle_db.filter(lambda p: p.name.startswith("pi") and len(p.name) == 3)

ParticleCollection({
  Particle(
    name='pi-',
    pid=-211,
    latex='\\pi^{-}',
    spin=0.0,
    mass=0.13957039000000002,
    width=2.5284e-17,
    charge=-1,
    isospin=Spin(1, -1),
    parity=-1,
    g_parity=-1,
  ),
  Particle(
    name='pi0',
    pid=111,
    latex='\\pi^{0}',
    spin=0.0,
    mass=0.1349768,
    width=7.81e-09,
    isospin=Spin(1, 0),
    parity=-1,
    c_parity=+1,
    g_parity=-1,
  ),
  Particle(
    name='pi+',
    pid=211,
    latex='\\pi^{+}',
    spin=0.0,
    mass=0.13957039000000002,
    width=2.5284e-17,
    charge=1,
    isospin=Spin(1, +1),
    parity=-1,
    g_parity=-1,
  ),
})

LaTeX representation¶

Particles also contain a latex tag. Here, we use ipython to render them nicely as mathematical symbols:

from IPython.display import Math

sigmas = particle_db.filter(
    lambda p: p.name.startswith("Sigma") and p.charmness == 1
)
Math(", ".join([p.latex for p in sigmas]))

\[\displaystyle \Sigma_{c}(2455)^{+}, \Sigma_{c}(2520)^{+}, \Sigma_{c}(2520)^{0}, \Sigma_{c}(2455)^{++}, \Sigma_{c}^{0}, \Sigma_{c}(2520)^{++}\]

Adding custom particle definitions through Python¶

A quick way to modify or overwrite particles, is through your Python script or notebook. Notice that the instances in the database are Particle instances:

N1650_plus = particle_db["N(1650)+"]
N1650_plus

Particle(
  name='N(1650)+',
  pid=32212,
  latex='N(1650)^{+}',
  spin=0.5,
  mass=1.65,
  width=0.125,
  charge=1,
  isospin=Spin(1/2, +1/2),
  baryon_number=1,
  parity=-1,
)

The instances in the database are immutable. Therefore, if you want to modify, say, the width, you have to create a new Particle instance from the particle you want to modify and add() it back to the database. You can do this with create_particle():

from qrules.particle import create_particle

new_N1650_plus = create_particle(
    template_particle=N1650_plus, name="Modified N(1650)+", width=0.2
)

particle_db.add(new_N1650_plus)
particle_db["Modified N(1650)+"].width

0.2

You often also want to add the antiparticle of the particle you modified to the database. Using create_antiparticle(), it is easy to create the corresponding antiparticle object.

from qrules.particle import create_antiparticle

new_N1650_minus = create_antiparticle(
    new_N1650_plus, new_name="Modified N(1650)-"
)

particle_db.add(new_N1650_minus)
particle_db["Modified N(1650)-"]

Particle(
  name='Modified N(1650)-',
  pid=-32212,
  latex='\\overline{N(1650)^{+}}',
  spin=0.5,
  mass=1.65,
  width=0.2,
  charge=-1,
  isospin=Spin(1/2, -1/2),
  baryon_number=-1,
  parity=+1,
)

When adding additional particles you may need for your research, it is easiest to work with an existing particle as template. Let’s say we want to study \(e^+e^-\) collisions of several energies:

energies_mev = {4180, 4220, 4420, 4600}
template_particle = particle_db["J/psi(1S)"]
for energy_mev in energies_mev:
    energy_gev = energy_mev / 1e3
    new_particle = create_particle(
        template_particle,
        name=f"EpEm ({energy_mev} MeV)",
        mass=energy_gev,
    )
    particle_db.add(new_particle)
len(particle_db)

particle_db.filter(lambda p: "EpEm" in p.name).names

['EpEm (4180 MeV)', 'EpEm (4220 MeV)', 'EpEm (4420 MeV)', 'EpEm (4600 MeV)']

Of course, it’s also possible to add any kind of custom Particle, as long as its quantum numbers comply with the gellmann_nishijima() rule:

from qrules.particle import Particle

custom = Particle(
    name="custom",
    pid=99999,
    latex=R"p_\mathrm{custom}",
    spin=1.0,
    mass=1,
    charge=1,
    isospin=(1.5, 0.5),
    charmness=1,
)
custom

Particle(
  name='custom',
  pid=99999,
  latex='p_\\mathrm{custom}',
  spin=1.0,
  mass=1.0,
  charge=1,
  isospin=Spin(3/2, +1/2),
  charmness=1,
)

particle_db += custom
len(particle_db)

Loading custom definitions from a YAML file¶

It’s also possible to add particles from a config file, with io.load(). Existing entries remain and if the imported file of particle definitions contains a particle with the same name, it is overwritten in the database.

It’s easiest to work with YAML. Here, we use the provided additional_particles.yml example file:

from qrules import io

particle_db += io.load("additional_particles.yml")

Writing to YAML¶

You can also dump the existing particle lists to YAML. You do this with the io.write() function.

io.write(instance=particle_db, filename="dumped_particle_list.yaml")

Note that the function write can dump any ParticleCollection to an output file, also a specific subset.

from qrules.particle import ParticleCollection

output = ParticleCollection()
output += particle_db["J/psi(1S)"]
output += particle_db.find(22)  # gamma
output += particle_db.filter(lambda p: p.name.startswith("f(0)"))
output += particle_db["pi0"]
output += particle_db["pi+"]
output += particle_db["pi-"]
output += particle_db["custom"]
io.write(output, "particle_list_selection.yml")
output.names

['custom',
 'f(0)(500)',
 'f(0)(980)',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'gamma',
 'J/psi(1S)',
 'pi0',
 'pi-',
 'pi+']

As a side note, qrules provides JSON schemas (reaction/particle-validation.json) to validate your particle list files (see also jsonschema.validate()). If you have installed qrules as an Editable installation and use VSCode, your YAML particle list are checked automatically in the GUI.

Visualize decay topologies¶

The io module allows you to convert StateTransitionGraph, Topology instances, and ProblemSets to DOT language with asdot(). You can visualize its output with third-party libraries, such as Graphviz. This is particularly useful after running find_solutions(), which produces a ReactionInfo object with a list of StateTransitionGraph instances (see Generate transitions).

Topologies¶

First of all, here are is an example of how to visualize a group of Topology instances. We use create_isobar_topologies() and create_n_body_topology() to create a few standard topologies.

import graphviz

import qrules
from qrules.topology import create_isobar_topologies, create_n_body_topology

topology = create_n_body_topology(2, 4)
graphviz.Source(qrules.io.asdot(topology, render_initial_state_id=True))

Note the IDs of the nodes is also rendered if there is more than node:

topologies = create_isobar_topologies(4)
graphviz.Source(qrules.io.asdot(topologies))

This can be turned on or off with the arguments of asdot():

topologies = create_isobar_topologies(3)
graphviz.Source(qrules.io.asdot(topologies, render_node=False))

asdot() provides other options as well:

topologies = create_isobar_topologies(5)
dot = qrules.io.asdot(
    topologies[0],
    render_final_state_id=False,
    render_resonance_id=True,
    render_node=False,
)
display(graphviz.Source(dot))

`ProblemSet`s¶

As noted in Generate transitions, the StateTransitionManager provides more control than the façade function generate_transitions(). One advantages, is that the StateTransitionManager first generates a set of ProblemSets with create_problem_sets() that you can further configure if you wish.

stm = qrules.StateTransitionManager(
    initial_state=["J/psi(1S)"],
    final_state=["K0", "Sigma+", "p~"],
    formalism="canonical-helicity",
)
problem_sets = stm.create_problem_sets()

Note that the output of create_problem_sets() is a dict with float values as keys (representing the interaction strength) and lists of ProblemSets as values.

sorted(problem_sets, reverse=True)

[3600.0, 60.0, 1.0, 0.006, 0.0001, 1e-08]

len(problem_sets[60.0])

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

`StateTransition`s¶

Here, we’ll visualize the allowed transitions for the decay \(\psi' \to \gamma\eta\eta\) as an example.

import qrules

reaction = qrules.generate_transitions(
    initial_state="psi(2S)",
    final_state=["gamma", "eta", "eta"],
    allowed_interaction_types="EM",
)

As noted in 3. Find solutions, the transitions contain all spin projection combinations (which is necessary for the ampform package). It is possible to convert all these solutions to DOT language with asdot(). To avoid visualizing all solutions, we just take a subset of the transitions:

dot = qrules.io.asdot(reaction.transitions[::50][:3])  # just some selection

This str of DOT language for the list of StateTransitionGraph instances can then be visualized with a third-party library, for instance, with graphviz.Source:

import graphviz

dot = qrules.io.asdot(
    reaction.transitions[::50][:3], render_node=False
)  # just some selection
graphviz.Source(dot)

You can also serialize the DOT string to file with io.write(). The file extension for a DOT file is .gv:

qrules.io.write(reaction, "decay_topologies_with_spin.gv")

Collapse graphs¶

Since this list of all possible spin projections transitions is rather long, it is often useful to use strip_spin=True or collapse_graphs=True to bundle comparable graphs. First, strip_spin=True allows one collapse (ignore) the spin projections (we again show a selection only):

dot = qrules.io.asdot(reaction.transitions[:3], strip_spin=True)
graphviz.Source(dot)

Note

By default, .asdot renders edge IDs, because they represent the (final) state IDs as well. In the example above, we switched this off.

If that list is still too much, there is collapse_graphs=True, which bundles all graphs with the same final state groupings:

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

Conservation rules¶

import attr
import graphviz

import qrules
from qrules.conservation_rules import (
    SpinEdgeInput,
    SpinNodeInput,
    parity_conservation,
    spin_conservation,
    spin_magnitude_conservation,
)
from qrules.quantum_numbers import Parity

QRules generates StateTransitionGraphs, populates them with quantum numbers (edge properties representing states and nodes properties representing interactions), then checks whether the generated StateTransitionGraphs comply with the rules formulated in the conservation_rules module.

The conservation_rules module can also be used separately. In this notebook, we will illustrate this by checking spin and parity conservation.

Parity conservation¶

See parity_conservation():

parity_conservation(
    ingoing_edge_qns=[Parity(-1)],
    outgoing_edge_qns=[Parity(+1), Parity(+1)],
    l_magnitude=1,
)

True

Spin conservation¶

spin_conservation() checks whether spin magnitude and spin projections are conserved. In addition, it checks whether the Clebsch-Gordan coefficients are non-zero, meaning that the coupled spins on the interaction nodes are valid as well.

No spin and angular momentum¶

spin_conservation(
    ingoing_spins=[
        SpinEdgeInput(0, 0),
    ],
    outgoing_spins=[
        SpinEdgeInput(0, 0),
        SpinEdgeInput(0, 0),
    ],
    interaction_qns=SpinNodeInput(
        l_magnitude=0,  # false if 1
        l_projection=0,
        s_magnitude=0,
        s_projection=0,
    ),
)

True

Non-zero example¶

spin_conservation(
    ingoing_spins=[
        SpinEdgeInput(1, 0),
    ],
    outgoing_spins=[
        SpinEdgeInput(1, +1),
        SpinEdgeInput(1, -1),
    ],
    interaction_qns=SpinNodeInput(
        l_magnitude=1,
        l_projection=0,
        s_magnitude=2,
        s_projection=0,
    ),
)

True

Example with a `StateTransition`¶

First, generate some StateTransitions with generate_transitions(), then select one of them:

reaction = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", "Sigma+", "p~"],
    allowed_interaction_types="strong",
    formalism="canonical",
)
transition = reaction.transitions[0]

Next, have a look at the edge and node properties, and use the underlying Topology to extract one of the node InteractionProperties with the surrounding states (these are tuples of a Particle and a float spin projection).

dot = qrules.io.asdot(transition, render_node=True)
display(graphviz.Source(dot))

dot = qrules.io.asdot(
    transition.topology,
    render_node=True,
    render_resonance_id=True,
    render_initial_state_id=True,
)
display(graphviz.Source(dot))

We select node \((0)\), which has incoming state ID \(-1\) and outgoing state IDs \(0\) and \(3\):

topology = transition.topology
node_id = 0
in_id, *_ = topology.get_edge_ids_ingoing_to_node(node_id)
out_id1, out_id2, *_ = topology.get_edge_ids_outgoing_from_node(node_id)

incoming_state = transition.states[in_id]
outgoing_state1 = transition.states[out_id1]
outgoing_state2 = transition.states[out_id2]
interaction = transition.interactions[node_id]

spin_magnitude_conservation(
    ingoing_spins=[
        SpinEdgeInput(
            spin_magnitude=incoming_state.particle.spin,
            spin_projection=incoming_state.spin_projection,
        )
    ],
    outgoing_spins=[
        SpinEdgeInput(
            spin_magnitude=outgoing_state1.particle.spin,
            spin_projection=outgoing_state1.spin_projection,
        ),
        SpinEdgeInput(
            spin_magnitude=outgoing_state2.particle.spin,
            spin_projection=outgoing_state2.spin_projection,
        ),
    ],
    interaction_qns=interaction,
)

True

Contrary to expectations, this transition does not conserve spin projection and therefore spin_conservation() returns False:

spin_conservation(
    ingoing_spins=[
        SpinEdgeInput(
            spin_magnitude=incoming_state.particle.spin,
            spin_projection=incoming_state.spin_projection,
        )
    ],
    outgoing_spins=[
        SpinEdgeInput(
            spin_magnitude=outgoing_state1.particle.spin,
            spin_projection=outgoing_state1.spin_projection,
        ),
        SpinEdgeInput(
            spin_magnitude=outgoing_state2.particle.spin,
            spin_projection=outgoing_state2.spin_projection,
        ),
    ],
    interaction_qns=interaction,
)

False

The reason is that AmpForm formulates the HelicityModel with the helicity formalism first and then uses a transformation to get the model in the canonical basis (see formulate_clebsch_gordan_coefficients()). The canonical basis does not conserve helicity (taken to be State.spin_projection).

Modifying `StateTransition`s¶

When checking conservation rules, you may want to modify the properties on the StateTransitions. However, a StateTransition is frozen, so it is not possible to modify its interactions and states. The only way around this is to create a new instance with attr.evolve().

First, we get the instance (in this case one of the InteractionProperties) and substitute its InteractionProperties.l_magnitude:

new_interaction = attr.evolve(transition.interactions[node_id], l_magnitude=2)
new_interaction

InteractionProperties(
  l_magnitude=2,
  l_projection=None,
  s_magnitude=1.0,
  s_projection=None,
  parity_prefactor=None,
)

We then again use attr.evolve() to substitute the StateTransition.interactions of the original StateTransition:

new_interaction_dict = dict(transition.interactions)  # make mutable
new_interaction_dict[node_id] = new_interaction
new_transition = attr.evolve(transition, interactions=new_interaction_dict)

dot = qrules.io.asdot(new_transition, render_node=True)
graphviz.Source(dot)

Custom topologies¶

As illustrated in Generate transitions, the StateTransitionManager offers you a bit more flexibility than the façade function generate_transitions() used in the main Usage page. In this notebook, we go one step further, by specifying a custom Topology via StateTransitionManager.topologies.

import graphviz

import qrules
from qrules import InteractionType, StateTransitionManager
from qrules.topology import Edge, Topology

2-to-2 topology¶

As a simple example, we start with a 2-to-2 scattering topology. We define it as follows:

topology = Topology(
    nodes=range(2),
    edges=enumerate(
        [
            Edge(None, 0),
            Edge(None, 0),
            Edge(1, None),
            Edge(1, None),
            Edge(0, 1),
        ],
        -2,
    ),
)

dot = qrules.io.asdot(
    topology,
    render_resonance_id=True,
    render_node=True,
    render_initial_state_id=True,
)
graphviz.Source(dot)

First, we construct a StateTransitionManager for the transition \(K^-K^+ \to \pi^+\pi^-\). The constructed Topology can then be inserted via its topologies attribute:

stm = StateTransitionManager(
    initial_state=["K-", "K+"],
    final_state=["pi-", "pi+"],
    formalism="canonical",
)
stm.set_allowed_interaction_types([InteractionType.STRONG, InteractionType.EM])
stm.topologies = (topology,)  # tuple is immutable

For the rest, the process is just the same as in Generate transitions:

problem_sets = stm.create_problem_sets()
reaction_kk = stm.find_solutions(problem_sets)

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

Warning

It is not yet possible to give the initial state a certain energy. So some collider process like \(e^-e^+\to\pi^+\pi\) does not result in a large number of resonances.

stm.initial_state = ["e-", "e+"]
problem_sets = stm.create_problem_sets()
reaction_ep = stm.find_solutions(problem_sets)

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

What can do at most, is switch off MassConservation, either through the constructor of the StateTransitionManager, or by modifying ProblemSet.

stm = StateTransitionManager(
    initial_state=["e-", "e+"],
    final_state=["pi-", "pi+"],
    formalism="canonical",
    mass_conservation_factor=None,
)
stm.set_allowed_interaction_types([InteractionType.STRONG, InteractionType.EM])
stm.topologies = [topology]
problem_sets = stm.create_problem_sets()
reaction_ep_no_mass = stm.find_solutions(problem_sets)

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

Bibliography¶

Tip

Download this bibliography as BibTeX here.

1: D. M. Asner. Dalitz Plot Analysis Formalism. In Review of Particle Physics: Volume I Reviews. January 2006. pdg.lbl.gov/2010/reviews/rpp2010-rev-dalitz-analysis-formalism.pdf.
2: S.-U. Chung et al. Partial wave analysis in 𝐾-matrix formalism. Annalen der Physik, 507(5):404–430, May 1995. doi:10.1002/andp.19955070504.

qrules¶

import qrules

A rule based system that facilitates particle reaction analysis.

QRules generates allowed particle transitions from a set of conservation rules and boundary conditions as specified by the user. The user boundary conditions for a particle reaction problem are for example the initial state, final state, and allowed interactions.

The core of qrules computes which transitions (represented by a StateTransitionGraph) are allowed between a certain initial and final state. Internally, the system propagates the quantum numbers defined by the particle module through the StateTransitionGraph, while satisfying the rules define by the conservation_rules module. See Generate transitions and Particle database.

Finally, the io module provides tools that can read and write the objects of this framework.

check_reaction_violations(initial_state: Union[str, Tuple[str, Sequence[float]], Sequence[Union[str, Tuple[str, Sequence[float]]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], mass_conservation_factor: Optional[float] = 3.0, particle_db: Optional[ParticleCollection] = None, max_angular_momentum: int = 1, max_spin_magnitude: float = 2.0) → Set[FrozenSet[str]][source]¶

Determine violated interaction rules for a given particle reaction.

Warning

This function only guarantees to find P, C and G parity violations, if it’s a two body decay. If all initial and final states have the C/G parity defined, then these violations are also determined correctly.

Parameters

initial_state – Shortform description of the initial state w/o spin projections.
final_state – Shortform description of the final state w/o spin projections.
mass_conservation_factor – Factor with which the width is multiplied when checking for MassConservation. Set to None in order to deactivate mass conservation.
particle_db (Optional) – Custom ParticleCollection object. Defaults to the ParticleCollection returned by load_pdg.
max_angular_momentum – Maximum angular momentum over which to generate \(LS\)-couplings.
max_spin_magnitude – Maximum spin magnitude over which to generate \(LS\)-couplings.

Returns

Set of least violating rules. The set can have multiple entries, as several quantum numbers can be violated. Each entry in the frozenset represents a group of rules that together violate all possible quantum number configurations.

Example

>>> import qrules
>>> qrules.check_reaction_violations(
...     initial_state="pi0",
...     final_state=["gamma", "gamma", "gamma"],
... )
{frozenset({'c_parity_conservation'})}

See also

Check allowed reactions

generate_transitions(initial_state: Union[str, Tuple[str, Sequence[float]], Sequence[Union[str, Tuple[str, Sequence[float]]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], allowed_intermediate_particles: Optional[List[str]] = None, allowed_interaction_types: Optional[Union[str, Iterable[str]]] = None, formalism: str = 'canonical-helicity', particle_db: Optional[ParticleCollection] = None, mass_conservation_factor: Optional[float] = 3.0, max_angular_momentum: int = 2, max_spin_magnitude: float = 2.0, topology_building: str = 'isobar', number_of_threads: Optional[int] = None) → ReactionInfo[source]¶

Generate allowed transitions between an initial and final state.

Serves as a facade to the StateTransitionManager (see Generate transitions).

Parameters

initial_state (list) – A list of particle names in the initial state. You can specify spin projections for these particles with a tuple, e.g. ("J/psi(1S)", [-1, 0, +1]). If spin projections are not specified, all projections are taken, so the example here would be equivalent to "J/psi(1S)".
final_state (list) – Same as initial_state, but for final state particles.
allowed_intermediate_particles (list, optional) – A list of particle states that you want to allow as intermediate states. This helps (1) filter out resonances and (2) speed up computation time.
allowed_interaction_types – Interaction types you want to consider. For instance, ["s", "em"] results in EM and STRONG and ["strong"] results in STRONG.
formalism (str, optional) – Formalism that you intend to use in the eventual amplitude model.
particle_db (ParticleCollection, optional) – The particles that you want to be involved in the reaction. Uses load_pdg by default. It’s better to use a subset for larger reactions, because of the computation times. This argument is especially useful when you want to use your own particle definitions (see Particle database).
mass_conservation_factor – Width factor that is taken into account for for the MassConservation rule.
max_angular_momentum – Maximum angular momentum over which to generate angular momenta.
max_spin_magnitude – Maximum spin magnitude over which to generate spins.
topology_building (str) –
Technique with which to build the Topology instances. Allowed values are:
- "isobar": Isobar model (each state decays into two states)
- "nbody": Use one central node and connect initial and final states to it
number_of_threads (int) – Number of cores with which to compute the allowed transitions. Defaults to all cores on the system.

An example (where, for illustrative purposes only, we specify all arguments) would be:

>>> import qrules
>>> reaction = qrules.generate_transitions(
...     initial_state="D0",
...     final_state=["K~0", "K+", "K-"],
...     allowed_intermediate_particles=["a(0)(980)", "a(2)(1320)-"],
...     allowed_interaction_types=["e", "w"],
...     formalism="helicity",
...     particle_db=qrules.load_pdg(),
...     topology_building="isobar",
... )
>>> len(reaction.transition_groups)
3
>>> len(reaction.transitions)
4

load_default_particles() → ParticleCollection[source]¶

Load the default particle list that comes with qrules.

Runs load_pdg and supplements its output definitions from the file additional_definitions.yml.

Submodules and Subpackages

io¶

import qrules.io

Serialization module for the qrules.

The io module provides tools to export or import objects from qrules to and from disk, so that they can be used by external packages, or just to store (cache) the state of the system.

asdict(instance: object) → dict[source]¶

asdot(instance: object, *, render_node: bool = False, render_final_state_id: bool = True, render_resonance_id: bool = False, render_initial_state_id: bool = False, strip_spin: bool = False, collapse_graphs: bool = False) → str[source]¶

Convert a object to a DOT language str.

Only works for objects that can be represented as a graph, particularly a StateTransitionGraph or a list of StateTransitionGraph instances.

Parameters

instance – the input object that is to be rendered as DOT (graphviz) language.
strip_spin – Normally, each StateTransitionGraph has a Particle with a spin projection on its edges. This option hides the projections, leaving only Particle names on edges.
collapse_graphs – Group all transitions by equivalent kinematic topology and combine all allowed particles on each edge.
render_node –
Whether or not to render node ID (in the case of a Topology) and/or node properties (in the case of a StateTransitionGraph). Meaning of the labels:
- \(P\): parity prefactor
- \(s\): tuple of coupled spin magnitude and its projection
- \(l\): tuple of angular momentum and its projection
See InteractionProperties for more info.
render_final_state_id – Add edge IDs for the final state edges.
render_resonance_id – Add edge IDs for the intermediate state edges.
render_initial_state_id – Add edge IDs for the initial state edges.

See also

Visualize decay topologies

fromdict(definition: dict) → object[source]¶

load(filename: str) → object[source]¶

write(instance: object, filename: str) → None[source]¶

argument_handling¶

import qrules.argument_handling

Handles argument handling for rules.

Responsibilities are the check of requirements for rules and the creation of the arguments from general graph property maps. The information is extracted from the type annotations of the rules.

class RuleArgumentHandler[source]¶

Welcome to QRules!¶

Installation¶

Usage¶

Main interface¶

Investigate intermediate resonances¶

Quantum number search¶

Check allowed reactions¶

Advanced¶

Generate transitions¶

1. Define the problem set¶

2. Prepare Problem Sets¶

3. Find solutions¶

4. Export generated transitions¶

Particle database¶

Loading the default database¶

Finding particles¶

LaTeX representation¶

Adding custom particle definitions through Python¶

Loading custom definitions from a YAML file¶

Writing to YAML¶

Visualize decay topologies¶

Topologies¶

ProblemSets¶

StateTransitions¶

Collapse graphs¶

Conservation rules¶

Parity conservation¶

Spin conservation¶

No spin and angular momentum¶

Non-zero example¶

Example with a StateTransition¶

Modifying StateTransitions¶

Custom topologies¶

2-to-2 topology¶

Bibliography¶

qrules¶

io¶

argument_handling¶

combinatorics¶

conservation_rules¶

particle¶

quantum_numbers¶

settings¶

solving¶

topology¶

transition¶

version¶

`ProblemSet`s¶

`StateTransition`s¶

Example with a `StateTransition`¶

Modifying `StateTransition`s¶