Welcome to QRules!#

10.5281/zenodo.5526360 Supported Python versions Google Colab Binder

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:

  1. State transition graphs

    A MutableTransition 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\).

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

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

  1. Preparation

    1.1. Build all possible topologies. A topology is represented by a MutableTransition, 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).

  2. 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#

PyPI package Conda package Supported Python versions

Quick installation#

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

python3 -m pip install qrules
conda install -c conda-forge qrules

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

Optionally, you can install the dependencies required for visualizing topologies with the following optional dependency syntax:

pip install qrules[viz]  # installs qrules with graphviz

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

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

Editable installation#

It is highly recommend to use the more dynamic ‘editable installation’. This allows you to:

For this, you first need to get the source code with Git:

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

Next, you install the project in editable mode with either Conda or pip. It’s recommended to use Conda, because this also pins the version of Python.

conda env create

This installs the project in a Conda environment following the definitions in environment.yml.

  1. [Recommended] Create a virtual environment with venv (see here).

  2. Install the project as an ‘editable installation’ with additional packages for the developer and all dependencies pinned through constraints files:

    python3 -m pip install -c .constraints/py3.x.txt -e .[dev]
    

See Updating for how to update the dependencies when new commits come in.

That’s all! Have a look at Usage to try out the package. You can also have a look at Help developing 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)
_images/755de53f4668ca9454b15cbe2773b3e45a31d8eb7242e21298655b3b151fae1b.svg

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!

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",
    max_angular_momentum=2,
)

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)
_images/d903399fda3cad542916bec74fafe7fc7e1bd6d8222648e5f723ce1d7970f73c.svg

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 396 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(0)(1500)',
 "f(2)'(1525)",
 'f(0)(1710)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(782)',
 'omega(1650)',
 'phi(1020)',
 'phi(1680)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)0']

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 192 solutions!
['b(1)(1235)0',
 'f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(0)(1370)',
 'f(0)(1500)',
 "f(2)'(1525)",
 'f(0)(1710)',
 'f(2)(1950)',
 '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 306 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(0)(1500)',
 "f(2)'(1525)",
 'f(0)(1710)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(782)',
 'omega(1650)',
 'phi(1020)',
 'phi(1680)',
 'rho(770)0',
 'rho(1450)0',
 'rho(1700)0']

Be aware, however, that the EM interaction is now available globally, 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)
Hide code cell source
dot = io.asdot(lc2pkpi_reaction, collapse_graphs=True)
graphviz.Source(dot)
_images/0383732071c513e119b53fb43f5199908633bcaa982e6a0534702e6942996ed0.svg

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
Hide code cell output
['f(0)(500)',
 'f(0)(980)',
 'f(2)(1270)',
 'f(0)(1370)',
 'f(0)(1500)',
 "f(2)'(1525)",
 'f(0)(1710)',
 'f(2)(1950)',
 '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) == 11

Now we have selected all amplitudes that involve f states:

Hide code cell source
dot = io.asdot(reaction, collapse_graphs=True, render_node=False)
graphviz.Source(dot)
_images/cf345014b6ff271082ecc7fb401106f7b420721a5ea0c14770f75ab63812b545.svg
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}}}

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: 530

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(tau)~',
 'nu(e)~',
 'nu(tau)',
 'nu(mu)~',
 'nu(e)',
 'nu(mu)',
 '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='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,
  ),
  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}^{0}, \Sigma_{c}(2455)^{++}, \Sigma_{c}(2520)^{++}, \Sigma_{c}(2520)^{0}, \Sigma_{c}(2455)^{+}, \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)
536
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)
537
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.validators.validate()). If you have installed qrules as an Editable installation and use VSCode, your YAML particle list are checked automatically in the GUI.

Visualize solutions#

The io module allows you to convert MutableTransition, 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 MutableTransition 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.

Hide code cell content
import graphviz

import qrules
from qrules.particle import Spin
from qrules.topology import create_isobar_topologies, create_n_body_topology
from qrules.transition import State
topology = create_n_body_topology(2, 4)
graphviz.Source(qrules.io.asdot(topology, render_initial_state_id=True))
_images/fde1d867e8ea71d86f44f28aa57dd3f21a0163db17447e6ed67d29fcb8a3d27e.svg

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))
_images/61a2589985401a6892b71ecc3a98085016762a2206bb4c5848453fa854cb2b59.svg

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))
_images/3412df546f05b64af047fa81bbc9c540ebcacf00d855bd9c807ca11feb7cf405.svg

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))
_images/c2108fa9a892b20b51867851cf6c74e9ed763d6e68624700a14cbd6c41b20a93.svg
ProblemSets#

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.

from qrules.settings import InteractionType

stm = qrules.StateTransitionManager(
    initial_state=["J/psi(1S)"],
    final_state=["K0", "Sigma+", "p~"],
    formalism="canonical-helicity",
)
stm.set_allowed_interaction_types([InteractionType.STRONG, InteractionType.EM])
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]
problem_set = problem_sets[60.0][0]
dot = qrules.io.asdot(problem_set, render_node=True)
graphviz.Source(dot)
_images/69c797cca4d6d887d103eaa560d51992d50e1d9ba88a8d6b59a5a3ded0115a16.svg
Quantum number solutions#

As noted in 3. Find solutions, a ProblemSet can be fed to StateTransitionManager.find_solutions() directly to get a ReactionInfo object. ReactionInfo is a final result that consists of Particles, but in the intermediate steps, QRules works with sets of quantum numbers. One can inspect these intermediate generated quantum numbers by using find_quantum_number_transitions() and inspecting is output. Note that the resulting object is again a dict with strengths as keys and a list of solution as values.

qn_solutions = stm.find_quantum_number_transitions(problem_sets)
{strength: len(values) for strength, values in qn_solutions.items()}
{3600.0: 36, 60.0: 72, 1.0: 36}

The list of solutions consist of a tuple of a QNProblemSet (compare ProblemSets) and a QNResult:

strong_qn_solutions = qn_solutions[3600.0]
qn_problem_set, qn_result = strong_qn_solutions[0]
Hide code cell source
dot = qrules.io.asdot(qn_problem_set, render_node=True)
graphviz.Source(dot)
_images/e76ea8ad1a0671c5b876284aac5296883d6ca2425fad2e93fef036233a0e8f41.svg
Hide code cell source
dot = qrules.io.asdot(qn_result, render_node=True)
graphviz.Source(dot)
_images/78dc5f11685c561833724b20c22a7f9bb21f1825bb3c9a3827683b54d190b7ef.svg
StateTransitions#

After finding the Quantum number solutions, QRules finds Particle definitions that match these quantum numbers. All these steps are hidden in the convenience functions StateTransitionManager.find_solutions() and generate_transitions(). In the following, 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 MutableTransition 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)
_images/9d4867b4d850c8d11280c2be9e8cdadc8807e445fd99b56fb0036001b95736eb.svg

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)
_images/f7dc60d4308b06a24eb27cde4aa9d67716fa33fdae420a3f8281ef534a2d019a.svg

or, with stripped node properties:

dot = qrules.io.asdot(reaction.transitions[:3], strip_spin=True, render_node=True)
graphviz.Source(dot)
_images/e38278546a79a8018b3efeccb393aa2bb0538842965f9dc997aaedb63ce4b7c1.svg

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)
_images/1820c015d11c8e710db3bb255ab74572954fc381c70267279521b0f6d484b932.svg
Other state renderings#

The convert() method makes it possible to convert the types of its states. This for instance allows us to only render the spin states on in a Transition:

spin_transitions = sorted(
    {
        t.convert(lambda s: Spin(s.particle.spin, s.spin_projection))
        for t in reaction.transitions
    }
)
some_selection = spin_transitions[::67][:3]
dot = qrules.io.asdot(some_selection, render_node=True)
graphviz.Source(dot)
_images/82c05c92eac40c1f52802e4a85930038fa97bc8dd690ce9ae0a4bdea03cc6483.svg

Or any other properties of a State:

Hide code cell source
def render_mass(state: State, digits: int = 3) -> str:
    mass = round(state.particle.mass, digits)
    width = round(state.particle.width, digits)
    if width == 0:
        return str(mass)
    return f"{mass}±{width}"


mass_transitions = sorted(
    {
        t.convert(
            state_converter=render_mass,
            interaction_converter=lambda _: None,
        )
        for t in reaction.transitions
    }
)
dot = qrules.io.asdot(mass_transitions[::10])
graphviz.Source(dot)
_images/141a4ad26eea2a322785b4d3aa86ff297ea6a1037de42bea3db0a32f9f5e94a7.svg
Styling#

The asdot() function also takes Graphviz attributes. These can be used to modify the layout of the whole figure. Examples are the size, color, and fontcolor. Edges and nodes can be styled with edge_style and node_style respectively:

dot = qrules.io.asdot(
    reaction.transitions[0],
    render_node=True,
    size=12,
    bgcolor="white",
    edge_style={
        "color": "red",
        "arrowhead": "open",
        "fontcolor": "blue",
        "fontsize": 25,
    },
    node_style={
        "color": "gray",
        "penwidth": 2,
        "shape": "ellipse",
        "style": "dashed",
    },
)
display(graphviz.Source(dot))
_images/e8df3d8eb37c61d4af3f94fc013cc410e79c2f764e293b1c4b2e171bdff93fae.svg

Conservation rules#

Hide code cell content
import attrs
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 MutableTransitions, populates them with quantum numbers (edge properties representing states and nodes properties representing interactions), then checks whether the generated MutableTransitions 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).

Hide code cell source
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))
_images/c9ed721c4d0c45c809aa88e07291e812583b658d08bf9451882c03bfbf5c9d69.svg_images/ee9b5ff812c4113f2d7d0b64be322f07c7df633d9241b3e72345756c3add590c.svg

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 StateTransitions#

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

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

new_interaction = attrs.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 attrs.evolve() to substitute the Transition.interactions of the original StateTransition:

new_interaction_dict = dict(transition.interactions)  # make mutable
new_interaction_dict[node_id] = new_interaction
new_transition = attrs.evolve(transition, interactions=new_interaction_dict)
Hide code cell source
dot = qrules.io.asdot(new_transition, render_node=True)
graphviz.Source(dot)
_images/c6d95b5fab1fb8e66b47d2e41de99cbe025d0134a330fb1b30322f5559897bfe.svg

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.

Hide code cell content
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,
    ),
)
Hide code cell source
dot = qrules.io.asdot(
    topology,
    render_resonance_id=True,
    render_node=True,
    render_initial_state_id=True,
)
graphviz.Source(dot)
_images/9f7ecac02305547942d78e88c13de9e5738933ad23db910c3e1023ddc940fee6.svg

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)
Hide code cell source
dot = qrules.io.asdot(reaction_kk, collapse_graphs=True)
graphviz.Source(dot)
_images/4e975ab6bdccb8c500869e18a6b70ea1d1f5b41f7941e2a23be0587aab9ea40e.svg

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)
Hide code cell source
dot = qrules.io.asdot(reaction_ep, collapse_graphs=True)
graphviz.Source(dot)
_images/0c0607fb9926da357c27d7fb258825a056a7822da7436e4750e48b61752f21bc.svg

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)
Hide code cell source
dot = qrules.io.asdot(reaction_ep_no_mass, collapse_graphs=True)
graphviz.Source(dot)
_images/e45d0323c1445706fa2817fb88e7f5f2cf5ffefa8d561b52b5497dd5e7928f1a.svg

LS-couplings#

The spin_conservation() rule is one of the more complicated checks in the conservation_rules module. It provides an implementation of \(LS\)-couplings, which is a procedure to determine which values for total angular momentum \(L\) and coupled spin \(S\) are allowed in an interaction node. In this notebook, we illustrate this procedure with the following decay chain as an example:

\[ J/\psi \to \Sigma^+ \bar\Sigma(1670)^-, \quad \bar\Sigma(1670)^- \rightarrow \bar p K^0. \]

In this decay chain, there are two decay nodes that we investigate separately. In addition, both decays are mediated interactions by the strong force, which means there is also parity conservation.

In the following derivations, the Particle.spin and Particle.parity values are of importance:

Hide code cell source
from IPython.display import Math

import qrules

PDG = qrules.load_pdg()
particle_names = [
    "J/psi(1S)",
    "Sigma+",
    "Sigma(1670)~-",
    "p~",
    "K0",
]
latex_expressions = []
for name in particle_names:
    particle = PDG[name]
    parity = "+" if particle.parity > 0 else "-"
    if particle.spin.is_integer():
        spin = int(particle.spin)
    else:
        nominator = int(particle.spin * 2)
        spin = Rf"\tfrac{{{nominator}}}{2}"
    latex_expressions.append(f"{particle.latex}={spin}^{parity}")
Math(R"\qquad ".join(latex_expressions))
\[\displaystyle J/\psi(1S)=1^-\qquad \Sigma^{+}=\tfrac{1}2^+\qquad \overline{\Sigma}(1670)^{-}=\tfrac{3}2^+\qquad \overline{p}=\tfrac{1}2^-\qquad K^{0}=0^-\]
Procedure#

Imagine we have a two-body decay of \(p_0\rightarrow p_1p_2\). We denote the Spin.magnitude of each particle \(p_i\) as \(s_i\) and their parity as \(\eta_i\). The values for \(L\) and \(S\) can now be determined as follows:

  1. Determine all values for \(S\) that satisfy \(\left| s_1-s_2 \right| \le S \le s_1+s_2\). The difference between each value for \(S\) has to integer, so \(S = \left| s_1-s_2 \right|, \left| s_1-s_2 \right|+1, \dots, s_1+s_2\).

  2. Determine all values for \(L\) that satisfy \(\left| L-S \right| \le s_0 \le L+S\), with \(L\) being a non-negative integer.

  3. If there is parity conservation, \(L\) has to satisfy an additional constraint: \(\eta_0 = \eta_1\cdot\eta_2\cdot(-1)^L\).

\(J/\psi \to \Sigma^+\bar\Sigma(1670)^-\)#

The spin and parity of each particle in the first transition can be summarized as \(1^-\to\frac{1}{2}^+\frac{3}{2}^+\). Following step 1 in the procedure, we get:

\[\begin{split} \begin{eqnarray} \left|s_{\Sigma^+} - s_{\bar\Sigma(1670)^-}\right| & \le S & \le s_{\Sigma^+} + s_{\bar\Sigma(1670)^-} \\ \left|\tfrac{1}{2}-\tfrac{3}{2}\right| & \le S & \le \tfrac{1}{2} + \tfrac{3}{2} \\ 1 & \le S & \le 2 \end{eqnarray} \end{split}\]
\[ \Rightarrow S=1 \quad \text{or} \quad S=2 \]

Next, we determine the allowed total angular momentum values \(L\) with step 2:

\[\begin{split} \begin{eqnarray} |L-S| & \le s_{J/\psi} & \le L+S \\ |L-S| & \le 1 & \le L+S \end{eqnarray} \end{split}\]
\[\begin{split} \Rightarrow \begin{cases} L=0,1,2 & \text{if} & S=1\\ L=1,2,3 & \text{if} & S=2. \end{cases} \end{split}\]

So in total, we have 6 \(LS\)-combinations:

\[ (L,S) = (0,1), (1,1), (2,1), (1,2), (2,2), (3,2). \]

This decay however goes via the strong force. This means that parity has to be conserved and we have to follow step 3:

\[\begin{split} \begin{eqnarray} \eta_{J/\psi} & = & \eta_{\Sigma^+} \cdot \eta_{\bar\Sigma(1670)^-} \cdot(-1)^L \\ (-1) & = & (+1)\cdot(+1)\cdot(-1)^L \\ (-1) & = & (-1)^{L}. \end{eqnarray} \end{split}\]

From this, we can easily see that only odd \(L\) values are possible, which leaves us with 3 \(LS\)-combinations:

\[ (L,S) = (1,1), (1,2), (3,2). \]
\(\bar \Sigma(1670)^-\to \bar pK^0\)#

The second part of the decay chain can be expressed as \(\frac{3}{2}^+ \to \frac{1}{2}^- 0^-\). Following step 1, we see:

\[\begin{split} \begin{eqnarray} |s_{\bar p} - s_{K^0}| & \le S & \le s_{\bar p} + s_{K^0} \\ \left|\tfrac{1}{2}-0 \right| & \le S & \le \tfrac{1}{2} + 0 \end{eqnarray} \end{split}\]
\[ \Rightarrow S = \tfrac{1}{2}. \]

This time, only one coupled spin value is allowed. That allows for the following values of \(L\):

\[\begin{split} \begin{eqnarray} |L-S| & \le s_0 & \le L+S \\ \left|L-\tfrac{1}{2}\right| & \le \tfrac{3}{2} & \le L+\tfrac{1}{2}. \end{eqnarray} \end{split}\]
\[ \Rightarrow L = 1,2 \]

By now, only two \(LS\)-combinations are possible:

\[ (L,S)=\left(1,\tfrac{1}{2}\right), \left(2,\tfrac{1}{2}\right). \]

This again is a strong interaction, which means we have to check for parity conservation.

\[\begin{split} \begin{eqnarray} \eta_{\bar \Sigma(1670)^-} & = & \eta_{\bar p}\cdot\eta_{K^0}\cdot(-1)^L\\ (+1) & = & (-1)\cdot(-1)\cdot(-1)^L\\ (+1) & = & (-1)^L. \end{eqnarray} \end{split}\]

Again, it is clear that only even \(L\)’s are allowed. This means that only one \(LS\)-combination is possible:

\[ (L,S)=\left(2,\tfrac{1}{2}\right) \]
Check with QRules#

Finally, let’s use generate_transitions() to check whether the allowed \(LS\)-couplings are found by qrules as well. Note that we have to increase the maximum angular momentum to find the \((L,S)=(3,2)\) combination as well.

Hide code cell source
import logging

import graphviz

LOGGER = logging.getLogger()
LOGGER.setLevel(logging.ERROR)

reaction = qrules.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", "Sigma+", "p~"],
    allowed_intermediate_particles=["Sigma(1670)"],
    allowed_interaction_types="strong",
    max_angular_momentum=3,
)
dot = qrules.io.asdot(reaction, render_node=True, strip_spin=True)
graphviz.Source(dot)
_images/fe42e3b12985d50dcaf2c014275fad29cc881c4a977b95efbf3ed26ff4d76a16.svg

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.