From 0112064af21b55a68a5c353d56517adca73c8543 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Tue, 15 Sep 2026 09:26:55 +0200
Subject: [PATCH 1/5] madspace: 2 -> 1 processes with PDFs, per-beam energies
 and PDFs

2 -> 1 (p p > z, p p > h): a collision Diagram may have one outgoing leg.
PhaseSpaceMapping integrates dPhi_1 = 2 pi delta(s_hat - m^2) against the
luminosity: s_hat is fixed at the mass and the single random number sets the
rapidity. Leptonic beams or a massless final state are refused.

run_card.toml [beam] now has ebeam1/ebeam2 and pdf1/pdf2. beam.e_cm and
beam.pdf stay readable and settable (both beams), and cards written with
e_cm/pdf are translated on read. The phase space is generated in the beams'
centre-of-mass frame (s = 4 E1 E2) and boosted by 0.5 ln(E1/E2) into the lab
frame before the cuts, so events are written, and eta cuts applied, in the
frame the user gave.

Asymmetric beams break the Integrand's mirror shortcut (mirroring accepted
lab events with unswapped PDFs). With mirror_beams the orientation is drawn
before the mapping, the mirror is taken in the beams' frame, and each leg
takes the PDF of the beam it comes from. Identical beams keep the old path.
The second beam's grid is passed to Integrand and SystematicsCalculator
(pdf_grid2 / nominal_pdf2); PDF member variations are dropped when the sets
differ, scale variations use both.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .../iolibs/template_files/mg7/gridpack.py     |   7 +-
 madgraph/iolibs/template_files/mg7/launch.py  |  99 ++++++++--
 .../iolibs/template_files/mg7/run_card.toml   |  11 +-
 madgraph/various/RunCardLO_to_MG7_mapping.md  |   4 +-
 madgraph/various/banner.py                    |  80 ++++++--
 .../include/madspace/driver/systematics.hpp   |   6 +-
 .../include/madspace/phasespace/integrand.hpp |  14 +-
 .../madspace/phasespace/phasespace.hpp        |  51 ++++-
 madspace/src/driver/systematics.cpp           |  77 +++++++-
 madspace/src/phasespace/integrand.cpp         |  73 +++++++-
 madspace/src/phasespace/phasespace.cpp        | 176 +++++++++++++++---
 madspace/src/phasespace/topology.cpp          |  14 +-
 madspace/src/python/madspace.cpp              |  30 ++-
 madspace/tests/test_beam_frame.py             | 144 ++++++++++++++
 tests/unit_tests/various/test_banner.py       |  54 +++++-
 15 files changed, 747 insertions(+), 93 deletions(-)
 create mode 100644 madspace/tests/test_beam_frame.py

diff --git a/madgraph/iolibs/template_files/mg7/gridpack.py b/madgraph/iolibs/template_files/mg7/gridpack.py
index da706c59af..0103caa37d 100644
--- a/madgraph/iolibs/template_files/mg7/gridpack.py
+++ b/madgraph/iolibs/template_files/mg7/gridpack.py
@@ -291,9 +291,12 @@ def load_systematics(run_card, backends=(), param_card_path=None):
     subproc_args = [
         ms.SubprocessSystArgs.from_json(json.dumps(a)) for a in data["subproc_args"]
     ]
-    nominal_pdf = None
+    nominal_pdf = nominal_pdf2 = None
     if config.has_pdf:
         nominal_pdf = ms.PdfGrid(_locate_pdf_file(data["nominal_grid_file"]))
+        # the second beam's set, recorded only when it differs from the first
+        if data.get("nominal_grid_file2"):
+            nominal_pdf2 = ms.PdfGrid(_locate_pdf_file(data["nominal_grid_file2"]))
     # The mu_R variations always need an alpha_s grid, so a gridpack that does
     # not name one cannot reweight anything -- drop the systematics rather than
     # die here (gridpacks written before the launcher recorded the file for a
@@ -336,7 +339,7 @@ def load_systematics(run_card, backends=(), param_card_path=None):
             ))
     systematics = ms.SystematicsCalculator(
         config, subproc_args, nominal_pdf, nominal_alpha_s,
-        context, matrix_elements, flavor_remap)
+        context, matrix_elements, flavor_remap, nominal_pdf2=nominal_pdf2)
     for warning in systematics.warnings:
         print("WARNING systematics: %s" % warning)
     return systematics
diff --git a/madgraph/iolibs/template_files/mg7/launch.py b/madgraph/iolibs/template_files/mg7/launch.py
index ed79f0b26e..1d2d872416 100644
--- a/madgraph/iolibs/template_files/mg7/launch.py
+++ b/madgraph/iolibs/template_files/mg7/launch.py
@@ -1,4 +1,5 @@
 import argparse
+import math
 import os
 import sys
 import time
@@ -269,6 +270,14 @@ def init_decay_mode(self) -> None:
         self.is_decay = incoming_counts == {1}
         if not self.is_decay:
             self.decaying_mass = None
+            single = [meta["outgoing"] for meta in self.subprocess_data
+                      if len(clean_pids(meta["outgoing"])) == 1]
+            if single and self.run_card["beam"]["leptonic"]:
+                raise ValueError(
+                    f"2 -> 1 processes ({single[0]} ...) need beams with a PDF "
+                    "(leptonic = false in [beam]): with fixed beam energies the "
+                    "one-particle phase space is a delta function in s_hat"
+                )
             return
 
         masses = {
@@ -499,7 +508,7 @@ def _no_pdf_message(self, pdf_set: str) -> str:
             found = ("no usable lhapdf-config was found, so the set could not "
                      "be downloaded automatically")
         return (
-            "LHAPDF set %r, requested by [beam] pdf in Cards/run_card.toml, "
+            "LHAPDF set %r, requested by [beam] pdf1/pdf2 in Cards/run_card.toml, "
             "was not found.\n"
             "  searched: %s\n"
             "  %s\n"
@@ -529,9 +538,19 @@ def init_beam(self) -> None:
             # No beams: the total energy is the decaying particle's mass, and
             # "leptonic" is what the mappings call "no parton luminosity".
             self.e_cm = self.decaying_mass
+            self.beam_energies = (float(self.e_cm), 0.0)
+            self.beam_rapidity = 0.0
+            self.asymmetric_beams = False
             self.leptonic = True
         else:
-            self.e_cm = beam_args["e_cm"]
+            ebeam1 = float(beam_args["ebeam1"])
+            ebeam2 = float(beam_args["ebeam2"])
+            self.beam_energies = (ebeam1, ebeam2)
+            # The phase space is generated in the beams' centre-of-mass frame
+            # (massless beams: s = 4 E1 E2) and boosted by this rapidity into
+            # the lab frame, where the events are written and the cuts applied.
+            self.e_cm = 2.0 * math.sqrt(ebeam1 * ebeam2)
+            self.beam_rapidity = 0.5 * math.log(ebeam1 / ebeam2)
             self.leptonic = beam_args["leptonic"]
 
         dynamical_scales = {
@@ -563,8 +582,10 @@ def init_beam(self) -> None:
                 fact_scale2=self.decaying_mass,
             )
             self.pdf_grid = None
+            self.pdf_grid2 = None
             self.lhapdf = None
             self.pdf_dir = None
+            self.pdf_dir2 = None
             self.fixed_alphas_info = self.write_fixed_alphas_info()
             self.alphas_grid = ms.AlphaSGrid(self.fixed_alphas_info)
             for context in self.contexts:
@@ -572,18 +593,40 @@ def init_beam(self) -> None:
             self.running_coupling = ms.RunningCoupling(self.alphas_grid)
             return
 
-        pdf_set = beam_args["pdf"]
+        pdf_set = beam_args["pdf1"]
         self.pdf_set = pdf_set
+        self.pdf_set2 = beam_args["pdf2"]
         self.fixed_alphas_info = None
         self.lhapdf = self.ensure_pdf_set(pdf_set)
         self.pdf_dir = self.lhapdf.find_set(pdf_set)
         if self.pdf_dir is None:
             raise RuntimeError(self._no_pdf_message(pdf_set))
         self.pdf_grid = ms.PdfGrid(os.path.join(self.pdf_dir, pdf_set, f"{pdf_set}_0000.dat"))
+        # A second grid only when the beams differ (and have a PDF at all);
+        # otherwise beam 2 shares the first one.
+        self.pdf_grid2 = None
+        self.pdf_dir2 = self.pdf_dir
+        if self.pdf_set2 != pdf_set and not self.leptonic:
+            self.pdf_dir2 = self.ensure_pdf_set(self.pdf_set2).find_set(self.pdf_set2)
+            if self.pdf_dir2 is None:
+                raise RuntimeError(self._no_pdf_message(self.pdf_set2))
+            self.pdf_grid2 = ms.PdfGrid(os.path.join(
+                self.pdf_dir2, self.pdf_set2, f"{self.pdf_set2}_0000.dat"))
+            logger.info("PDF sets: %s (beam 1), %s (beam 2); alpha_s from %s",
+                        pdf_set, self.pdf_set2, pdf_set)
         self.alphas_grid = ms.AlphaSGrid(os.path.join(self.pdf_dir, pdf_set, f"{pdf_set}.info"))
         for context in self.contexts:
             self.pdf_grid.initialize_globals(context)
+            if self.pdf_grid2 is not None:
+                # the prefix ms.Integrand reads the second beam's grid from
+                # (Integrand::pdf2_prefix)
+                self.pdf_grid2.initialize_globals(context, "beam2")
             self.alphas_grid.initialize_globals(context)
+        # Beams that differ in energy or PDF are not mirror images: a flavor
+        # standing for both orientations then needs the orientation drawn
+        # before the phase-space mapping (ms.PhaseSpaceMapping mirror_beams).
+        self.asymmetric_beams = (
+            self.beam_rapidity != 0.0 or self.pdf_grid2 is not None)
         self.running_coupling = ms.RunningCoupling(self.alphas_grid)
 
     # ------------------------------------------------------------------
@@ -761,7 +804,14 @@ def build_systematics(self):
             config.nominal_lhaid = info["SetIndex"]
             config.nominal_error_type = info["ErrorType"]
             config.nominal_description = info["SetDesc"]
-            config.pdf_members = self.resolve_pdf_variations()
+            if self.pdf_grid2 is None:
+                config.pdf_members = self.resolve_pdf_variations()
+            else:
+                logger.warning(
+                    "systematics: the beams use different PDF sets (%s, %s), so "
+                    "no PDF variation is computed; the scale variations use both "
+                    "sets", self.pdf_set, self.pdf_set2)
+                config.pdf_members = []
         else:
             config.pdf_members = []
         args = self.build_systematics_args()
@@ -773,7 +823,8 @@ def build_systematics(self):
             args, self.systematics_context)
         self.systematics = ms.SystematicsCalculator(
             config, args, self.pdf_grid, self.alphas_grid,
-            self.systematics_context, matrix_elements, flavor_remap)
+            self.systematics_context, matrix_elements, flavor_remap,
+            nominal_pdf2=self.pdf_grid2 if config.has_pdf else None)
         for warning in self.systematics.warnings:
             logger.warning("systematics: %s", warning)
         self.systematics_data = {
@@ -782,6 +833,9 @@ def build_systematics(self):
             "nominal_grid_file": os.path.join(
                 self.pdf_dir, self.pdf_set, f"{self.pdf_set}_0000.dat")
                 if config.has_pdf else None,
+            "nominal_grid_file2": os.path.join(
+                self.pdf_dir2, self.pdf_set2, f"{self.pdf_set2}_0000.dat")
+                if config.has_pdf and self.pdf_grid2 is not None else None,
             "nominal_info_file": info_path,
             # matrix elements re-evaluated for the mixed-order subprocesses
             "me_paths": [meta["me_path"] for meta in self.subprocess_data],
@@ -1510,7 +1564,6 @@ def _beam_info(self):
         `incoming` holds the *partonic* initial state (e.g. gluons), not the
         beam particle, so hadronic beams are protons (2212); leptonic beams are
         the incoming leptons themselves."""
-        half_e = float(self.e_cm) / 2.
         if self.is_decay:
             # No beams. LHE has no way to say that, so report the decaying
             # particle at rest as a single "beam"; the second slot is empty.
@@ -1519,22 +1572,23 @@ def _beam_info(self):
             return [pdg, 0], [float(self.e_cm), 0.0]
         if not self.leptonic:
             # hadronic collider: proton beams (p-pbar is not distinguished)
-            return [2212, 2212], [half_e, half_e]
+            return [2212, 2212], list(self.beam_energies)
         data = self.subprocess_data[0]
         incoming = data["incoming"]
         flavor0 = data["flavors"][0]["options"][0]
         init_pdgs = flavor0[:len(incoming)]
         beam_pdgs = [pdg if abs(code) in (81, 82) else code
                      for code, pdg in zip(incoming, init_pdgs)]
-        return beam_pdgs, [half_e, half_e]
+        return beam_pdgs, list(self.beam_energies)
 
-    def _lhapdf_id(self):
-        """Central LHAPDF id of the beam PDF set (read from its .info SetIndex),
-        or -1 for a leptonic beam (no PDF)."""
+    def _lhapdf_id(self, beam=1):
+        """Central LHAPDF id of the PDF set of `beam` (read from its .info
+        SetIndex), or -1 for a leptonic beam (no PDF)."""
         if self.leptonic or not self.pdf_dir:
             return -1
-        pdf_set = self.run_card["beam"]["pdf"]
-        info = os.path.join(self.pdf_dir, pdf_set, "%s.info" % pdf_set)
+        pdf_set = self.pdf_set if beam == 1 else self.pdf_set2
+        pdf_dir = self.pdf_dir if beam == 1 else self.pdf_dir2
+        info = os.path.join(pdf_dir, pdf_set, "%s.info" % pdf_set)
         try:
             for line in open(info):
                 if line.strip().startswith("SetIndex:"):
@@ -1548,7 +1602,6 @@ def build_lhe_meta(self):
         run_card.toml (<MG7RunCard>) headers plus the beam/PDF/cross-section info
         needed by downstream tools (systematics, MadSpin, ...)."""
         beam_pdgs, energies = self._beam_info()
-        lhaid = self._lhapdf_id()
         pdf_group = -1 if self.leptonic else 0
         status = self.event_generator.status()
         xsec, err = status.mean, status.error
@@ -1571,7 +1624,7 @@ def build_lhe_meta(self):
             beam1_pdg_id=beam_pdgs[0], beam2_pdg_id=beam_pdgs[1],
             beam1_energy=energies[0], beam2_energy=energies[1],
             beam1_pdf_authors=pdf_group, beam2_pdf_authors=pdf_group,
-            beam1_pdf_id=lhaid, beam2_pdf_id=lhaid,
+            beam1_pdf_id=self._lhapdf_id(1), beam2_pdf_id=self._lhapdf_id(2),
             weight_mode=3,
             # positional: the pybind arg name for max_weight is non-kwarg-safe
             processes=[ms.LHEProcess(xsec, err, xsec, 1)],
@@ -2034,6 +2087,8 @@ def build_multichannel_phasespace(self) -> PhaseSpace:
                     invariant_power=self.process.run_card["phasespace"]["invariant_power"],
                     permutations=chan_permutations,
                     leptonic=self.process.leptonic,
+                    beam_rapidity=self.process.beam_rapidity,
+                    mirror_beams=self.process.asymmetric_beams,
                 )
                 prefix = f"subproc{self.subproc_id}.channel{channel_id}"
                 if topo_count > 1:
@@ -2094,6 +2149,8 @@ def build_flat_phasespace(self) -> PhaseSpace:
             mode=self.t_channel_mode(self.process.run_card["phasespace"]["flat_mode"]),
             cuts=self.cuts,
             leptonic=self.process.leptonic,
+            beam_rapidity=self.process.beam_rapidity,
+            mirror_beams=self.process.asymmetric_beams,
         )
         prefix = f"subproc{self.subproc_id}.flat"
         discrete_sym, discrete_flavor = self.build_discrete(
@@ -2627,7 +2684,8 @@ def build_integrands(
                 flavor_diff_xs_indices,
                 flavor_subproc_indices,
                 flavor_per_subproc_remap,
-                madnis_args["compressed_channel_weight_count"]
+                madnis_args["compressed_channel_weight_count"],
+                None if self.process.leptonic else self.process.pdf_grid2,
             ))
         #print(integrands[1].function())
         #for i in integrands: print(i.function())
@@ -3030,6 +3088,15 @@ def do_set(self, line, *args, **kwargs):
                         self.modified_card.add("run")
                         return
 
+                    # quantities the card stores per beam: set both beams
+                    if rest and nlow in ("e_cm", "beam.e_cm", "beam.pdf"):
+                        value = rest if nlow == "beam.pdf" else run_card.evaluate(rest, masses)
+                        run_card.set(nlow, value, user=True)
+                        logger.info("set %s of the run_card.toml to %s (both beams)",
+                                    nlow, run_card["beam." + nlow.split(".")[-1]])
+                        self.modified_card.add("run")
+                        return
+
                     # legacy madevent run_card names -> mg7 "section.key", so a
                     # madevent-style launch script ("set nevents 500", "set
                     # use_syst F", "set bwcutoff 10", ...) edits the mg7
diff --git a/madgraph/iolibs/template_files/mg7/run_card.toml b/madgraph/iolibs/template_files/mg7/run_card.toml
index 327842f8d6..95a10fc2c7 100644
--- a/madgraph/iolibs/template_files/mg7/run_card.toml
+++ b/madgraph/iolibs/template_files/mg7/run_card.toml
@@ -34,9 +34,16 @@ include_madspace = %(gridpack.include_madspace)s
 include_madspace_source = %(gridpack.include_madspace_source)s
 
 [beam]
-e_cm = %(beam.e_cm)s
+# energy of each beam in GeV, beam 1 moving along +z. The events are written in
+# this (lab) frame and the eta cuts apply in it; the collision energy is
+# 2 sqrt(ebeam1 ebeam2). "set e_cm X" sets both beams to X/2.
+ebeam1 = %(beam.ebeam1)s
+ebeam2 = %(beam.ebeam2)s
 leptonic = %(beam.leptonic)s
-pdf = %(beam.pdf)s
+# LHAPDF set of each beam. "set beam.pdf X" sets both. PDF member variations
+# in [systematics] need the same set on both beams.
+pdf1 = %(beam.pdf1)s
+pdf2 = %(beam.pdf2)s
 fixed_ren_scale = %(beam.fixed_ren_scale)s
 fixed_fact_scale = %(beam.fixed_fact_scale)s
 ren_scale = %(beam.ren_scale)s
diff --git a/madgraph/various/RunCardLO_to_MG7_mapping.md b/madgraph/various/RunCardLO_to_MG7_mapping.md
index 5e9c50954d..ffa8115789 100644
--- a/madgraph/various/RunCardLO_to_MG7_mapping.md
+++ b/madgraph/various/RunCardLO_to_MG7_mapping.md
@@ -24,7 +24,7 @@ Legend for the mapping:
 
 | LO parameter(s) | MG7 target | class | notes |
 |---|---|---|---|
-| `ebeam1`, `ebeam2` | `beam.e_cm` | [~] | `e_cm = ebeam1 + ebeam2` for a head-on collider. Asymmetric beams keep e_cm but lose the boost. |
+| `ebeam1`, `ebeam2` | `beam.ebeam1`, `beam.ebeam2` | [x] | One energy per beam; asymmetric beams keep their boost (events written and eta cuts applied in the lab frame). `beam.e_cm = 2 sqrt(ebeam1 ebeam2)` is derived. |
 | `lpp1`, `lpp2` | `beam.leptonic` | [!] | `leptonic = |lpp| in {3,4}` (e/µ). MG7 only stores hadronic-vs-leptonic; it cannot express antiproton vs proton (`lpp=±1`), elastic photon (`lpp=2`), EVA (`lpp=±3/±4` variants), plugin (`lpp=9`) or the no-PDF fixed-energy beam (`lpp=0`). |
 | `polbeam1`, `polbeam2` | — | [x] | beam polarization: no MG7 field. |
 | `nb_proton1/2`, `nb_neutron1/2`, `mass_ion1/2` | — | [x] | heavy-ion beams: no MG7 field. |
@@ -137,7 +137,7 @@ simplified_channel_count,decays}`, all of `[madnis]`.
 ## Summary / recommendation
 
 **Safe to port now (the converter should cover these):**
-- Beams: `ebeam1/2`→`e_cm`, `lpp`→`leptonic` (common case).
+- Beams: `ebeam1/2`→`ebeam1/2`, `lpp`→`leptonic` (common case).
 - PDF: `lhaid`/`pdlabel`→`pdf` (with a label/id → name table).
 - Scales: `fixed_ren_scale`, `scale`→`ren_scale`, `dsqrt_q2fact1/2`→`fact_scale1/2`,
   `dynamical_scale_choice` (int→str table), `fixed_fac_scale`.
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index c946245dcf..5dd1dc0f05 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -6451,7 +6451,7 @@ class RunCardMG7(RunCard):
                     'gev': 1.0, 'tev': 1e3, 'pev': 1e6}
     # typed parameters carrying an energy dimension (stored in GeV); a value
     # like "13 TeV" is converted to 13000.0 when assigned to one of these.
-    energy_params = ('beam.e_cm', 'beam.ren_scale',
+    energy_params = ('beam.ebeam1', 'beam.ebeam2', 'beam.e_cm', 'beam.ren_scale',
                      'beam.fact_scale1', 'beam.fact_scale2')
     # name -> pdg for masses that can be referenced in expressions (e.g. a scale
     # set to "mz/2"); resolved against the param_card by get_mass_shortcuts.
@@ -6556,7 +6556,11 @@ def default_setup(self):
         self.add_toml_param('gridpack', 'include_madspace_source', False)
 
         # ----------------------------- [beam] -------------------------
-        self.add_toml_param('beam', 'e_cm', 13000.0)
+        # One energy per beam, beam 1 along +z. The events are written in this
+        # lab frame; beam.e_cm = 2 sqrt(ebeam1 ebeam2) is derived from them
+        # (and setting it sets both beams, see set()).
+        self.add_toml_param('beam', 'ebeam1', 6500.0)
+        self.add_toml_param('beam', 'ebeam2', 6500.0)
         self.add_toml_param('beam', 'leptonic', False)
         # NNPDF4.0 LO, 5-flavour scheme, alpha_s(M_Z) = 0.118 -- the same set
         # the legacy LO run_card now defaults to (lhaid 331900). It carries a
@@ -6566,7 +6570,9 @@ def default_setup(self):
         # NNPDF4.0 has no 4-flavour LO counterpart, so unlike the NLO card there
         # is no scheme-dependent choice here: this set is used whatever the
         # b-quark treatment.
-        self.add_toml_param('beam', 'pdf', "NNPDF40_lo_as_01180")
+        # one set per beam; beam.pdf reads/sets both at once
+        self.add_toml_param('beam', 'pdf1', "NNPDF40_lo_as_01180")
+        self.add_toml_param('beam', 'pdf2', "NNPDF40_lo_as_01180")
         # Default to the dynamical scale set by dynamical_scale_choice below
         # (half_transverse_mass, i.e. HT/2) rather than to the fixed ren_scale
         # / fact_scale values. Those fixed values are kept as the fallback used
@@ -6746,8 +6752,23 @@ def __getitem__(self, name):
                 return TOMLSectionView(self, lname)
             if lname in self.legacy_compat_keys:
                 return self._legacy_compat(lname)
+            if lname in self.derived_beam_keys:
+                return self._derived_beam_value(lname)
         return super(RunCardMG7, self).__getitem__(name)
 
+    # [beam] quantities kept readable (and settable, see set()) although the
+    # card stores them per beam
+    derived_beam_keys = ('beam.e_cm', 'beam.pdf')
+
+    def _derived_beam_value(self, lname):
+        if lname == 'beam.e_cm':
+            return 2. * (float(self['beam.ebeam1']) * float(self['beam.ebeam2'])) ** 0.5
+        pdf1, pdf2 = self['beam.pdf1'], self['beam.pdf2']
+        if pdf1 != pdf2:
+            raise KeyError("beam.pdf is ambiguous: the beams use different PDF "
+                           "sets (%s, %s), read beam.pdf1 / beam.pdf2" % (pdf1, pdf2))
+        return pdf1
+
     get = __getitem__
 
     def __setitem__(self, name, value, *args, **opts):
@@ -6803,7 +6824,6 @@ def __setitem__(self, name, value, *args, **opts):
     def _legacy_compat(self, key):
         beam = self['beam']
         leptonic = bool(beam['leptonic'])
-        half_e = float(beam['e_cm']) / 2.
         if key == 'nevents':
             return self['generation']['events']
         if key in ('pdlabel', 'pdlabel1', 'pdlabel2'):
@@ -6812,7 +6832,7 @@ def _legacy_compat(self, key):
         if key in ('lpp1', 'lpp2'):
             return 0 if leptonic else 1
         if key in ('ebeam1', 'ebeam2'):
-            return half_e
+            return float(beam[key])
         if key in ('nb_proton1', 'nb_proton2'):
             return 0 if leptonic else 1
         if key in ('nb_neutron1', 'nb_neutron2'):
@@ -6852,11 +6872,11 @@ def _legacy_compat(self, key):
             return -2            # -2: reuse iseed for the python RNG
         raise KeyError(key)
 
-    def get_lhapdf_id(self):
-        """Central LHAPDF id of the beam PDF set, read from the set's .info
-        (SetIndex) so it does not depend on lhapdf being able to resolve the
-        set by name; falls back to lhapdf.getPDFSet, then 0."""
-        name = self['beam']['pdf']
+    def get_lhapdf_id(self, beam=1):
+        """Central LHAPDF id of the PDF set of `beam` (1 or 2), read from the
+        set's .info (SetIndex) so it does not depend on lhapdf being able to
+        resolve the set by name; falls back to lhapdf.getPDFSet, then 0."""
+        name = self['beam']['pdf%d' % beam]
         if not name:
             return 0
         search = []
@@ -6925,7 +6945,10 @@ def read_data(self, data, unknown_warning=True):
             elif sl in self.toml_sections:
                 for key, value in content.items():
                     internal = '%s.%s' % (sl, key.lower())
-                    if internal in self:
+                    if internal in self.derived_beam_keys:
+                        # a card written before the per-beam energies/PDFs
+                        self.set(internal, value, user=True)
+                    elif internal in self:
                         self.set(internal, value, user=True)
                     else:
                         if unknown_warning:
@@ -6967,6 +6990,27 @@ def set(self, name, value, *args, **opts):
         across all sections."""
         if isinstance(name, str):
             lname = name.lower()
+            if lname == 'e_cm':
+                lname = name = 'beam.e_cm'
+            if lname == 'pdf':
+                logger.warning("Ambiguous key 'pdf' — use beam.pdf (both beams), "
+                               "beam.pdf1, beam.pdf2 or systematics.pdf")
+                return
+            if lname == 'beam.e_cm':
+                value = self.parse_energy(value)
+                if isinstance(value, str):
+                    if value.strip().lower().startswith('scan'):
+                        raise InvalidRunCard(
+                            "e_cm cannot be scanned: scan beam.ebeam1 and "
+                            "beam.ebeam2 together instead (scan1:[...] on both)")
+                    value = self.format_variable(value, float, name='e_cm')
+                for key in ('beam.ebeam1', 'beam.ebeam2'):
+                    super(RunCardMG7, self).set(key, float(value) / 2., *args, **opts)
+                return
+            if lname == 'beam.pdf':
+                for key in ('beam.pdf1', 'beam.pdf2'):
+                    super(RunCardMG7, self).set(key, value, *args, **opts)
+                return
             if '.' not in lname:
                 matches = ['%s.%s' % (sec, lname)
                            for sec, keys in self.toml_sections.items()
@@ -7099,7 +7143,7 @@ def set_collider(self, name, arg, masses=None):
             val = self.format_variable(val, float, name=name)
         if not had_unit and name in ('lhc', 'lcc'):
             val = val * 1000.0  # unit-less hadron collider energy is in TeV
-        self.set('beam.e_cm', float(val), user=True)
+        self.set('beam.e_cm', float(val), user=True)  # both beams at val/2
         self.set('beam.leptonic', name in ('lep', 'ilc'), user=True)
         return float(val)
 
@@ -7125,6 +7169,10 @@ def check_validity(self):
         """Minimal consistency checks for the TOML run_card."""
         if self['generation']['survey_min_iters'] > self['generation']['survey_max_iters']:
             raise InvalidRunCard("survey_min_iters can not be larger than survey_max_iters")
+        for beam in (1, 2):
+            energy = self['beam']['ebeam%d' % beam]
+            if isinstance(energy, (int, float)) and energy <= 0:
+                raise InvalidRunCard("beam.ebeam%d must be positive (got %s)" % (beam, energy))
 
         # 'device' is list-valued and accepts a "<type>:<index>" syntax, so the
         # generic 'allowed' machinery cannot check it on its own.
@@ -7275,7 +7323,8 @@ def create_default_for_process(self, proc_characteristic, history, proc_def):
         if beam_id and beam_id <= leptons:
             # lepton collider: no PDF, lower default energy, no hadronic cuts
             self['beam']['leptonic'] = True
-            self['beam']['e_cm'] = 1000.0
+            self['beam']['ebeam1'] = 500.0
+            self['beam']['ebeam2'] = 500.0
             self['beam']['fixed_fact_scale'] = True
             self.remove_jet_cuts()
         elif beam_id and not (beam_id & hadronic):
@@ -7309,6 +7358,8 @@ def remove_jet_cuts(self):
         'nevents': 'generation.events',
         'gridpack': 'gridpack.save_gridpack',
         'fixed_ren_scale': 'beam.fixed_ren_scale',
+        'ebeam1': 'beam.ebeam1',
+        'ebeam2': 'beam.ebeam2',
         'scale': 'beam.ren_scale',
         'dsqrt_q2fact1': 'beam.fact_scale1',
         'dsqrt_q2fact2': 'beam.fact_scale2',
@@ -7435,7 +7486,6 @@ def nondefault(name):
                 mg7.set(mg7key, lo[loname])
 
         # --- beams ---
-        mg7.set('beam.e_cm', float(lo['ebeam1']) + float(lo['ebeam2']))
         lpps = [lo['lpp1'], lo['lpp2']]
         is_lep = lambda l: isinstance(l, int) and abs(l) in (3, 4)
         # lpp 0 = fixed-energy beam with no PDF, typically a lepton collider
@@ -7447,7 +7497,7 @@ def nondefault(name):
             elif isinstance(lpp, int) and lpp < 0:
                 dropped.append('lpp%d=%s (antiparticle beam: mg7 keeps the particle PDF)' % (i, lpp))
             elif lpp == 0:
-                dropped.append('lpp%d=0 (no-PDF fixed-energy beam: check beam.leptonic/beam.pdf)' % i)
+                dropped.append('lpp%d=0 (no-PDF fixed-energy beam: check beam.leptonic/beam.pdf1/beam.pdf2)' % i)
 
         # --- scales ---
         if 'fixed_fac_scale' in lo:
diff --git a/madspace/include/madspace/driver/systematics.hpp b/madspace/include/madspace/driver/systematics.hpp
index 3f817da963..5d2c139a97 100644
--- a/madspace/include/madspace/driver/systematics.hpp
+++ b/madspace/include/madspace/driver/systematics.hpp
@@ -107,7 +107,10 @@ class SystematicsCalculator {
         const std::optional<AlphaSGrid>& nominal_alpha_s,
         ContextPtr context = nullptr,
         const std::vector<std::optional<MatrixElement>>& matrix_elements = {},
-        const nested_vector2<me_int_t>& me_flavor_remap = {}
+        const nested_vector2<me_int_t>& me_flavor_remap = {},
+        // the second beam's PDF when it differs from nominal_pdf (the first
+        // beam's); PDF member variations are then dropped with a warning
+        const std::optional<PdfGrid>& nominal_pdf2 = std::nullopt
     );
 
     const SystematicsConfig& config() const { return _config; }
@@ -179,6 +182,7 @@ class SystematicsCalculator {
     ContextPtr _context;
     std::string _prefix;
     std::optional<PdfEvaluator> _nominal_pdf;   // has_pdf only
+    std::optional<PdfEvaluator> _nominal_pdf2;  // second beam, if its set differs
     std::vector<PdfEvaluator> _member_pdfs;     // one per member
     std::vector<RuntimePtr> _alpha_s_runtimes;  // index 0: nominal set
     std::vector<std::optional<MatrixElementData>> _matrix_elements;
diff --git a/madspace/include/madspace/phasespace/integrand.hpp b/madspace/include/madspace/phasespace/integrand.hpp
index 40f9fe8365..f9d93dadad 100644
--- a/madspace/include/madspace/phasespace/integrand.hpp
+++ b/madspace/include/madspace/phasespace/integrand.hpp
@@ -67,8 +67,13 @@ class Integrand : public FunctionGenerator {
         const std::vector<std::size_t>& flavor_diff_xs_indices = {},
         const std::vector<std::size_t>& flavor_subproc_indices = {},
         const std::vector<std::size_t>& flavor_per_subproc_remap = {},
-        std::size_t compressed_channel_weight_count = 50
+        std::size_t compressed_channel_weight_count = 50,
+        // PDF of the second beam when it differs from pdf_grid (which is then
+        // the first beam's). Its globals must be initialized with the prefix
+        // Integrand::pdf2_prefix.
+        const std::optional<PdfGrid>& pdf_grid2 = std::nullopt
     );
+    inline static const std::string pdf2_prefix = "beam2";
     std::size_t particle_count() const { return _mapping.particle_count(); }
     bool madnis_training() const { return _madnis_training; }
     std::optional<std::string> vegas_grid_name() const {
@@ -126,6 +131,9 @@ class Integrand : public FunctionGenerator {
     AdaptiveDiscrete _discrete_flavor;
     nested_vector2<me_int_t> _pid_options;
     std::array<std::optional<PartonDensity>, 2> _pdfs;
+    // with different PDFs per beam: each leg's density from the other beam,
+    // for the events whose orientation is mirrored
+    std::array<std::optional<PartonDensity>, 2> _pdfs_swapped;
     std::array<std::vector<me_int_t>, 2> _pdf_indices;
     std::optional<RunningCoupling> _running_coupling;
     std::optional<EnergyScale> _energy_scale;
@@ -148,7 +156,11 @@ class Integrand : public FunctionGenerator {
     std::vector<me_int_t> _flavor_remap;
     std::vector<double> _flavor_factors;
     std::vector<me_int_t> _flavor_mirror;
+    std::vector<double> _flavor_mirror_factors;
     bool _has_mirror;
+    // Asymmetric beams: the orientation is drawn before the phase-space
+    // mapping (mirror_index condition) instead of mirroring accepted events.
+    bool _mirror_before_cuts;
     NamedVector<Type> _channel_part_ret_types;
     std::vector<me_int_t> _flavor_diff_xs_indices;
     std::vector<me_int_t> _flavor_subproc_indices;
diff --git a/madspace/include/madspace/phasespace/phasespace.hpp b/madspace/include/madspace/phasespace/phasespace.hpp
index 2851627607..073a06fa82 100644
--- a/madspace/include/madspace/phasespace/phasespace.hpp
+++ b/madspace/include/madspace/phasespace/phasespace.hpp
@@ -25,7 +25,9 @@ class PhaseSpaceMapping : public Mapping {
         TChannelMode t_channel_mode = propagator,
         const std::optional<Cuts>& cuts = std::nullopt,
         const std::vector<std::vector<std::size_t>>& permutations = {},
-        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt
+        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt,
+        double beam_rapidity = 0.,
+        bool mirror_beams = false
     );
 
     PhaseSpaceMapping(
@@ -35,13 +37,37 @@ class PhaseSpaceMapping : public Mapping {
         double invariant_power = 0.8,
         TChannelMode mode = rambo,
         const std::optional<Cuts>& cuts = std::nullopt,
-        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt
+        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt,
+        double beam_rapidity = 0.,
+        bool mirror_beams = false
     );
 
+    // mirror_beams: add a "mirror_index" condition (after permutation_index).
+    // For index 1 the event is rotated by pi about the x axis in the beams'
+    // centre-of-mass frame -- leg 1 then comes from beam 2 and leg 2 from beam
+    // 1 -- before the boost into the lab frame and the cuts. With identical
+    // beams that is a rotation of the lab event as well, which is why the
+    // Integrand can otherwise mirror accepted events after the cuts.
+    //
+    // beam_rapidity: rapidity of the beams' centre-of-mass frame in the frame
+    // the momenta are returned in (the lab frame), 0.5 ln(E1 / E2) for beam
+    // energies E1 (along +z) and E2. The phase space is generated in the beams'
+    // centre-of-mass frame, of energy cm_energy = 2 sqrt(E1 E2), and boosted
+    // into the lab frame before the cuts, so rapidity cuts act in the lab.
+
     // A 1 -> n decay and a leptonic (fixed-s) 2 -> n collision both have 3n-4
     // degrees of freedom; a hadronic 2 -> n adds the beam momentum fractions,
     // i.e. one further sampled invariant on top of s_hat -> 3n-2.
+    // A hadronic 2 -> 1 collision keeps a single random number, the rapidity
+    // of the produced particle; s_hat is fixed at its mass.
     static std::size_t random_dim_for(const Topology& topology, bool leptonic) {
+        if (topology.outgoing_masses().size() == 1 &&
+            (leptonic || topology.is_decay())) {
+            throw std::invalid_argument(
+                "a 2 -> 1 process needs beams with a PDF: with fixed beam "
+                "energies s_hat is fixed and there is nothing to integrate"
+            );
+        }
         return 3 * topology.outgoing_masses().size() -
             ((leptonic || topology.is_decay()) ? 4 : 2);
     }
@@ -51,6 +77,8 @@ class PhaseSpaceMapping : public Mapping {
         return _topology.outgoing_masses().size() + _topology.incoming_masses().size();
     }
     std::size_t channel_count() const { return _permutations.size(); }
+    double beam_rapidity() const { return _beam_rapidity; }
+    bool mirror_beams() const { return _mirror_beams; }
 
 private:
     Result build_forward_impl(
@@ -63,13 +91,32 @@ class PhaseSpaceMapping : public Mapping {
         const NamedVector<Value>& inputs,
         const NamedVector<Value>& conditions
     ) const override;
+    // partonic centre-of-mass frame <-> lab frame (partonic boost, mirror,
+    // beam_rapidity boost)
+    Value to_lab(
+        FunctionBuilder& fb,
+        Value momenta,
+        Value x1,
+        Value x2,
+        const NamedVector<Value>& conditions
+    ) const;
+    Value from_lab(
+        FunctionBuilder& fb,
+        Value momenta,
+        Value x1,
+        Value x2,
+        const NamedVector<Value>& conditions
+    ) const;
 
     Topology _topology;
     Cuts _cuts;
     double _pi_factors;
     double _sqrt_s_lab;
+    double _beam_rapidity;
+    bool _mirror_beams;
     bool _leptonic;
     bool _map_luminosity;
+    bool _two_to_one;
     std::size_t _n_discrete;
     std::vector<Invariant> _s_invariants;
     std::variant<
diff --git a/madspace/src/driver/systematics.cpp b/madspace/src/driver/systematics.cpp
index 2d9cd96108..ce8f54259f 100644
--- a/madspace/src/driver/systematics.cpp
+++ b/madspace/src/driver/systematics.cpp
@@ -173,7 +173,8 @@ SystematicsCalculator::SystematicsCalculator(
     const std::optional<AlphaSGrid>& nominal_alpha_s,
     ContextPtr context,
     const std::vector<std::optional<MatrixElement>>& matrix_elements,
-    const nested_vector2<me_int_t>& me_flavor_remap
+    const nested_vector2<me_int_t>& me_flavor_remap,
+    const std::optional<PdfGrid>& nominal_pdf2
 ) :
     _config(config),
     _subproc_args(subproc_args),
@@ -222,6 +223,10 @@ SystematicsCalculator::SystematicsCalculator(
     std::sort(pids.begin(), pids.end());
     if (_config.has_pdf) {
         _nominal_pdf = make_pdf_evaluator(nominal_pdf.value(), pids, "nominal", 0);
+        if (nominal_pdf2) {
+            _nominal_pdf2 =
+                make_pdf_evaluator(nominal_pdf2.value(), pids, "nominal2", 0);
+        }
     }
 
     // matrix elements for the subprocesses with mixed alpha_s powers
@@ -276,6 +281,16 @@ SystematicsCalculator::SystematicsCalculator(
             ));
             continue;
         }
+        if (_nominal_pdf2) {
+            // a member of one set cannot stand in for two different sets
+            _warnings.push_back(std::format(
+                "PDF variation {} member {} ignored: the two beams use different "
+                "PDF sets",
+                spec.set_name,
+                spec.member
+            ));
+            continue;
+        }
         bool same_set = spec.set_lhaid == _config.nominal_lhaid;
         if (!same_set && !_mur_supported) {
             _warnings.push_back(std::format(
@@ -628,6 +643,12 @@ SystematicsCalculator::reweight_info(EventBuffer& buffer, std::size_t event_inde
         .fact_scale1 = has_beam1 ? event.fact_scale1().value() : 0.,
         .fact_scale2 = has_beam2 ? event.fact_scale2().value() : 0.,
     };
+    // beam 1 moves along +z: leg 1 of a mirrored event belongs to beam 2
+    if (has_beam1 && has_beam2 && buffer.particle(event_index, 0).pz() < 0.) {
+        std::swap(info.pdg1, info.pdg2);
+        std::swap(info.x1, info.x2);
+        std::swap(info.fact_scale1, info.fact_scale2);
+    }
     return info;
 }
 
@@ -654,6 +675,9 @@ void SystematicsCalculator::compute(
         int subproc, qcd_power;
         me_int_t slot1, slot2;
         bool use_me;
+        // leg 1 comes from beam 2 (mirrored event); only tracked when the
+        // beams have different PDFs
+        bool swapped;
         std::array<double, 5> dyn_scale;
     };
     std::vector<EventInput> inputs(count);
@@ -677,6 +701,7 @@ void SystematicsCalculator::compute(
         in.muf2 = has_beam2 ? event.fact_scale2().value() : 0.;
         in.slot1 = in.slot2 = 0;
         in.nominal_product = 1.;
+        in.swapped = _nominal_pdf2 && buffer.particle(i, 0).pz() < 0.;
         if (use_pdf) {
             auto& pdgs = args.beam_pdgs.at(event.flavor_index());
             in.slot1 = pid_slot(pdgs.at(0));
@@ -696,20 +721,57 @@ void SystematicsCalculator::compute(
     auto scale_of = [&](const EventInput& in, int dyn, double generated) {
         return dyn == -1 ? generated : in.dyn_scale[dyn];
     };
+    // x f(x, q) at the given points, point n evaluated with `second` if
+    // on_second[n] and with `first` otherwise. One batched call when both are
+    // the same grid.
+    auto evaluate_beams = [&](const PdfEvaluator& first,
+                              const PdfEvaluator& second,
+                              const std::vector<double>& x,
+                              const std::vector<double>& q,
+                              const std::vector<me_int_t>& slots,
+                              const std::vector<char>& on_second) {
+        if (&first == &second) {
+            return evaluate_pdf(first, x, q, slots);
+        }
+        std::array<std::vector<double>, 2> grid_x, grid_q;
+        std::array<std::vector<me_int_t>, 2> grid_slots;
+        for (std::size_t n = 0; n < x.size(); ++n) {
+            grid_x[on_second[n]].push_back(x[n]);
+            grid_q[on_second[n]].push_back(q[n]);
+            grid_slots[on_second[n]].push_back(slots[n]);
+        }
+        std::array<std::vector<double>, 2> grid_values{
+            evaluate_pdf(first, grid_x[0], grid_q[0], grid_slots[0]),
+            evaluate_pdf(second, grid_x[1], grid_q[1], grid_slots[1]),
+        };
+        std::vector<double> values(x.size());
+        std::array<std::size_t, 2> next{0, 0};
+        for (std::size_t n = 0; n < x.size(); ++n) {
+            values[n] = grid_values[on_second[n]][next[on_second[n]]++];
+        }
+        return values;
+    };
+    const PdfEvaluator* nominal_pdf2 =
+        _nominal_pdf2 ? &_nominal_pdf2.value() : (_nominal_pdf ? &_nominal_pdf.value() : nullptr);
 
     // Nominal PDF product when the events do not carry it: one batched call
     if (use_pdf && !has_partial) {
         std::vector<double> x, q;
         std::vector<me_int_t> slots;
+        std::vector<char> on_second;
         for (auto& in : inputs) {
             if (has_beam1) {
                 x.push_back(in.x1), q.push_back(in.muf1), slots.push_back(in.slot1);
+                on_second.push_back(in.swapped);
             }
             if (has_beam2) {
                 x.push_back(in.x2), q.push_back(in.muf2), slots.push_back(in.slot2);
+                on_second.push_back(!in.swapped);
             }
         }
-        auto values = evaluate_pdf(_nominal_pdf.value(), x, q, slots);
+        auto values = evaluate_beams(
+            _nominal_pdf.value(), *nominal_pdf2, x, q, slots, on_second
+        );
         for (std::size_t i = 0, n = 0; i < count; ++i) {
             double product = 1.;
             if (has_beam1) {
@@ -727,9 +789,11 @@ void SystematicsCalculator::compute(
     std::vector<std::vector<double>> r_pdf(count, std::vector<double>(var_count, 1.));
     if (use_pdf) {
         auto evaluate_variations = [&](const PdfEvaluator& evaluator,
+                                       const PdfEvaluator& evaluator2,
                                        const std::vector<std::size_t>& var_indices) {
             std::vector<double> x, q;
             std::vector<me_int_t> slots;
+            std::vector<char> on_second;
             for (auto& in : inputs) {
                 for (std::size_t k : var_indices) {
                     auto& var = _variations.at(k);
@@ -737,15 +801,18 @@ void SystematicsCalculator::compute(
                         x.push_back(in.x1);
                         q.push_back(var.muf * scale_of(in, var.dyn, in.muf1));
                         slots.push_back(in.slot1);
+                        on_second.push_back(in.swapped);
                     }
                     if (has_beam2) {
                         x.push_back(in.x2);
                         q.push_back(var.muf * scale_of(in, var.dyn, in.muf2));
                         slots.push_back(in.slot2);
+                        on_second.push_back(!in.swapped);
                     }
                 }
             }
-            auto values = evaluate_pdf(evaluator, x, q, slots);
+            auto values =
+                evaluate_beams(evaluator, evaluator2, x, q, slots, on_second);
             for (std::size_t i = 0, n = 0; i < count; ++i) {
                 for (std::size_t k : var_indices) {
                     double product = 1.;
@@ -777,11 +844,11 @@ void SystematicsCalculator::compute(
             }
         }
         if (!nominal_grid_vars.empty()) {
-            evaluate_variations(_nominal_pdf.value(), nominal_grid_vars);
+            evaluate_variations(_nominal_pdf.value(), *nominal_pdf2, nominal_grid_vars);
         }
         for (std::size_t m = 0; m < _member_pdfs.size(); ++m) {
             if (!member_vars[m].empty()) {
-                evaluate_variations(_member_pdfs[m], member_vars[m]);
+                evaluate_variations(_member_pdfs[m], _member_pdfs[m], member_vars[m]);
             }
         }
     }
diff --git a/madspace/src/phasespace/integrand.cpp b/madspace/src/phasespace/integrand.cpp
index b2e74cb3dc..08f592ec9d 100644
--- a/madspace/src/phasespace/integrand.cpp
+++ b/madspace/src/phasespace/integrand.cpp
@@ -62,7 +62,8 @@ Integrand::Integrand(
     const std::vector<std::size_t>& flavor_diff_xs_indices,
     const std::vector<std::size_t>& flavor_subproc_indices,
     const std::vector<std::size_t>& flavor_per_subproc_remap,
-    std::size_t compressed_channel_weight_count
+    std::size_t compressed_channel_weight_count,
+    const std::optional<PdfGrid>& pdf_grid2
 ) :
     FunctionGenerator(
         "Integrand",
@@ -210,7 +211,21 @@ Integrand::Integrand(
                     std::distance(pids.begin(), pids.find(option.at(i)))
                 );
             }
-            _pdfs.at(i) = PartonDensity(pdf_grid.value(), {pids.begin(), pids.end()});
+            bool own_grid = i == 1 && pdf_grid2;
+            _pdfs.at(i) = PartonDensity(
+                own_grid ? pdf_grid2.value() : pdf_grid.value(),
+                {pids.begin(), pids.end()},
+                false,
+                own_grid ? pdf2_prefix : ""
+            );
+            if (pdf_grid2) {
+                _pdfs_swapped.at(i) = PartonDensity(
+                    own_grid ? pdf_grid.value() : pdf_grid2.value(),
+                    {pids.begin(), pids.end()},
+                    false,
+                    own_grid ? "" : pdf2_prefix
+                );
+            }
         }
     }
 
@@ -239,8 +254,20 @@ Integrand::Integrand(
     _has_mirror = false;
     for (bool mirror : flavor_mirror) {
         _flavor_mirror.push_back(mirror ? 2 : 1);
+        _flavor_mirror_factors.push_back(mirror ? 2. : 1.);
         _has_mirror |= mirror;
     }
+    // Mirroring an accepted lab-frame event swaps the beams only if they are
+    // identical: otherwise the mirrored event has another boost, other cuts
+    // and other PDFs, so its orientation must be chosen before the mapping.
+    if (_has_mirror && !mapping.mirror_beams() &&
+        (mapping.beam_rapidity() != 0. || pdf_grid2)) {
+        throw std::invalid_argument(
+            "mirrored flavors with asymmetric beams need a PhaseSpaceMapping "
+            "built with mirror_beams = true"
+        );
+    }
+    _mirror_before_cuts = _has_mirror && mapping.mirror_beams();
     _channel_part_ret_types = compute_channel_part_ret_types();
 }
 
@@ -301,7 +328,7 @@ NamedVector<Type> Integrand::compute_channel_part_ret_types() const {
     // outputs after cuts
     ret.push_back("indices_acc", Type(DataType::dt_int, acc_batch_size, {}));
     ret.push_back("momenta_acc", acc_four_vec_array(particle_count));
-    if (_has_mirror) {
+    if (_has_mirror && !_mirror_before_cuts) {
         if (!_madnis_training) {
             ret.push_back("momenta_mirror_acc", acc_four_vec_array(particle_count));
         }
@@ -427,6 +454,20 @@ NamedVector<Value> Integrand::build_channel_part(
         chan_index_in_group = fb.full({static_cast<me_int_t>(0), batch_size_val});
     }
 
+    // Orientation of the event, for a mapping with a mirror_index condition.
+    // Drawn uniformly; the flavor-dependent factor follows the flavor sampling.
+    Value mirror_index;
+    if (_mapping.mirror_beams()) {
+        if (_mirror_before_cuts) {
+            auto [index, mirror_det] =
+                fb.sample_discrete(mirror_random, static_cast<me_int_t>(2));
+            mirror_index = index;
+        } else {
+            mirror_index = fb.full({static_cast<me_int_t>(0), batch_size_val});
+        }
+        mapping_conditions.push_back(mirror_index);
+    }
+
     // Apply phase space mapping
     auto mapping_result = _mapping.build_forward(fb, {latent}, mapping_conditions);
     weights_before_cuts.push_back(mapping_result["det"]);
@@ -463,6 +504,18 @@ NamedVector<Value> Integrand::build_channel_part(
                         .value()
                         .build_function(fb, {x_acc.at(i), scales.at(i + 1)})
                         .at(0);
+                if (_mirror_before_cuts && _pdfs_swapped.at(i)) {
+                    // leg i of a mirrored event comes from the other beam
+                    auto pdf_swapped =
+                        _pdfs_swapped.at(i)
+                            .value()
+                            .build_function(fb, {x_acc.at(i), scales.at(i + 1)})
+                            .at(0);
+                    pdf = fb.gather_vector(
+                        fb.batch_gather(indices_acc, mirror_index),
+                        fb.stack({pdf, pdf_swapped})
+                    );
+                }
                 pdf_results.at(i) = pdf;
                 pdf_priors.push_back(fb.select(pdf, _pdf_indices.at(i)));
             }
@@ -544,8 +597,16 @@ NamedVector<Value> Integrand::build_channel_part(
         }
     }
 
+    if (_mirror_before_cuts) {
+        // Both orientations were sampled with probability 1/2. A flavor that
+        // stands for itself and its mirror image collects both (factor 2);
+        // for any other flavor the two orientations cover the same phase
+        // space (factor 1).
+        weights_after_cuts.push_back(fb.gather(flavor_id, _flavor_mirror_factors));
+    }
+
     Value momenta_mirror_acc, mirror_id_acc;
-    if (_has_mirror) {
+    if (_has_mirror && !_mirror_before_cuts) {
         Value option_count;
         if (std::all_of(
                 _flavor_mirror.begin(), _flavor_mirror.end(), [](me_int_t mirror) {
@@ -593,7 +654,7 @@ NamedVector<Value> Integrand::build_channel_part(
     // outputs after cuts
     out.push_back("indices_acc", indices_acc);
     out.push_back("momenta_acc", momenta_acc);
-    if (_has_mirror) {
+    if (_has_mirror && !_mirror_before_cuts) {
         if (!_madnis_training) {
             out.push_back("momenta_mirror_acc", momenta_mirror_acc);
         }
@@ -898,7 +959,7 @@ NamedVector<Value> Integrand::build_common_part(
         }
     } else {
         outputs.push_back("weight", optional_cut(weight));
-        if (_has_mirror) {
+        if (_has_mirror && !_mirror_before_cuts) {
             outputs.push_back(
                 "momenta",
                 scatter_or_drop(args.at("momenta"), args.at("momenta_mirror_acc"))
diff --git a/madspace/src/phasespace/phasespace.cpp b/madspace/src/phasespace/phasespace.cpp
index da86671962..709bc863bb 100644
--- a/madspace/src/phasespace/phasespace.cpp
+++ b/madspace/src/phasespace/phasespace.cpp
@@ -2,6 +2,7 @@
 #include "madspace/constants.hpp"
 #include "madspace/util.hpp"
 #include <algorithm>
+#include <cmath>
 
 using namespace madspace;
 
@@ -126,7 +127,9 @@ PhaseSpaceMapping::PhaseSpaceMapping(
     TChannelMode t_channel_mode,
     const std::optional<Cuts>& cuts,
     const std::vector<std::vector<std::size_t>>& permutations,
-    const std::optional<std::vector<std::size_t>>& color_order
+    const std::optional<std::vector<std::size_t>>& color_order,
+    double beam_rapidity,
+    bool mirror_beams
 ) :
     Mapping(
         "PhaseSpaceMapping",
@@ -154,9 +157,16 @@ PhaseSpaceMapping::PhaseSpaceMapping(
           )},
          {"x1", batch_float},
          {"x2", batch_float}},
-        permutations.size() > 1
-            ? NamedVector<Type>{{"permutation_index", batch_int}}
-            : NamedVector<Type>{}
+        [&] {
+            NamedVector<Type> conditions;
+            if (permutations.size() > 1) {
+                conditions.push_back("permutation_index", batch_int);
+            }
+            if (mirror_beams) {
+                conditions.push_back("mirror_index", batch_int);
+            }
+            return conditions;
+        }()
     ),
     _topology(topology),
     _cuts(cuts.value_or(Cuts(
@@ -166,6 +176,8 @@ PhaseSpaceMapping::PhaseSpaceMapping(
         std::pow(2 * PI, 4 - 3 * static_cast<int>(topology.outgoing_masses().size()))
     ),
     _sqrt_s_lab(cm_energy),
+    _beam_rapidity(beam_rapidity),
+    _mirror_beams(mirror_beams),
     _leptonic(leptonic),
     // A decay has no beams to sample momentum fractions for: the root
     // virtuality is fixed at the decaying particle's mass (passed as
@@ -175,7 +187,17 @@ PhaseSpaceMapping::PhaseSpaceMapping(
         (_topology.t_propagator_count() == 0 ||
          t_channel_mode != PhaseSpaceMapping::chili)
     ),
+    _two_to_one(topology.outgoing_masses().size() == 1),
     _t_mapping(std::monostate{}) {
+    if ((_beam_rapidity != 0. || _mirror_beams) && _topology.is_decay()) {
+        throw std::invalid_argument("a decay has no beams to give a rapidity");
+    }
+    if (_two_to_one && _topology.outgoing_masses().at(0) <= 0.) {
+        // s_hat = m^2 = 0 is outside the luminosity (tau = 0)
+        throw std::invalid_argument(
+            "a 2 -> 1 process needs a massive final-state particle"
+        );
+    }
     bool has_t_channel = _topology.t_propagator_count() > 0;
     struct DecayInfo {
         double m_min, pt_min, eta_max;
@@ -330,7 +352,9 @@ PhaseSpaceMapping::PhaseSpaceMapping(
         std::vector<double> eta_max, pt_min;
         for (std::size_t index : topology.decays().at(0).child_indices) {
             auto& info = decay_info.at(index);
-            eta_max.push_back(info.eta_max);
+            // eta cuts hold in the lab frame, the t-channel is generated in the
+            // beams' centre-of-mass frame: |y_cm| <= |eta_lab| + |beam_rapidity|
+            eta_max.push_back(info.eta_max + std::abs(beam_rapidity));
             pt_min.push_back(info.pt_min);
         }
         if (t_channel_mode == PhaseSpaceMapping::chili) {
@@ -378,8 +402,14 @@ PhaseSpaceMapping::PhaseSpaceMapping(
                     }
                     m_inv_co.at(a).at(b) =
                         m_inv_full.at(child_to_out.at(a)).at(child_to_out.at(b));
-                    dr_co.at(a).at(b) =
-                        dr_full.at(child_to_out.at(a)).at(child_to_out.at(b));
+                    // Delta R uses pseudorapidities, which a longitudinal boost
+                    // does not shift uniformly for massive particles: a lab
+                    // cut is no bound in the beams' frame. The cut itself is
+                    // still applied to the lab momenta.
+                    if (beam_rapidity == 0.) {
+                        dr_co.at(a).at(b) =
+                            dr_full.at(child_to_out.at(a)).at(child_to_out.at(b));
+                    }
                 }
             }
             _t_mapping = ColorOrderedMapping(
@@ -416,12 +446,26 @@ PhaseSpaceMapping::PhaseSpaceMapping(
     double invariant_power,
     TChannelMode mode,
     const std::optional<Cuts>& cuts,
-    const std::optional<std::vector<std::size_t>>& color_order
+    const std::optional<std::vector<std::size_t>>& color_order,
+    double beam_rapidity,
+    bool mirror_beams
 ) :
     PhaseSpaceMapping(
         Topology([&] {
-            if (external_masses.size() < 4) {
-                throw std::invalid_argument("The number of masses must be at least 4");
+            if (external_masses.size() < 3) {
+                throw std::invalid_argument("The number of masses must be at least 3");
+            }
+            if (external_masses.size() == 3) {
+                return Diagram(
+                    {external_masses.at(0), external_masses.at(1)},
+                    {external_masses.at(2)},
+                    {},
+                    {{
+                        {Diagram::incoming, 0},
+                        {Diagram::incoming, 1},
+                        {Diagram::outgoing, 0},
+                    }}
+                );
             }
             std::vector<Diagram::Vertex> vertices;
             auto n_out = external_masses.size() - 2;
@@ -455,7 +499,9 @@ PhaseSpaceMapping::PhaseSpaceMapping(
         mode,
         cuts,
         {},
-        color_order
+        color_order,
+        beam_rapidity,
+        mirror_beams
     ) {}
 
 Mapping::Result PhaseSpaceMapping::build_forward_impl(
@@ -479,6 +525,39 @@ Mapping::Result PhaseSpaceMapping::build_forward_impl(
     ValueVec dets{_pi_factors};
     Value x1 = 1.0, x2 = 1.0;
 
+    // permute, boost into the lab frame and apply the cuts
+    auto finish = [&](Value p_ext_stack) -> Result {
+        if (_permutations.size() > 1) {
+            p_ext_stack =
+                fb.permute_momenta(p_ext_stack, _permutations, conditions.at(0));
+        } else if (_permutations.size() == 1 &&
+                   !std::is_sorted(
+                       _permutations.at(0).begin(), _permutations.at(0).end()
+                   )) {
+            p_ext_stack = fb.permute_momenta(
+                p_ext_stack, _permutations, static_cast<me_int_t>(0)
+            );
+        }
+        auto p_ext_lab = to_lab(fb, p_ext_stack, x1, x2, conditions);
+        dets.push_back(_cuts.build_function(fb, {p_ext_lab}).at(0));
+        auto ps_weight = fb.cut_unphysical(fb.product(dets), p_ext_lab, x1, x2);
+        return {{{"momenta", p_ext_lab}, {"x1", x1}, {"x2", x2}}, ps_weight};
+    };
+
+    if (_two_to_one) {
+        // dPhi_1 = 2 pi delta(s_hat - m^2) (the 2 pi is _pi_factors): s_hat is
+        // fixed and the one random number left is the rapidity, through the
+        // split of tau = m^2 / s between x1 and x2.
+        double mass = _topology.outgoing_masses().at(0);
+        auto [x1_new, x2_new, det_x] =
+            fb.r_to_x1x2(next_random(), mass * mass, _sqrt_s_lab * _sqrt_s_lab);
+        x1 = x1_new;
+        x2 = x2_new;
+        dets.push_back(det_x);
+        auto [p1, p2] = fb.com_p_in(mass);
+        return finish(fb.stack({p1, p2, fb.add(p1, p2)}));
+    }
+
     // initialize masses and square masses
     std::vector<DecayData> decay_data(
         _topology.decays().begin(), _topology.decays().end()
@@ -621,24 +700,62 @@ Mapping::Result PhaseSpaceMapping::build_forward_impl(
     for (std::size_t decay_index : _topology.outgoing_indices()) {
         p_ext.push_back(decay_data.at(decay_index).momentum.value());
     }
-    auto p_ext_stack = fb.stack(p_ext);
+    return finish(fb.stack(p_ext));
+}
 
-    // permute momenta if permutations are given
-    if (_permutations.size() > 1) {
-        p_ext_stack = fb.permute_momenta(p_ext_stack, _permutations, conditions.at(0));
-    } else if (_permutations.size() == 1 &&
-               !std::is_sorted(
-                   _permutations.at(0).begin(), _permutations.at(0).end()
-               )) {
-        p_ext_stack =
-            fb.permute_momenta(p_ext_stack, _permutations, static_cast<me_int_t>(0));
+Value PhaseSpaceMapping::to_lab(
+    FunctionBuilder& fb,
+    Value momenta,
+    Value x1,
+    Value x2,
+    const NamedVector<Value>& conditions
+) const {
+    // boost_beam boosts by the rapidity 0.5 ln(x1 / x2). Without a mirror the
+    // partonic and beam boosts are one: scaling x1 by exp(y) and x2 by
+    // exp(-y) adds the beam rapidity y.
+    double exp_plus = std::exp(_beam_rapidity), exp_minus = std::exp(-_beam_rapidity);
+    if (_map_luminosity && !_mirror_beams && _beam_rapidity != 0.) {
+        return fb.boost_beam(momenta, fb.mul(x1, exp_plus), fb.mul(x2, exp_minus));
     }
+    if (_map_luminosity) {
+        momenta = fb.boost_beam(momenta, x1, x2);
+    }
+    if (_mirror_beams) {
+        momenta = fb.mirror_momenta(
+            momenta, conditions.at(_permutations.size() > 1 ? 1 : 0)
+        );
+    }
+    if (_beam_rapidity != 0.) {
+        momenta = fb.boost_beam(momenta, exp_plus, exp_minus);
+    }
+    return momenta;
+}
 
-    // boost into correct frame and apply cuts
-    auto p_ext_lab = _map_luminosity ? fb.boost_beam(p_ext_stack, x1, x2) : p_ext_stack;
-    dets.push_back(_cuts.build_function(fb, {p_ext_lab}).at(0));
-    auto ps_weight = fb.cut_unphysical(fb.product(dets), p_ext_lab, x1, x2);
-    return {{{"momenta", p_ext_lab}, {"x1", x1}, {"x2", x2}}, ps_weight};
+Value PhaseSpaceMapping::from_lab(
+    FunctionBuilder& fb,
+    Value momenta,
+    Value x1,
+    Value x2,
+    const NamedVector<Value>& conditions
+) const {
+    double exp_plus = std::exp(_beam_rapidity), exp_minus = std::exp(-_beam_rapidity);
+    if (_map_luminosity && !_mirror_beams && _beam_rapidity != 0.) {
+        return fb.boost_beam_inverse(
+            momenta, fb.mul(x1, exp_plus), fb.mul(x2, exp_minus)
+        );
+    }
+    if (_beam_rapidity != 0.) {
+        momenta = fb.boost_beam_inverse(momenta, exp_plus, exp_minus);
+    }
+    if (_mirror_beams) {
+        momenta = fb.mirror_momenta(
+            momenta, conditions.at(_permutations.size() > 1 ? 1 : 0)
+        );
+    }
+    if (_map_luminosity) {
+        momenta = fb.boost_beam_inverse(momenta, x1, x2);
+    }
+    return momenta;
 }
 
 Mapping::Result PhaseSpaceMapping::build_inverse_impl(
@@ -647,8 +764,11 @@ Mapping::Result PhaseSpaceMapping::build_inverse_impl(
     const NamedVector<Value>& conditions
 ) const {
     Value p_ext_lab = inputs.at(0), x1 = inputs.at(1), x2 = inputs.at(2);
-    Value p_ext_stack =
-        _map_luminosity ? fb.boost_beam_inverse(p_ext_lab, x1, x2) : p_ext_lab;
+    if (_two_to_one) {
+        auto [r, det_x] = fb.x1x2_to_r(x1, x2, _sqrt_s_lab * _sqrt_s_lab);
+        return {{{"random", fb.stack({r})}}, fb.mul(det_x, 1. / _pi_factors)};
+    }
+    Value p_ext_stack = from_lab(fb, p_ext_lab, x1, x2, conditions);
 
     // permute momenta if permutations are given
     if (_permutations.size() > 1) {
diff --git a/madspace/src/phasespace/topology.cpp b/madspace/src/phasespace/topology.cpp
index 25f875ee39..56303597ae 100644
--- a/madspace/src/phasespace/topology.cpp
+++ b/madspace/src/phasespace/topology.cpp
@@ -245,9 +245,15 @@ Diagram::Diagram(
             "Diagram must have one incoming particle (a decay) or two (a collision)"
         );
     }
-    if (outgoing_masses.size() < 2) {
+    // A collision may produce a single particle (p p > h): its one-particle
+    // phase space is a delta function in s_hat, which PhaseSpaceMapping
+    // integrates against the beam momentum fractions. A decay into one
+    // particle has nothing to integrate at all.
+    if (outgoing_masses.size() < (incoming_masses.size() == 2 ? 1 : 2)) {
         throw std::invalid_argument(
-            "Diagram must have at least two outgoing particles"
+            incoming_masses.size() == 2
+                ? "Diagram must have at least one outgoing particle"
+                : "A decay diagram must have at least two outgoing particles"
         );
     }
 
@@ -490,7 +496,9 @@ Topology::propagator_momentum_terms(bool only_decays) const {
     }
     for (auto& decay : std::views::reverse(_decays)) {
         if (decay.index == 0) {
-            if (_t_integration_order.size() == 0) {
+            // In a 2 -> 1 topology the root is the outgoing particle itself:
+            // its invariant is fixed at its mass, so it is no propagator.
+            if (_t_integration_order.size() == 0 && decay.child_indices.size() != 0) {
                 std::vector<int> factors(n_ext);
                 for (std::size_t i = 0; i < n_in; ++i) {
                     factors.at(i) = 1;
diff --git a/madspace/src/python/madspace.cpp b/madspace/src/python/madspace.cpp
index 244c949708..f95aea2430 100644
--- a/madspace/src/python/madspace.cpp
+++ b/madspace/src/python/madspace.cpp
@@ -719,7 +719,9 @@ PYBIND11_MODULE(_madspace_py, m) {
                 PhaseSpaceMapping::TChannelMode,
                 const std::optional<Cuts>&,
                 const nested_vector2<std::size_t>&,
-                const std::optional<std::vector<std::size_t>>&>(),
+                const std::optional<std::vector<std::size_t>>&,
+                double,
+                bool>(),
             py::arg("topology"),
             py::arg("cm_energy"),
             py::arg("leptonic") = false,
@@ -727,7 +729,9 @@ PYBIND11_MODULE(_madspace_py, m) {
             py::arg("t_channel_mode") = PhaseSpaceMapping::propagator,
             py::arg("cuts") = std::nullopt,
             py::arg("permutations") = std::vector<Topology>{},
-            py::arg("color_order") = std::nullopt
+            py::arg("color_order") = std::nullopt,
+            py::arg("beam_rapidity") = 0.,
+            py::arg("mirror_beams") = false
         )
         .def(
             py::init<
@@ -737,16 +741,22 @@ PYBIND11_MODULE(_madspace_py, m) {
                 double,
                 PhaseSpaceMapping::TChannelMode,
                 std::optional<Cuts>,
-                const std::optional<std::vector<std::size_t>>&>(),
+                const std::optional<std::vector<std::size_t>>&,
+                double,
+                bool>(),
             py::arg("masses"),
             py::arg("cm_energy"),
             py::arg("leptonic") = false,
             py::arg("invariant_power") = 0.8,
             py::arg("mode") = PhaseSpaceMapping::rambo,
             py::arg("cuts") = std::nullopt,
-            py::arg("color_order") = std::nullopt
+            py::arg("color_order") = std::nullopt,
+            py::arg("beam_rapidity") = 0.,
+            py::arg("mirror_beams") = false
         )
         .def("random_dim", &PhaseSpaceMapping::random_dim)
+        .def("beam_rapidity", &PhaseSpaceMapping::beam_rapidity)
+        .def("mirror_beams", &PhaseSpaceMapping::mirror_beams)
         .def("discrete_dim", &PhaseSpaceMapping::discrete_dim)
         .def("particle_count", &PhaseSpaceMapping::particle_count)
         .def("channel_count", &PhaseSpaceMapping::channel_count);
@@ -1266,7 +1276,8 @@ PYBIND11_MODULE(_madspace_py, m) {
                 const std::vector<std::size_t>&,
                 const std::vector<std::size_t>&,
                 const std::vector<std::size_t>&,
-                std::size_t>(),
+                std::size_t,
+                const std::optional<PdfGrid>&>(),
             py::arg("mapping"),
             py::arg("diff_xs"),
             py::arg("adaptive_map") = std::monostate{},
@@ -1294,7 +1305,8 @@ PYBIND11_MODULE(_madspace_py, m) {
             py::arg("flavor_diff_xs_indices") = std::vector<std::size_t>{},
             py::arg("flavor_subproc_indices") = std::vector<std::size_t>{},
             py::arg("flavor_per_subproc_remap") = std::vector<std::size_t>{},
-            py::arg("compressed_channel_weight_count") = 50
+            py::arg("compressed_channel_weight_count") = 50,
+            py::arg("pdf_grid2") = std::nullopt
         )
         .def("particle_count", &Integrand::particle_count)
         .def("madnis_training", &Integrand::madnis_training)
@@ -1915,14 +1927,16 @@ PYBIND11_MODULE(_madspace_py, m) {
                 const std::optional<AlphaSGrid>&,
                 ContextPtr,
                 const std::vector<std::optional<MatrixElement>>&,
-                const nested_vector2<me_int_t>&>(),
+                const nested_vector2<me_int_t>&,
+                const std::optional<PdfGrid>&>(),
             py::arg("config"),
             py::arg("subproc_args"),
             py::arg("nominal_pdf") = std::nullopt,
             py::arg("nominal_alpha_s") = std::nullopt,
             py::arg("context") = nullptr,
             py::arg("matrix_elements") = std::vector<std::optional<MatrixElement>>{},
-            py::arg("me_flavor_remap") = nested_vector2<me_int_t>{}
+            py::arg("me_flavor_remap") = nested_vector2<me_int_t>{},
+            py::arg("nominal_pdf2") = std::nullopt
         )
         .def_property_readonly("config", &SystematicsCalculator::config)
         .def_property_readonly(
diff --git a/madspace/tests/test_beam_frame.py b/madspace/tests/test_beam_frame.py
new file mode 100644
index 0000000000..8a82ecac54
--- /dev/null
+++ b/madspace/tests/test_beam_frame.py
@@ -0,0 +1,144 @@
+"""2 -> 1 phase space and asymmetric beams (PhaseSpaceMapping beam_rapidity)."""
+
+import math
+
+import numpy as np
+import pytest
+from pytest import approx
+
+import madspace as ms
+
+N = 2000
+
+
+@pytest.fixture
+def rng():
+    return np.random.default_rng(4321)
+
+
+def boost_z(p, y):
+    q = p.copy()
+    q[..., 0] = p[..., 0] * math.cosh(y) + p[..., 3] * math.sinh(y)
+    q[..., 3] = p[..., 3] * math.cosh(y) + p[..., 0] * math.sinh(y)
+    return q
+
+
+def two_to_one_topology(mass):
+    return ms.Topology(ms.Diagram([0.0, 0.0], [mass], [], [["i0", "i1", "o0"]]))
+
+
+def test_two_to_one(rng):
+    mass, sqrt_s = 91.188, 13000.0
+    mapping = ms.PhaseSpaceMapping(two_to_one_topology(mass), sqrt_s)
+    assert mapping.random_dim() == 1
+    r = rng.random((N, 1))
+    p, x1, x2, det = map(np.asarray, mapping.map_forward([r], []))
+
+    tau = mass**2 / sqrt_s**2
+    assert x1 * x2 == approx(tau)
+    assert p[:, 0] + p[:, 1] == approx(p[:, 2], abs=1e-9)
+    assert p[:, 0, 0] == approx(x1 * sqrt_s / 2)
+    assert p[:, 1, 0] == approx(x2 * sqrt_s / 2)
+    mass2 = p[:, 2, 0] ** 2 - np.sum(p[:, 2, 1:] ** 2, axis=1)
+    assert mass2 == approx(mass**2)
+    # dPhi_1 = 2 pi delta(s_hat - m^2), dx1 dx2 = |ln tau| / s ds_hat dr
+    assert det == approx(2 * math.pi * abs(math.log(tau)) / sqrt_s**2)
+
+    r_inv, det_inv = mapping.map_inverse([p, x1, x2], [])
+    assert np.asarray(r_inv) == approx(r)
+    assert np.asarray(det_inv) * det == approx(1.0)
+
+    flat = ms.PhaseSpaceMapping([0.0, 0.0, mass], sqrt_s)
+    assert np.asarray(flat.map_forward([r], [])[0]) == approx(p)
+
+
+def test_two_to_one_needs_pdf():
+    with pytest.raises(ValueError):
+        ms.PhaseSpaceMapping(two_to_one_topology(125.0), 13000.0, leptonic=True)
+    with pytest.raises(ValueError):
+        ms.PhaseSpaceMapping(two_to_one_topology(0.0), 13000.0)
+
+
+@pytest.mark.parametrize("leptonic", [False, True], ids=["hadronic", "leptonic"])
+@pytest.mark.parametrize("mode", ["rambo", "propagator"])
+def test_beam_rapidity(rng, mode, leptonic):
+    e1, e2 = 7000.0, 4000.0
+    e_cm = 2 * math.sqrt(e1 * e2)
+    y0 = 0.5 * math.log(e1 / e2)
+    masses = [0.0, 0.0, 0.0, 0.0, 173.0]
+    kwargs = dict(leptonic=leptonic, mode=getattr(ms.PhaseSpaceMapping, mode))
+    sym = ms.PhaseSpaceMapping(masses, e_cm, **kwargs)
+    lab = ms.PhaseSpaceMapping(masses, e_cm, beam_rapidity=y0, **kwargs)
+    r = rng.random((N, sym.random_dim()))
+    p_sym, x1_sym, x2_sym, det_sym = map(np.asarray, sym.map_forward([r], []))
+    p_lab, x1, x2, det = map(np.asarray, lab.map_forward([r], []))
+    x1, x2 = np.broadcast_to(x1, (N,)), np.broadcast_to(x2, (N,))
+    ok = det_sym > 0
+
+    # same point, same weight: the lab frame is a longitudinal boost away
+    assert det == approx(det_sym)
+    assert p_lab[ok] == approx(boost_z(p_sym[ok], y0), abs=1e-8)
+    # the incoming partons carry the momentum fractions of their own beam
+    assert p_lab[ok, 0, 0] == approx(x1[ok] * e1)
+    assert p_lab[ok, 1, 0] == approx(x2[ok] * e2)
+
+    r_inv, det_inv = lab.map_inverse([p_lab, x1, x2], [])
+    assert np.asarray(r_inv)[ok] == approx(r[ok], abs=1e-6)
+    assert np.asarray(det_inv)[ok] * det[ok] == approx(1.0, rel=1e-5)
+
+
+def test_beam_rapidity_cut_in_lab_frame(rng):
+    """An eta cut is applied to the boosted (lab) momenta."""
+    e1, e2 = 7000.0, 1000.0
+    e_cm = 2 * math.sqrt(e1 * e2)
+    y0 = 0.5 * math.log(e1 / e2)
+    pids = [2, -2, 11, -11]
+    cuts = ms.Cuts([
+        ms.CutItem(observable=ms.Observable(pids, "eta", [[11, -11]]), min=0.0)
+    ])
+    mapping = ms.PhaseSpaceMapping(
+        [0.0] * 4, e_cm, cuts=cuts, beam_rapidity=y0, mode=ms.PhaseSpaceMapping.rambo
+    )
+    free = ms.PhaseSpaceMapping(
+        [0.0] * 4, e_cm, beam_rapidity=y0, mode=ms.PhaseSpaceMapping.rambo
+    )
+    r = rng.random((N, mapping.random_dim()))
+    p, _, _, det = map(np.asarray, mapping.map_forward([r], []))
+    p_free, _, _, det_free = map(np.asarray, free.map_forward([r], []))
+    eta = np.arctanh(p_free[:, 2:, 3] / np.linalg.norm(p_free[:, 2:, 1:], axis=2))
+    passed = (eta > 0).all(axis=1) & (det_free > 0)
+    assert passed.any() and not passed.all()
+    assert ((det > 0) == passed).all()
+
+
+@pytest.mark.parametrize("pdf_asymmetric_only", [False, True], ids=["energies", "pdfs"])
+def test_mirror_beams(rng, pdf_asymmetric_only):
+    """mirror_index 1: leg 1 from beam 2, i.e. rotated in the beams' frame"""
+    e1, e2 = (6500.0, 6500.0) if pdf_asymmetric_only else (7000.0, 4000.0)
+    e_cm = 2 * math.sqrt(e1 * e2)
+    y0 = 0.5 * math.log(e1 / e2)
+    masses = [0.0, 0.0, 0.0, 0.0]
+    plain = ms.PhaseSpaceMapping(masses, e_cm, beam_rapidity=y0)
+    mirror = ms.PhaseSpaceMapping(masses, e_cm, beam_rapidity=y0, mirror_beams=True)
+    assert mirror.mirror_beams() and not plain.mirror_beams()
+    r = rng.random((N, plain.random_dim()))
+    p, x1, x2, det = map(np.asarray, plain.map_forward([r], []))
+    zeros, ones = np.zeros(N, dtype=np.int32), np.ones(N, dtype=np.int32)
+    p0, x1_0, x2_0, det0 = map(np.asarray, mirror.map_forward([r], [zeros]))
+    p1, x1_1, x2_1, det1 = map(np.asarray, mirror.map_forward([r], [ones]))
+    ok = det > 0
+    assert p0[ok] == approx(p[ok], abs=1e-8)
+    assert det1 == approx(det)
+    assert x1_1 == approx(x1)
+    assert x2_1 == approx(x2)
+    # leg 1 now moves along -z with its fraction of beam 2, leg 2 along +z
+    assert (p1[ok, 0, 3] < 0).all()
+    assert p1[ok, 0, 0] == approx(x1[ok] * e2)
+    assert (p1[ok, 1, 3] > 0).all()
+    assert p1[ok, 1, 0] == approx(x2[ok] * e1)
+    # which is the beams'-frame rotation of the unmirrored event
+    rot = boost_z(p[ok], -y0)
+    rot[..., 2:] *= -1
+    assert p1[ok] == approx(boost_z(rot, y0), abs=1e-8)
+    r_inv, det_inv = mirror.map_inverse([p1, x1_1, x2_1], [ones])
+    assert np.asarray(r_inv)[ok] == approx(r[ok], abs=1e-6)
diff --git a/tests/unit_tests/various/test_banner.py b/tests/unit_tests/various/test_banner.py
index b7a78e3d6a..7948a3aa57 100755
--- a/tests/unit_tests/various/test_banner.py
+++ b/tests/unit_tests/various/test_banner.py
@@ -24,6 +24,7 @@
 import os
 import models
 import io
+import re
 import sys
 from madgraph import MG5DIR
 
@@ -1487,6 +1488,49 @@ def test_collider_shortcuts(self):
         self.assertEqual(rc['beam']['e_cm'], 250.0)
         self.assertIs(rc['beam']['leptonic'], True)
 
+    def test_per_beam_energy_and_pdf(self):
+        """ebeam1/ebeam2 and pdf1/pdf2 are per beam; e_cm and pdf are derived
+        from them and setting either sets both beams"""
+        rc = bannermod.RunCardMG7()
+        rc.set('beam.ebeam1', '7 TeV', user=True)
+        rc.set('beam.ebeam2', 4000, user=True)
+        self.assertEqual(rc['beam']['ebeam1'], 7000.0)
+        self.assertAlmostEqual(rc['beam']['e_cm'], 2 * (7000.0 * 4000.0) ** 0.5)
+        self.assertEqual(rc['ebeam2'], 4000.0)  # legacy run_card.dat name
+        rc.set('e_cm', '14 TeV', user=True)
+        self.assertEqual((rc['beam']['ebeam1'], rc['beam']['ebeam2']), (7000.0, 7000.0))
+        rc.set('beam.pdf', 'CT18LO', user=True)
+        self.assertEqual((rc['beam']['pdf1'], rc['beam']['pdf2']), ('CT18LO', 'CT18LO'))
+        self.assertEqual(rc['beam']['pdf'], 'CT18LO')
+        rc.set('beam.pdf2', 'nNNPDF30_nlo_as_0118_A208_Z82', user=True)
+        self.assertRaises(KeyError, lambda: rc['beam']['pdf'])
+        # a bare "pdf" could be the beam or the systematics one: left untouched
+        rc.set('pdf', 'NNPDF31_lo_as_0118', user=True)
+        self.assertEqual(rc['beam']['pdf1'], 'CT18LO')
+        self.assertEqual(rc['systematics']['pdf'], ['errorset'])
+        with self.assertRaises(bannermod.InvalidRunCard):
+            rc.set('beam.e_cm', 'scan:[7000, 13000]', user=True)
+        rc.set('beam.ebeam2', 0, user=True)
+        self.assertRaises(bannermod.InvalidRunCard, rc.check_validity)
+
+    def test_read_card_with_e_cm_and_pdf(self):
+        """a card written before the per-beam keys still reads"""
+        rc = bannermod.RunCardMG7()
+        out = io.StringIO()
+        rc.write(out, template=self.template)
+        text = out.getvalue()
+        text = text.replace('ebeam1 = 6500.0\nebeam2 = 6500.0', 'e_cm = 1000.0')
+        text = re.sub(r'pdf1 = .*\npdf2 = .*', 'pdf = "NNPDF23_lo_as_0130_qed"', text)
+        self.assertNotIn('ebeam1 =', text)
+        self.assertNotIn('pdf1 =', text)
+        rc2 = bannermod.RunCardMG7(text)
+        self.assertEqual((rc2['beam']['ebeam1'], rc2['beam']['ebeam2']), (500.0, 500.0))
+        self.assertEqual(rc2['beam']['pdf2'], 'NNPDF23_lo_as_0130_qed')
+        out = io.StringIO()
+        rc2.write(out, template=self.template)
+        self.assertIn('ebeam2 = 500.0', out.getvalue())
+        self.assertNotIn('e_cm =', out.getvalue())
+
     def test_fixed_scale_and_remove_cuts(self):
         """fixed_scale sets all scales; remove_all_cut clears cuts"""
         rc = bannermod.RunCardMG7()
@@ -1506,7 +1550,9 @@ def test_run_card_scan_expansion(self):
         out = io.StringIO()
         rc.write(out, template=self.template)
         text = out.getvalue()
-        text = text.replace('e_cm = 13000.0', 'e_cm = "scan:[7000, 13000]"')
+        # the collision energy is scanned through both beams, moving together
+        text = text.replace('ebeam1 = 6500.0', 'ebeam1 = "scan2:[3500, 6500]"')
+        text = text.replace('ebeam2 = 6500.0', 'ebeam2 = "scan2:[3500, 6500]"')
         text = text.replace('ren_scale = 91.188', 'ren_scale = "scan1:[91.0, 172.0]"')
         text = text.replace('fact_scale1 = 91.188', 'fact_scale1 = "scan1:[45.5, 86.0]"')
         # valid TOML, and detected as a scan
@@ -1566,6 +1612,7 @@ def test_from_LO_conversion(self):
         mg7, dropped = bannermod.RunCardMG7.from_LO(lo, warn=False)
         # ported
         self.assertEqual(mg7['beam']['e_cm'], 13000.0)
+        self.assertEqual(mg7['beam']['ebeam1'], 6500.0)
         self.assertEqual(mg7['generation']['events'], 25000)
         self.assertEqual(mg7['beam']['dynamical_scale_choice'], 'half_transverse_mass')
         self.assertEqual(mg7['phasespace']['sde_strategy'], 'denominators')
@@ -1610,7 +1657,10 @@ def test_write_is_valid_toml(self):
         rc.write(out, template=self.template)
         data = tomllib.loads(out.getvalue())
         self.assertEqual(data['run']['output_format'], 'lhe')
-        self.assertEqual(data['beam']['e_cm'], 13000.0)
+        self.assertEqual(data['beam']['ebeam1'], 6500.0)
+        self.assertEqual(data['beam']['ebeam2'], 6500.0)
+        self.assertEqual(data['beam']['pdf1'], data['beam']['pdf2'])
+        self.assertNotIn('e_cm', data['beam'])
         self.assertIs(data['vegas']['enable'], True)
         self.assertEqual(data['multiparticles']['photon'], [22])
         self.assertEqual(data['cuts']['sqrt_s'], {'min': 0.0})

From 2a297a7941ba2d65144c7953dcd4ee10999c2698 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Tue, 15 Sep 2026 11:53:36 +0200
Subject: [PATCH 2/5] check_xsec_processes: pin the reference PDF on both beams

The run card now stores pdf1/pdf2, so the `pdf = "..."` substitution
matched nothing and the reference runs silently used the default
NNPDF4.0 set instead of NNPDF23_lo_as_0130_qed (3jets, dy1j, ttx1j and
ttx_decay off by 2-43%). Pin pdf1 and pdf2 (or a legacy pdf line) and fail
loudly when nothing was pinned.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 tests/acceptance_tests/test_check_xsec_processes_mg7.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tests/acceptance_tests/test_check_xsec_processes_mg7.py b/tests/acceptance_tests/test_check_xsec_processes_mg7.py
index 4c79c62686..eb8795df2a 100644
--- a/tests/acceptance_tests/test_check_xsec_processes_mg7.py
+++ b/tests/acceptance_tests/test_check_xsec_processes_mg7.py
@@ -148,7 +148,11 @@ def _edit_run_card(toml_path, events, disable_jet_cuts):
     t = open(toml_path).read()
     t = t.replace('fixed_ren_scale = false', 'fixed_ren_scale = true')
     t = t.replace('fixed_fact_scale = false', 'fixed_fact_scale = true')
-    t = re.sub(r'(?m)^pdf = ".*"$', 'pdf = "%s"' % _REFERENCE_PDF, t)
+    # one set per beam (pdf1/pdf2); a card from before the split has pdf
+    t, pinned = re.subn(r'(?m)^(pdf[12]?) = ".*"$',
+                        r'\1 = "%s"' % _REFERENCE_PDF, t)
+    if not pinned:
+        raise RuntimeError('no [beam] pdf entry to pin in %s' % toml_path)
     t = re.sub(r'events = \d+', 'events = %d' % events, t)
     if disable_jet_cuts:
         # jet cuts must be disabled for the hadronic tt~ decay processes to

From b44c4a78f015d55f4c0f84fd3f26b88d0465c332 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Thu, 17 Sep 2026 22:57:23 +0200
Subject: [PATCH 3/5] madspace: the mirrored orientation is the one the cuts
 and the ME see

The initial-state mirror covers the beam-swapped copy of a subprocess by
rotating the event by pi about x (py, pz -> -py, -pz, kernel_mirror_momenta),
which moves each leg onto the other beam. With mirror_beams the mapping draws
the orientation first, so the rotation happens before the boost into the lab
frame and before the cuts. Without it the Integrand keeps the cheaper path:
generate, cut, then mirror the accepted event. Everything downstream of that
mirror then sees an orientation nothing upstream ever looked at.

Two places were reading the wrong side of it.

The matrix element. Integrand::build_common_part kept feeding the cross
section momenta_acc, the momenta as generated, while the event was written
mirrored. It now takes the mirrored set when there is one, which is what
madevent does (it flips the momentum array that then goes to both the matrix
element and the event record, super_auto_dsig_group_v4.inc "Flip momenta
(rotate around x axis)") and what SystematicsCalculator::matrix_elements
already assumed, since it reads the momenta back out of the event buffer.
This is the same one-line change as "madspace: evaluate the matrix element on
the momenta the event is written with" on the mirror-me-momenta branch,
carried onto the post-cut path that only exists here; on the mirror_beams
path momenta_acc is already the written orientation and nothing changes.
Feeding the unmirrored momenta only ever worked because |M|^2 is invariant
under the rotation -- true for any Lorentz invariant matrix element, false
for a polarised one evaluated in a frame that holds the polarised particle at
rest, where HELAS quantises the spin along the frame z axis (48% on
g q > z{+} q in the Z rest frame, 23 sigma between p p > z{+} j and
p p > z{-} j end to end).

The cuts. Nothing stated that the cuts have to be blind to the mirror, and
nothing enforced it; it happened to hold because pt, |eta|, delta_r, mass and
sqrt_s all are. Observable::mirror_invariant classifies each observable (and
the one it is ordered by -- sorting by a quantity that flips picks a different
particle), Cuts::non_mirror_invariant_cuts names the offenders, and Integrand
refuses the post-cut mirror when any cut is one, next to the asymmetric-beams
check it already had. mg7 does not hit that refusal: launch.py now turns
mirror_beams on for a subprocess whose cuts are not mirror invariant, the same
switch asymmetric beams already trigger, and logs why.

This is not latent. The run_card [cuts] table takes any observable name, and
the template even advertises "eta" -- p p > e+ e- with lepton-eta.min = 0.0:

  post-cut mirror   276.3 +- 2.1 pb, 1017 of 2000 events (50.9%) written
                    with a lepton at eta < 0, i.e. outside the cut
  orientation first 233.9 +- 2.1 pb, no event outside the cut

and the second is right: lepton-eta.max = 0.0, which a pp collider's symmetry
forces to be equal, gives 232.9 +- 2.2 pb (0.3 sigma). The old number is 18%
high because it integrates one beam orientation twice, and a single ordering
of a valence-sea pair is forward-skewed.

Nothing on the invariant-cut path moves: default cuts keep the post-cut mirror
(no switch logged) and g _quark > e+ e- _quark gives 59.0349 +- 0.0490 pb
against the 59.0235 +- 0.0287 reference, 0.2 sigma.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 madgraph/iolibs/template_files/mg7/launch.py  |  33 ++-
 madspace/include/madspace/phasespace/cuts.hpp |   8 +
 .../madspace/phasespace/observable.hpp        |   6 +
 .../madspace/phasespace/phasespace.hpp        |   1 +
 madspace/src/phasespace/cuts.cpp              |  12 +
 madspace/src/phasespace/integrand.cpp         |  64 +++-
 madspace/src/phasespace/observable.cpp        |  39 +++
 madspace/src/python/madspace.cpp              |   7 +
 madspace/tests/test_beam_frame.py             | 278 ++++++++++++++++++
 9 files changed, 439 insertions(+), 9 deletions(-)

diff --git a/madgraph/iolibs/template_files/mg7/launch.py b/madgraph/iolibs/template_files/mg7/launch.py
index 1d2d872416..89be700dfb 100644
--- a/madgraph/iolibs/template_files/mg7/launch.py
+++ b/madgraph/iolibs/template_files/mg7/launch.py
@@ -2025,6 +2025,24 @@ def __init__(
             if len(self.process.cut_data) > 0
             else None
         )
+        # Mirroring an accepted event after the cuts hands the event writer an
+        # orientation the cuts never saw. That reproduces the same sample only
+        # if no cut can tell the two orientations apart -- true of pt, |eta|,
+        # delta_r and any invariant, false as soon as one is on a signed
+        # rapidity, eta, phi or pz. Draw the orientation before the mapping in
+        # that case, as asymmetric beams already do; ms.Integrand refuses the
+        # other way round.
+        self.mirror_beams = self.process.asymmetric_beams
+        if (self.cuts is not None and not self.process.is_decay
+                and not self.mirror_beams and self.has_mirrored_flavors()):
+            not_invariant = self.cuts.non_mirror_invariant_cuts()
+            if not_invariant:
+                logger.info(
+                    "subprocess %d: the cut(s) %s change under the initial-state "
+                    "mirror, so the beam orientation is drawn before the phase-"
+                    "space mapping rather than after the cuts",
+                    self.subproc_id, ", ".join(not_invariant))
+                self.mirror_beams = True
         self.histograms = (
             ms.ObservableHistograms([
                 ms.HistItem(
@@ -2088,7 +2106,7 @@ def build_multichannel_phasespace(self) -> PhaseSpace:
                     permutations=chan_permutations,
                     leptonic=self.process.leptonic,
                     beam_rapidity=self.process.beam_rapidity,
-                    mirror_beams=self.process.asymmetric_beams,
+                    mirror_beams=self.mirror_beams,
                 )
                 prefix = f"subproc{self.subproc_id}.channel{channel_id}"
                 if topo_count > 1:
@@ -2142,6 +2160,17 @@ def build_multichannel_phasespace(self) -> PhaseSpace:
             subchan_weights=subchan_weights,
         )
 
+    def has_mirrored_flavors(self) -> bool:
+        """Does any flavor of this subprocess stand for both beam orientations?
+        Resolved the same way build_integrands resolves flavor_mirror."""
+        for flav in self.meta["flavors"]:
+            if self.unmerged_meta is not None:
+                subproc_index = self.meta["subprocesses"][flav["subprocess"]]
+                flav = self.unmerged_meta[subproc_index]["flavors"][flav["flavor"]]
+            if flav["mirror"]:
+                return True
+        return False
+
     def build_flat_phasespace(self) -> PhaseSpace:
         mapping = ms.PhaseSpaceMapping(
             self.incoming_masses + self.outgoing_masses,
@@ -2150,7 +2179,7 @@ def build_flat_phasespace(self) -> PhaseSpace:
             cuts=self.cuts,
             leptonic=self.process.leptonic,
             beam_rapidity=self.process.beam_rapidity,
-            mirror_beams=self.process.asymmetric_beams,
+            mirror_beams=self.mirror_beams,
         )
         prefix = f"subproc{self.subproc_id}.flat"
         discrete_sym, discrete_flavor = self.build_discrete(
diff --git a/madspace/include/madspace/phasespace/cuts.hpp b/madspace/include/madspace/phasespace/cuts.hpp
index 7d8347fc75..72ba52d49a 100644
--- a/madspace/include/madspace/phasespace/cuts.hpp
+++ b/madspace/include/madspace/phasespace/cuts.hpp
@@ -5,6 +5,7 @@
 #include "madspace/phasespace/observable.hpp"
 
 #include <functional>
+#include <string>
 #include <utility>
 #include <vector>
 
@@ -22,6 +23,13 @@ class Cuts : public FunctionGenerator {
 
     Cuts(const std::vector<CutItem>& cut_data);
     Cuts(std::size_t particle_count);
+    // Names of the configured cuts that are not invariant under the
+    // initial-state mirror (py, pz -> -py, -pz); empty if every cut is.
+    // Mirroring an accepted event after the cuts only reproduces the mirrored
+    // half of the initial state if the cuts cannot tell the two orientations
+    // apart, since the event that gets written is the mirrored one.
+    std::vector<std::string> non_mirror_invariant_cuts() const;
+    bool mirror_invariant() const { return non_mirror_invariant_cuts().empty(); }
     double sqrt_s_min() const;
     std::vector<double> eta_max() const;
     std::vector<double> pt_min() const;
diff --git a/madspace/include/madspace/phasespace/observable.hpp b/madspace/include/madspace/phasespace/observable.hpp
index 96b291f65f..dbbe2eefb6 100644
--- a/madspace/include/madspace/phasespace/observable.hpp
+++ b/madspace/include/madspace/phasespace/observable.hpp
@@ -49,6 +49,12 @@ class Observable : public FunctionGenerator {
         const std::string& name = ""
     );
     ObservableOption observable() const { return _observable; }
+    // Unchanged by the initial-state mirror, the rotation by pi about x
+    // (py, pz -> -py, -pz) that moves each leg onto the other beam? Both the
+    // observable itself and the one it is ordered by have to be: sorting by a
+    // quantity that flips picks a different particle out of the event.
+    static bool mirror_invariant(ObservableOption observable);
+    bool mirror_invariant() const;
     const nested_vector2<me_int_t>& indices() const { return _indices; }
     bool sum_momenta() const { return _sum_momenta; }
     std::vector<std::size_t> simple_observable_indices() const {
diff --git a/madspace/include/madspace/phasespace/phasespace.hpp b/madspace/include/madspace/phasespace/phasespace.hpp
index 073a06fa82..ab9f8ce004 100644
--- a/madspace/include/madspace/phasespace/phasespace.hpp
+++ b/madspace/include/madspace/phasespace/phasespace.hpp
@@ -79,6 +79,7 @@ class PhaseSpaceMapping : public Mapping {
     std::size_t channel_count() const { return _permutations.size(); }
     double beam_rapidity() const { return _beam_rapidity; }
     bool mirror_beams() const { return _mirror_beams; }
+    const Cuts& cuts() const { return _cuts; }
 
 private:
     Result build_forward_impl(
diff --git a/madspace/src/phasespace/cuts.cpp b/madspace/src/phasespace/cuts.cpp
index 2c3fd03a84..222c352133 100644
--- a/madspace/src/phasespace/cuts.cpp
+++ b/madspace/src/phasespace/cuts.cpp
@@ -39,6 +39,18 @@ Cuts::build_function_impl(FunctionBuilder& fb, const NamedVector<Value>& args) c
     return {{"mask", fb.product(weights)}};
 }
 
+std::vector<std::string> Cuts::non_mirror_invariant_cuts() const {
+    std::vector<std::string> names;
+    for (auto& item : _cut_data) {
+        if (item.observable.mirror_invariant()) {
+            continue;
+        }
+        auto name = item.observable.name();
+        names.push_back(name.empty() ? "(unnamed cut)" : name);
+    }
+    return names;
+}
+
 double Cuts::sqrt_s_min() const {
     double sqrt_s_min = 0.;
     for (auto& item : _cut_data) {
diff --git a/madspace/src/phasespace/integrand.cpp b/madspace/src/phasespace/integrand.cpp
index 08f592ec9d..7651479bd4 100644
--- a/madspace/src/phasespace/integrand.cpp
+++ b/madspace/src/phasespace/integrand.cpp
@@ -2,6 +2,7 @@
 
 #include "madspace/util.hpp"
 
+#include <format>
 #include <set>
 
 using namespace madspace;
@@ -260,12 +261,35 @@ Integrand::Integrand(
     // Mirroring an accepted lab-frame event swaps the beams only if they are
     // identical: otherwise the mirrored event has another boost, other cuts
     // and other PDFs, so its orientation must be chosen before the mapping.
-    if (_has_mirror && !mapping.mirror_beams() &&
-        (mapping.beam_rapidity() != 0. || pdf_grid2)) {
-        throw std::invalid_argument(
-            "mirrored flavors with asymmetric beams need a PhaseSpaceMapping "
-            "built with mirror_beams = true"
-        );
+    if (_has_mirror && !mapping.mirror_beams()) {
+        if (mapping.beam_rapidity() != 0. || pdf_grid2) {
+            throw std::invalid_argument(
+                "mirrored flavors with asymmetric beams need a PhaseSpaceMapping "
+                "built with mirror_beams = true"
+            );
+        }
+        // Even with identical beams, mirroring after the cuts hands the event
+        // writer an orientation the cuts never saw. That is only the same
+        // sample if no cut can tell the two orientations apart -- which is
+        // true of every cut on an invariant, a pt or an |eta|, and false as
+        // soon as one is on a signed rapidity, eta, phi or pz.
+        auto bad_cuts = mapping.cuts().non_mirror_invariant_cuts();
+        if (!bad_cuts.empty()) {
+            std::string names;
+            for (auto& name : bad_cuts) {
+                if (!names.empty()) {
+                    names += ", ";
+                }
+                names += name;
+            }
+            throw std::invalid_argument(std::format(
+                "the cut(s) {} are not invariant under the initial-state mirror "
+                "(py, pz -> -py, -pz), so they have to be applied to the mirrored "
+                "event: mirrored flavors then need a PhaseSpaceMapping built with "
+                "mirror_beams = true, which draws the orientation before the cuts",
+                names
+            ));
+        }
     }
     _mirror_before_cuts = _has_mirror && mapping.mirror_beams();
     _channel_part_ret_types = compute_channel_part_ret_types();
@@ -698,6 +722,32 @@ NamedVector<Value> Integrand::build_common_part(
     auto flavor_id = args.at("flavor_id");
     auto batch_size_val = fb.batch_size({args.at("weight_before_cuts")});
 
+    // The matrix element is evaluated on the momenta the event is written
+    // with, which for the mirrored half of a beam-swapped subprocess is the
+    // mirrored set (mirror_momenta: py, pz -> -py, -pz, the rotation by pi
+    // about x that moves each leg onto the other beam). This is what madevent
+    // does too -- it flips the momentum array that then goes to both the matrix
+    // element and the event record (super_auto_dsig_group_v4.inc, "Flip momenta
+    // (rotate around x axis)") -- and it is what the systematics reweighting
+    // here already assumes, since that reads the momenta back out of the event
+    // buffer.
+    //
+    // Feeding the unmirrored momenta instead only ever worked because |M|^2 is
+    // invariant under that rotation. A polarised matrix element is not, once it
+    // is evaluated in a frame that holds the polarised particle at rest: HELAS
+    // quantises such a particle along the *frame* z axis (the pp == 0 branch of
+    // vxxxxx) rather than along its own momentum, the mirror flips that axis,
+    // and the + and - states swap. Measured at 48% on g q > z{+} q evaluated in
+    // the Z rest frame.
+    //
+    // momenta_mirror_acc is only produced on the post-cut path. With
+    // mirror_beams the mirror is inside the mapping, before the boost into the
+    // lab frame and before the cuts, so momenta_acc is already the orientation
+    // the event is written with and there is nothing to pick.
+    auto momenta_me = args.index_map().contains("momenta_mirror_acc")
+        ? args.at("momenta_mirror_acc")
+        : momenta_acc;
+
     auto scatter_or_drop = [&](Value default_value, Value value) -> Value {
         if (_drop_cuts_and_rescale) {
             return value;
@@ -741,7 +791,7 @@ NamedVector<Value> Integrand::build_common_part(
 
     // Evaluate differential cross section
     ValueVec xs_args{
-        momenta_acc,
+        momenta_me,
         _flavor_remap.size() > 0 ? fb.gather_int(flavor_id, _flavor_remap) : flavor_id,
     };
     xs_args.push_back(x1_acc);
diff --git a/madspace/src/phasespace/observable.cpp b/madspace/src/phasespace/observable.cpp
index a3c777b38c..1f031fb41a 100644
--- a/madspace/src/phasespace/observable.cpp
+++ b/madspace/src/phasespace/observable.cpp
@@ -229,6 +229,45 @@ std::tuple<nested_vector2<me_int_t>, nested_vector2<me_int_t>, Type> build_indic
 
 } // namespace
 
+bool Observable::mirror_invariant(ObservableOption observable) {
+    // The initial-state mirror is the rotation by pi about the x axis
+    // (E, px, py, pz) -> (E, px, -py, -pz). Every case is listed so that a new
+    // observable has to declare which side it falls on.
+    switch (observable) {
+    case obs_py:        // -> -py
+    case obs_pz:        // -> -pz
+    case obs_phi:       // atan2(py, px) -> -phi
+    case obs_theta:     // -> pi - theta
+    case obs_y:         // -> -y
+    case obs_eta:       // -> -eta
+    case obs_delta_eta: // signed difference of two flipped etas -> -delta_eta
+    case obs_delta_phi: // signed difference of two flipped phis -> -delta_phi
+        return false;
+    case obs_e:
+    case obs_px:
+    case obs_mass:
+    case obs_pt:
+    case obs_p_mag:
+    case obs_y_abs:
+    case obs_eta_abs:
+    case obs_delta_r: // sqrt(delta_eta^2 + delta_phi^2): both flip together
+    case obs_pair_mass:
+    case obs_sqrt_s:
+        return true;
+    }
+}
+
+bool Observable::mirror_invariant() const {
+    // An observable that matched no particle is the constant 0, whatever the
+    // orientation: a cut naming a particle this process does not have imposes
+    // nothing, and must not stand in the way of the post-cut mirror.
+    if (not_found()) {
+        return true;
+    }
+    return mirror_invariant(_observable) &&
+        (!_order_observable || mirror_invariant(_order_observable.value()));
+}
+
 const std::vector<int> Observable::jet_pids{1, 2, 3, 4, -1, -2, -3, -4, 21};
 const std::vector<int> Observable::bottom_pids{-5, 5};
 const std::vector<int> Observable::lepton_pids{11, 13, 15, -11, -13, -15};
diff --git a/madspace/src/python/madspace.cpp b/madspace/src/python/madspace.cpp
index f95aea2430..c7355a33cf 100644
--- a/madspace/src/python/madspace.cpp
+++ b/madspace/src/python/madspace.cpp
@@ -585,6 +585,10 @@ PYBIND11_MODULE(_madspace_py, m) {
            py::arg("ignore_incoming") = true,
            py::arg("name") = ""
     )
+        .def(
+            "mirror_invariant",
+            py::overload_cast<>(&Observable::mirror_invariant, py::const_)
+        )
         .def_readonly_static("jet_pids", &Observable::jet_pids)
         .def_readonly_static("bottom_pids", &Observable::bottom_pids)
         .def_readonly_static("lepton_pids", &Observable::lepton_pids)
@@ -614,6 +618,8 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def_readonly("mode", &Cuts::CutItem::mode);
     cuts.def(py::init<const std::vector<Cuts::CutItem>&>(), py::arg("cut_data"))
         .def(py::init<std::size_t>(), py::arg("particle_count"))
+        .def("non_mirror_invariant_cuts", &Cuts::non_mirror_invariant_cuts)
+        .def("mirror_invariant", &Cuts::mirror_invariant)
         .def("sqrt_s_min", &Cuts::sqrt_s_min)
         .def("eta_max", &Cuts::eta_max)
         .def("pt_min", &Cuts::pt_min)
@@ -757,6 +763,7 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def("random_dim", &PhaseSpaceMapping::random_dim)
         .def("beam_rapidity", &PhaseSpaceMapping::beam_rapidity)
         .def("mirror_beams", &PhaseSpaceMapping::mirror_beams)
+        .def("cuts", &PhaseSpaceMapping::cuts)
         .def("discrete_dim", &PhaseSpaceMapping::discrete_dim)
         .def("particle_count", &PhaseSpaceMapping::particle_count)
         .def("channel_count", &PhaseSpaceMapping::channel_count);
diff --git a/madspace/tests/test_beam_frame.py b/madspace/tests/test_beam_frame.py
index 8a82ecac54..53b5465668 100644
--- a/madspace/tests/test_beam_frame.py
+++ b/madspace/tests/test_beam_frame.py
@@ -142,3 +142,281 @@ def test_mirror_beams(rng, pdf_asymmetric_only):
     assert p1[ok] == approx(boost_z(rot, y0), abs=1e-8)
     r_inv, det_inv = mirror.map_inverse([p1, x1_1, x2_1], [ones])
     assert np.asarray(r_inv)[ok] == approx(r[ok], abs=1e-6)
+
+
+# --- the mirror and the cuts -------------------------------------------------
+#
+# The initial-state mirror is the rotation by pi about x, (E, px, py, pz) ->
+# (E, px, -py, -pz), which moves each leg onto the other beam. With mirror_beams
+# the mapping draws the orientation first, so the cuts and the written event see
+# the same momenta. Without it the Integrand mirrors accepted events afterwards,
+# which only reproduces the same sample if no cut can tell the two orientations
+# apart; Integrand refuses the combination otherwise.
+
+
+def mirror(p):
+    """Rotate by pi about the x axis, as kernel_mirror_momenta does."""
+    q = p.copy()
+    q[..., 2:] *= -1
+    return q
+
+
+def eta(p):
+    return np.arctanh(p[..., 3] / np.linalg.norm(p[..., 1:], axis=-1))
+
+
+PIDS = [2, -2, 11, -11]
+
+
+def lepton_cut(observable, **bounds):
+    return ms.Cuts([
+        ms.CutItem(
+            observable=ms.Observable(
+                PIDS, observable, [[11, -11]], name=f"lepton-{observable}"
+            ),
+            **bounds,
+        )
+    ])
+
+
+@pytest.mark.parametrize(
+    "observable, invariant",
+    [
+        ("e", True),
+        ("px", True),
+        ("py", False),
+        ("pz", False),
+        ("mass", True),
+        ("pt", True),
+        ("p_mag", True),
+        ("phi", False),
+        ("theta", False),
+        ("y", False),
+        ("y_abs", True),
+        ("eta", False),
+        ("eta_abs", True),
+    ],
+)
+def test_observable_mirror_invariance(observable, invariant):
+    obs = ms.Observable(PIDS, observable, [[11, -11]])
+    assert obs.mirror_invariant() == invariant
+
+
+@pytest.mark.parametrize(
+    "observable, invariant",
+    [("delta_eta", False), ("delta_phi", False), ("delta_r", True),
+     ("pair_mass", True)],
+)
+def test_pairwise_observable_mirror_invariance(observable, invariant):
+    obs = ms.Observable(PIDS, observable, [[11, -11], [11, -11]])
+    assert obs.mirror_invariant() == invariant
+
+
+def test_sqrt_s_is_mirror_invariant():
+    assert ms.Observable(PIDS, "sqrt_s", []).mirror_invariant()
+
+
+def test_ordering_observable_decides_mirror_invariance():
+    """Sorting by a quantity that flips picks another particle out of the event."""
+    kwargs = dict(select_pids=[[11, -11]], order_indices=[1])
+    by_pt = ms.Observable(PIDS, "pt", order_observable="pt", **kwargs)
+    by_eta = ms.Observable(PIDS, "pt", order_observable="eta", **kwargs)
+    assert by_pt.mirror_invariant()
+    assert not by_eta.mirror_invariant()
+
+
+def test_cut_on_absent_particle_is_mirror_invariant():
+    """It matched nothing, so it is the constant 0 whatever the orientation."""
+    obs = ms.Observable(PIDS, "eta", [[5, -5]])
+    assert obs.mirror_invariant()
+    assert ms.Cuts([ms.CutItem(observable=obs, min=0.0)]).mirror_invariant()
+
+
+def test_cuts_report_the_offending_cuts_by_name():
+    cuts = ms.Cuts([
+        ms.CutItem(
+            observable=ms.Observable(PIDS, "pt", [[11, -11]], name="lepton-pt"),
+            min=10.0,
+        ),
+        ms.CutItem(
+            observable=ms.Observable(PIDS, "eta", [[11, -11]], name="lepton-eta"),
+            min=0.0,
+        ),
+        ms.CutItem(
+            observable=ms.Observable(PIDS, "y", [[11, -11]], name="lepton-y"),
+            min=0.0,
+        ),
+    ])
+    assert not cuts.mirror_invariant()
+    assert cuts.non_mirror_invariant_cuts() == ["lepton-eta", "lepton-y"]
+    mapping = ms.PhaseSpaceMapping([0.0] * 4, 13000.0, cuts=cuts)
+    assert mapping.cuts().non_mirror_invariant_cuts() == ["lepton-eta", "lepton-y"]
+
+
+@pytest.mark.parametrize(
+    "observable, bounds, separates",
+    [("eta", dict(min=0.0), True), ("eta_abs", dict(max=1.0), False)],
+    ids=["eta", "eta_abs"],
+)
+def test_mirror_beams_cuts_the_written_orientation(rng, observable, bounds, separates):
+    """With mirror_beams the cuts act on the momenta the mapping returns.
+
+    A signed eta cut tells the two orientations apart -- the configuration the
+    post-cut mirror may not be used for -- while |eta| cannot, which is why the
+    post-cut path is sound for the cuts mg7 ships.
+    """
+    e_cm, mode = 13000.0, ms.PhaseSpaceMapping.rambo
+    cut = ms.PhaseSpaceMapping(
+        [0.0] * 4, e_cm, cuts=lepton_cut(observable, **bounds),
+        mirror_beams=True, mode=mode
+    )
+    free = ms.PhaseSpaceMapping([0.0] * 4, e_cm, mirror_beams=True, mode=mode)
+    r = rng.random((N, cut.random_dim()))
+    accepted, momenta, physical = {}, {}, {}
+    for index in (0, 1):
+        condition = np.full(N, index, dtype=np.int32)
+        accepted[index] = np.asarray(cut.map_forward([r], [condition])[3]) > 0
+        p, _, _, det = map(np.asarray, free.map_forward([r], [condition]))
+        momenta[index], physical[index] = p, det > 0
+
+    ok = physical[0]
+    assert ok.any()
+    # orientation 1 is the pi rotation of orientation 0
+    assert momenta[1][ok] == approx(mirror(momenta[0][ok]), abs=1e-7)
+    # and both are cut on the momenta that come out, not on the other orientation
+    for index in (0, 1):
+        obs = eta(momenta[index][:, 2:])
+        if observable == "eta_abs":
+            passed = (np.abs(obs) < bounds["max"]).all(axis=1)
+        else:
+            passed = (obs > bounds["min"]).all(axis=1)
+        assert (accepted[index][ok] == (passed & physical[index])[ok]).all()
+    # which matters only if the cut can tell the two orientations apart
+    differ = (accepted[0][ok] != accepted[1][ok]).any()
+    assert differ == separates
+
+
+# --- the mirror and the matrix element ---------------------------------------
+
+
+@pytest.fixture(scope="module")
+def running_coupling(tmp_path_factory):
+    """A RunningCoupling over a hand-written alpha_s grid: the Integrand needs
+    one to build, and nothing here depends on the values."""
+    info = tmp_path_factory.mktemp("alphas") / "grid.info"
+    info.write_text(
+        "AlphaS_Qs: [1.0, 10.0, 100.0, 1000.0, 10000.0]\n"
+        "AlphaS_Vals: [0.4, 0.2, 0.12, 0.09, 0.08]\n"
+    )
+    return ms.RunningCoupling(ms.AlphaSGrid(str(info)))
+
+
+PID_OPTIONS = [[2, -2], [-2, 2]]
+
+
+def build_integrand(running_coupling, cuts=None, mirror_beams=False,
+                    flavor_mirror=(True, True)):
+    """A 2 -> 2 integrand over a beam-swapped pair of flavors. The matrix
+    element is never called; only the compute graph is built."""
+    e_cm = 13000.0
+    me = ms.MatrixElement(
+        matrix_element_index=0,
+        particle_count=4,
+        inputs=[
+            ms.MatrixElement.momenta_in,
+            ms.MatrixElement.flavor_in,
+            ms.MatrixElement.alpha_s_in,
+        ],
+        outputs=[
+            ms.MatrixElement.matrix_element_out,
+            ms.MatrixElement.diagram_amp2_out,
+            ms.MatrixElement.color_index_out,
+            ms.MatrixElement.helicity_index_out,
+            ms.MatrixElement.diagram_index_out,
+        ],
+        diagram_count=2,
+        sample_random_inputs=True,
+    )
+    diff_xs = ms.DifferentialCrossSection(
+        matrix_element=me,
+        cm_energy=e_cm,
+        running_coupling=running_coupling,
+        energy_scale=ms.CachedScale(),
+        pid_options=PID_OPTIONS,
+    )
+    mapping = ms.PhaseSpaceMapping(
+        [0.0] * 4, e_cm, cuts=cuts, mirror_beams=mirror_beams
+    )
+    return ms.Integrand(
+        mapping=mapping,
+        diff_xs=[diff_xs],
+        pid_options=PID_OPTIONS,
+        running_coupling=running_coupling,
+        energy_scale=ms.EnergyScale(4),
+        flavor_mirror=list(flavor_mirror),
+        channel_indices=[0],
+    )
+
+
+def test_post_cut_mirror_accepts_mirror_invariant_cuts(running_coupling):
+    build_integrand(running_coupling, cuts=lepton_cut("eta_abs", max=2.5))
+
+
+def test_post_cut_mirror_refuses_a_cut_that_is_not_mirror_invariant(running_coupling):
+    """The event is written mirrored, so a cut the mirror changes would be
+    applied to an orientation nobody keeps."""
+    with pytest.raises(ValueError, match="lepton-eta"):
+        build_integrand(running_coupling, cuts=lepton_cut("eta", min=0.0))
+
+
+def test_mirror_beams_allows_any_cut(running_coupling):
+    """Drawing the orientation before the mapping puts the cuts on the written
+    event, so the restriction does not apply."""
+    build_integrand(
+        running_coupling, cuts=lepton_cut("eta", min=0.0), mirror_beams=True
+    )
+
+
+def test_unmirrored_flavors_allow_any_cut(running_coupling):
+    build_integrand(
+        running_coupling, cuts=lepton_cut("eta", min=0.0), flavor_mirror=(False, False)
+    )
+
+
+def instruction_names(function):
+    return [instruction.instruction.name for instruction in function.instructions]
+
+
+def find_instruction(function, name):
+    names = instruction_names(function)
+    assert names.count(name) == 1, f"expected one {name}, got {names.count(name)}"
+    return function.instructions[names.index(name)]
+
+
+def test_matrix_element_sees_the_written_momenta(running_coupling):
+    """On the post-cut path the matrix element is evaluated on the mirrored
+    momenta, the ones the event is written with -- not on the momenta as
+    generated. |M|^2 is invariant under the rotation for a Lorentz invariant
+    matrix element, but not for a polarised one evaluated in a frame that holds
+    the polarised particle at rest."""
+    function = build_integrand(running_coupling).function()
+    names = instruction_names(function)
+    # the orientation is drawn after the cuts here
+    assert names.index("mirror_momenta") > names.index("nonzero")
+    mirrored = find_instruction(function, "mirror_momenta").outputs[0]
+    matrix_element = find_instruction(function, "matrix_element")
+    assert any(
+        str(value) == str(mirrored) for value in matrix_element.inputs
+    ), "the matrix element is not evaluated on the mirrored momenta"
+
+
+def test_mirror_beams_needs_no_second_mirror(running_coupling):
+    """With mirror_beams the mapping already returns the written orientation,
+    so there is nothing left for the Integrand to pick between."""
+    function = build_integrand(running_coupling, mirror_beams=True).function()
+    names = instruction_names(function)
+    # the only mirror is the one inside the mapping, before the cuts
+    assert names.index("mirror_momenta") < names.index("nonzero")
+    mirrored = find_instruction(function, "mirror_momenta").outputs[0]
+    matrix_element = find_instruction(function, "matrix_element")
+    assert not any(str(value) == str(mirrored) for value in matrix_element.inputs)

From cf8583fff84062c512b5e89b687b35d089a0db48 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Fri, 18 Sep 2026 11:05:10 +0200
Subject: [PATCH 4/5] madspace: tidy the mirror-invariance API after the main
 merge

Follow-up to the merge of origin/main (12e2e1602):

- Observable::mirror_invariant(ObservableOption), the static per-option
  overload, becomes option_mirror_invariant() in observable.cpp's anonymous
  namespace. Nothing outside that file called it, and as an overload with its
  own doc comment it gave generate_docstrings two entries for one name
  (Observable::mirror_invariant / #2), which is why the Python method was
  bound without a docstring. It now gets pydoc::doc like every other binding.
- clang-format on the Integrand error message and black/isort on
  test_beam_frame.py, which the madspace pre-commit hooks run.

pydoc keys checked against docstrings generated with doxygen (ENABLE_DOCS=ON
only runs on pushes to main): all 613 doc() calls resolve. madspace suite
1605 passed, 11 skipped, 4 failed (scipy not installed, as before the merge).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .../madspace/phasespace/observable.hpp        | 11 ++-
 madspace/src/phasespace/integrand.cpp         | 16 ++--
 madspace/src/phasespace/observable.cpp        | 60 +++++++--------
 madspace/src/python/madspace.cpp              |  3 +-
 madspace/tests/test_beam_frame.py             | 74 +++++++++++--------
 5 files changed, 89 insertions(+), 75 deletions(-)

diff --git a/madspace/include/madspace/phasespace/observable.hpp b/madspace/include/madspace/phasespace/observable.hpp
index e5c4b26c6d..d95b81b0e8 100644
--- a/madspace/include/madspace/phasespace/observable.hpp
+++ b/madspace/include/madspace/phasespace/observable.hpp
@@ -100,13 +100,12 @@ class Observable : public FunctionGenerator {
     );
     /// The kinematic quantity this observable computes.
     ObservableOption observable() const { return _observable; }
-    /// Whether @p observable is unchanged by the initial-state mirror, the
+    /// Whether this observable is unchanged by the initial-state mirror, the
     /// rotation by pi about x (py, pz -> -py, -pz) that moves each leg onto
-    /// the other beam.
-    static bool mirror_invariant(ObservableOption observable);
-    /// Whether this observable is mirror invariant. Both the observable itself
-    /// and the one it is ordered by have to be: sorting by a quantity that
-    /// flips picks a different particle out of the event.
+    /// the other beam. Both the observable itself and the one it is ordered by
+    /// have to be: sorting by a quantity that flips picks a different particle
+    /// out of the event. An observable that matched no particle (see
+    /// not_found()) is the constant 0 and counts as invariant.
     bool mirror_invariant() const;
     /// Per-selection lists of particle indices the observable is evaluated on.
     const nested_vector2<me_int_t>& indices() const { return _indices; }
diff --git a/madspace/src/phasespace/integrand.cpp b/madspace/src/phasespace/integrand.cpp
index 7651479bd4..a97e8493e8 100644
--- a/madspace/src/phasespace/integrand.cpp
+++ b/madspace/src/phasespace/integrand.cpp
@@ -282,13 +282,15 @@ Integrand::Integrand(
                 }
                 names += name;
             }
-            throw std::invalid_argument(std::format(
-                "the cut(s) {} are not invariant under the initial-state mirror "
-                "(py, pz -> -py, -pz), so they have to be applied to the mirrored "
-                "event: mirrored flavors then need a PhaseSpaceMapping built with "
-                "mirror_beams = true, which draws the orientation before the cuts",
-                names
-            ));
+            throw std::invalid_argument(
+                std::format(
+                    "the cut(s) {} are not invariant under the initial-state mirror "
+                    "(py, pz -> -py, -pz), so they have to be applied to the mirrored "
+                    "event: mirrored flavors then need a PhaseSpaceMapping built with "
+                    "mirror_beams = true, which draws the orientation before the cuts",
+                    names
+                )
+            );
         }
     }
     _mirror_before_cuts = _has_mirror && mapping.mirror_beams();
diff --git a/madspace/src/phasespace/observable.cpp b/madspace/src/phasespace/observable.cpp
index 1f031fb41a..deb800788b 100644
--- a/madspace/src/phasespace/observable.cpp
+++ b/madspace/src/phasespace/observable.cpp
@@ -33,6 +33,34 @@ int observable_type(Observable::ObservableOption observable) {
     }
 }
 
+bool option_mirror_invariant(Observable::ObservableOption observable) {
+    // The initial-state mirror is the rotation by pi about the x axis
+    // (E, px, py, pz) -> (E, px, -py, -pz). Every case is listed so that a new
+    // observable has to declare which side it falls on.
+    switch (observable) {
+    case Observable::obs_py:        // -> -py
+    case Observable::obs_pz:        // -> -pz
+    case Observable::obs_phi:       // atan2(py, px) -> -phi
+    case Observable::obs_theta:     // -> pi - theta
+    case Observable::obs_y:         // -> -y
+    case Observable::obs_eta:       // -> -eta
+    case Observable::obs_delta_eta: // difference of two flipped etas -> -delta_eta
+    case Observable::obs_delta_phi: // difference of two flipped phis -> -delta_phi
+        return false;
+    case Observable::obs_e:
+    case Observable::obs_px:
+    case Observable::obs_mass:
+    case Observable::obs_pt:
+    case Observable::obs_p_mag:
+    case Observable::obs_y_abs:
+    case Observable::obs_eta_abs:
+    case Observable::obs_delta_r: // sqrt(delta_eta^2 + delta_phi^2): both flip
+    case Observable::obs_pair_mass:
+    case Observable::obs_sqrt_s:
+        return true;
+    }
+}
+
 Value build_observable(
     FunctionBuilder& fb,
     Observable::ObservableOption observable,
@@ -229,34 +257,6 @@ std::tuple<nested_vector2<me_int_t>, nested_vector2<me_int_t>, Type> build_indic
 
 } // namespace
 
-bool Observable::mirror_invariant(ObservableOption observable) {
-    // The initial-state mirror is the rotation by pi about the x axis
-    // (E, px, py, pz) -> (E, px, -py, -pz). Every case is listed so that a new
-    // observable has to declare which side it falls on.
-    switch (observable) {
-    case obs_py:        // -> -py
-    case obs_pz:        // -> -pz
-    case obs_phi:       // atan2(py, px) -> -phi
-    case obs_theta:     // -> pi - theta
-    case obs_y:         // -> -y
-    case obs_eta:       // -> -eta
-    case obs_delta_eta: // signed difference of two flipped etas -> -delta_eta
-    case obs_delta_phi: // signed difference of two flipped phis -> -delta_phi
-        return false;
-    case obs_e:
-    case obs_px:
-    case obs_mass:
-    case obs_pt:
-    case obs_p_mag:
-    case obs_y_abs:
-    case obs_eta_abs:
-    case obs_delta_r: // sqrt(delta_eta^2 + delta_phi^2): both flip together
-    case obs_pair_mass:
-    case obs_sqrt_s:
-        return true;
-    }
-}
-
 bool Observable::mirror_invariant() const {
     // An observable that matched no particle is the constant 0, whatever the
     // orientation: a cut naming a particle this process does not have imposes
@@ -264,8 +264,8 @@ bool Observable::mirror_invariant() const {
     if (not_found()) {
         return true;
     }
-    return mirror_invariant(_observable) &&
-        (!_order_observable || mirror_invariant(_order_observable.value()));
+    return option_mirror_invariant(_observable) &&
+        (!_order_observable || option_mirror_invariant(_order_observable.value()));
 }
 
 const std::vector<int> Observable::jet_pids{1, 2, 3, 4, -1, -2, -3, -4, 21};
diff --git a/madspace/src/python/madspace.cpp b/madspace/src/python/madspace.cpp
index f07f5d375c..41d71451cf 100644
--- a/madspace/src/python/madspace.cpp
+++ b/madspace/src/python/madspace.cpp
@@ -889,7 +889,8 @@ PYBIND11_MODULE(_madspace_py, m) {
     )
         .def(
             "mirror_invariant",
-            py::overload_cast<>(&Observable::mirror_invariant, py::const_)
+            &Observable::mirror_invariant,
+            pydoc::doc("Observable::mirror_invariant")
         )
         .def_readonly_static("jet_pids", &Observable::jet_pids)
         .def_readonly_static("bottom_pids", &Observable::bottom_pids)
diff --git a/madspace/tests/test_beam_frame.py b/madspace/tests/test_beam_frame.py
index 53b5465668..497c4cfcaa 100644
--- a/madspace/tests/test_beam_frame.py
+++ b/madspace/tests/test_beam_frame.py
@@ -93,9 +93,9 @@ def test_beam_rapidity_cut_in_lab_frame(rng):
     e_cm = 2 * math.sqrt(e1 * e2)
     y0 = 0.5 * math.log(e1 / e2)
     pids = [2, -2, 11, -11]
-    cuts = ms.Cuts([
-        ms.CutItem(observable=ms.Observable(pids, "eta", [[11, -11]]), min=0.0)
-    ])
+    cuts = ms.Cuts(
+        [ms.CutItem(observable=ms.Observable(pids, "eta", [[11, -11]]), min=0.0)]
+    )
     mapping = ms.PhaseSpaceMapping(
         [0.0] * 4, e_cm, cuts=cuts, beam_rapidity=y0, mode=ms.PhaseSpaceMapping.rambo
     )
@@ -169,14 +169,16 @@ def eta(p):
 
 
 def lepton_cut(observable, **bounds):
-    return ms.Cuts([
-        ms.CutItem(
-            observable=ms.Observable(
-                PIDS, observable, [[11, -11]], name=f"lepton-{observable}"
-            ),
-            **bounds,
-        )
-    ])
+    return ms.Cuts(
+        [
+            ms.CutItem(
+                observable=ms.Observable(
+                    PIDS, observable, [[11, -11]], name=f"lepton-{observable}"
+                ),
+                **bounds,
+            )
+        ]
+    )
 
 
 @pytest.mark.parametrize(
@@ -204,8 +206,12 @@ def test_observable_mirror_invariance(observable, invariant):
 
 @pytest.mark.parametrize(
     "observable, invariant",
-    [("delta_eta", False), ("delta_phi", False), ("delta_r", True),
-     ("pair_mass", True)],
+    [
+        ("delta_eta", False),
+        ("delta_phi", False),
+        ("delta_r", True),
+        ("pair_mass", True),
+    ],
 )
 def test_pairwise_observable_mirror_invariance(observable, invariant):
     obs = ms.Observable(PIDS, observable, [[11, -11], [11, -11]])
@@ -233,20 +239,22 @@ def test_cut_on_absent_particle_is_mirror_invariant():
 
 
 def test_cuts_report_the_offending_cuts_by_name():
-    cuts = ms.Cuts([
-        ms.CutItem(
-            observable=ms.Observable(PIDS, "pt", [[11, -11]], name="lepton-pt"),
-            min=10.0,
-        ),
-        ms.CutItem(
-            observable=ms.Observable(PIDS, "eta", [[11, -11]], name="lepton-eta"),
-            min=0.0,
-        ),
-        ms.CutItem(
-            observable=ms.Observable(PIDS, "y", [[11, -11]], name="lepton-y"),
-            min=0.0,
-        ),
-    ])
+    cuts = ms.Cuts(
+        [
+            ms.CutItem(
+                observable=ms.Observable(PIDS, "pt", [[11, -11]], name="lepton-pt"),
+                min=10.0,
+            ),
+            ms.CutItem(
+                observable=ms.Observable(PIDS, "eta", [[11, -11]], name="lepton-eta"),
+                min=0.0,
+            ),
+            ms.CutItem(
+                observable=ms.Observable(PIDS, "y", [[11, -11]], name="lepton-y"),
+                min=0.0,
+            ),
+        ]
+    )
     assert not cuts.mirror_invariant()
     assert cuts.non_mirror_invariant_cuts() == ["lepton-eta", "lepton-y"]
     mapping = ms.PhaseSpaceMapping([0.0] * 4, 13000.0, cuts=cuts)
@@ -267,8 +275,11 @@ def test_mirror_beams_cuts_the_written_orientation(rng, observable, bounds, sepa
     """
     e_cm, mode = 13000.0, ms.PhaseSpaceMapping.rambo
     cut = ms.PhaseSpaceMapping(
-        [0.0] * 4, e_cm, cuts=lepton_cut(observable, **bounds),
-        mirror_beams=True, mode=mode
+        [0.0] * 4,
+        e_cm,
+        cuts=lepton_cut(observable, **bounds),
+        mirror_beams=True,
+        mode=mode,
     )
     free = ms.PhaseSpaceMapping([0.0] * 4, e_cm, mirror_beams=True, mode=mode)
     r = rng.random((N, cut.random_dim()))
@@ -314,8 +325,9 @@ def running_coupling(tmp_path_factory):
 PID_OPTIONS = [[2, -2], [-2, 2]]
 
 
-def build_integrand(running_coupling, cuts=None, mirror_beams=False,
-                    flavor_mirror=(True, True)):
+def build_integrand(
+    running_coupling, cuts=None, mirror_beams=False, flavor_mirror=(True, True)
+):
     """A 2 -> 2 integrand over a beam-swapped pair of flavors. The matrix
     element is never called; only the compute graph is built."""
     e_cm = 13000.0

From 1120a165c013eea80b8aa1f788b8321b0b75116f Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Sat, 10 Oct 2026 17:01:28 +0200
Subject: [PATCH 5/5] mg7 launch: madevent beam names onto the per-beam [beam]
 entries

At the launch question, set ebeam/ebeam1/ebeam2, lpp/lpp1/lpp2,
pdlabel/pdlabel1/pdlabel2 and lhaid were "invalid set command" (except the
ebeam ones, which reached _LO_SCALAR_MAP) and the run kept the defaults.
RunCardMG7.set_lo_beam_param now maps them, a trailing 1/2 acting on that
beam and the other names on both:

- ebeam*: beam.ebeam1/2, with units and mass expressions.
- pdlabel*/lhaid: beam.pdf1/2 through the label and lhaid tables of the
  run_card.dat conversion; an lhaid outside the table is looked up in the
  pdfsets.index of the LHAPDF data paths. 'lhapdf' alone changes nothing,
  'none' on both beams is beam.leptonic = true.
- lpp*: the single beam.leptonic switch (0, +-3, +-4 lepton, +-1 hadron),
  with a warning when the two beams were given different types.

A value mg7 cannot represent is reported and leaves the card alone.

from_LO keeps different pdlabel1/pdlabel2 per beam instead of reporting
them as dropped.
---
 docs/source/madevent_to_mg7.rst               |   8 +-
 madgraph/iolibs/template_files/mg7/launch.py  |  13 ++
 madgraph/various/banner.py                    | 140 ++++++++++++++++--
 tests/unit_tests/interface/test_mg7_launch.py |  20 +++
 tests/unit_tests/various/test_banner.py       |  52 +++++++
 5 files changed, 216 insertions(+), 17 deletions(-)

diff --git a/docs/source/madevent_to_mg7.rst b/docs/source/madevent_to_mg7.rst
index 8b91e21166..74fd850fce 100644
--- a/docs/source/madevent_to_mg7.rst
+++ b/docs/source/madevent_to_mg7.rst
@@ -15,7 +15,9 @@ existing scripts keep working. These are ``nevents``, ``gridpack``, ``fixed_ren_
 ``scale``, ``dsqrt_q2fact1``, ``dsqrt_q2fact2``, ``scalefact``, ``bwcutoff``, ``cut_decays``
 and ``use_syst``. The shortcuts ``set lhc``, ``set lep``, ``set fixed_scale``,
 ``set no_parton_cut`` and ``set dynamical_scale_choice HT/n`` work as before, and so does
-``set iseed``. MadEvent draws a random seed for ``iseed = 0``. MG7 does the same for
+``set iseed``. The beam settings ``ebeam``, ``ebeam1``, ``ebeam2``, ``lpp``, ``lpp1``,
+``lpp2``, ``pdlabel``, ``pdlabel1``, ``pdlabel2`` and ``lhaid`` are translated as in the
+table below. MadEvent draws a random seed for ``iseed = 0``. MG7 does the same for
 ``seed = -1`` and treats ``0`` as an ordinary fixed seed, so ``set iseed 0`` is translated.
 All other settings must use their new names.
 
@@ -43,8 +45,8 @@ Beams and PDFs
         ``"NNPDF40_lo_as_01180"``. ``set beam.pdf`` sets both beams.
     * - ``pdlabel1``, ``pdlabel2``
       - ``beam.pdf1``, ``beam.pdf2``
-      - The conversion of a ``run_card.dat`` uses the same set for both beams. Set
-        ``beam.pdf1`` and ``beam.pdf2`` by hand to use different ones.
+      - One set per beam. A beam without PDF exists only for a lepton collider, as
+        ``beam.leptonic``.
 
 Scales
 ------
diff --git a/madgraph/iolibs/template_files/mg7/launch.py b/madgraph/iolibs/template_files/mg7/launch.py
index 8b15bc4348..a8bfedc526 100644
--- a/madgraph/iolibs/template_files/mg7/launch.py
+++ b/madgraph/iolibs/template_files/mg7/launch.py
@@ -3645,6 +3645,19 @@ def do_set(self, line, *args, **kwargs):
                         self.modified_card.add("run")
                         return
 
+                    # madevent beam settings (ebeam, lpp, pdlabel, lhaid, per
+                    # beam or for both) onto the per-beam [beam] entries
+                    if rest and nlow in run_card.lo_beam_params:
+                        data_paths = lhapdf_paths().data_paths if nlow == "lhaid" else ()
+                        try:
+                            done = run_card.set_lo_beam_param(nlow, rest, masses, data_paths)
+                        except _banner_mod.InvalidRunCard as error:
+                            logger.warning("ignoring 'set %s %s': %s", name, rest, error)
+                            return
+                        logger.info("set %s %s: %s", nlow, rest, done)
+                        self.modified_card.add("run")
+                        return
+
                     # quantities the card stores per beam: set both beams
                     if rest and nlow in ("e_cm", "beam.e_cm", "beam.pdf"):
                         value = rest if nlow == "beam.pdf" else run_card.evaluate(rest, masses)
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index 4de7c55e4b..a1aaa5c59d 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -7427,6 +7427,85 @@ def set_collider(self, name, arg, masses=None):
         self.set('beam.leptonic', name in ('lep', 'ilc'), user=True)
         return float(val)
 
+    # madevent beam settings understood by set_lo_beam_param
+    lo_beam_params = ('ebeam', 'ebeam1', 'ebeam2', 'lpp', 'lpp1', 'lpp2',
+                      'pdlabel', 'pdlabel1', 'pdlabel2', 'lhaid')
+
+    def set_lo_beam_param(self, name, value, masses=None, data_paths=()):
+        """``set <name> <value>`` with a madevent beam setting (lo_beam_params);
+        returns a description of what was set. A name ending in 1 or 2 acts
+        on that beam, the others on both.
+
+        * ebeam: the beam energies, with units and expressions as for any
+          energy.
+        * pdlabel / lhaid: the PDF set, through the label and id tables of
+          the run_card.dat conversion (an lhaid also through pdfsets.index).
+          'lhapdf' alone changes nothing, since the set comes from lhaid;
+          'none' on both beams is a lepton collider.
+        * lpp: mg7 has one switch for both beams, beam.leptonic, true for
+          0 (no PDF) and +-3/+-4 (e/mu), false for +-1.
+
+        Raises InvalidRunCard for a value mg7 cannot represent."""
+        name = name.lower()
+        beams = (name[-1],) if name[-1] in '12' else ('1', '2')
+        which = 'beam ' + beams[0] if len(beams) == 1 else 'both beams'
+        value = str(value).strip()
+
+        if name.startswith('ebeam'):
+            energy = self.evaluate(value, masses)
+            if isinstance(energy, str):
+                energy = self.format_variable(energy, float, name=name)
+            for beam in beams:
+                self.set('beam.ebeam' + beam, float(energy), user=True)
+            return 'beam energy of %s: %s GeV' % (which, float(energy))
+
+        if name.startswith('lpp'):
+            try:
+                lpp = int(value)
+            except ValueError:
+                raise InvalidRunCard('%s = %s: not an integer' % (name, value))
+            if abs(lpp) not in (0, 1, 3, 4):
+                raise InvalidRunCard('%s = %s: beam type not representable in mg7'
+                                     % (name, lpp))
+            leptonic = abs(lpp) != 1
+            beam_types = getattr(self, '_lo_beam_leptonic', {})
+            beam_types.update((beam, leptonic) for beam in beams)
+            self._lo_beam_leptonic = beam_types
+            if len(set(beam_types.values())) > 1:
+                logger.warning('mg7 cannot collide a lepton beam with a hadron '
+                               'beam: beam.leptonic = %s now applies to both',
+                               leptonic)
+            if lpp == -1:
+                logger.warning('%s = -1: mg7 has no antiproton beam and keeps '
+                               'the proton PDF', name)
+            self.set('beam.leptonic', leptonic, user=True)
+            return 'beam.leptonic = %s' % leptonic
+
+        if name == 'lhaid':
+            label, lhaid = 'lhapdf', value.split()[0] if value else value
+            try:
+                lhaid = int(lhaid)
+            except ValueError:
+                raise InvalidRunCard('lhaid = %s: not an integer' % value)
+        else:
+            label, lhaid = value, None
+            if label.lower() == 'lhapdf':
+                return ('unchanged: mg7 always uses LHAPDF, the set is chosen '
+                        'with lhaid or beam.pdf1/pdf2')
+            if label.lower() in ('none', 'no'):
+                if len(beams) == 1:
+                    raise InvalidRunCard('%s = none: only a lepton collider has '
+                                         'no PDF in mg7, on both beams' % name)
+                self.set('beam.leptonic', True, user=True)
+                return 'beam.leptonic = True'
+        pdf = self._resolve_pdf_label(label, lhaid, [], data_paths)
+        if pdf is None:
+            raise InvalidRunCard('%s = %s: unknown %s' % (
+                name, value, 'LHAPDF id' if name == 'lhaid' else 'PDF label'))
+        for beam in beams:
+            self.set('beam.pdf' + beam, pdf, user=True)
+        return 'PDF of %s: %s' % (which, pdf)
+
     def set_fixed_scale(self, value, masses=None):
         """``set fixed_scale X``: switch to fixed scales and set them all to X
         (X may use units, arithmetic or a mass reference, e.g. ``mz``)."""
@@ -8110,27 +8189,60 @@ def remove_jet_cuts(self):
     }
 
     @classmethod
-    def _resolve_pdf(cls, lo, dropped):
-        """Map the LO pdlabel/lhaid to an LHAPDF set name (or None to keep the
-        MG7 default, recording the reason in ``dropped``)."""
-        pdlabel = str(lo['pdlabel']).lower() if 'pdlabel' in lo else ''
-        if str(lo['pdlabel1']).lower() != str(lo['pdlabel2']).lower():
-            dropped.append('pdlabel1/pdlabel2 differ (mg7 uses a single pdf)')
+    def _lhaid_to_pdf_name(cls, lhaid, data_paths=()):
+        """LHAPDF set name of the set whose central member is ``lhaid``: the
+        common ids of _LO_LHAID_MAP, else the pdfsets.index of ``data_paths``
+        or $LHAPDF_DATA_PATH; None when neither knows it."""
+        if lhaid in cls._LO_LHAID_MAP:
+            return cls._LO_LHAID_MAP[lhaid]
+        search = list(data_paths)
+        if os.environ.get('LHAPDF_DATA_PATH'):
+            search += os.environ['LHAPDF_DATA_PATH'].split(os.pathsep)
+        for base in search:
+            try:
+                with open(os.path.join(base, 'pdfsets.index')) as fsock:
+                    for line in fsock:
+                        parts = line.split()
+                        if len(parts) >= 2 and parts[0] == str(lhaid):
+                            return parts[1]
+            except OSError:
+                continue
+        return None
+
+    @classmethod
+    def _resolve_pdf_label(cls, pdlabel, lhaid, dropped, data_paths=()):
+        """Map one LO pdlabel (and the lhaid for 'lhapdf') to an LHAPDF set
+        name, or None to keep the MG7 default, recording the reason in
+        ``dropped`` (nothing for 'none': that beam has no PDF)."""
+        pdlabel = str(pdlabel).lower()
         if pdlabel in ('none', '', 'no'):
             return None
         if pdlabel == 'lhapdf':
-            lhaid = lo['lhaid']
             if isinstance(lhaid, list):
                 lhaid = lhaid[0]
-            if lhaid in cls._LO_LHAID_MAP:
-                return cls._LO_LHAID_MAP[lhaid]
-            dropped.append('lhaid=%s (unknown LHAPDF id; keeping mg7 default pdf)' % lhaid)
-            return None
+            name = cls._lhaid_to_pdf_name(lhaid, data_paths)
+            if name is None:
+                dropped.append('lhaid=%s (unknown LHAPDF id; keeping mg7 default pdf)' % lhaid)
+            return name
         if pdlabel in cls._LO_PDF_LABEL_MAP:
             return cls._LO_PDF_LABEL_MAP[pdlabel]
         dropped.append('pdlabel=%s (unknown PDF label; keeping mg7 default pdf)' % pdlabel)
         return None
 
+    @classmethod
+    def _resolve_pdf(cls, lo, dropped):
+        """Map the LO pdlabel/lhaid to an LHAPDF set name per beam, (pdf1,
+        pdf2), None keeping the MG7 default of that beam. pdlabel1/pdlabel2
+        are only read when they differ: pdlabel holds the common value."""
+        lhaid = lo['lhaid'] if 'lhaid' in lo else None
+        labels = [str(lo['pdlabel%d' % i]).lower() if 'pdlabel%d' % i in lo else ''
+                  for i in (1, 2)]
+        if labels[0] == labels[1]:
+            pdlabel = lo['pdlabel'] if 'pdlabel' in lo else ''
+            name = cls._resolve_pdf_label(pdlabel, lhaid, dropped)
+            return name, name
+        return tuple(cls._resolve_pdf_label(label, lhaid, dropped) for label in labels)
+
     @classmethod
     def from_LO(cls, lo, warn=True):
         """Build a :class:`RunCardMG7` reproducing, as far as possible, the setup
@@ -8199,9 +8311,9 @@ def nondefault(name):
             mg7.set('run.me_frame', list(lo['me_frame']) or [1, 2])
 
         # --- PDF ---
-        pdf_name = cls._resolve_pdf(lo, dropped)
-        if pdf_name:
-            mg7.set('beam.pdf', pdf_name)
+        for beam, pdf_name in zip((1, 2), cls._resolve_pdf(lo, dropped)):
+            if pdf_name:
+                mg7.set('beam.pdf%d' % beam, pdf_name)
 
         # --- maxjetflavor -> jet multiparticle ---
         if 'maxjetflavor' in lo:
diff --git a/tests/unit_tests/interface/test_mg7_launch.py b/tests/unit_tests/interface/test_mg7_launch.py
index b5874b043a..9b5180e2ee 100644
--- a/tests/unit_tests/interface/test_mg7_launch.py
+++ b/tests/unit_tests/interface/test_mg7_launch.py
@@ -676,6 +676,26 @@ def test_set_histograms_rejects_anything_else(self):
         self.assertEqual(dict(obj.run_card['histograms']), before)
         self.assertNotIn('run', obj.modified_card)
 
+    def test_set_madevent_beam_names(self):
+        """set ebeam1/ebeam/pdlabel2/lpp, as a madevent launch script writes
+        them, edit the per-beam [beam] entries"""
+        obj = self.histogram_selector()
+        obj.do_set('ebeam 6.8 TeV')
+        obj.do_set('ebeam1 7 TeV')
+        obj.do_set('pdlabel2 cteq6l1')
+        obj.do_set('lpp 1')
+        beam = obj.run_card['beam']
+        self.assertEqual((beam['ebeam1'], beam['ebeam2']), (7000.0, 6800.0))
+        self.assertEqual(beam['pdf2'], 'cteq6l1')
+        self.assertNotEqual(beam['pdf1'], 'cteq6l1')
+        self.assertFalse(beam['leptonic'])
+        self.assertIn('run', obj.modified_card)
+        # a value mg7 cannot represent leaves the card alone
+        obj.modified_card.clear()
+        obj.do_set('pdlabel1 none')
+        self.assertNotEqual(obj.run_card['beam']['pdf1'], 'none')
+        self.assertNotIn('run', obj.modified_card)
+
     def test_no_post_processing_keeps_the_npy_output(self):
         self.assertEqual(self.output_format({}), 'compact_npy')
         self.assertEqual(self.output_format(None), 'compact_npy')
diff --git a/tests/unit_tests/various/test_banner.py b/tests/unit_tests/various/test_banner.py
index 73ca16b3e6..65b73fecf2 100755
--- a/tests/unit_tests/various/test_banner.py
+++ b/tests/unit_tests/various/test_banner.py
@@ -1648,6 +1648,44 @@ def test_per_beam_energy_and_pdf(self):
         rc.set('beam.ebeam2', 0, user=True)
         self.assertRaises(bannermod.InvalidRunCard, rc.check_validity)
 
+    def test_set_lo_beam_param(self):
+        """madevent beam settings (set ebeam/lpp/pdlabel/lhaid at the launch
+        question) go to the per-beam entries"""
+        rc = bannermod.RunCardMG7()
+        beams = lambda key: (rc['beam'][key + '1'], rc['beam'][key + '2'])
+        rc.set_lo_beam_param('ebeam', '6.8 TeV')
+        self.assertEqual(beams('ebeam'), (6800.0, 6800.0))
+        rc.set_lo_beam_param('ebeam2', 'mz*40', masses={'mz': 91.188})
+        self.assertEqual(beams('ebeam'), (6800.0, 3647.52))
+        rc.set_lo_beam_param('pdlabel', 'nn23lo1')
+        self.assertEqual(beams('pdf'), ('NNPDF23_lo_as_0130_qed',) * 2)
+        rc.set_lo_beam_param('pdlabel1', 'cteq6l1')
+        self.assertEqual(beams('pdf'), ('cteq6l1', 'NNPDF23_lo_as_0130_qed'))
+        # the set is chosen by lhaid: 'lhapdf' alone changes nothing
+        rc.set_lo_beam_param('pdlabel', 'lhapdf')
+        self.assertEqual(beams('pdf'), ('cteq6l1', 'NNPDF23_lo_as_0130_qed'))
+        # an id outside the built-in table, through the LHAPDF index
+        index_dir = tempfile.mkdtemp()
+        try:
+            with open(os.path.join(index_dir, 'pdfsets.index'), 'w') as fsock:
+                fsock.write('244800 NNPDF23_nlo_as_0119_qed 5\n')
+            rc.set_lo_beam_param('lhaid', '244800', data_paths=[index_dir])
+        finally:
+            shutil.rmtree(index_dir)
+        self.assertEqual(beams('pdf'), ('NNPDF23_nlo_as_0119_qed',) * 2)
+        for name, value in (('lhaid', '999999'), ('pdlabel2', 'isronlyll'),
+                            ('pdlabel2', 'none'), ('lpp1', '2'), ('lpp', 'p')):
+            self.assertRaises(bannermod.InvalidRunCard,
+                              rc.set_lo_beam_param, name, value)
+        self.assertEqual(beams('pdf'), ('NNPDF23_nlo_as_0119_qed',) * 2)
+        # one lepton/hadron switch for both beams
+        rc.set_lo_beam_param('lpp', '0')
+        self.assertTrue(rc['beam']['leptonic'])
+        rc.set_lo_beam_param('lpp1', '1')
+        self.assertFalse(rc['beam']['leptonic'])
+        rc.set_lo_beam_param('pdlabel', 'none')
+        self.assertTrue(rc['beam']['leptonic'])
+
     def test_read_card_with_e_cm_and_pdf(self):
         """a card written before the per-beam keys still reads"""
         rc = bannermod.RunCardMG7()
@@ -1777,6 +1815,20 @@ def test_from_LO_conversion(self):
         mg7.write(buf, template=self.template)
         tomllib.loads(buf.getvalue())
 
+    def test_from_LO_conversion_per_beam(self):
+        """asymmetric beam energies and PDFs are kept per beam"""
+        lo = bannermod.RunCardLO()
+        lo['ebeam1'] = 7000
+        lo['ebeam2'] = 4000
+        lo['pdlabel1'] = 'nn23lo1'
+        lo['pdlabel2'] = 'cteq6l1'
+        mg7, dropped = bannermod.RunCardMG7.from_LO(lo, warn=False)
+        self.assertEqual((mg7['beam']['ebeam1'], mg7['beam']['ebeam2']), (7000.0, 4000.0))
+        self.assertAlmostEqual(mg7['beam']['e_cm'], 2 * (7000.0 * 4000.0) ** 0.5)
+        self.assertEqual((mg7['beam']['pdf1'], mg7['beam']['pdf2']),
+                         ('NNPDF23_lo_as_0130_qed', 'cteq6l1'))
+        self.assertNotIn('pdlabel', ' '.join(dropped))
+
     def test_from_LO_conversion_of_composite_cuts(self):
         """ordered, HT, summed-momentum, energy and per-pdg LO cuts are ported"""
         lo = bannermod.RunCardLO()