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
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\).
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
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).
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#
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:
exactly pin all dependencies to a specific version, so that your work is reproducible.
edit the source code of the framework and help improving it.
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
.
[Recommended] Create a virtual environment with
venv
(see here).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)
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!
See also
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
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,
)
State Transition Manager
The StateTransitionManager
(STM) is the main user interface class of qrules
. The boundary conditions of your physics problem, such as the initial state, final state, formalism type, etc., are defined through this interface.
create_problem_sets()
of the STM creates all problem sets ― using the boundary conditions of theStateTransitionManager
instance. In total 3 steps are performed. The creation of reaction topologies. The creation ofInitialFacts
, based on a topology and the initial and final state information. And finally the solving settings such as the conservation laws and quantum number domains to use at which point of the topology.By default, all three interaction types (
EM
,STRONG
, andWEAK
) are used in the preparation stage. However, it is also possible to choose the allowed interaction types globally viaset_allowed_interaction_types()
.After the preparation step, you can modify the problem sets returned by
create_problem_sets()
to your liking. Since the output of this function contains quite a lot of information,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 instanceallowed_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 ProblemSet
s represent, you can use asdot()
and Graphviz to visualize one of them (see Visualize solutions):
import graphviz
from qrules import io
some_problem_set = problem_sets[60.0][0]
dot = io.asdot(some_problem_set, render_node=True)
graphviz.Source(dot)
Each ProblemSet
provides a mapping of initial_facts
that represent the initial and final states with spin projections. The nodes and edges in between these initial_facts
are still to be generated. This will be done from the provided solving_settings
. There are two mechanisms there:
One the one hand, the
EdgeSettings.qn_domains
andNodeSettings.qn_domains
contained in thesolving_settings
define the domain over which quantum number sets can be generated.On the other, the
EdgeSettings.rule_priorities
andNodeSettings.rule_priorities
insolving_settings
define whichconservation_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 MutableTransition
s 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']
About the number of solutions
The “number of transitions
” is the total number of allowed MutableTransition
instances that the StateTransitionManager
has found. This also includes all allowed spin projection combinations. In this channel, we for example consider a \(J/\psi\) with spin projection \(\pm1\) that decays into a \(\gamma\) with spin projection \(\pm1\), which already gives us four possibilities.
On the other hand, the intermediate state names that was extracted with ReactionInfo.get_intermediate_particles()
, is just a set
of the state names on the intermediate edges of the list of transitions
, regardless of spin projection.
Now we have a lot of solutions that are actually heavily suppressed (they involve two weak decays).
In general, you can modify the ProblemSet
s returned by create_problem_sets()
directly, but the STM also comes with functionality to globally choose the allowed interaction types. So, go ahead and disable the EM
and 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)
Show code cell source
dot = io.asdot(lc2pkpi_reaction, collapse_graphs=True)
graphviz.Source(dot)
Great! Now we selected only the strongest contributions. Be aware, though, that there are more effects that can suppress certain decays, like small branching ratios. In this example, the initial state \(J/\Psi\) can decay into \(\pi^0 + \rho^0\) or \(\pi^0 + \omega\).
decay |
branching ratio |
---|---|
\(\omega\to\gamma+\pi^0\) |
0.0828 |
\(\rho^0\to\gamma+\pi^0\) |
0.0006 |
Unfortunately, the \(\rho^0\) mainly decays into \(\pi^0+\pi^0\), not \(\gamma+\pi^0\) and is therefore suppressed. This information is currently not known to qrules
, but it is possible to hand qrules
a list of allowed intermediate states,
stm.set_allowed_intermediate_particles(["f(0)", "f(2)"])
reaction = stm.find_solutions(problem_sets)
reaction.get_intermediate_particles().names
Show code cell output
['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:
Show code cell source
dot = io.asdot(reaction, collapse_graphs=True, render_node=False)
graphviz.Source(dot)
See also
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’ Particle
s 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#
Particle
s 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]))
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 ProblemSet
s 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.
Show 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))
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 ProblemSet
s 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 list
s of ProblemSet
s 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)
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 Particle
s, 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]
Show code cell source
dot = qrules.io.asdot(qn_problem_set, render_node=True)
graphviz.Source(dot)
Show code cell source
dot = qrules.io.asdot(qn_result, render_node=True)
graphviz.Source(dot)
StateTransition
s#
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)
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)
or, with stripped node properties:
dot = qrules.io.asdot(reaction.transitions[:3], strip_spin=True, render_node=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)
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)
Or any other properties of a State
:
Show 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)
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))
Conservation rules#
Show 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 MutableTransition
s, populates them with quantum numbers (edge properties representing states and nodes properties representing interactions), then checks whether the generated MutableTransition
s 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#
parity_conservation(
ingoing_edge_qns=[Parity(-1)],
outgoing_edge_qns=[Parity(+1), Parity(+1)],
l_magnitude=1,
)
True
Spin conservation#
See also
spin_conservation()
, tests/unit/conservation_rules/test_spin.py
, PDG2020, §Quark Model, and these lecture notes by Curtis Meyer.
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 StateTransition
s 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 tuple
s of a Particle
and a float
spin projection).
Show 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))
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 StateTransition
s. 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)
Show code cell source
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
.
Show 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,
),
)
Show code cell source
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)
Show code cell source
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)
Show code cell source
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)
Show code cell source
dot = qrules.io.asdot(reaction_ep_no_mass, collapse_graphs=True)
graphviz.Source(dot)
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:
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:
Show 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))
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:
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\).
Determine all values for \(L\) that satisfy \(\left| L-S \right| \le s_0 \le L+S\), with \(L\) being a non-negative integer.
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:
Next, we determine the allowed total angular momentum values \(L\) with step 2:
So in total, we have 6 \(LS\)-combinations:
This decay however goes via the strong force. This means that parity has to be conserved and we have to follow step 3:
From this, we can easily see that only odd \(L\) values are possible, which leaves us with 3 \(LS\)-combinations:
\(\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:
This time, only one coupled spin value is allowed. That allows for the following values of \(L\):
By now, only two \(LS\)-combinations are possible:
This again is a strong interaction, which means we have to check for parity conservation.
Again, it is clear that only even \(L\)’s are allowed. This means that only one \(LS\)-combination is possible:
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.
Show 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)
Bibliography#
Tip
Download this bibliography as BibTeX here
.
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.
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.