Vectores de Lorentz, IDs de partículas del PDG, agrupamiento de jets, ¡ay!

Vectores de Lorentz¶

Siguiendo con la filosofía de “muchos paquetes pequeños”, los vectores 2D/3D/de Lorentz los maneja un paquete llamado Vector. Aquí es donde puedes encontrar cálculos como deltaR y transformaciones de coordenadas.

import vector

uno = vector.obj(px=1, py=0, pz=0)
dos = vector.obj(px=0, py=1, pz=1)

uno + dos

MomentumObject3D(px=1, py=1, pz=1)

uno.deltaR(dos)

np.float64(1.8011719796199461)

uno.to_rhophieta()

MomentumObject3D(pt=1.0, phi=0.0, eta=0.0)

dos.to_rhophieta()

MomentumObject3D(pt=1.0, phi=1.5707963267948966, eta=0.881373587019543)

Para integrarse con el resto del ecosistema, Vector debe ser una librería orientada a arrays. Los arrays de vectores 2D/3D/de Lorentz se procesan en bloque.

MomentumNumpy2D, MomentumNumpy3D y MomentumNumpy4D son subtipos de arrays de NumPy: los arrays de NumPy se pueden convertir a estos tipos y obtener todas las funciones vectoriales.

import skhep_testdata, uproot
import awkward as ak
import vector

tree = uproot.open(skhep_testdata.data_path("uproot-Zmumu.root"))["events"]

uno = ak.to_numpy(tree.arrays(filter_name=["E1", "p[xyz]1"]))
dos = ak.to_numpy(tree.arrays(filter_name=["E2", "p[xyz]2"]))

uno.dtype.names = ("E", "px", "py", "pz")
dos.dtype.names = ("E", "px", "py", "pz")

uno = uno.view(vector.MomentumNumpy4D)
dos = dos.view(vector.MomentumNumpy4D)

uno + dos

MomentumNumpy4D([( -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')])

uno.deltaR(dos)

array([3.10530402, 3.10550819, 3.10547199, ..., 2.81363587, 2.81359933,
       2.81354694], shape=(2304,))

uno.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')])

dos.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')])

Después de llamar a vector.register_awkward(), "Momentum2D", "Momentum3D" y "Momentum4D" son nombres de registros que Awkward Array reconocerá para obtener todas las funciones vectoriales.

vector.register_awkward()

tree = uproot.open(skhep_testdata.data_path("uproot-HZZ.root"))["events"]

array = tree.arrays(filter_name=["Muon_E", "Muon_P[xyz]"])

muones = 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(muones, 2))

mu1 + mu2

mu1.deltaR(mu2)

muones.to_rhophieta()

Propiedades de las partículas e identificadores del PDG¶

La librería Particle proporciona todas las masas de las partículas, anchos de desintegración y más, tomados del PDG. Además, contiene una serie de herramientas para consultar programáticamente las propiedades de las partículas y usar varios esquemas de identificación.

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>]

boson_z = particle.Particle.from_name("Z0")
boson_z.mass / GeV, boson_z.width / GeV

(91.1879, 2.4955)

print(boson_z.describe())

Name: Z0             ID: 23           Latex: $Z^{0}$
Mass  = 91187.9 ± 2.0 MeV
Width = 2495.5 ± 2.3 MeV
Q (charge)        = 0       J (total angular) = 1        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)

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.12 ± 0.06 MeV>,
 <Particle: name="D(s1)(2536)-", pdgid=-10433, mass=2535.12 ± 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
L              0
S              0
Z              None
abspid         211
charge         1
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)

Agrupamiento de jets¶

En una colisión pp de alta energía, por ejemplo, se produce una lluvia de hadrones que se agrupan en “jets” de partículas, y a este método/proceso se le llama agrupamiento de jets (jet clustering). El algoritmo de agrupamiento anti-kt es uno de esos algoritmos usados para combinar en jets las partículas/hadrones que están cerca entre sí.

Algunas personas necesitan hacer agrupamiento de jets a nivel de análisis. El paquete fastjet permite hacerlo un array (Awkward) a la vez.

import skhep_testdata, uproot
import awkward as ak
import particle
from hepunits import GeV
import vector

vector.register_awkward()

picodst = uproot.open(
    "https://zenodo.org/records/21777191/files/pythia_ppZee_run17emb.picoDst.root:PicoDst"
)
px, py, pz = ak.unzip(
    picodst.arrays(filter_name=["Track/Track.mPMomentum[XYZ]"], entry_stop=100)
)

masa_probable = particle.Particle.from_name("pi+").mass / GeV

pseudojets = ak.zip(
    {"px": px, "py": py, "pz": pz, "mass": masa_probable}, with_name="Momentum4D"
)
buenos_pseudojets = pseudojets[pseudojets.pt > 0.1]

import fastjet

jetdef = fastjet.JetDefinition(fastjet.antikt_algorithm, 1.0)

clusterseq = fastjet.ClusterSequence(buenos_pseudojets, jetdef)
clusterseq.inclusive_jets()

ak.num(buenos_pseudojets), ak.num(clusterseq.inclusive_jets())

(<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'>)

#--------------------------------------------------------------------------
#                         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.
#--------------------------------------------------------------------------

Este paquete fastjet usa Vector para obtener las transformaciones de coordenadas y todos los métodos de vectores de Lorentz a los que estás acostumbrado/a en los objetos pseudo-jet. Usamos Particle para imputar la masa de partículas de las que solo tenemos información a nivel de traza.

¿Ves cómo se van acumulando todas las piezas?