User Guide

This package provides a common interface to work with reactions and decay trees for several kinds of elements (string-based, PDG, …). It is therefore good separating the part corresponding to the treatment of the reactions and decays to that corresponding to the underlying elements.

Reactions and decays

Reactions and decays can be acessed through the reactions.make_reaction and reactions.make_decay functions, respectively. The syntax for both is very similar, and can be consulted in the syntax section.

import reactions
r = reactions.reaction('A B -> C D')

This reaction corresponds to two reactants A and B that generate two products C and D. By default, the elements of this reaction contain simply a string. We can create more complex chains of reactions:

r = reactions.reaction('A B -> {C D -> {E {F -> G H} -> G H I J}}')

Note that we have used curly braces to define the nested reactions. We can now define a function to explore the reaction and print the corresponding tree. This package provides the function reactions.is_element to distinguish when the node corresponds to an element and when it corresponds to a reaction or a decay.

def element_info(e, attributes):
    """ Build a string with the information of an element """
    attributes = attributes or ('name',)
    return ', '.join(f'{a}={getattr(e, a)}' for a in attributes)

def print_reaction(reaction, indent=0, attributes=None):
    """ Display a reaction on a simple way """
    prefix = indent * ' '
    print(f'{prefix}reactants:')
    for n in reaction.reactants:
        if reactions.is_element(n):
            print(f'{prefix} - {element_info(n, attributes)}')
        else:
            print_reaction(n, indent=indent + 2)
    print(f'{prefix}products:')
    for n in reaction.products:
        if reactions.is_element(n):
            print(f'{prefix} - {element_info(n, attributes)}')
        else:
            print_reaction(n, indent=indent + 2)

r = reactions.reaction('A B -> {C D -> {E {F -> G H} -> G H I J}}')

print_reaction(r)

reactants:
 - name=A
 - name=B
products:
  reactants:
   - name=C
   - name=D
  products:
    reactants:
     - name=E
      reactants:
       - name=F
      products:
       - name=G
       - name=H
    products:
     - name=G
     - name=H
     - name=I
     - name=J

A similar function can be defined for decays, that are composed by a head and a set of products:

def print_decay(decay, indent=0, attributes=None):
    """ Display a reaction on a simple way """
    prefix = indent * ' '
    print(f'{prefix}- {element_info(decay.head, attributes)}')
    print(f'{prefix}  products:')
    for n in decay.products:
        if reactions.is_element(n):
            print(f'{prefix}  - {element_info(n, attributes)}')
        else:
            print_decay(n, indent=indent + 2)

d = reactions.decay('A -> B {C -> D E} F')

print_decay(d)

- name=A
  products:
  - name=B
  - name=C
    products:
    - name=D
    - name=E
  - name=F

The comparison between reactions and decays is not sensitive to the orther specified in the reactants or the products, thus these expressions are all true:

assert(reactions.reaction('A B -> C D') == reactions.reaction('B A -> C D'))
assert(reactions.reaction('A B -> C D') == reactions.reaction('A B -> D C'))
assert(reactions.reaction('A B -> C D') == reactions.reaction('B A -> D C'))
assert(reactions.decay('A -> B C') == reactions.decay('A -> C B'))

However, note that we can not compare reactions with decays, and the comparison between objects must be done for the same underlying type. The kind of element can be specified on construction using the kind keyword argument, as can be seen in the following.

String elements

This is the default kind of element used when constructing reactions and decays. It has just one property, a string.

B = reactions.string_element('B')
assert(B.name == 'B')
r1 = reactions.reaction('A B -> C D')
assert(r.reactants[1] == B)
r2 = reactions.reaction('A B -> C D', kind='string')
assert(r.reactants[1] == B)
assert(r1 == r2)

This element can be used for simple operations, but is not useful for scientific applications.

PDG elements

This kind of elements are based on the information from the Particle Data Group (PDG). Their construction is dones through the reactions.pdg_database class, that acts as a service. This class is provided as a singleton, in such a way that any instance depending storing information of the PDG will access it through this class.

Constructing PDG elements

There are two ways of building this kind of elements: through the database or through the reactions.pdg_element constructor.

z0_db = reactions.pdg_database('Z0')
z0_el = reactions.pdg_element('Z0')
assert(z0_db == z0_el)

We can also access the elements using the PDG identification number.

z0_db = reactions.pdg_database(310)
z0_el = reactions.pdg_element(310)
assert(z0_db == z0_el)

The PDG elements contain information about the name and PDG ID (unique for each element), three times the charge, mass and width with their lower and upper errors, and whether the element is self charge-conjugate or not. Reactions and decays can be built with this particles providing the pdg value to the kind keyword argument.

decay = reactions.decay('Z0 -> mu+ mu-', kind='pdg')
print_decay(decay, attributes=['name', 'pdg_id', 'three_charge', 'mass', 'width', 'is_self_cc'])

- name=Z0, pdg_id=23, three_charge=0, mass=91.1876, width=2.4952, is_self_cc=True
  products:
  - name=mu+, pdg_id=-13, three_charge=3, mass=0.1056583745, width=2.9959837e-19, is_self_cc=False
  - name=mu-, pdg_id=13, three_charge=-3, mass=0.1056583745, width=2.9959837e-19, is_self_cc=False

The values of the mass and width depend on whether these have been measured by the experiments, so for certain particles this information is missing (also for their errors). To check if the information is available you can check whether the returned value is None.

assert(reactions.pdg_element('H0').width is None)

The full table used to construct this kind of elements can be consulted here.

Registering new PDG elements

It is possible to register custom elements in the database for later use.

# register the element from its initialization arguments
reactions.pdg_database.register_element('A', 99999999, 0, None, None, True)
A = reactions.pdg_database('A')

# directly register the element
B = reactions.pdg_element('B', 99999998, 0, None, None, True)
reactions.pdg_database.register_element(B)

print(element_info(A, attributes=['name', 'pdg_id', 'three_charge', 'is_self_cc']))
print(element_info(B, attributes=['name', 'pdg_id', 'three_charge', 'is_self_cc']))

name=A, pdg_id=99999999, three_charge=0, is_self_cc=True
name=B, pdg_id=99999998, three_charge=0, is_self_cc=True

There is one single condition that need to be satisfied in order to register an element, and is that none of the elements registered in the PDG database must have the same name or PDG ID.

Changing the energy units

It is possible to change the energy units used by the reactions.pdg_element classes with the use of the reactions.pdg_system_of_units object. This class is another singleton which determines the units to be used by all the PDG-related object. The PDG uses GeV units by default. If you want to change them, you simply need to provide the new units as a string.

z0_mass_gev = reactions.pdg_database('Z0').mass
reactions.pdg_system_of_units.set_energy_units('MeV')
z0_mass_mev = reactions.pdg_database('Z0').mass
assert(abs(z0_mass_gev - z0_mass_mev * 1e-3) < 1e-12)

Using the cache

The PDG database allows to use a cache, loading all the elements in memory to boost the access to them. You can do this through the reactions.pdg_database.enable_cache function. The cache can later be disabled using reactions.pdg_database.disable_cache. Note that this will not remove the elements registered by the user. If you wish to remove all elements you must call to reactions.pdg_database.clear_cache.

reactions.pdg_database('Z0') # this is taken from the cache
reactions.pdg_database.enable_cache() # load all particles
reactions.pdg_database('Z0') # this is taken from the cache
reactions.pdg_database.register_element("Z0'", 99999997, 0, None, None, True)
reactions.pdg_database('Z0') # this is taken from the cache

reactions.pdg_database.disable_cache()

reactions.pdg_database("Z0'") # our element is still there

reactions.pdg_database.clear_cache() # remove all elements in the cache

try:
    reactions.pdg_database("Z0'") # this will fail
    raise RutimeError('Should have failed before')
except reactions.LookupError as e:
    pass