Lorentz vectors¶
In keeping with the “many small packages” philosophy, 2D/3D/Lorentz vectors are handled by a package named Vector. This is where you can find calculations like deltaR and coordinate transformations.
import vector
one = vector.obj(px=1, py=0, pz=0)
two = vector.obj(px=0, py=1, pz=1)
one + twoMomentumObject3D(px=1, py=1, pz=1)one.deltaR(two)np.float64(1.8011719796199461)one.to_rhophieta()MomentumObject3D(pt=1.0, phi=0.0, eta=0.0)two.to_rhophieta()MomentumObject3D(pt=1.0, phi=1.5707963267948966, eta=0.881373587019543)To fit in with the rest of the ecosystem, Vector must be an array-oriented library. Arrays of 2D/3D/Lorentz vectors are processed in bulk.
MomentumNumpy2D, MomentumNumpy3D, MomentumNumpy4D are NumPy array subtypes: NumPy arrays can be cast to these types and get all the vector functions.
import skhep_testdata, uproot
import awkward as ak
import vector
tree = uproot.open(skhep_testdata.data_path("uproot-Zmumu.root"))["events"]
one = ak.to_numpy(tree.arrays(filter_name=["E1", "p[xyz]1"]))
two = ak.to_numpy(tree.arrays(filter_name=["E2", "p[xyz]2"]))
one.dtype.names = ("E", "px", "py", "pz")
two.dtype.names = ("E", "px", "py", "pz")
one = one.view(vector.MomentumNumpy4D)
two = two.view(vector.MomentumNumpy4D)
one + twoMomentumNumpy4D([( -7.0508504 , 1.31371932, -116.3919462 , 142.82374098),
( -6.07723788, 0.86288157, -117.74020834, 144.54679534),
( -5.76527367, 0.72893471, -117.22250173, 143.92770728), ...,
(-35.664423 , -24.9064416 , -226.76744871, 250.05025691),
(-36.41664408, -25.19899466, -228.38003444, 251.853268 ),
(-36.30874217, -25.19705013, -228.65597631, 252.14934978)],
shape=(2304,), dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8'), ('t', '<f8')])one.deltaR(two)array([3.10530402, 3.10550819, 3.10547199, ..., 2.81363587, 2.81359933,
2.81354694], shape=(2304,))one.to_rhophieta()MomentumNumpy3D([(44.7322, 2.74126 , -1.21769), (38.8311, -0.440873 , -1.05139),
(38.8311, -0.440873 , -1.05139), ...,
(32.3997, 0.0370275, -1.57044), (32.3997, 0.0370275, -1.57044),
(32.5076, 0.0369644, -1.57078)],
shape=(2304,), dtype=[('rho', '<f8'), ('phi', '<f8'), ('eta', '<f8')])two.to_rhophieta()MomentumNumpy3D([(37.7582, -0.441078, -1.05138), (44.7322, 2.74126 , -1.21769),
(44.3927, 2.7413 , -1.21776), ...,
(72.8781, -2.77524 , -1.4827 ), (73.6852, -2.77519 , -1.48227),
(73.6852, -2.77519 , -1.48227)],
shape=(2304,), dtype=[('rho', '<f8'), ('phi', '<f8'), ('eta', '<f8')])After vector.register_awkward() is called, "Momentum2D", "Momentum3D", "Momentum4D" are record names that Awkward Array will recognize to get all the vector functions.
vector.register_awkward()
tree = uproot.open(skhep_testdata.data_path("uproot-HZZ.root"))["events"]
array = tree.arrays(filter_name=["Muon_E", "Muon_P[xyz]"])
muons = ak.zip(
{"px": array.Muon_Px, "py": array.Muon_Py, "pz": array.Muon_Pz, "E": array.Muon_E},
with_name="Momentum4D",
)
mu1, mu2 = ak.unzip(ak.combinations(muons, 2))
mu1 + mu2Loading...
mu1.deltaR(mu2)Loading...
muons.to_rhophieta()Loading...
Particle properties and PDG identifiers¶
The Particle library provides all of the particle masses, decay widths and more from the PDG. It further contains a series of tools to programmatically query particle properties and use several identification schemes.
import particle
from hepunits import GeV
particle.Particle.findall("pi")[<Particle: name="pi0", pdgid=111, mass=134.9768 ± 0.0005 MeV>,
<Particle: name="pi+", pdgid=211, mass=139.57039 ± 0.00018 MeV>,
<Particle: name="pi-", pdgid=-211, mass=139.57039 ± 0.00018 MeV>,
<Particle: name="pi(2)(1670)0", pdgid=10115, mass=1670.6 + 2.9 - 1.2 MeV>,
<Particle: name="pi(2)(1670)+", pdgid=10215, mass=1670.6 + 2.9 - 1.2 MeV>,
<Particle: name="pi(2)(1670)-", pdgid=-10215, mass=1670.6 + 2.9 - 1.2 MeV>,
<Particle: name="pi(1300)0", pdgid=100111, mass=1300 ± 100 MeV>,
<Particle: name="pi(1300)+", pdgid=100211, mass=1300 ± 100 MeV>,
<Particle: name="pi(1300)-", pdgid=-100211, mass=1300 ± 100 MeV>,
<Particle: name="pi(1)(1400)0", pdgid=9000113, mass=None>,
<Particle: name="pi(1)(1400)+", pdgid=9000213, mass=None>,
<Particle: name="pi(1)(1400)-", pdgid=-9000213, mass=None>,
<Particle: name="pi(1800)0", pdgid=9010111, mass=1810 + 9 - 11 MeV>,
<Particle: name="pi(1)(1600)0", pdgid=9010113, mass=1645 + 40 - 17 MeV>,
<Particle: name="pi(1800)+", pdgid=9010211, mass=1810 + 9 - 11 MeV>,
<Particle: name="pi(1800)-", pdgid=-9010211, mass=1810 + 9 - 11 MeV>,
<Particle: name="pi(1)(1600)+", pdgid=9010213, mass=1645 + 40 - 17 MeV>,
<Particle: name="pi(1)(1600)-", pdgid=-9010213, mass=1645 + 40 - 17 MeV>]z_boson = particle.Particle.from_name("Z0")
z_boson.mass / GeV, z_boson.width / GeV(91.188, 2.4955)print(z_boson.describe())Name: Z0 ID: 23 Latex: $Z^{0}$
Mass = 91188.0 ± 2.0 MeV
Width = 2495.5 ± 2.3 MeV
Q (charge) = 0 J (total angular) = 1.0 P (space parity) = None
C (charge parity) = None I (isospin) = None G (G-parity) = None
Antiparticle name: Z0 (antiparticle status: Same)
particle.Particle.from_pdgid(111)Loading...
particle.Particle.findall(
lambda p: p.pdgid.is_meson and p.pdgid.has_strange and p.pdgid.has_charm
)[<Particle: name="D(s)+", pdgid=431, mass=1968.35 ± 0.07 MeV>,
<Particle: name="D(s)-", pdgid=-431, mass=1968.35 ± 0.07 MeV>,
<Particle: name="D(s)*+", pdgid=433, mass=2112.2 ± 0.4 MeV>,
<Particle: name="D(s)*-", pdgid=-433, mass=2112.2 ± 0.4 MeV>,
<Particle: name="D(s2)*(2573)+", pdgid=435, mass=2569.1 ± 0.8 MeV>,
<Particle: name="D(s2)*(2573)-", pdgid=-435, mass=2569.1 ± 0.8 MeV>,
<Particle: name="D(s0)*(2317)+", pdgid=10431, mass=2317.8 ± 0.5 MeV>,
<Particle: name="D(s0)*(2317)-", pdgid=-10431, mass=2317.8 ± 0.5 MeV>,
<Particle: name="D(s1)(2536)+", pdgid=10433, mass=2535.11 ± 0.06 MeV>,
<Particle: name="D(s1)(2536)-", pdgid=-10433, mass=2535.11 ± 0.06 MeV>,
<Particle: name="D(s1)(2460)+", pdgid=20433, mass=2459.5 ± 0.6 MeV>,
<Particle: name="D(s1)(2460)-", pdgid=-20433, mass=2459.5 ± 0.6 MeV>]print(particle.PDGID(211).info())A None
J 0.0
L 0
S 0
Z None
abspid 211
charge 1.0
has_bottom False
has_charm False
has_down True
has_fundamental_anti False
has_strange False
has_top False
has_up True
is_Qball False
is_Rhadron False
is_SUSY False
is_baryon False
is_diquark False
is_dyon False
is_excited_quark_or_lepton False
is_gauge_boson_or_higgs False
is_generator_specific False
is_hadron True
is_lepton False
is_meson True
is_nucleus False
is_pentaquark False
is_quark False
is_sm_gauge_boson_or_higgs False
is_sm_lepton False
is_sm_quark False
is_special_particle False
is_technicolor False
is_valid True
j_spin 1
l_spin 1
s_spin 1
three_charge 3
pdgid = particle.Corsika7ID(11).to_pdgid()
particle.Particle.from_pdgid(pdgid)Loading...
Jet clustering¶
In a high-energy pp collision, for instance, a spray of hadrons is produced which is clustered into `jets’ of particles and this method/process is called jet-clustering. The anti-kt jet clustering algorithm is one such algorithm used to combine particles/hadrons that are close to each other into jets.
Some people need to do jet-clustering at the analysis level. The fastjet package makes it possible to do that an (Awkward) array at a time.
import skhep_testdata, uproot
import awkward as ak
import particle
from hepunits import GeV
import vector
vector.register_awkward()
picodst = uproot.open(
"https://pivarski-princeton.s3.amazonaws.com/pythia_ppZee_run17emb.picoDst.root:PicoDst"
)
px, py, pz = ak.unzip(
picodst.arrays(filter_name=["Track/Track.mPMomentum[XYZ]"], entry_stop=100)
)
probable_mass = particle.Particle.from_name("pi+").mass / GeV
pseudojets = ak.zip(
{"px": px, "py": py, "pz": pz, "mass": probable_mass}, with_name="Momentum4D"
)
good_pseudojets = pseudojets[pseudojets.pt > 0.1]
import fastjet
jetdef = fastjet.JetDefinition(fastjet.antikt_algorithm, 1.0)
clusterseq = fastjet.ClusterSequence(good_pseudojets, jetdef)
clusterseq.inclusive_jets()
ak.num(good_pseudojets), ak.num(clusterseq.inclusive_jets())#--------------------------------------------------------------------------
# FastJet release 3.5.1
# M. Cacciari, G.P. Salam and G. Soyez
# A software package for jet finding and analysis at colliders
# https://fastjet.fr
#
# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package
# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210].
#
# FastJet is provided without warranty under the GNU GPL v2 or higher.
# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code,
# CGAL and 3rd party plugin jet algorithms. See COPYING file for details.
#--------------------------------------------------------------------------
(<Array [24, 19, 13, 36, 39, 48, ..., 57, 41, 35, 34, 20, 93] type='100 * int64'>,
<Array [8, 6, 6, 7, 8, 9, 8, 7, ..., 9, 8, 7, 6, 7, 4, 11] type='100 * int64'>)This fastjet package uses Vector to get coordinate transformations and all the Lorentz vector methods you’re accustomed to having in pseudo-jet objects. I used Particle to impute the mass of particles with only track-level information.
See how all the pieces accumulate?