The Particle & DecayLanguage packages
Eduardo Rodrigues
University of Liverpool
PyHEP 2022 Workshop
Online, 12-16 September 2022
Quick intro to the following packages¶
Particle- PDG particle data, MC identification codes, and more.DecayLanguage- Decay files (notably for EvtGen), universal description of decay chains.
We will be using also the little package hepunits for the HEP system of units.
Further reading - intermezzo on hepunits¶
The package hepunits collects the most commonly used units and constants in the
HEP System of Units, which are not the same as the international system of units (aka SI units).
The HEP system of units is based on the following:
| Quantity | Name | Unit |
|---|---|---|
| Length | millimeter | mm |
| Time | nanosecond | ns |
| Energy | Mega electron Volt | MeV |
| Positron charge | eplus | |
| Temperature | kelvin | K |
| Amount of substance | mole | mol |
| Luminous intensity | candela | cd |
| Plane angle | radian | rad |
| Solid angle | steradian | sr |
Note: no need to make use of sophisticated packages (e.g. as in AstroPy) since we basically never need to change systems of units (we never use ergs as energy, for example ;-)).
Basic usage is straightforward, though it may be confusing at first. Remember, all variables are written wrt to the units:
from hepunits import mm, ns, MeV, eplus, GeV, kelvin, mol, cd, rad, sr
mm == ns == MeV == eplus == kelvin == mol == cd == rad == sr == 1
GeV == 1000*MeV
Add two quantities with different length units:
from hepunits import units as u
1*u.meter + 5*u.cm
The result is in the unit of length, meaning mm!
PDG particle data, MC identification codes
Package motivation - particle data¶
- The PDG provides a downloadable table of particle masses, widths, charges and Monte Carlo particle ID numbers (PDG IDs).
- Most recent file here.
- It also provided an experimental file with extended information
(spin, quark content, P and C parities, etc.) until 2008 only, see here (not widely known!).
- But anyone wanting to use these data, the only readily available,
has to parse the file programmatically.
- Why not make a Python package to deal with all these data, for everyone?
Package motivation - MC identification codes¶
- The C++ HepPID and HepPDT libraries provide functions for processing particle ID codes
in the standard particle (aka PDG) numbering scheme.
- Different event generators may have their separate set of particle IDs: Geant3, etc.
- Again, why not make a package providing all functionality/conversions, Python-ically, for everyone?
Pythonic interface to¶
- Particle Data Group (PDG) particle data table.
- Particle MC identification codes, with inter-MC converters.
- With various extra goodies.
Package, in short¶
- Particle - loads extended PDG data tables and implements search and manipulations / display.
- PDGID - find out as much as possible from the PDG ID number. No table lookup.
- Converters for MC IDs used in Pythia and Geant.
- Basic usage via the command line.
- Fexible / advanced usage programmatically.
1. PDGID class and MC ID classes¶
- Classes
PDGID,PythiaID,Geant3ID. - Converters in module
particle.converters:Geant2PDGIDBiMap, etc.
PDG IDs module overview¶
- Process and query PDG IDs, and more – no look-up table needed.
- Current version of package reflects the latest version of the HepPID & HepPDT utility functions defined in the C++ HepPID and HepPDT versions 3.04.01
- It contains more functionality than that available in the C++ code … and minor fixes too.
- Definition of a PDGID class, PDG ID literals,
and set of standalone HepPID functions to query PDG IDs (is_meson, has_bottom, j_spin, charge, etc.).
- All PDGID class functions are available standalone.
PDGID class¶
- Wrapper class
PDGIDfor PDG IDs. - Behaves like an int, with extra goodies.
- Large spectrum of properties and methods, with a Pythonic interface, and yet more!
from particle import PDGID
pid = PDGID(211) # a pi+
pid
You will know if you create an invalid PDG ID ...
PDGID(99999999)
Class properties are mirroed from standalone functions.
from particle.pdgid import is_meson
pid.is_meson, is_meson(pid)
To print all PDGID properties:
print(pid.info())
Further reading - MC ID classes and converters¶
- Classes for MC IDs used in Pythia and Geant3:
PythiaIDandGeant3ID. - ID converters in module
particle.converters:Geant2PDGIDBiMap, etc.
from particle import PythiaID, Geant3ID
pyid = PythiaID(10221)
pyid.to_pdgid()
Conversions are directly available via mapping classes.
E.g., bi-directional map Pythia ID - PDG ID:
from particle.converters import Pythia2PDGIDBiMap
Pythia2PDGIDBiMap[PDGID(9010221)], Pythia2PDGIDBiMap[PythiaID(10221)]
2. Particle class¶
There are various ways to create a particle. The often used method is via its PDG ID.
from particle import Particle
Particle.from_pdgid(211)
Searching
Simple and natural API to deal with the PDG particle data table,
with powerful 1-line search and look-up utilities!
Particle.findall(…)– search a list of candidates.Particle.finditer(…)– search for a particle, returning an iterator of candidates.Search methods that can query any particle property!
Particle.findall('J/psi')
You can specify search terms as keywords - any particle property:
next(Particle.finditer(latex_name=r'\phi(1020)')) # first item in iterator of particles
You can directly check the numeric charge:
Particle.findall('pi', charge=-1)
Or use a lambda function for the ultimate in generality! For example, to find all the neutral particles with a bottom quark between 5.2 and 5.3 GeV:
from hepunits import GeV, s # Units are good. Use them.
Particle.findall(lambda p:
p.pdgid.has_bottom
and p.charge==0
and 5.2*GeV < p.mass < 5.3*GeV
)
Another lambda function example: You can use the width or the lifetime:
Particle.findall(lambda p: p.lifetime > 1000*s)
Trivially find all pseudoscalar charm mesons:
from particle import SpinType
Particle.findall(lambda p: p.pdgid.is_meson and p.pdgid.has_charm and p.spin_type==SpinType.PseudoScalar)
Display
Nice display in Jupyter notebooks, as well as str and repr support:
p = Particle.from_pdgid(-415)
p
print(p)
print(repr(p))
Full descriptions:
print(p.describe())
You may find LaTeX or HTML to be more useful in your program; both are supported:
print(p.latex_name, '\n', p.html_name)
Particle properties
You can do things to particles, like invert them:
~p
There are a plethora of properties you can access:
# p.<tab>
3. Further reading - literals¶
They provide a handy way to manipulate things with human-readable names!
Particle defines literals for most common particles, with easily recognisable names.
- Literals are dynamically generated on import for both
PDGIDandParticleclasses.
PDGID literals
from particle.pdgid import literals as lid
lid.phi_1020
Particle literals
from particle import literals as lp
lp.phi_1020
4. Data files, stored in particle/data/¶
PDG particle data files
- Original PDG data files, which are in a fixed-width format - simply for bookkeeping and reference.
- Code rather uses “digested forms” of these, produced within
Particle, stored as CSV, for optimised querying. - Latest PDG data (2020) used by default.
- Advanced usage: user can load older PDG tables, load a “user table” with new particles, append to default table.
Other data files
- CSV file for mapping of PDG IDs to particle LaTeX names.
Dump table contents
The package provides the 2 methods Particle.to_dict(...) and Particle.to_list(...), which make it easy to dump (selected) particle properties in an easy way. No need to dig into the package installation directory to inspect the particle data table ;-).
Tabular output can be formatted with the powerful package tabulate, for example (other similar libraries exist).
from tabulate import tabulate
# help(Particle.to_dict)
fields = ['pdgid', 'pdg_name', 'mass', 'mass_upper', 'mass_lower', 'three_charge']
query_as_dict = Particle.to_dict(exclusive_fields=fields, n_rows=10)
print(tabulate(query_as_dict, headers='keys'))
Be fancy - table with all pseudoscalar charm hadrons, in reStructuredText format:
fields = ['pdgid', 'name', 'evtgen_name', 'mass', 'mass_upper', 'mass_lower', 'three_charge']
query_as_dict = Particle.to_dict(filter_fn=lambda p: p.pdgid.is_meson and p.pdgid.has_charm and p.spin_type==SpinType.PseudoScalar,
exclusive_fields=fields)
print(tabulate(query_as_dict, headers='keys', tablefmt='rst'))
Notebook-friendly HTML is just as easy:
query_as_dict = Particle.to_dict(filter_fn=lambda p: p.pdgid.is_meson and p.pdgid.has_charm and p.spin_type==SpinType.PseudoScalar,
exclusive_fields=['pdgid', 'pdg_name', 'html_name'])
tabulate(query_as_dict, headers='keys', tablefmt='html')
5. Further reading - advanced usage¶
You can:
- Extend or replace the default particle data table in
Particle. - Adjust properties for a particle.
- Make custom particles.
I hope you enjoy the package :-).
Decay files, universal description of decay chains
Package motivation¶
- Any experiment uses event generators which, among many things, need to describe particle decay chains.
- Programs such as EvtGen rely on so-called .dec decay files.
- Many experiments need decay data files.
- Why not make a Python package to deal with decay files, for everyone?
Package, in short¶
DecayLanguage is designed for the manipulation of decay structures in Python
- Decay file parsers:
- Read .dec DecFiles, such as EvtGen decay files typically used in Flavour Physics experiments.
- Tools to parse decay files and programmatically manipulate them, query, display and visualise information.
- Descriptions and parsing built atop the Lark parser.
- Ability to describe decay-tree-like structures.
- Amplitude Analysis decay language:
- Tools to translate decay amplitude models from AmpGen to GooFit, and manipulate them.
- Input based on AmpGen generator, output format for GooFit C++ program.
1. Decay files¶
*Master file” DECAY.DEC¶
Gigantic file defining decay modes for all relevant particles, including decay model specifications. LHCb uses one. Belle II as well, and other collaborations.
User .dec files¶
- Needed to produce specific simulation samples.
- Typically contain a single decay chain (except if defining inclusive samples).
Example user decay file:
# Decay file for [B_c+ -> (B_s0 -> K+ K-) pi+]ccAlias B_c+sig B_c+ Alias B_c-sig B_c- ChargeConj B_c+sig B_c-sig Alias MyB_s0 B_s0 Alias Myanti-B_s0 anti-B_s0 ChargeConj MyB_s0 Myanti-B_s0
Decay B_c+sig 1.000 MyB_s0 pi+ PHOTOS PHSP; Enddecay CDecay B_c-sig
Decay MyB_s0 1.000 K+ K- SSD_CP 20.e12 0.1 1.0 0.04 9.6 -0.8 8.4 -0.6; Enddecay CDecay Myanti-B_s0
2. Decay file parsing¶
- Parsing should be simple and (reasonably) fast!
- Expert users can configure parser choice and settings, etc.
After parsing, many queries are possible!
from decaylanguage import DecFileParser
The LHCb "master decay file"¶
It's a big file! ~ 500 particle decays defined, thousands of decay modes, over 11k lines in total.
dfp = DecFileParser('data/DECAY_LHCB.DEC')
%%time
dfp.parse()
dfp
Let's parse and play with a small decay file:
with open('data/Dst.dec') as f:
print(f.read())
dfp_Dst = DecFileParser('data/Dst.dec')
dfp_Dst
dfp_Dst.parse()
dfp_Dst
It can be handy to parse from a multi-line string rather than a file:
s = """
# Decay file for [B_c+ -> (B_s0 -> K+ K-) pi+]cc
Alias B_c+sig B_c+
Alias B_c-sig B_c-
ChargeConj B_c+sig B_c-sig
Alias MyB_s0 B_s0
Alias Myanti-B_s0 anti-B_s0
ChargeConj MyB_s0 Myanti-B_s0
Decay B_c+sig
1.000 MyB_s0 pi+ PHOTOS PHSP;
Enddecay
CDecay B_c-sig
Decay MyB_s0
1.000 K+ K- SSD_CP 20.e12 0.1 1.0 0.04 9.6 -0.8 8.4 -0.6;
Enddecay
CDecay Myanti-B_s0
"""
dfp = DecFileParser.from_string(s)
dfp.parse()
dfp
Decay file information¶
There are several methods such as print_*, dict_*. Play around ...
For example:
dfp_Dst.list_decay_mother_names()
# dfp_Dst.print_decay_modes('D*+')
# Info such as charge conjugates
# dfp.dict_charge_conjugates()
3. Display of decay chains¶
The parser can provide a simple dict representation of any decay chain found in the input decay file(s). Being generic and simple, that is what is used as input information for the viewer class (see below).
dc = dfp_Dst.build_decay_chains('D+')
dc
from decaylanguage import DecayChainViewer
Let's visualise, as the above can quickly become unreadable to the human eye ...
DecayChainViewer(dc)
DecayChainViewer(dfp_Dst.build_decay_chains('D*+')) # Warning: decay chains can quickly become large !
dc = dfp_Dst.build_decay_chains('D*+', stable_particles=['D+', 'D0', 'pi0'])
DecayChainViewer(dc)
Charge conjugation¶
dc_cc = dfp_Dst.build_decay_chains('D*-', stable_particles=['D-', 'anti-D0', 'pi0'])
DecayChainViewer(dc_cc)
Want to save a graph?¶
Try for example
dcv = DecayChainViewer(...)
dcv.graph.render(filename='test', format='png', view=False, cleanup=True)
This saves the graph as "test.png", does not open the default application to view the file, and deletes the temporary .dot file.
Comment: if you reckon the command is non-trivial, understand that DecayLanguage merely delegates to graphviz, to minimise (package) couplings.
# DecayChainViewer(dc).graph.render(filename='test', format='png', view=False, cleanup=True)
Further reading - parsing several files¶
Typically useful when the user decay file needs information from the master decay file.
s = u"""
Alias MyXic+ Xi_c+
Alias MyantiXic- anti-Xi_c-
ChargeConj MyXic+ MyantiXic-
Decay Xi_cc+sig
1.000 MyXic+ pi- pi+ PHSP;
Enddecay
CDecay anti-Xi_cc-sig
Decay MyXic+
1.000 p+ K- pi+ PHSP;
Enddecay
CDecay MyantiXic-
End
"""
dfp = DecFileParser.from_string(s)
dfp.parse()
dfp
Note the subtletly: 3, not 4 decays, are found! This is because the file contains no statement
ChargeConj anti-Xi_cc-sigXi_cc+sig, hence the parser cannot know to which particle (matching Decay statement) the charge-conjugate decay of anti-Xi_cc-sig relates to (code does not rely on position of statements to guess ;-)).
d = dfp.build_decay_chains('Xi_cc+sig')
DecayChainViewer(d)
As said in the warning, the information provided is not enough for the anti-Xi_cc-sig to make sense:
from decaylanguage.dec.dec import DecayNotFound
try:
d = dfp.build_decay_chains('anti-Xi_cc-sig')
except DecayNotFound:
print("Decays of particle 'anti-Xi_cc-sig' not found in .dec file!")
But the missing information is easily providing parsing two files simultaneously ...! (Any number of files is allowed.)
from tempfile import NamedTemporaryFile
with NamedTemporaryFile(delete=False) as tf:
tf.write(s.encode('utf-8'))
dfp = DecFileParser(tf.name, 'data/DECAY_LHCB.DEC')
dfp.parse()
dc = dfp.build_decay_chains('Xi_cc+sig')
DecayChainViewer(dc)
dc_cc = dfp.build_decay_chains('anti-Xi_cc-sig')
DecayChainViewer(dc_cc)
4. Representation of decay chains¶
The universal (and digital) representation of decay chains is of interest well outside the context of decay file parsing!
Building blocks¶
- A daughters list - list of final-state particles.
- A decay mode - typically a branching fraction and a list of final-state particles (may also contain any metadata such as decay model and optional decay-model parameters, as defined for example in .dec decay files).
- A decay chain - can be seen as a mother particle and a list of decay modes.
from decaylanguage.decay.decay import DaughtersDict, DecayMode, DecayChain
Daughters list (actually a Counter dictionary, internally):
# Constructor from a dictionary
dd = DaughtersDict({'K+': 1, 'K-': 2, 'pi+': 1, 'pi0': 1})
# Constructor from a list of particle names
dd = DaughtersDict(['K+', 'K-', 'K-', 'pi+', 'pi0'])
# Constructor from a string representing the final state
dd = DaughtersDict('K+ K- pi0')
dd
Decay Modes¶
# A 'default' and hence empty, decay mode
dm = DecayMode()
# Decay mode with minimal input information
dd = DaughtersDict('K+ K-')
dm = DecayMode(0.5, dd)
# Decay mode with decay model information and user metadata
dm = DecayMode(0.2551, # branching fraction
'pi- pi0 nu_tau', # final-state particles
model='TAUHADNU', # decay model
model_params=[-0.108, 0.775, 0.149, 1.364, 0.400], # decay-model parameters
study='toy', year=2019 # user metadata
)
dm
print(dm.describe())
Various manipulations are available:
dm = DecayMode.from_pdgids(0.5, [321, -321])
dm
dm = DecayMode(1.0, 'K+ K+ pi-')
dm.charge_conjugate()
Decay chains¶
dm1 = DecayMode(0.0124, 'K_S0 pi0', model='PHSP')
dm2 = DecayMode(0.692, 'pi+ pi-')
dm3 = DecayMode(0.98823, 'gamma gamma')
dc = DecayChain('D0', {'D0':dm1, 'K_S0':dm2, 'pi0':dm3})
dc
dc.decays
Flatten the decay chain, i.e. replace all intermediate, decaying particles, with their final states:
- The BF is now the visible BF, i.e. the product of all involved decay BFs
dc.flatten()
Of course you can still just as easily visualise decays defined via this DecayChain class:
DecayChainViewer(dc.to_dict())