Reading LHE Files

A Les Houches Event file is the hand-off between a matrix-element generator and a parton shower: a plain-text record of hard-scattering events, one block each, with just enough information for the next program to continue the story. This page explains what every number in the file means, why the format looks the way it does, and how to write a reader that will not fall over on a real 10 GB MadGraph output.

01 · Three documents, one accord

The format grew in three steps, and the field names you will meet are inherited from the first of them.

• The Les Houches Accord (2001), hep-ph/0109068, was written at the Les Houches workshop to let matrix-element codes talk to Pythia and Herwig. It defines two Fortran common blocks: HEPRUP (the run: beams, PDFs, process list, cross sections) and HEPEUP (one event: particles, momenta, colours). Every variable name below — IDBMUP, XWGTUP, ICOLUP — comes from here; the trailing UP stands for “user process”.
• The Les Houches Event File, LHEF 1.0 (2006), hep-ph/0609017, put those common blocks into a file so the two programs no longer had to be linked together. The layout is XML-like: <init> holds HEPRUP, each <event> holds one HEPEUP, and the numbers inside are whitespace-separated in the order the common blocks list them.
• LHEF 3.0 (2014), arXiv:1405.1067, added what NLO and multi-leg generators needed: named alternative weights (<initrwgt>, <rwgt>), per-event scale information, attributes on the <event> tag, and a versioned <generator> tag. Everything is backwards-compatible: a 1.0 reader that ignores what it does not understand still gets the events right.

The current specification of the numbering scheme for particle identities is the PDG Monte Carlo Particle Numbering Scheme, summarized on the Newcomers page.

02 · Anatomy of a file

Here is a complete, valid, two-event file. It is small enough to read in full, and every construction you will meet in a real file appears at least once: a header with an LHEF 3 weight definition, an init block listing two processes, one Drell–Yan event with an intermediate resonance, and one gluon-fusion top-pair event with a negative weight and a non-trivial colour flow. Tags are shown in red, weight lines in green.

<LesHouchesEvents version="3.0">
<header>
<!-- a hand-written sample for lhe_reader.py; two events, two processes -->
<MGVersion>
3.5.x
</MGVersion>
<initrwgt>
<weightgroup name="scale_variation" combine="envelope">
<weight id="1"> muR=1.0 muF=1.0 </weight>
<weight id="2"> muR=2.0 muF=2.0 </weight>
<weight id="3"> muR=0.5 muF=0.5 </weight>
</weightgroup>
</initrwgt>
</header>
<init>
  2212  2212  6.500000e+03  6.500000e+03  0  0  303600  303600  -4  2
  1.9376400e+03  2.100000e+00  1.9376400e+01  1
  8.3200000e+02  1.500000e+00  1.9376400e+01  2
</init>
<event npNLO="0">
  5   1  1.9376400000e+01  9.1651513900e+01  7.5467e-03  1.2440e-01
        2  -1    0    0  501    0 +0.0000000000e+00 +0.0000000000e+00 +6.0000000000e+01 +6.0000000000e+01 0.0000000000e+00 0.0000e+00 9.0
       -2  -1    0    0    0  501 +0.0000000000e+00 +0.0000000000e+00 -3.5000000000e+01 +3.5000000000e+01 0.0000000000e+00 0.0000e+00 9.0
       23   2    1    2    0    0 +0.0000000000e+00 +0.0000000000e+00 +2.5000000000e+01 +9.5000000000e+01 9.1651513900e+01 0.0000e+00 9.0
       13   1    3    3    0    0 +3.1236264540e+01 +2.6309942670e+01 +3.4045758500e+01 +5.3169936450e+01 1.0566000000e-01 0.0000e+00 -1.0
      -13   1    3    3    0    0 -3.1236264540e+01 -2.6309942670e+01 -9.0457584970e+00 +4.1830063550e+01 1.0566000000e-01 0.0000e+00 1.0
<rwgt>
<wgt id='1'> +1.9376400e+01 </wgt>
<wgt id='2'> +1.8801000e+01 </wgt>
<wgt id='3'> +2.0155000e+01 </wgt>
</rwgt>
</event>
<event npNLO="1">
  4   1 -1.9376400000e+01  3.6200000000e+02  7.5467e-03  1.0620e-01
       21  -1    0    0  501  502 +0.0000000000e+00 +0.0000000000e+00 +4.0000000000e+02 +4.0000000000e+02 0.0000000000e+00 0.0000e+00 9.0
       21  -1    0    0  502  503 +0.0000000000e+00 +0.0000000000e+00 -2.5000000000e+02 +2.5000000000e+02 0.0000000000e+00 0.0000e+00 9.0
        6   1    1    2  501    0 +6.4386824175e+01 -2.3192793698e+02 -3.8213734472e+01 +2.9887375358e+02 1.7300000000e+02 0.0000e+00 1.0
       -6   1    1    2    0  503 -6.4386824175e+01 +2.3192793698e+02 +1.8821373447e+02 +3.5112624642e+02 1.7300000000e+02 0.0000e+00 -1.0
<mgrwt>
<rscale>  0 0.36200000E+03</rscale>
</mgrwt>
<rwgt>
<wgt id='1'> -1.9376400e+01 </wgt>
<wgt id='2'> -1.7100000e+01 </wgt>
<wgt id='3'> -2.2300000e+01 </wgt>
</rwgt>
</event>
</LesHouchesEvents>

Momenta in GeV, cross sections in pb, lifetimes in mm. Numbers may be in C or Fortran notation (1.0E+00 or 1.0D+00).

Read from the top: the file is one <LesHouchesEvents> element. Inside it comes an optional <header> (free text: generator banners, run cards, weight definitions), then exactly one <init>, then any number of <event> blocks, each of which is a header line followed by one line per particle, followed by optional extra tagged blocks. That is the whole grammar.

03 · The init block

The first line of <init> is HEPRUP's run-level scalars, in this order:

Then NPRUP lines, one per subprocess:

Note what is not here: the number of events. An LHE file never says how many events it contains. You count them.

04 · The event block

The first line inside <event> is the event-level part of HEPEUP:

Then NUP particle lines, each with thirteen numbers:

In the Drell–Yan event, line 3 is the Z: its momentum is the sum of the incoming quark momenta (0, 0, 25, 95), its stored mass is 91.65 GeV, and its two daughters both point back to it through MOTHUP = 3 3. The muon lines carry SPINUP = −1 and +1: definite helicities, as a parity-violating Z decay should.

05 · Status codes and mothers

Mother indices point to earlier lines. MOTHUP(1) = MOTHUP(2) = 0 marks the incoming partons; outgoing particles produced directly in the hard scattering usually give 1 2 (both incomers); decay products give the resonance line. When MOTHUP(1) < MOTHUP(2) the pair is a range of mothers, not two individual ones.

06 · Colour flow

Showers do not evolve SU(3) amplitudes; they evolve colour lines in the large-N_c limit. So every coloured particle carries two integer labels: ICOLUP(1), the colour line it carries, and ICOLUP(2), the anticolour line. Quarks have only a colour, antiquarks only an anticolour, gluons one of each, colourless particles 0 0. Labels are arbitrary integers (convention: start at 501) and mean nothing across events; what matters is that each label appears exactly twice so that the lines close.

In the top-pair event of the sample:

gluon 1   colour 501  anticolour 502
gluon 2   colour 502  anticolour 503
top       colour 501
antitop              anticolour 503

Line 501 runs from gluon 1 to the top; line 502 connects the two gluons; line 503 runs from gluon 2 to the antitop. A colour line on an incoming particle is the same object as an anticolour line on an outgoing one — crossing swaps the roles — which is why a validity check has to flip colour and anticolour for status −1 entries before demanding that every label balances. Which of the several possible colour flows a given event carries was chosen probabilistically by the generator, weighted by the colour-ordered amplitudes; that choice is what determines where the shower radiates and how the hadronization strings are drawn.

07 · Weights and cross sections

IDWTUP in the init block fixes what XWGTUP means. Its sign says whether negative weights can occur, and its magnitude says how the shower should treat the events:

Two of these matter in practice. IDWTUP = 3 is what a leading-order MadGraph run produces: every event is worth the same, and the cross section is read from init. IDWTUP = −4 is the NLO convention (MadGraph5_aMC@NLO, POWHEG, Sherpa): weights are in pb, sum to Nσ, and a fraction of them are negative because the subtraction terms are distributed among the events. For a weighted sample the physics number is

σ = (1/N) ∑_i w_i, δσ = √( Var(w) / N ),

and a histogram bin contains the sum of the weights of events that landed in it, with a statistical error given by the square root of the sum of squared weights. The effective sample size is

N_eff = (∑w_i)² / ∑w_i² ≈ N (1 − 2f)²

for a fraction f of negative events of equal magnitude. At f = 0.25 you have paid for four events and got one; at f = 0.4 you have paid for twenty-five. That arithmetic is why negative-weight reduction is a research topic rather than a bookkeeping detail.

08 · LHEF 3: weights, scales, attributes

Everything in this section can be ignored by a reader that only wants kinematics, but it is where the systematics live.

Named weights

The header declares weights inside <initrwgt>, grouped in <weightgroup> elements whose name and combine attributes say what the group is (scale variation, PDF replicas, α_s variations) and how to combine it (envelope, hessian, replica…). Each <weight id="…"> carries a free-text description — muR=2.0 muF=1.0, or an LHAPDF member number. Every event then repeats those ids in a <rwgt> block with a number each. The ids are strings, not integers: MadGraph writes 1, 2, 3, POWHEG writes 1001, 1002, Sherpa writes names like MUR2_MUF1_PDF303600. Match on the string.

Event attributes

The <event> tag may carry attributes. The most common are npLO and npNLO: the number of extra partons this event contributes at LO or NLO, which multi-jet merging (FxFx, MEPS@NLO, UNLOPS) needs to decide how to treat the event in the shower. Generators are free to add their own.

Scales

A <scales> tag can carry per-parton shower starting scales (pt_start_1="…") plus muf, mur, mups; when present it overrides SCALUP. This is how POWHEG communicates the p_T veto scale to the shower.

Generator-specific blocks

MadGraph writes <mgrwt> (with the information needed to reweight to a new model) and, in the header, <MGGenerationInfo>, <MGRunCard>, <slha>. POWHEG writes #pdf, #scale and #rwgt comment lines inside the event. Pythia's LHEF3 writer adds <generator name="…" version="…"> to the header. Keep the raw text of anything you do not parse; you may need it downstream.

09 · The parsing logic

The tempting approach is to feed the file to an XML parser. Do not. Real LHE files are not well-formed XML: the header contains unescaped <, & and > from run cards and banners; POWHEG puts # comment lines inside events; some generators emit attributes with single quotes; and a strict parser has to read the entire document before it returns the first event, which for a 10 GB file means it never returns. Every production reader — Pythia's LHAupLHEF, Rivet, HepMC3's ReaderLHEF, pylhe — is a line-oriented state machine. This is the whole algorithm:

1. Open in text mode, transparently ungzipping if the name ends in .gz. Never read the whole file into memory.
2. Walk lines until <init>, collecting everything between <header> and </header> as raw text and, within it, the lines between <initrwgt> and </initrwgt> for the weight names.
3. Parse <init>: strip # comments, drop blank lines, split the first line into ten fields, then read NPRUP process lines. Keep any further lines as opaque extras.
4. Now yield events lazily. For each line: if you are outside an event and see <event, record its attributes and start buffering; if inside and you see </event>, parse the buffer and yield the result.
5. To parse an event buffer, first peel off tagged sub-blocks: lines between <rwgt> and </rwgt> become the weight dictionary; any other line beginning with < is stored verbatim. What remains is numeric: the first line is the six-field event header, the next NUP lines are particles with thirteen fields each. Replace Fortran D exponents before converting to float.
6. Everything else — blank lines, </LesHouchesEvents>, stray text — is ignored.

Two design decisions make the reader robust. Reading NUP particle lines by count rather than “until the next tag” means extra material after the particles cannot be mistaken for a fourteen-column particle. And returning the events as a generator means memory use is constant: you process a terabyte the same way you process a kilobyte.

10 · A streaming reader in Python

The implementation below has no dependencies. read_lhe(file) returns the parsed init block and an iterator of Event objects; check_event is the validator of section 11; cross_section applies the IDWTUP rules of section 07. Run it as a script on any file to get a summary.

#!/usr/bin/env python3
"""lhe_reader.py — a dependency-free, streaming reader for Les Houches Event files.

Handles LHEF 1.0 (hep-ph/0609017) and the LHEF 3.0 additions (arXiv:1405.1067):
<init>, <event> blocks, per-event <rwgt> weights, <initrwgt> weight definitions,
event attributes, and gzip-compressed files. It deliberately does NOT use an XML
parser: real LHE files carry unescaped '<', '&' and raw generator banners in the
<header>, so a strict parser will choke. A line-oriented state machine is what
every production reader (Pythia's LHAup, Rivet, pylhe) actually does.

Usage:
    from lhe_reader import read_lhe, check_event
    with open("events.lhe") as f:          # or gzip.open(..., "rt")
        init, events = read_lhe(f)
        for ev in events:
            print(ev.weight, [p.pid for p in ev.particles if p.status == 1])
"""
from __future__ import annotations
import gzip, math, re, sys
from dataclasses import dataclass, field
from typing import Dict, Iterator, List, Optional, TextIO, Tuple


# ------------------------------------------------------------------ records
@dataclass
class Particle:
    pid: int            # IDUP     PDG Monte Carlo numbering scheme
    status: int         # ISTUP    -1 in, +1 out, 2 resonance, -2 spacelike, 3 doc, -9 beam
    mother1: int        # MOTHUP(1)  1-based index into the event's particle list (0 = none)
    mother2: int        # MOTHUP(2)
    color: int          # ICOLUP(1)  colour line label (0 = none)
    anticolor: int      # ICOLUP(2)  anticolour line label
    px: float; py: float; pz: float; e: float; m: float   # PUP(1..5), GeV
    lifetime: float     # VTIMUP   c*tau, mm
    spin: float         # SPINUP   helicity / cos(angle to mother); 9 = unknown

    @property
    def pt(self) -> float: return math.hypot(self.px, self.py)
    @property
    def p(self) -> float: return math.sqrt(self.px**2 + self.py**2 + self.pz**2)
    @property
    def eta(self) -> float:
        p = self.p
        return 0.5 * math.log((p + self.pz) / (p - self.pz)) if p > abs(self.pz) else math.copysign(float("inf"), self.pz)
    @property
    def phi(self) -> float: return math.atan2(self.py, self.px)
    @property
    def mass_from_p4(self) -> float:
        m2 = self.e**2 - self.px**2 - self.py**2 - self.pz**2
        return math.copysign(math.sqrt(abs(m2)), m2)


@dataclass
class Event:
    nup: int            # NUP      number of particle entries
    process: int        # IDPRUP   which LPRUP process this event belongs to
    weight: float       # XWGTUP   event weight (meaning fixed by IDWTUP)
    scale: float        # SCALUP   factorisation/shower starting scale, GeV
    alpha_qed: float    # AQEDUP
    alpha_qcd: float    # AQCDUP
    particles: List[Particle] = field(default_factory=list)
    weights: Dict[str, float] = field(default_factory=dict)   # LHEF3 <rwgt><wgt id=...>
    attributes: Dict[str, str] = field(default_factory=dict)  # attributes on the <event> tag
    extra: List[str] = field(default_factory=list)            # any other lines inside <event>

    def final_state(self) -> List[Particle]:
        return [p for p in self.particles if p.status == 1]

    def incoming(self) -> List[Particle]:
        return [p for p in self.particles if p.status == -1]

    def p4_sum(self, parts) -> Tuple[float, float, float, float]:
        return (sum(p.px for p in parts), sum(p.py for p in parts),
                sum(p.pz for p in parts), sum(p.e for p in parts))


@dataclass
class ProcessInfo:
    xsec: float         # XSECUP  pb
    xsec_err: float     # XERRUP  pb
    xmax: float         # XMAXUP  maximum weight
    lprup: int          # LPRUP   process label used by IDPRUP


@dataclass
class Init:
    beam_pids: Tuple[int, int]        # IDBMUP
    beam_energies: Tuple[float, float] # EBMUP, GeV
    pdf_groups: Tuple[int, int]        # PDFGUP
    pdf_sets: Tuple[int, int]          # PDFSUP  (LHAPDF ids)
    weight_strategy: int               # IDWTUP  ±1..±4 (see docs)
    processes: List[ProcessInfo]       # NPRUP entries
    weight_names: Dict[str, str] = field(default_factory=dict)   # LHEF3 <initrwgt> id -> description
    extra: List[str] = field(default_factory=list)
    header: str = ""                   # raw text of the <header> block


# ------------------------------------------------------------------ parsing
_TAG_ATTR = re.compile(r'(\w+)\s*=\s*"([^"]*)"')


def _strip_comment(line: str) -> str:
    """Everything after '#' is a comment in the numeric blocks."""
    i = line.find("#")
    return line if i < 0 else line[:i]


def _parse_init(lines: List[str], weight_names: Dict[str, str], header: str) -> Init:
    numeric = [l for l in (_strip_comment(x).strip() for x in lines) if l]
    f = numeric[0].split()
    if len(f) < 10:
        raise ValueError(f"malformed <init> first line: {numeric[0]!r}")
    nprup = int(f[9])
    procs = []
    for l in numeric[1:1 + nprup]:
        x = l.split()
        procs.append(ProcessInfo(float(x[0]), float(x[1]), float(x[2]), int(x[3])))
    return Init(beam_pids=(int(f[0]), int(f[1])),
                beam_energies=(float(f[2]), float(f[3])),
                pdf_groups=(int(f[4]), int(f[5])),
                pdf_sets=(int(f[6]), int(f[7])),
                weight_strategy=int(f[8]),
                processes=procs,
                weight_names=weight_names,
                extra=numeric[1 + nprup:],
                header=header)


def _parse_event(lines: List[str], attrs: Dict[str, str]) -> Event:
    """lines = everything between <event ...> and </event>."""
    # 1. peel off LHEF3 <rwgt> ... </rwgt> and remember any other tagged blocks verbatim
    weights: Dict[str, float] = {}
    body: List[str] = []
    extra: List[str] = []
    in_rwgt = False
    for raw in lines:
        s = raw.strip()
        if not s:
            continue
        if s.startswith("<rwgt"):
            in_rwgt = True; continue
        if s.startswith("</rwgt"):
            in_rwgt = False; continue
        if in_rwgt:
            m = re.match(r"<wgt\s+id\s*=\s*['\"]([^'\"]+)['\"]\s*>\s*([-+0-9.eEdD]+)\s*</wgt>", s)
            if m:
                weights[m.group(1)] = float(m.group(2).replace("D", "E").replace("d", "e"))
            continue
        if s.startswith("<"):            # <mgrwt>, <scales>, <clustering>, ... keep, don't parse
            extra.append(s); continue
        body.append(_strip_comment(s))
    body = [b for b in body if b.strip()]
    if not body:
        raise ValueError("empty <event> block")
    h = body[0].split()
    nup = int(h[0])
    ev = Event(nup=nup, process=int(h[1]), weight=float(h[2]), scale=float(h[3]),
               alpha_qed=float(h[4]), alpha_qcd=float(h[5]), weights=weights,
               attributes=attrs, extra=extra)
    if len(body) < 1 + nup:
        raise ValueError(f"event declares NUP={nup} but has {len(body)-1} particle lines")
    for l in body[1:1 + nup]:
        x = l.replace("D", "E").replace("d", "e").split()   # Fortran 1.0D+00 style
        ev.particles.append(Particle(
            pid=int(x[0]), status=int(x[1]), mother1=int(x[2]), mother2=int(x[3]),
            color=int(x[4]), anticolor=int(x[5]),
            px=float(x[6]), py=float(x[7]), pz=float(x[8]), e=float(x[9]), m=float(x[10]),
            lifetime=float(x[11]), spin=float(x[12])))
    return ev


def _parse_initrwgt(lines: List[str]) -> Dict[str, str]:
    names: Dict[str, str] = {}
    group = ""
    for s in (l.strip() for l in lines):
        if s.startswith("<weightgroup"):
            m = re.search(r'name\s*=\s*["\']([^"\']*)["\']', s)
            group = m.group(1) if m else ""
        m = re.match(r"<weight\s+id\s*=\s*['\"]([^'\"]+)['\"][^>]*>(.*?)</weight>", s)
        if m:
            names[m.group(1)] = (group + " : " if group else "") + m.group(2).strip()
    return names


def read_lhe(f: TextIO) -> Tuple[Init, Iterator[Event]]:
    """Return (init, iterator over events). The iterator streams: a 20 GB file is fine."""
    header_lines: List[str] = []
    init_lines: List[str] = []
    initrwgt_lines: List[str] = []
    state = "top"
    line_no = 0
    # ---- phase 1: read up to the end of <init> (header + init are always first)
    for raw in f:
        line_no += 1
        s = raw.strip()
        if state == "top":
            if s.startswith("<header"):
                state = "header"
            elif s.startswith("<init"):
                state = "init"
        elif state == "header":
            if s.startswith("</header"):
                state = "top"
            else:
                header_lines.append(raw)
                if s.startswith("<initrwgt"):
                    state = "initrwgt"
        elif state == "initrwgt":
            header_lines.append(raw)
            if s.startswith("</initrwgt"):
                state = "header"
            else:
                initrwgt_lines.append(raw)
        elif state == "init":
            if s.startswith("</init"):
                break
            init_lines.append(raw)
    else:
        raise ValueError("no <init> block found — is this an LHE file?")
    init = _parse_init(init_lines, _parse_initrwgt(initrwgt_lines), "".join(header_lines))

    # ---- phase 2: a generator that yields events one at a time
    def events() -> Iterator[Event]:
        buf: List[str] = []
        attrs: Dict[str, str] = {}
        inside = False
        for raw in f:
            s = raw.strip()
            if not inside:
                if s.startswith("<event"):
                    inside = True
                    attrs = dict(_TAG_ATTR.findall(s))
                    buf = []
                # anything else at top level (blank lines, </LesHouchesEvents>) is ignored
            else:
                if s.startswith("</event"):
                    inside = False
                    yield _parse_event(buf, attrs)
                else:
                    buf.append(raw)
    return init, events()


def open_lhe(path: str) -> TextIO:
    """Transparently open plain or gzipped files in text mode."""
    return gzip.open(path, "rt") if path.endswith(".gz") else open(path, "r")


# ------------------------------------------------------------------ checks
def check_event(ev: Event, tol: float = 1e-3) -> List[str]:
    """Sanity checks that catch most broken files. Returns a list of complaints (empty = fine)."""
    out: List[str] = []
    n = len(ev.particles)
    if n != ev.nup:
        out.append(f"NUP={ev.nup} but {n} particles read")
    for i, p in enumerate(ev.particles, 1):
        if p.status not in (-9, -2, -1, 1, 2, 3):
            out.append(f"particle {i}: unknown status {p.status}")
        for mo in (p.mother1, p.mother2):
            if mo < 0 or mo > n:
                out.append(f"particle {i}: mother index {mo} out of range 1..{n}")
            elif mo >= i and mo != 0:
                out.append(f"particle {i}: mother {mo} does not precede it")
        # on-shell consistency: the stored mass should match E^2 - p^2
        m_calc = p.mass_from_p4
        if abs(m_calc - p.m) > tol * max(1.0, p.e):
            out.append(f"particle {i} (pid {p.pid}): m={p.m:.6g} but sqrt(E^2-p^2)={m_calc:.6g}")
    # four-momentum conservation: incoming == outgoing (status 1 only; resonances are intermediate)
    pin, pout = ev.p4_sum(ev.incoming()), ev.p4_sum(ev.final_state())
    for k, (a, b) in zip("px py pz E".split(), zip(pin, pout)):
        if abs(a - b) > tol * max(1.0, abs(a), abs(b)):
            out.append(f"{k} not conserved: in {a:.6g}, out {b:.6g}")
    # colour flow: every colour label must appear once as colour and once as anticolour
    # among the 'external' legs (incoming legs count with colour/anticolour swapped)
    seen: Dict[int, int] = {}
    for p in ev.particles:
        if p.status not in (-1, 1):
            continue
        c, a = (p.anticolor, p.color) if p.status == -1 else (p.color, p.anticolor)
        if c: seen[c] = seen.get(c, 0) + 1
        if a: seen[a] = seen.get(a, 0) - 1
    for lbl, bal in seen.items():
        if bal != 0:
            out.append(f"colour line {lbl} unbalanced ({bal:+d})")
    return out


def cross_section(init: Init, weights: List[float]) -> Tuple[float, float]:
    """Total cross section (pb) implied by a sample, and its MC statistical error.

    IDWTUP = ±3 : unweighted, |w| = 1 -> sigma is XSECUP from <init>; only the sign matters.
    IDWTUP = ±4 : weighted, sum(w)/N = sigma (pb) -> sigma = mean weight.
    IDWTUP = ±1,±2 : weights are only meaningful relative to XMAXUP; use XSECUP.
    """
    n = len(weights)
    strat = abs(init.weight_strategy)
    if strat == 4:
        mean = sum(weights) / n
        var = sum((w - mean) ** 2 for w in weights) / (n - 1) if n > 1 else 0.0
        return mean, math.sqrt(var / n)
    sigma = sum(p.xsec for p in init.processes)
    if strat == 3:
        # sign-only weights: sigma * <sign> with the same fractional error
        mean_sign = sum(math.copysign(1.0, w) for w in weights) / n
        return sigma * mean_sign, sigma / math.sqrt(n)
    return sigma, math.sqrt(sum(p.xsec_err ** 2 for p in init.processes))


# ------------------------------------------------------------------ CLI
if __name__ == "__main__":
    if len(sys.argv) < 2:
        print(__doc__); sys.exit(0)
    with open_lhe(sys.argv[1]) as fh:
        init, evs = read_lhe(fh)
        print(f"beams {init.beam_pids} at {init.beam_energies} GeV, PDF sets {init.pdf_sets}, IDWTUP={init.weight_strategy}")
        for p in init.processes:
            print(f"  process {p.lprup}: sigma = {p.xsec:.6g} +- {p.xsec_err:.2g} pb")
        if init.weight_names:
            print(f"  {len(init.weight_names)} named weights, e.g. {next(iter(init.weight_names.items()))}")
        n = neg = bad = 0
        ws: List[float] = []
        for ev in evs:
            n += 1; ws.append(ev.weight); neg += ev.weight < 0
            problems = check_event(ev)
            if problems and bad < 5:
                print(f"event {n}: " + "; ".join(problems))
            bad += bool(problems)
        sig, err = cross_section(init, ws)
        print(f"{n} events, {neg} with negative weight ({100*neg/max(n,1):.1f}%), {bad} failed checks")
        print(f"sigma from weights = {sig:.6g} +- {err:.2g} pb")

Using it

from lhe_reader import open_lhe, read_lhe, check_event, cross_section
import math

with open_lhe("unweighted_events.lhe.gz") as f:
    init, events = read_lhe(f)
    print(init.beam_energies, init.pdf_sets, init.weight_strategy)

    masses, weights, neg = [], [], 0
    for ev in events:
        mus = [p for p in ev.final_state() if abs(p.pid) == 13]
        if len(mus) == 2:
            px, py, pz, e = ev.p4_sum(mus)
            masses.append(math.sqrt(e*e - px*px - py*py - pz*pz))   # m(mu mu)
        weights.append(ev.weight)
        neg += ev.weight < 0
        # scale envelope from LHEF3 weights, if present
        if ev.weights:
            central = ev.weights.get("1", ev.weight)
            env = max(abs(w - central) for w in ev.weights.values())

print(f"sigma = {cross_section(init, weights)} pb, negative fraction {neg/len(weights):.3f}")

If you would rather not maintain your own, pylhe (Scikit-HEP) does the same job with the same design, returns Awkward arrays, and draws the event graph; pylhe.read_lhe_with_attributes() is the LHEF 3-aware entry point. Inside MadGraph, madgraph/various/lhe_parser.py is what the reweighting module uses. In C++, the header LHEF3.h shipped with Pythia (Leif Lönnblad's reference reader) is the canonical implementation.

11 · Validating what you read

Files break in predictable ways — a truncated write, a generator patch that changed a column, a merging script that mixed conventions — and the checks that catch them are cheap. check_event in the reader runs all of these:

• Particle count: the number of particle lines equals NUP.
• Status codes are in the allowed set; exactly two entries have status −1.
• Mother indices lie in 1…NUP and point to earlier lines.
• On-shell consistency: |√(E² − p²) − m| is small compared to E. Generators write the mass they used; if it disagrees with the four-momentum, one of them has been modified.
• Four-momentum conservation between incoming (status −1) and outgoing (status +1). Resonances are excluded from the sum; they are already represented by their daughters.
• Colour balance: after swapping colour and anticolour on incoming legs, every non-zero label appears once as colour and once as anticolour among external particles.

At the file level, count the events, confirm the negative-weight fraction is what the generator log claims, and check that the mean weight reproduces XSECUP for IDWTUP = ±4 samples. A file that fails these is not worth showering.

12 · Pitfalls, tools and further reading

• Units. GeV, pb, mm. Some older tools quote cross sections in mb (Pythia internally does); convert at the boundary and nowhere else.
• Order of the momentum components is (p_x, p_y, p_z, E, m). HepMC and ROOT's TLorentzVector also use (px, py, pz, E), but generator-internal formats and some textbooks put E first.
• Never derive the mass from E and p for light particles; use column 11. For a 1 TeV electron the difference is at the sixth digit and the stored value is the one the generator meant.
• Resonances are informational. Do not add status-2 entries to a momentum sum, and do not assume they are present: some generators write only the final state.
• Beams versus incoming partons. The proton beams (status −9) are almost never written; the two status −1 entries are partons after PDF sampling and carry only a fraction of the beam energy, which is why p_z does not sum to zero.
• Weight ids are strings, and the central weight is not always id "1". Read <initrwgt> to know which is which.
• Event attributes are strings too. npNLO="1" needs an int().
• Do not merge files by concatenation unless their init blocks are identical; the cross sections and IDWTUP must agree. MadGraph's merge.pl and pylhe's utilities do the bookkeeping.
• Showering an LHE file: in Pythia 8 set Beams:frameType = 4 and Beams:LHEF = file.lhe; in Herwig use LesHouchesFileReader; Sherpa reads its own and others' through EVENT_INPUT. Each honours SCALUP or the <scales> tag as the shower starting scale.

Further reading

• E. Boos et al., Generic user process interface for event generators, hep-ph/0109068 — the original accord.
• J. Alwall et al., A standard format for Les Houches Event Files, hep-ph/0609017 — LHEF 1.0.
• J. R. Andersen et al., Les Houches 2013: Physics at TeV Colliders: Standard Model Working Group Report, arXiv:1405.1067, section on LHEF 3.0 — weights, scales and attributes.
• A. Buckley et al., General-purpose event generators for LHC physics, arXiv:1101.2599 — what happens to the event after the shower reads it.
• The PDG Monte Carlo Particle Numbering Scheme for the identities in column 1.

Found an error? Write to saptaparnab@smu.edu. Back to the Newcomers page.