10: Scikit-HEP

The Scikit-HEP project is a well maintained collection of HEP packages in Python that are useful for analysis. It contains tools ranging from plotting helpers, PDG lookup, DecayLanguage converters over high-performance histogramming libraries and ROOT I/O up to likelihood fitting and statistics libraries.

This is a minimal overview over the packages that are available and have not yet been used in the other tutorials.

formulate - converting expressions

The package formulate converts between different style of expressions, notably between ROOT and numexpr. The former was used in the previous notebooks to apply cuts and create new variables with Pandas DataFrames query and eval methods.

The standard usage is with a from_style and to_style, but it allows for more as it parses the string such as the variables used can be extracted.

import formulate

momentum = formulate.from_root('TMath::Sqrt(X_PX**2 + X_PY**2 + X_PZ**2)')
momentum

Expression<SQRT>(Expression<ADD>(Expression<POW>(Variable(X_PX), UnnamedConstant(2)), Expression<POW>(Variable(X_PY), UnnamedConstant(2)), Expression<POW>(Variable(X_PZ), UnnamedConstant(2))))

momentum.to_numexpr()  # as used in Pandas eval/query

'sqrt(((X_PX ** 2) + (X_PY ** 2) + (X_PZ ** 2)))'

momentum.to_root()  # as used in ROOT

'TMath::Sqrt(((X_PX ** 2) + (X_PY ** 2) + (X_PZ ** 2)))'

Particle

particle provides an easy access to particle properties such as mass and decay width from PDG. One example is to get the actual particle in a Particle object which holds all the relevant information.

# Particle

from particle import Particle

piplus = Particle.from_pdgid(211)

Here we can access all the properties we wish to such as mass and decay width.

piplus.mass

139.57039

piplus.charge

1.0

piplus.width

2.5284e-14

piplus.name

'pi+'

It also serves as a data-base like structure and we can do searches. Let’s find all neutral beauty hadrons.

Particle.findall(lambda p: p.pdgid.has_bottom and p.charge==0)

[<Particle: name="B0", pdgid=511, mass=5279.66 ± 0.12 MeV>,
 <Particle: name="B~0", pdgid=-511, mass=5279.66 ± 0.12 MeV>,
 <Particle: name="B*0", pdgid=513, mass=5324.71 ± 0.21 MeV>,
 <Particle: name="B*~0", pdgid=-513, mass=5324.71 ± 0.21 MeV>,
 <Particle: name="B(2)*(5747)0", pdgid=515, mass=5739.5 ± 0.7 MeV>,
 <Particle: name="B(2)*(5747)~0", pdgid=-515, mass=5739.5 ± 0.7 MeV>,
 <Particle: name="B(s)0", pdgid=531, mass=5366.92 ± 0.10 MeV>,
 <Particle: name="B(s)~0", pdgid=-531, mass=5366.92 ± 0.10 MeV>,
 <Particle: name="B(s)*0", pdgid=533, mass=5415.4 + 1.8 - 1.5 MeV>,
 <Particle: name="B(s)*~0", pdgid=-533, mass=5415.4 + 1.8 - 1.5 MeV>,
 <Particle: name="B(s2)*(5840)0", pdgid=535, mass=5839.86 ± 0.12 MeV>,
 <Particle: name="B(s2)*(5840)~0", pdgid=-535, mass=5839.86 ± 0.12 MeV>,
 <Particle: name="eta(b)(1S)", pdgid=551, mass=9398.7 ± 2.0 MeV>,
 <Particle: name="Upsilon(1S)", pdgid=553, mass=9460.40 ± 0.10 MeV>,
 <Particle: name="chi(b2)(1P)", pdgid=555, mass=9912.2 ± 0.4 MeV>,
 <Particle: name="Lambda(b)0", pdgid=5122, mass=5619.60 ± 0.17 MeV>,
 <Particle: name="Lambda(b)~0", pdgid=-5122, mass=5619.60 ± 0.17 MeV>,
 <Particle: name="Xi(b)0", pdgid=5232, mass=5791.9 ± 0.5 MeV>,
 <Particle: name="Xi(b)~0", pdgid=-5232, mass=5791.9 ± 0.5 MeV>,
 <Particle: name="chi(b0)(1P)", pdgid=10551, mass=9859.4 ± 0.5 MeV>,
 <Particle: name="h(b)(1P)", pdgid=10553, mass=9899.3 ± 0.8 MeV>,
 <Particle: name="chi(b1)(1P)", pdgid=20553, mass=9892.8 ± 0.4 MeV>,
 <Particle: name="Upsilon(2)(1D)", pdgid=20555, mass=10163.7 ± 1.4 MeV>,
 <Particle: name="Upsilon(2S)", pdgid=100553, mass=10023.4 ± 0.5 MeV>,
 <Particle: name="chi(b2)(2P)", pdgid=100555, mass=10268.6 ± 0.5 MeV>,
 <Particle: name="chi(b0)(2P)", pdgid=110551, mass=10232.5 ± 0.6 MeV>,
 <Particle: name="h(b)(2P)", pdgid=110553, mass=10259.8 ± 1.2 MeV>,
 <Particle: name="chi(b1)(2P)", pdgid=120553, mass=10255.5 ± 0.5 MeV>,
 <Particle: name="Upsilon(3S)", pdgid=200553, mass=10355.1 ± 0.5 MeV>,
 <Particle: name="chi(b2)(3P)", pdgid=200555, mass=10524.0 ± 0.8 MeV>,
 <Particle: name="chi(b1)(3P)", pdgid=220553, mass=10513.4 ± 0.7 MeV>,
 <Particle: name="Upsilon(4S)", pdgid=300553, mass=10579.4 ± 1.2 MeV>,
 <Particle: name="Upsilon(10860)", pdgid=9000553, mass=10885.2 + 2.6 - 1.6 MeV>,
 <Particle: name="Upsilon(11020)", pdgid=9010553, mass=11000 ± 4 MeV>]

hepunits

Ever hardcoded a physics constant? Or needed to manually convert units? That’s what the neat hepunits package is for.

It contains constants as well as units that all get converted into the standard units (e.g. MeV) which allows for calculations without the need to explicitly convert by just adding the unit; furthermore it greatly increases readability

from hepunits.constants import c_light

c_light

299.792458

import hepunits as u  # u for units

150 * u.MeV + 1.1 * u.GeV  # result in MeV

1250.0

Vector

Manipulating and calculation vector quantities can be done using the vector library. It provides functionality to load, combine and work with vectors.

import vector

# NumPy-like arguments (literally passed through to NumPy)
vector.array([
    (1.1, 2.1), (1.2, 2.2), (1.3, 2.3), (1.4, 2.4), (1.5, 2.5)
], dtype=[("x", float), ("y", float)])

# Pandas-like arguments (dict from names to column arrays)
vector.array({"x": [1.1, 1.2, 1.3, 1.4, 1.5], "y": [2.1, 2.2, 2.3, 2.4, 2.5]})

# As with objects, the coordinate system and dimension is taken from the names of the fields.
vec1 = vector.array({
    "x": [1.1, 1.2, 1.3, 1.4, 1.5],
    "y": [2.1, 2.2, 2.3, 2.4, 2.5],
    "z": [3.1, 3.2, 3.3, 3.4, 3.5],
    "t": [4.1, 4.2, 4.3, 4.4, 4.5],
})

vector.array({
    "pt": [1.1, 1.2, 1.3, 1.4, 1.5],
    "phi": [2.1, 2.2, 2.3, 2.4, 2.5],
    "eta": [3.1, 3.2, 3.3, 3.4, 3.5],
    "M": [4.1, 4.2, 4.3, 4.4, 4.5],
})

MomentumNumpy4D([(1.1, 2.1, 3.1, 4.1), (1.2, 2.2, 3.2, 4.2), (1.3, 2.3, 3.3, 4.3),
                 (1.4, 2.4, 3.4, 4.4), (1.5, 2.5, 3.5, 4.5)],
                dtype=[('rho', '<f8'), ('phi', '<f8'), ('eta', '<f8'), ('tau', '<f8')])

Vector properties

Vectors can have properties and provide quantities as attributes

vector.obj(x=1, y=2, z=3).theta
vector.obj(x=1, y=2, z=3).eta

1.1035868415601453

vector.obj(x=3, y=4, z=-2, t=10)   # Cartesian 4D vector

VectorObject4D(x=3, y=4, z=-2, t=10)