madspace: always draw the beam orientation before the mapping

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.
 
@@ -30,19 +32,21 @@ Beams and PDFs
       - ``run_card.toml``
       - Notes
     * - ``ebeam1``, ``ebeam2``
-      - ``beam.e_cm``
-      - The sum of the two energies. Only symmetric beams are supported.
+      - ``beam.ebeam1``, ``beam.ebeam2``
+      - One energy per beam. With different energies the events are written, and the
+        :math:`\eta` cuts applied, in the lab frame.
     * - ``lpp1``, ``lpp2``
       - ``beam.leptonic``
       - ``false`` for proton beams and ``true`` for lepton beams without PDFs. Other beam
         types, such as elastic photons or antiprotons, are not supported.
     * - ``pdlabel``, ``lhaid``
-      - ``beam.pdf``
+      - ``beam.pdf1``, ``beam.pdf2``
       - The name of the LHAPDF set instead of its id. For example, ``lhaid = 331900`` becomes
-        ``"NNPDF40_lo_as_01180"``.
+        ``"NNPDF40_lo_as_01180"``. ``set beam.pdf`` sets both beams.
     * - ``pdlabel1``, ``pdlabel2``
-      - ``beam.pdf``
-      - MG7 uses one PDF set for both beams.
+      - ``beam.pdf1``, ``beam.pdf2``
+      - One set per beam. A beam without PDF exists only for a lepton collider, as
+        ``beam.leptonic``.
 
 Scales
 ------
@@ -212,8 +216,8 @@ The following MadEvent features have no MG7 equivalent. Settings for them are ig
 * Beam polarization and heavy-ion or equivalent-photon beams: ``polbeam1``, ``polbeam2``,
   ``nb_proton1``, ``nb_proton2``, ``nb_neutron1``, ``nb_neutron2``, ``mass_ion1``,
   ``mass_ion2``, ``ievo_eva``, ``evaorder``, ``eva_xcut``.
-* Different PDFs or fixed factorization scales for the two beams, and the scale
-  parameters ``mue_over_ref``, ``mue_ref_fixed`` and ``fixed_extra_scale``.
+* Fixed factorization scales for the two beams, and the scale parameters
+  ``mue_over_ref``, ``mue_ref_fixed`` and ``fixed_extra_scale``.
 * Cuts that are not a bound on a single observable: ``cutuse = 1``, which accepts an event if
   any of the ordered jet cuts passes, and the vector-boson-fusion cuts ``xetamin`` and
   ``deltaeta``, which require the two hardest jets to lie in opposite hemispheres.

docs/source/run_card.rst

@@ -118,17 +118,22 @@ See :doc:`gridpacks` for details.
 [beam]
 ------
 
-``e_cm`` (``13000.0``)
-    Center-of-mass energy in GeV. Units are accepted, for example ``"13 TeV"``.
+``ebeam1`` (``6500.0``), ``ebeam2`` (``6500.0``)
+    Energy of each beam in GeV, beam 1 moving along :math:`+z`. Units are accepted, for
+    example ``"7 TeV"``. The collision energy is :math:`2\sqrt{E_1 E_2}`. With different
+    energies the events are written, and the :math:`\eta` cuts applied, in this lab frame.
+    ``e_cm`` can still be read and set: ``set e_cm 13 TeV`` sets both beams to 6.5 TeV.
 
 ``leptonic`` (``false``)
-    Treat the beams as leptons: no PDFs are used and the partonic energy equals ``e_cm``.
-    This is set automatically for lepton-collider processes, together with ``e_cm = 1000``
-    and no jet cuts.
-
-``pdf`` (``"NNPDF40_lo_as_01180"``)
-    Name of the LHAPDF set. MG7 tries to download it if it is not installed. The default is
-    the 5-flavor NNPDF 4.0 LO set, which has 100 error members for the PDF variations.
+    Treat the beams as leptons: no PDFs are used and the partonic energy equals the
+    collision energy. This is set automatically for lepton-collider processes, together
+    with 500 GeV beams and no jet cuts.
+
+``pdf1`` (``"NNPDF40_lo_as_01180"``), ``pdf2`` (``"NNPDF40_lo_as_01180"``)
+    Name of the LHAPDF set of each beam. ``set beam.pdf X`` sets both. MG7 tries to
+    download a set that is not installed. The default is the 5-flavor NNPDF 4.0 LO set,
+    which has 100 error members for the PDF variations. With different sets on the two
+    beams the PDF member variations of ``[systematics]`` are dropped.
 
 ``fixed_ren_scale`` (``false``), ``fixed_fact_scale`` (``false``)
     Use the fixed values below instead of the dynamical scale choice. The two switches act

madgraph/iolibs/template_files/mg7/gridpack.py

@@ -346,9 +346,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
@@ -397,7 +400,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

madgraph/iolibs/template_files/mg7/launch.py

@@ -1,5 +1,6 @@
 import argparse
 import gc
+import math
 import os
 import sys
 import time
@@ -399,6 +400,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 = {
@@ -638,7 +647,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"
@@ -668,9 +677,18 @@ 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.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 = {
@@ -704,26 +722,45 @@ def init_beam(self) -> None:
                 scale_factor=1.,
             )
             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:
                 self.alphas_grid.initialize_globals(context)
             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)
         self.running_coupling = ms.RunningCoupling(self.alphas_grid)
 
@@ -904,7 +941,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()
@@ -916,7 +960,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 = {
@@ -925,6 +970,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],
@@ -1773,7 +1821,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.
@@ -1782,22 +1829,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:"):
@@ -1816,7 +1864,6 @@ def build_lhe_meta(self, systematics=None):
         building meta for a writer that bypasses that injection (e.g.
         write_lhe_header)."""
         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
@@ -1844,7 +1891,7 @@ def build_lhe_meta(self, systematics=None):
             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=weight_mode,
             # positional: the pybind arg name for max_weight is non-kwarg-safe
             processes=[ms.LHEProcess(xsec, err, max_weight, 1)],
@@ -2423,6 +2470,9 @@ def cut_pids(observable_kwargs):
             if len(self.process.cut_data) > 0
             else None
         )
+        # A flavor standing for both beam orientations draws the orientation
+        # before the phase-space mapping, so boost, cuts and PDFs all see it.
+        self.mirror_beams = not self.process.is_decay and self.has_mirrored_flavors()
         # the integration histograms are functions of the momenta; the weight
         # distribution is a property of the final event sample, so it is only
         # filled by MadgraphProcess.build_event_histograms at combine time
@@ -2501,6 +2551,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.mirror_beams,
                 )
                 prefix = f"subproc{self.subproc_id}.channel{channel_id}"
                 if topo_count > 1:
@@ -2554,13 +2606,26 @@ 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,
             self.process.e_cm,
             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.mirror_beams,
         )
         prefix = f"subproc{self.subproc_id}.flat"
         discrete_sym, discrete_flavor = self.build_discrete(
@@ -3171,7 +3236,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())
@@ -3611,6 +3677,28 @@ 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)
+                        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