Basic IO, event inspection, and visualization
We load an event from a file using pyhepmc.open and demonstrate how it renders in Jupyter (this requires graphviz).
[1]:
import pyhepmc
# pyhepmc.open can read most HepMC formats using auto-detection
with pyhepmc.open("../../tests/pythia6_simple.dat") as f:
event = f.read()
event
[1]:
The rendered SVG has tooltips that provide extra information on each particle and vertex.
Not all events have beam particles. The Sibyll-2.1 stack does not contain them.
[2]:
with pyhepmc.open("../../tests/sibyll21.dat") as f:
event2 = f.read()
event2
[2]:
pyhepmc.open also supports iteration like the built-in open. This is the default way to read many events from a file. This file only contains one event, however.
[3]:
with pyhepmc.open("../../tests/sibyll21.dat") as f:
for event in f:
print(event)
________________________________________________________________________
GenEvent: #0
Momentum units: GEV Position units: MM
Entries in this event: 7 vertices, 23 particles, 0 weights.
Position offset: 0, 0, 0, 0
GenParticle Legend
ID PDG ID ( px, py, pz, E ) Stat ProdVtx
________________________________________________________________________
Vtx: -1 stat: 0 (X,cT): 0
I: 1 2224 -8.92e-02,+1.35e-01,+2.03e+00,+2.38e+00 2 0
O: 8 2212 -1.10e-01,+9.29e-02,+1.17e+00,+1.51e+00 1 -1
9 211 +2.11e-02,+4.21e-02,+8.64e-01,+8.76e-01 1 -1
Vtx: -2 stat: 0 (X,cT): 0
I: 2 111 +8.95e-02,-4.49e-01,+3.42e+00,+3.45e+00 2 0
O: 10 22 -1.59e-02,-2.60e-01,+1.94e+00,+1.96e+00 1 -2
11 22 +1.05e-01,-1.89e-01,+1.48e+00,+1.49e+00 1 -2
Vtx: -3 stat: 0 (X,cT): 0
I: 3 331 +1.50e-01,+4.60e-01,+2.04e+00,+2.31e+00 2 0
O: 12 211 +1.30e-02,+9.59e-02,+2.26e-01,+2.83e-01 1 -3
13 -211 +9.56e-02,+6.02e-02,+1.22e-01,+2.18e-01 1 -3
14 221 +4.13e-02,+3.04e-01,+1.69e+00,+1.81e+00 2 -3
Vtx: -4 stat: 0 (X,cT): 0
I: 4 221 -1.62e-01,-2.44e-01,+1.95e+00,+2.05e+00 2 0
O: 15 22 -1.35e-01,+4.44e-02,+1.53e+00,+1.54e+00 1 -4
16 22 -2.70e-02,-2.89e-01,+4.21e-01,+5.12e-01 1 -4
Vtx: -5 stat: 0 (X,cT): 0
I: 6 111 -2.49e-02,+3.66e-01,-1.89e+00,+1.93e+00 2 0
O: 17 22 -4.50e-02,+2.62e-01,-1.59e+00,+1.61e+00 1 -5
18 22 +2.01e-02,+1.04e-01,-3.07e-01,+3.25e-01 1 -5
Vtx: -6 stat: 0 (X,cT): 0
I: 14 221 +4.13e-02,+3.04e-01,+1.69e+00,+1.81e+00 2 -3
O: 19 211 +1.62e-02,+1.12e-01,+3.07e-01,+3.57e-01 1 -6
20 -211 -5.02e-02,+7.78e-02,+2.68e-01,+3.17e-01 1 -6
21 111 +7.53e-02,+1.14e-01,+1.11e+00,+1.13e+00 2 -6
Vtx: -7 stat: 0 (X,cT): 0
I: 21 111 +7.53e-02,+1.14e-01,+1.11e+00,+1.13e+00 2 -6
O: 22 22 +6.89e-02,-3.50e-03,+5.46e-01,+5.51e-01 1 -7
23 22 +6.39e-03,+1.18e-01,+5.69e-01,+5.82e-01 1 -7
________________________________________________________________________
Coming back to visualization: you can customize it with pyhepmc.view.to_dot.
[4]:
from pyhepmc.view import to_dot
# see help(pyhepmc.view.to_dot) for details
g = to_dot(event2, size=(4, 4), color_hadron="red")
g
[4]:
You can also save these images like this. Many common formats are supported, see docs.
[5]:
from pyhepmc.view import savefig
savefig(event2, "event.svg")
savefig(event2, "event.png")
savefig(event2, "event.pdf")
Let’s look at a few particles and vertices.
[6]:
# first five particles in the record
event.particles[:5]
[6]:
[GenParticle(FourVector(-0.0892, 0.135, 2.03, 2.38), mass=1.231, pid=2224, status=2),
GenParticle(FourVector(0.0895, -0.449, 3.42, 3.45), mass=0.13497, pid=111, status=2),
GenParticle(FourVector(0.15, 0.46, 2.04, 2.31), mass=0.9575, pid=331, status=2),
GenParticle(FourVector(-0.162, -0.244, 1.95, 2.05), mass=0.5488, pid=221, status=2),
GenParticle(FourVector(-0.0529, 0.253, -0.0519, 0.299), mass=0.13957, pid=-211, status=1)]
[7]:
# final state particles sorted by energy
[p for p in sorted(event.particles, key=lambda p:p.momentum.e) if p.status == 1]
[7]:
[GenParticle(FourVector(0.0956, 0.0602, 0.122, 0.218), mass=0.13957, pid=-211, status=1),
GenParticle(FourVector(0.013, 0.0959, 0.226, 0.283), mass=0.13957, pid=211, status=1),
GenParticle(FourVector(-0.0529, 0.253, -0.0519, 0.299), mass=0.13957, pid=-211, status=1),
GenParticle(FourVector(-0.0502, 0.0778, 0.268, 0.317), mass=0.13957, pid=-211, status=1),
GenParticle(FourVector(0.0201, 0.104, -0.307, 0.325), mass=0, pid=22, status=1),
GenParticle(FourVector(0.0162, 0.112, 0.307, 0.357), mass=0.13957, pid=211, status=1),
GenParticle(FourVector(-0.027, -0.289, 0.421, 0.512), mass=0, pid=22, status=1),
GenParticle(FourVector(0.0689, -0.0035, 0.546, 0.551), mass=0, pid=22, status=1),
GenParticle(FourVector(0.00639, 0.118, 0.569, 0.582), mass=0, pid=22, status=1),
GenParticle(FourVector(0.0211, 0.0421, 0.864, 0.876), mass=0.13957, pid=211, status=1),
GenParticle(FourVector(0.105, -0.189, 1.48, 1.49), mass=0, pid=22, status=1),
GenParticle(FourVector(-0.11, 0.0929, 1.17, 1.51), mass=0.93827, pid=2212, status=1),
GenParticle(FourVector(-0.135, 0.0444, 1.53, 1.54), mass=0, pid=22, status=1),
GenParticle(FourVector(-0.045, 0.262, -1.59, 1.61), mass=0, pid=22, status=1),
GenParticle(FourVector(-0.0159, -0.26, 1.94, 1.96), mass=0, pid=22, status=1),
GenParticle(FourVector(0.0898, -0.521, -7.5, 7.58), mass=0.93827, pid=2212, status=1)]
All particles which produce photons.
[8]:
def has_photon_children(p):
return any(q.pid == 22 for q in p.children)
[p for p in event.particles if has_photon_children(p)]
[8]:
[GenParticle(FourVector(0.0895, -0.449, 3.42, 3.45), mass=0.13497, pid=111, status=2),
GenParticle(FourVector(-0.162, -0.244, 1.95, 2.05), mass=0.5488, pid=221, status=2),
GenParticle(FourVector(-0.0249, 0.366, -1.89, 1.93), mass=0.13497, pid=111, status=2),
GenParticle(FourVector(0.0753, 0.114, 1.11, 1.13), mass=0.13497, pid=111, status=2)]
This event has no attached metadata. HepMC3 supports metadata attached to individual particles and vertices, events, and the whole run. We use this to add information about the generator, Pythia-6.4.28. We add this information to GenRunInfo and mark all vertices that have photon children as “sparkly”.
Then, we save the modified event.
[9]:
# run info is empty
event.run_info
[9]:
GenRunInfo(tools=[ToolInfo(name='SIBYLL', version='2.1', description='')], weight_names=[], attributes={})
[10]:
event.run_info.tools = [pyhepmc.GenRunInfo.ToolInfo("Pythia", "6.4.28", "Legacy Pythia 6 event")]
# this equivalent shortcut is also supported:
# event.run_info.tools = [("Pythia", "6.4.28", "Legacy Pythia 6 event")]
for v in event.vertices:
if any(q.pid == 22 for q in v.particles_out):
v.attributes["sparkly"] = True
[11]:
# save event with reduced precision
with pyhepmc.open("modified_event.dat", "w", precision=3) as f:
f.write(event)
The HepMC3 format is human-readable ASCII.
[12]:
with open("modified_event.dat") as f:
print(f.read().strip())
HepMC::Version 3.02.05
HepMC::Asciiv3-START_EVENT_LISTING
T Pythia\|6.4.28\|Legacy Pythia 6 event
E 0 7 23
U GEV MM
A -7 sparkly 1
A -5 sparkly 1
A -4 sparkly 1
A -2 sparkly 1
P 1 0 2224 -8.921e-02 1.349e-01 2.034e+00 2.383e+00 1.231e+00 2
P 2 0 111 8.946e-02 -4.486e-01 3.418e+00 3.452e+00 1.350e-01 2
P 3 0 331 1.499e-01 4.600e-01 2.036e+00 2.307e+00 9.575e-01 2
P 4 0 221 -1.621e-01 -2.442e-01 1.955e+00 2.051e+00 5.488e-01 2
P 5 0 -211 -5.295e-02 2.535e-01 -5.190e-02 2.987e-01 1.396e-01 1
P 6 0 111 -2.490e-02 3.658e-01 -1.892e+00 1.932e+00 1.350e-01 2
P 7 0 2212 8.976e-02 -5.213e-01 -7.499e+00 7.576e+00 9.383e-01 1
P 8 1 2212 -1.103e-01 9.285e-02 1.171e+00 1.507e+00 9.383e-01 1
P 9 1 211 2.107e-02 4.207e-02 8.636e-01 8.760e-01 1.396e-01 1
P 10 2 22 -1.591e-02 -2.596e-01 1.942e+00 1.960e+00 0.000e+00 1
P 11 2 22 1.054e-01 -1.892e-01 1.477e+00 1.493e+00 0.000e+00 1
P 12 3 211 1.301e-02 9.594e-02 2.256e-01 2.833e-01 1.396e-01 1
P 13 3 -211 9.560e-02 6.017e-02 1.217e-01 2.180e-01 1.396e-01 1
P 14 3 221 4.130e-02 3.039e-01 1.689e+00 1.806e+00 5.488e-01 2
P 15 4 22 -1.351e-01 4.437e-02 1.533e+00 1.540e+00 0.000e+00 1
P 16 4 22 -2.698e-02 -2.886e-01 4.215e-01 5.115e-01 0.000e+00 1
P 17 6 22 -4.505e-02 2.616e-01 -1.586e+00 1.608e+00 0.000e+00 1
P 18 6 22 2.014e-02 1.042e-01 -3.067e-01 3.245e-01 0.000e+00 1
P 19 14 211 1.624e-02 1.118e-01 3.071e-01 3.567e-01 1.396e-01 1
P 20 14 -211 -5.020e-02 7.775e-02 2.676e-01 3.167e-01 1.396e-01 1
P 21 14 111 7.526e-02 1.144e-01 1.114e+00 1.132e+00 1.350e-01 2
P 22 21 22 6.889e-02 -3.499e-03 5.456e-01 5.509e-01 0.000e+00 1
P 23 21 22 6.393e-03 1.179e-01 5.686e-01 5.818e-01 0.000e+00 1
HepMC::Asciiv3-END_EVENT_LISTING
HepMC2 support
The HepMC3 format is used by default, but some programs still expects HepMC2 format. pyhepmc.open detects the format automatically when reading data. If you want to write in HepMC2, you can specify that with the format argument.
[18]:
# save event with reduced precision
with pyhepmc.open("modified_event.hepmc2", "w", format="hepmc2", precision=3) as f:
f.write(event)
with open("modified_event.hepmc2") as f:
print(f.read())
WARNING::WriterAsciiHepMC2::WriterAsciiHepMC2: HepMC2 IO_GenEvent format is outdated. Please use HepMC3 Asciiv3 format instead.
HepMC::Version 3.02.05
HepMC::IO_GenEvent-START_EVENT_LISTING
E 0 0 0.000000e+00 0.000000e+00 0.000000e+00 0 0 7 10000 10000 0 0
U GEV MM
V -1 0 0 0 0 0 1 2 0
P 10001 2224 -8.921e-02 1.349e-01 2.034e+00 2.383e+00 1.231e+00 2 0 0 -1 0
P 10002 2212 -1.103e-01 9.285e-02 1.171e+00 1.507e+00 9.383e-01 1 0 0 0 0
P 10003 211 2.107e-02 4.207e-02 8.636e-01 8.760e-01 1.396e-01 1 0 0 0 0
V -2 0 0 0 0 0 1 2 0
P 10004 111 8.946e-02 -4.486e-01 3.418e+00 3.452e+00 1.350e-01 2 0 0 -2 0
P 10005 22 -1.591e-02 -2.596e-01 1.942e+00 1.960e+00 0.000e+00 1 0 0 0 0
P 10006 22 1.054e-01 -1.892e-01 1.477e+00 1.493e+00 0.000e+00 1 0 0 0 0
V -3 0 0 0 0 0 1 3 0
P 10007 331 1.499e-01 4.600e-01 2.036e+00 2.307e+00 9.575e-01 2 0 0 -3 0
P 10008 211 1.301e-02 9.594e-02 2.256e-01 2.833e-01 1.396e-01 1 0 0 0 0
P 10009 -211 9.560e-02 6.017e-02 1.217e-01 2.180e-01 1.396e-01 1 0 0 0 0
P 10010 221 4.130e-02 3.039e-01 1.689e+00 1.806e+00 5.488e-01 2 0 0 -6 0
V -4 0 0 0 0 0 1 2 0
P 10011 221 -1.621e-01 -2.442e-01 1.955e+00 2.051e+00 5.488e-01 2 0 0 -4 0
P 10012 22 -1.351e-01 4.437e-02 1.533e+00 1.540e+00 0.000e+00 1 0 0 0 0
P 10013 22 -2.698e-02 -2.886e-01 4.215e-01 5.115e-01 0.000e+00 1 0 0 0 0
V -5 0 0 0 0 0 1 2 0
P 10014 111 -2.490e-02 3.658e-01 -1.892e+00 1.932e+00 1.350e-01 2 0 0 -5 0
P 10015 22 -4.505e-02 2.616e-01 -1.586e+00 1.608e+00 0.000e+00 1 0 0 0 0
P 10016 22 2.014e-02 1.042e-01 -3.067e-01 3.245e-01 0.000e+00 1 0 0 0 0
V -6 0 0 0 0 0 0 3 0
P 10017 211 1.624e-02 1.118e-01 3.071e-01 3.567e-01 1.396e-01 1 0 0 0 0
P 10018 -211 -5.020e-02 7.775e-02 2.676e-01 3.167e-01 1.396e-01 1 0 0 0 0
P 10019 111 7.526e-02 1.144e-01 1.114e+00 1.132e+00 1.350e-01 2 0 0 -7 0
V -7 0 0 0 0 0 0 2 0
P 10020 22 6.889e-02 -3.499e-03 5.456e-01 5.509e-01 0.000e+00 1 0 0 0 0
P 10021 22 6.393e-03 1.179e-01 5.686e-01 5.818e-01 0.000e+00 1 0 0 0 0
HepMC::IO_GenEvent-END_EVENT_LISTING
Using attributes
Gotcha: If you deserialize attributes, you need to cast them to the right types yourself, because the HepMC3 C++ library does not store the type information. This is a limitation of the current HepMC3 format, not a limitation in pyhepmc. We try to work around this issue as good as possible, but a fix would require a change to HepMC3 that was rejected by the HepMC3 developers.
[13]:
with pyhepmc.open("modified_event.dat") as f:
event = f.read()
for v in event.vertices:
if v.attributes:
print(v.attributes)
<AttributesView>{'sparkly': <UnparsedAttribute>('1')}
<AttributesView>{'sparkly': <UnparsedAttribute>('1')}
<AttributesView>{'sparkly': <UnparsedAttribute>('1')}
<AttributesView>{'sparkly': <UnparsedAttribute>('1')}
As you can see, the deserialized attribute is restored as a UnparsedAttribute. To turn it into a boolean, call its .astype method. The method name should remind you of numpy.
[14]:
# you need to remember what name corresponds to what type;
# sparkly was a boolean
key = "sparkly"
for v in event.vertices:
if key in v.attributes:
val = v.attributes[key]
b = val.astype(bool)
print(f"key: {key}, value: {b}")
key: sparkly, value: True
key: sparkly, value: True
key: sparkly, value: True
key: sparkly, value: True
You only have to do this once (although it does not hurt to do it more than once), since .astype also replaces the UnparsedAttribute with the typed attribute as a side effect. Next time you iterate over the attributes, they appear normally.
[15]:
for v in event.vertices:
if v.attributes:
print(v.attributes)
<AttributesView>{'sparkly': True}
<AttributesView>{'sparkly': True}
<AttributesView>{'sparkly': True}
<AttributesView>{'sparkly': True}