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)

NuBase elements

This kind of elements are based on the information of the NuBase database. The construction of elements is handled through the reactions.nubase_database and reactions.nubase_element objects, with a similar implementation to PDG elements.

Constructing NuBase elements

Similarly to PDG elements, there are two ways of building NuBase elements: through the database or through the reactions.nubase_element constructor.

proton_db_by_name = reactions.nubase_database('1H')
proton_el_by_name = reactions.nubase_element('1H')
assert(proton_db_by_name == proton_el_by_name)
proton_db_by_id = reactions.nubase_database(1001000)
proton_el_by_id = reactions.nubase_element(1001000)
assert(proton_db_by_id == proton_el_by_id)

assert(proton_db_by_id == proton_el_by_name)

A NuBase element contains information about its name, ID, whether the nucleus is stable or not, whether it is a ground state and the information about the mass excess and half-life with their corresponding errors, and whether these where obtained from systematics or not.

print(reactions.nubase_database('1H'))

reactions.nubase_element(name="1H", nubase_id=1001000, atomic_number=1, mass_number=1, mass_excess_and_error_with_tag=(value=7288.97, error=1.300000e-05, tag=False), is_stable=True, half_life_and_error_with_tag=None, is_ground_state=True)

The mass excess and half-life information can be missing, in which case the returned value associated to that quantity is None:

assert(reactions.nubase_element('1H').half_life is None)
assert(reactions.nubase_element('76Cu(m)').mass_excess is None)

Reactions and decays can be created providing nubase as the kind to reactions.reaction and reactions.decay.

decay = reactions.decay('1n -> 1H e-', kind='nubase')
print_decay(decay, attributes=['name', 'nubase_id', 'is_stable'])

- name=1n, nubase_id=1000000, is_stable=False
  products:
  - name=1H, nubase_id=1001000, is_stable=True
  - name=e-, nubase_id=1, is_stable=True

The full list of NuBase elements can be consulted here.

Registering new NuBase elements

# register the element from its initialization arguments
reactions.nubase_database.register_element('999Un', 999999000, 999, 999, None, True, None, True)
un999 = reactions.nubase_database('999Un')

# directly register the element
un998 = reactions.nubase_element('998Un', 999998000, 999, 998, None, False, None, True)
reactions.nubase_database.register_element(un998)

print(element_info(un999, attributes=['name', 'nubase_id', 'is_ground_state']))
print(element_info(un998, attributes=['name', 'nubase_id', 'is_ground_state']))

name=999Un, nubase_id=999999000, is_ground_state=True
name=998Un, nubase_id=999998000, is_ground_state=True

Newly registered elements must not have a similar name or ID to any existing element.

Changing the units

NuBase elements contain information about the mass excess, expressed with energy units (keV by default) and the half-life, expressed with time unit (seconds by default). It is possible to change both through the reactions.nubase_system_of_units singleton.

proton_mass_excess_kev = reactions.nubase_database('1H').mass_excess
reactions.nubase_system_of_units.set_energy_units('eV')
proton_mass_excess_ev = reactions.nubase_database('1H').mass_excess
assert(abs(proton_mass_excess_kev - proton_mass_excess_ev * 1e-3) < 1e-12)

neutron_half_life_sec = reactions.nubase_database('1n').half_life
reactions.nubase_system_of_units.set_time_units('ms')
neutron_half_life_ms = reactions.nubase_database('1n').half_life
assert(abs(neutron_half_life_sec - neutron_half_life_ms * 1e-3) < 1e-12)

Using the database cache

The PDG and NuBase databases 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 and reactions.nubase_database.enable_cache functions. The cache can later be disabled using reactions.pdg_database.disable_cache and reactions.nubase_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 or reactions.nubase_database.clear_cache. The following example shows how to use it for PDG elements.

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