diff --git a/docs/source/madevent_to_mg7.rst b/docs/source/madevent_to_mg7.rst index 7419a4f45..74fd850fc 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. @@ -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. diff --git a/docs/source/run_card.rst b/docs/source/run_card.rst index 9893e6473..e32aa4c3c 100644 --- a/docs/source/run_card.rst +++ b/docs/source/run_card.rst @@ -116,17 +116,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 diff --git a/madgraph/iolibs/template_files/mg7/gridpack.py b/madgraph/iolibs/template_files/mg7/gridpack.py index e0842f702..8029cfed3 100644 --- a/madgraph/iolibs/template_files/mg7/gridpack.py +++ b/madgraph/iolibs/template_files/mg7/gridpack.py @@ -342,9 +342,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 @@ -393,7 +396,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 2aa49beed..a8bfedc52 100644 --- a/madgraph/iolibs/template_files/mg7/launch.py +++ b/madgraph/iolibs/template_files/mg7/launch.py @@ -1,5 +1,6 @@ import argparse import gc +import math import os import sys import time @@ -348,6 +349,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 = { @@ -587,7 +596,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" @@ -617,9 +626,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 = { @@ -653,8 +672,10 @@ 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: @@ -662,18 +683,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) # ------------------------------------------------------------------ @@ -853,7 +896,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() @@ -865,7 +915,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 = { @@ -874,6 +925,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], @@ -1716,7 +1770,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. @@ -1725,22 +1778,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:"): @@ -1759,7 +1813,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 @@ -1786,7 +1839,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=3, # positional: the pybind arg name for max_weight is non-kwarg-safe processes=[ms.LHEProcess(xsec, err, xsec, 1)], @@ -2365,6 +2418,24 @@ def cut_pids(observable_kwargs): 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 # 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 @@ -2443,6 +2514,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: @@ -2496,6 +2569,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, @@ -2503,6 +2587,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.mirror_beams, ) prefix = f"subproc{self.subproc_id}.flat" discrete_sym, discrete_flavor = self.build_discrete( @@ -3112,7 +3198,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()) @@ -3558,6 +3645,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 diff --git a/madgraph/iolibs/template_files/mg7/run_card.toml b/madgraph/iolibs/template_files/mg7/run_card.toml index 05b92450e..1dccd5332 100644 --- a/madgraph/iolibs/template_files/mg7/run_card.toml +++ b/madgraph/iolibs/template_files/mg7/run_card.toml @@ -45,9 +45,13 @@ include_madspace = %(gridpack.include_madspace)s include_madspace_source = %(gridpack.include_madspace_source)s [beam] -e_cm = %(beam.e_cm)s +# beam 1 moves along +z; events and eta cuts are in this lab frame +ebeam1 = %(beam.ebeam1)s +ebeam2 = %(beam.ebeam2)s leptonic = %(beam.leptonic)s -pdf = %(beam.pdf)s +# "set beam.pdf X" sets 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/banner.py b/madgraph/various/banner.py index 407316f4b..a1aaa5c59 100755 --- a/madgraph/various/banner.py +++ b/madgraph/various/banner.py @@ -6616,7 +6616,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. @@ -6750,12 +6750,17 @@ def default_setup(self): comment="include the madspace source code") # ----------------------------- [beam] ------------------------- - self.add_toml_param('beam', 'e_cm', 13000.0, - comment="center-of-mass energy in GeV; units are accepted, e.g. \"13 TeV\"") + self.add_toml_param('beam', 'ebeam1', 6500.0, + comment="energy of beam 1 (moving along +z) in GeV; units are accepted. " + "Events and eta cuts are in this lab frame; set e_cm X sets both beams to X/2") + self.add_toml_param('beam', 'ebeam2', 6500.0, + comment="energy of beam 2 (moving along -z) in GeV; units are accepted") self.add_toml_param('beam', 'leptonic', False, comment="lepton collider: no PDFs, beam energy fixed") - self.add_toml_param('beam', 'pdf', "NNPDF40_lo_as_01180", - comment="LHAPDF set name") + self.add_toml_param('beam', 'pdf1', "NNPDF40_lo_as_01180", + comment="LHAPDF set name of beam 1; set beam.pdf X sets both beams") + self.add_toml_param('beam', 'pdf2', "NNPDF40_lo_as_01180", + comment="LHAPDF set name of beam 2; PDF member variations need the same set on both beams") self.add_toml_param('beam', 'fixed_ren_scale', False, comment="use ren_scale instead of the dynamical scale") self.add_toml_param('beam', 'fixed_fact_scale', False, @@ -6996,8 +7001,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): @@ -7066,7 +7086,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'): @@ -7075,7 +7094,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'): @@ -7124,11 +7143,11 @@ def _legacy_compat(self, key): return float(self['phasespace']['bw_cutoff']) 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 = [] @@ -7206,7 +7225,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: @@ -7248,6 +7270,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() @@ -7380,10 +7423,89 @@ 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) + # 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 `` 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``).""" @@ -7406,6 +7528,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)) beam = self['beam'] if float(beam['scale_factor']) <= 0.: @@ -7578,7 +7704,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): @@ -7946,6 +8073,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', 'scalefact': 'beam.scale_factor', 'scale': 'beam.ren_scale', 'dsqrt_q2fact1': 'beam.fact_scale1', @@ -8060,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 @@ -8103,7 +8265,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 @@ -8115,7 +8276,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: @@ -8150,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/madspace/include/madspace/driver/systematics.hpp b/madspace/include/madspace/driver/systematics.hpp index 4238ec86f..59431c661 100644 --- a/madspace/include/madspace/driver/systematics.hpp +++ b/madspace/include/madspace/driver/systematics.hpp @@ -154,6 +154,10 @@ class SystematicsCalculator { * @param me_flavor_remap `me_flavor_remap[subprocess][flavor_index]` is * the flavor passed to that subprocess's matrix * element. + * @param nominal_pdf2 PDF of the second beam when its set differs from + * @p nominal_pdf (which is then the first beam's). + * PDF member variations are then dropped with a + * warning; the scale variations use both sets. */ SystematicsCalculator( const SystematicsConfig& config, @@ -162,7 +166,8 @@ class SystematicsCalculator { const std::optional& nominal_alpha_s, ContextPtr context = nullptr, const std::vector>& matrix_elements = {}, - const nested_vector2& me_flavor_remap = {} + const nested_vector2& me_flavor_remap = {}, + const std::optional& nominal_pdf2 = std::nullopt ); /// The configuration passed to the constructor. @@ -257,6 +262,7 @@ class SystematicsCalculator { ContextPtr _context; std::string _prefix; std::optional _nominal_pdf; // has_pdf only + std::optional _nominal_pdf2; // second beam, if its set differs std::vector _member_pdfs; // one per member std::vector _alpha_s_runtimes; // index 0: nominal set std::vector> _matrix_elements; diff --git a/madspace/include/madspace/phasespace/cuts.hpp b/madspace/include/madspace/phasespace/cuts.hpp index f3ff571ac..08ada6866 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 +#include #include #include @@ -62,6 +63,14 @@ class Cuts : public FunctionGenerator { /// Build a pass-through mask with no cuts. /// @param particle_count Number of external particles. 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 non_mirror_invariant_cuts() const; + /// Whether every configured cut is invariant under the initial-state mirror. + bool mirror_invariant() const { return non_mirror_invariant_cuts().empty(); } /// Largest required partonic center-of-mass energy, or 0 if unconstrained. double sqrt_s_min() const; /// Per-outgoing-particle maximum pseudorapidity (infinity where inactive). diff --git a/madspace/include/madspace/phasespace/integrand.hpp b/madspace/include/madspace/phasespace/integrand.hpp index 0ca15b8ca..3449ad1c3 100644 --- a/madspace/include/madspace/phasespace/integrand.hpp +++ b/madspace/include/madspace/phasespace/integrand.hpp @@ -126,6 +126,11 @@ class Integrand : public FunctionGenerator { * @param flavor_per_subproc_remap Per-sub-process flavour remap. * @param compressed_channel_weight_count Channel-weight entries kept per * event. + * @param pdf_grid2 PDF of the second beam when it + * differs from @p pdf_grid (which + * is then the first beam's). Its + * globals must be registered with + * the prefix @ref pdf2_prefix. */ Integrand( const PhaseSpaceMapping& mapping, @@ -155,8 +160,12 @@ class Integrand : public FunctionGenerator { const std::vector& flavor_diff_xs_indices = {}, const std::vector& flavor_subproc_indices = {}, const std::vector& flavor_per_subproc_remap = {}, - std::size_t compressed_channel_weight_count = 50 + std::size_t compressed_channel_weight_count = 50, + const std::optional& pdf_grid2 = std::nullopt ); + /// Globals prefix under which the second beam's PDF grid is registered + /// (PdfGrid::initialize_globals) when `pdf_grid2` is given. + inline static const std::string pdf2_prefix = "beam2"; /// Total number of particles, incoming and outgoing. std::size_t particle_count() const { return _mapping.particle_count(); } /// Whether the integrand emits the @ref MadnisLoss training quantities. @@ -232,6 +241,9 @@ class Integrand : public FunctionGenerator { AdaptiveDiscrete _discrete_flavor; nested_vector2 _pid_options; std::array, 2> _pdfs; + // with different PDFs per beam: each leg's density from the other beam, + // for the events whose orientation is mirrored + std::array, 2> _pdfs_swapped; std::array, 2> _pdf_indices; std::optional _running_coupling; std::optional _energy_scale; @@ -254,7 +266,11 @@ class Integrand : public FunctionGenerator { std::vector _flavor_remap; std::vector _flavor_factors; std::vector _flavor_mirror; + std::vector _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 _channel_part_ret_types; std::vector _flavor_diff_xs_indices; std::vector _flavor_subproc_indices; diff --git a/madspace/include/madspace/phasespace/observable.hpp b/madspace/include/madspace/phasespace/observable.hpp index d8e5fd95b..8877c2317 100644 --- a/madspace/include/madspace/phasespace/observable.hpp +++ b/madspace/include/madspace/phasespace/observable.hpp @@ -104,6 +104,13 @@ class Observable : public FunctionGenerator { ); /// The kinematic quantity this observable computes. ObservableOption observable() const { return _observable; } + /// 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. 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& indices() const { return _indices; } /// Whether the tuple momenta are summed before evaluation. diff --git a/madspace/include/madspace/phasespace/phasespace.hpp b/madspace/include/madspace/phasespace/phasespace.hpp index 45d78d977..750cbe24f 100644 --- a/madspace/include/madspace/phasespace/phasespace.hpp +++ b/madspace/include/madspace/phasespace/phasespace.hpp @@ -47,6 +47,9 @@ namespace madspace { * **Conditions** * - `permutation_index` – `int`, shape `(batch,)` – which permutation channel * to use. Present only when more than one permutation is given. + * - `mirror_index` – `int`, shape `(batch,)` – 1 to take leg 1 from beam 2 and + * leg 2 from beam 1. Present only when `mirror_beams` is set; it follows + * `permutation_index`. * * **Outputs** * - `momenta` – `float`, shape `(batch, n_out + 2, 4)` – all momenta (incoming @@ -85,6 +88,18 @@ class PhaseSpaceMapping : public Mapping { * More than one enables the `permutation_index` * condition. * @param color_order Color ordering for @ref TChannelMode::color_ordered. + * @param beam_rapidity Rapidity of the beams' centre-of-mass frame in the + * lab frame the momenta are returned in: + * 0.5 ln(E1 / E2) for beam energies E1 (along +z) + * and E2, with @p cm_energy = 2 sqrt(E1 E2). The + * phase space is generated in the beams' frame and + * boosted into the lab before the cuts, so rapidity + * cuts act in the lab. + * @param mirror_beams Add the `mirror_index` condition: for index 1 the + * event is rotated by pi about the x axis in the + * beams' frame (leg 1 from beam 2, leg 2 from beam 1) + * before the lab boost and the cuts. With identical + * beams that is a rotation of the lab event too. */ PhaseSpaceMapping( const Topology& topology, @@ -94,7 +109,9 @@ class PhaseSpaceMapping : public Mapping { TChannelMode t_channel_mode = propagator, const std::optional& cuts = std::nullopt, const std::vector>& permutations = {}, - const std::optional>& color_order = std::nullopt + const std::optional>& color_order = std::nullopt, + double beam_rapidity = 0., + bool mirror_beams = false ); /** @@ -109,6 +126,18 @@ class PhaseSpaceMapping : public Mapping { * @ref TChannelMode::rambo). * @param cuts Fiducial cuts. See @ref Cuts. * @param color_order Color ordering for @ref TChannelMode::color_ordered. + * @param beam_rapidity Rapidity of the beams' centre-of-mass frame in the + * lab frame the momenta are returned in: + * 0.5 ln(E1 / E2) for beam energies E1 (along +z) + * and E2, with @p cm_energy = 2 sqrt(E1 E2). The + * phase space is generated in the beams' frame and + * boosted into the lab before the cuts, so rapidity + * cuts act in the lab. + * @param mirror_beams Add the `mirror_index` condition: for index 1 the + * event is rotated by pi about the x axis in the + * beams' frame (leg 1 from beam 2, leg 2 from beam 1) + * before the lab boost and the cuts. With identical + * beams that is a rotation of the lab event too. */ PhaseSpaceMapping( const std::vector& external_masses, @@ -117,7 +146,9 @@ class PhaseSpaceMapping : public Mapping { double invariant_power = 0.8, TChannelMode mode = rambo, const std::optional& cuts = std::nullopt, - const std::optional>& color_order = std::nullopt + const std::optional>& color_order = std::nullopt, + double beam_rapidity = 0., + bool mirror_beams = false ); /** @@ -126,12 +157,22 @@ class PhaseSpaceMapping : public Mapping { * * 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, giving 3n-2. + * i.e. one further sampled invariant on top of s_hat, giving 3n-2. A + * hadronic 2 -> 1 collision keeps a single random number, the rapidity of + * the produced particle, since s_hat is fixed at its mass. * * @param topology The decay topology. * @param leptonic If true, skip the PDF convolution. + * @throws std::invalid_argument for a 2 -> 1 topology without PDFs. */ 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); } @@ -145,6 +186,12 @@ class PhaseSpaceMapping : public Mapping { } /// Number of permutation channels. std::size_t channel_count() const { return _permutations.size(); } + /// Rapidity of the beams' centre-of-mass frame in the lab frame. + double beam_rapidity() const { return _beam_rapidity; } + /// Whether the mapping takes the `mirror_index` condition. + bool mirror_beams() const { return _mirror_beams; } + /// The cuts, applied to the lab-frame momenta. + const Cuts& cuts() const { return _cuts; } /// True when the cuts exclude the whole region this mapping samples: a /// pair mass cut puts a propagator's floor at or above the upper end of /// its on-shell window. Every point then fails the cuts, so the channel @@ -162,15 +209,35 @@ class PhaseSpaceMapping : public Mapping { const NamedVector& inputs, const NamedVector& 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& conditions + ) const; + Value from_lab( + FunctionBuilder& fb, + Value momenta, + Value x1, + Value x2, + const NamedVector& conditions + ) const; Topology _topology; Cuts _cuts; double _pi_factors; double _sqrt_s_lab; + double _beam_rapidity; + bool _mirror_beams; bool _leptonic; bool _map_luminosity; - // bound on the absolute rapidity of the partonic system, log(x1 / x2) / 2, - // implied by the cuts; negative when they imply none + bool _two_to_one; + // bound on the absolute rapidity of the partonic system in the beams' + // centre-of-mass frame, log(x1 / x2) / 2, implied by the cuts; negative + // when they imply none double _y_max_lab = -1.; bool _empty = false; std::size_t _n_discrete; diff --git a/madspace/src/driver/systematics.cpp b/madspace/src/driver/systematics.cpp index 9021fe154..9fd22e05e 100644 --- a/madspace/src/driver/systematics.cpp +++ b/madspace/src/driver/systematics.cpp @@ -173,7 +173,8 @@ SystematicsCalculator::SystematicsCalculator( const std::optional& nominal_alpha_s, ContextPtr context, const std::vector>& matrix_elements, - const nested_vector2& me_flavor_remap + const nested_vector2& me_flavor_remap, + const std::optional& 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 dyn_scale; }; std::vector 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)); @@ -697,20 +722,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& x, + const std::vector& q, + const std::vector& slots, + const std::vector& on_second) { + if (&first == &second) { + return evaluate_pdf(first, x, q, slots); + } + std::array, 2> grid_x, grid_q; + std::array, 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, 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 values(x.size()); + std::array 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 x, q; std::vector slots; + std::vector 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) { @@ -728,9 +790,11 @@ void SystematicsCalculator::compute( std::vector> r_pdf(count, std::vector(var_count, 1.)); if (use_pdf) { auto evaluate_variations = [&](const PdfEvaluator& evaluator, + const PdfEvaluator& evaluator2, const std::vector& var_indices) { std::vector x, q; std::vector slots; + std::vector on_second; for (auto& in : inputs) { for (std::size_t k : var_indices) { auto& var = _variations.at(k); @@ -738,15 +802,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.; @@ -778,11 +845,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/cuts.cpp b/madspace/src/phasespace/cuts.cpp index 69795e581..c5a6ae77c 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& args) c return {{"mask", fb.product(weights)}}; } +std::vector Cuts::non_mirror_invariant_cuts() const { + std::vector 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; +} + namespace { // The accessors below hand bounds to the phase-space mappings, which apply each diff --git a/madspace/src/phasespace/integrand.cpp b/madspace/src/phasespace/integrand.cpp index 379c3bbda..a97e8493e 100644 --- a/madspace/src/phasespace/integrand.cpp +++ b/madspace/src/phasespace/integrand.cpp @@ -2,6 +2,7 @@ #include "madspace/util.hpp" +#include #include using namespace madspace; @@ -62,7 +63,8 @@ Integrand::Integrand( const std::vector& flavor_diff_xs_indices, const std::vector& flavor_subproc_indices, const std::vector& flavor_per_subproc_remap, - std::size_t compressed_channel_weight_count + std::size_t compressed_channel_weight_count, + const std::optional& pdf_grid2 ) : FunctionGenerator( "Integrand", @@ -210,7 +212,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 +255,45 @@ 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()) { + 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(); } @@ -301,7 +354,7 @@ NamedVector 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 +480,20 @@ NamedVector Integrand::build_channel_part( chan_index_in_group = fb.full({static_cast(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(2)); + mirror_index = index; + } else { + mirror_index = fb.full({static_cast(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 +530,18 @@ NamedVector 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 +623,16 @@ NamedVector 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 +680,7 @@ NamedVector 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); } @@ -654,6 +741,11 @@ NamedVector Integrand::build_common_part( // 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; @@ -919,7 +1011,7 @@ NamedVector 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/observable.cpp b/madspace/src/phasespace/observable.cpp index dc4bdbc87..6049ca181 100644 --- a/madspace/src/phasespace/observable.cpp +++ b/madspace/src/phasespace/observable.cpp @@ -34,6 +34,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, @@ -256,6 +284,17 @@ std::tuple, nested_vector2, Type> build_indic } // namespace +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 option_mirror_invariant(_observable) && + (!_order_observable || option_mirror_invariant(_order_observable.value())); +} + const std::vector Observable::jet_pids{1, 2, 3, 4, -1, -2, -3, -4, 21}; const std::vector Observable::bottom_pids{-5, 5}; const std::vector Observable::lepton_pids{11, 13, 15, -11, -13, -15}; diff --git a/madspace/src/phasespace/phasespace.cpp b/madspace/src/phasespace/phasespace.cpp index 41202b82b..6a51a3650 100644 --- a/madspace/src/phasespace/phasespace.cpp +++ b/madspace/src/phasespace/phasespace.cpp @@ -258,7 +258,9 @@ PhaseSpaceMapping::PhaseSpaceMapping( TChannelMode t_channel_mode, const std::optional& cuts, const std::vector>& permutations, - const std::optional>& color_order + const std::optional>& color_order, + double beam_rapidity, + bool mirror_beams ) : Mapping( "PhaseSpaceMapping", @@ -286,9 +288,16 @@ PhaseSpaceMapping::PhaseSpaceMapping( )}, {"x1", batch_float}, {"x2", batch_float}}, - permutations.size() > 1 - ? NamedVector{{"permutation_index", batch_int}} - : NamedVector{} + [&] { + NamedVector 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( @@ -298,6 +307,8 @@ PhaseSpaceMapping::PhaseSpaceMapping( std::pow(2 * PI, 4 - 3 * static_cast(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 @@ -307,7 +318,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; @@ -438,6 +459,10 @@ PhaseSpaceMapping::PhaseSpaceMapping( // the rapidity of a sum of momenta, tanh(y) = sum_i m_T,i sinh(y_i) / // sum_i m_T,i cosh(y_i), is a weighted mean of the tanh(y_i), so it never // exceeds the largest |y_i|. This holds for massive particles too. + // The eta cuts hold in the lab frame, and the bounds are read in the + // beams' centre-of-mass frame, where the momenta are generated: + // |y_cm| <= |eta_lab| + |beam_rapidity|. The mirror only flips the sign of + // y there, so the bound holds for both orientations. std::vector node_y_max(node_leaves.size(), -1.); for (std::size_t d = 0; d < node_leaves.size(); ++d) { if (node_leaves.at(d).empty()) { @@ -450,7 +475,7 @@ PhaseSpaceMapping::PhaseSpaceMapping( y_max = -1.; break; } - y_max = std::max(y_max, eta_max); + y_max = std::max(y_max, eta_max + std::abs(beam_rapidity)); } node_y_max.at(d) = y_max; } @@ -676,7 +701,9 @@ PhaseSpaceMapping::PhaseSpaceMapping( std::vector eta_max, pt_min, y_max; 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); // rapidity bound of the child, composites included (node_y_max) y_max.push_back(node_y_max.at(index)); @@ -726,8 +753,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)); + } } } // The blocks of the chain that scatter the two beams (central 2->2, @@ -771,12 +804,26 @@ PhaseSpaceMapping::PhaseSpaceMapping( double invariant_power, TChannelMode mode, const std::optional& cuts, - const std::optional>& color_order + const std::optional>& 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 vertices; auto n_out = external_masses.size() - 2; @@ -810,7 +857,9 @@ PhaseSpaceMapping::PhaseSpaceMapping( mode, cuts, {}, - color_order + color_order, + beam_rapidity, + mirror_beams ) {} Mapping::Result PhaseSpaceMapping::build_forward_impl( @@ -834,6 +883,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(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 decay_data( _topology.decays().begin(), _topology.decays().end() @@ -995,26 +1077,68 @@ 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(0)); +Value PhaseSpaceMapping::to_lab( + FunctionBuilder& fb, + Value momenta, + Value x1, + Value x2, + const NamedVector& 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. + Value masses(lab_masses(_topology, _permutations)); + 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, masses, fb.mul(x1, exp_plus), fb.mul(x2, exp_minus) + ); + } + if (_map_luminosity) { + momenta = fb.boost_beam(momenta, masses, x1, x2); + } + if (_mirror_beams) { + // a rotation: the masses of the legs are unchanged by it + momenta = fb.mirror_momenta( + momenta, conditions.at(_permutations.size() > 1 ? 1 : 0) + ); + } + if (_beam_rapidity != 0.) { + momenta = fb.boost_beam(momenta, masses, 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, Value(lab_masses(_topology, _permutations)), 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& conditions +) const { + Value masses(lab_masses(_topology, _permutations)); + 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, masses, fb.mul(x1, exp_plus), fb.mul(x2, exp_minus) + ); + } + if (_beam_rapidity != 0.) { + momenta = fb.boost_beam_inverse(momenta, masses, exp_plus, exp_minus); + } + if (_mirror_beams) { + // a rotation: the masses of the legs are unchanged by it + momenta = fb.mirror_momenta( + momenta, conditions.at(_permutations.size() > 1 ? 1 : 0) + ); + } + if (_map_luminosity) { + momenta = fb.boost_beam_inverse(momenta, masses, x1, x2); + } + return momenta; } Mapping::Result PhaseSpaceMapping::build_inverse_impl( @@ -1023,11 +1147,11 @@ Mapping::Result PhaseSpaceMapping::build_inverse_impl( const NamedVector& 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, Value(lab_masses(_topology, _permutations)), 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 f5a3b3c6a..b1216c983 100644 --- a/madspace/src/phasespace/topology.cpp +++ b/madspace/src/phasespace/topology.cpp @@ -246,9 +246,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" ); } @@ -491,7 +497,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 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 b1fede50d..132e569a8 100644 --- a/madspace/src/python/madspace.cpp +++ b/madspace/src/python/madspace.cpp @@ -949,6 +949,11 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("name") = "", pydoc::doc("Observable::Observable") ) + .def( + "mirror_invariant", + &Observable::mirror_invariant, + pydoc::doc("Observable::mirror_invariant") + ) .def_readonly_static("jet_pids", &Observable::jet_pids) .def_readonly_static("bottom_pids", &Observable::bottom_pids) .def_readonly_static("lepton_pids", &Observable::lepton_pids) @@ -1000,6 +1005,16 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("particle_count"), pydoc::doc("Cuts::Cuts#2") ) + .def( + "non_mirror_invariant_cuts", + &Cuts::non_mirror_invariant_cuts, + pydoc::doc("Cuts::non_mirror_invariant_cuts") + ) + .def( + "mirror_invariant", + &Cuts::mirror_invariant, + pydoc::doc("Cuts::mirror_invariant") + ) .def("sqrt_s_min", &Cuts::sqrt_s_min, pydoc::doc("Cuts::sqrt_s_min")) .def("eta_max", &Cuts::eta_max, pydoc::doc("Cuts::eta_max")) .def("pt_min", &Cuts::pt_min, pydoc::doc("Cuts::pt_min")) @@ -1264,7 +1279,9 @@ PYBIND11_MODULE(_madspace_py, m) { PhaseSpaceMapping::TChannelMode, const std::optional&, const nested_vector2&, - const std::optional>&>(), + const std::optional>&, + double, + bool>(), py::arg("topology"), py::arg("cm_energy"), py::arg("leptonic") = false, @@ -1273,6 +1290,8 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("cuts") = std::nullopt, py::arg("permutations") = nested_vector2{}, py::arg("color_order") = std::nullopt, + py::arg("beam_rapidity") = 0., + py::arg("mirror_beams") = false, pydoc::doc("PhaseSpaceMapping::PhaseSpaceMapping") ) .def( @@ -1283,7 +1302,9 @@ PYBIND11_MODULE(_madspace_py, m) { double, PhaseSpaceMapping::TChannelMode, std::optional, - const std::optional>&>(), + const std::optional>&, + double, + bool>(), py::arg("external_masses"), py::arg("cm_energy"), py::arg("leptonic") = false, @@ -1291,6 +1312,8 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("mode") = PhaseSpaceMapping::rambo, py::arg("cuts") = std::nullopt, py::arg("color_order") = std::nullopt, + py::arg("beam_rapidity") = 0., + py::arg("mirror_beams") = false, pydoc::doc("PhaseSpaceMapping::PhaseSpaceMapping#2") ) .def( @@ -1313,6 +1336,17 @@ PYBIND11_MODULE(_madspace_py, m) { &PhaseSpaceMapping::channel_count, pydoc::doc("PhaseSpaceMapping::channel_count") ) + .def( + "beam_rapidity", + &PhaseSpaceMapping::beam_rapidity, + pydoc::doc("PhaseSpaceMapping::beam_rapidity") + ) + .def( + "mirror_beams", + &PhaseSpaceMapping::mirror_beams, + pydoc::doc("PhaseSpaceMapping::mirror_beams") + ) + .def("cuts", &PhaseSpaceMapping::cuts, pydoc::doc("PhaseSpaceMapping::cuts")) .def("empty", &PhaseSpaceMapping::empty, pydoc::doc("PhaseSpaceMapping::empty")); py::classh( @@ -2063,7 +2097,8 @@ PYBIND11_MODULE(_madspace_py, m) { const std::vector&, const std::vector&, const std::vector&, - std::size_t>(), + std::size_t, + const std::optional&>(), py::arg("mapping"), py::arg("diff_xs"), py::arg("adaptive_map") = std::monostate{}, @@ -2092,6 +2127,7 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("flavor_subproc_indices") = std::vector{}, py::arg("flavor_per_subproc_remap") = std::vector{}, py::arg("compressed_channel_weight_count") = 50, + py::arg("pdf_grid2") = std::nullopt, pydoc::doc("Integrand::Integrand") ) .def( @@ -3307,7 +3343,8 @@ PYBIND11_MODULE(_madspace_py, m) { const std::optional&, ContextPtr, const std::vector>&, - const nested_vector2&>(), + const nested_vector2&, + const std::optional&>(), py::arg("config"), py::arg("subproc_args"), py::arg("nominal_pdf") = std::nullopt, @@ -3315,6 +3352,7 @@ PYBIND11_MODULE(_madspace_py, m) { py::arg("context") = nullptr, py::arg("matrix_elements") = std::vector>{}, py::arg("me_flavor_remap") = nested_vector2{}, + py::arg("nominal_pdf2") = std::nullopt, pydoc::doc("SystematicsCalculator::SystematicsCalculator") ) .def_property_readonly( diff --git a/madspace/tests/test_beam_frame.py b/madspace/tests/test_beam_frame.py new file mode 100644 index 000000000..5a0ce1656 --- /dev/null +++ b/madspace/tests/test_beam_frame.py @@ -0,0 +1,496 @@ +"""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) + + +def s_channel_two_body(): + return ms.Topology( + ms.Diagram( + [0.0, 0.0], + [0.0, 0.0], + [ms.Propagator(0.0, 0.0)], + [["i0", "i1", "p0"], ["p0", "o0", "o1"]], + ) + ) + + +def t_channel_two_body(): + return ms.Topology( + ms.Diagram( + [0.0, 0.0], + [0.0, 0.0], + [ms.Propagator(0.0, 0.0)], + [["i0", "o0", "p0"], ["p0", "i1", "o1"]], + ) + ) + + +@pytest.mark.parametrize("mirror_index", [None, 0, 1], ids=["plain", "mirror0", "mirror1"]) +@pytest.mark.parametrize("topology", [s_channel_two_body, t_channel_two_body]) +def test_rapidity_windows_hold_in_the_lab_frame(rng, topology, mirror_index): + """The pt and |eta| cuts narrow the sampled Y = log(x1 / x2) / 2 and the + angles of the first decay or scattering (test_cut_rapidity_windows). Those + are generated in the beams' frame while the cuts act in the lab, so every + event that passes the cuts in the lab must still lie inside the windows: + it maps back into the unit cube.""" + 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] + leptons = [[11, -11]] + cuts = ms.Cuts( + [ + ms.CutItem(observable=ms.Observable(pids, "pt", leptons), min=10.0), + ms.CutItem(observable=ms.Observable(pids, "eta_abs", leptons), max=2.5), + ] + ) + kwargs = dict(beam_rapidity=y0, mirror_beams=mirror_index is not None) + free = ms.PhaseSpaceMapping(topology(), e_cm, **kwargs) + windowed = ms.PhaseSpaceMapping(topology(), e_cm, cuts=cuts, **kwargs) + n = 200_000 + conditions = ( + [] if mirror_index is None else [np.full(n, mirror_index, dtype=np.int32)] + ) + r = rng.random((n, free.random_dim())) + p, x1, x2, det = map(np.asarray, free.map_forward([r], conditions)) + x1, x2 = np.broadcast_to(x1, (n,)), np.broadcast_to(x2, (n,)) + keep = (det > 0) & (np.asarray(cuts(p)).reshape(-1) > 0.5) + # boosted far enough that the lab cut is not the beams'-frame one + assert (np.abs(0.5 * np.log(x1 / x2))[keep] > 2.5).any() + conditions = [c[keep] for c in conditions] + r_back, _ = windowed.map_inverse([p[keep], x1[keep], x2[keep]], conditions) + r_back = np.asarray(r_back) + assert np.all(np.isfinite(r_back)) + assert r_back.min() > -1e-6 and r_back.max() < 1 + 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, + ) + r = rng.random((N, cut.random_dim())) + accepted, momenta = {}, {} + for index in (0, 1): + condition = np.full(N, index, dtype=np.int32) + # the momenta of the cut mapping itself: an |eta| cut also narrows the + # sampled rapidity and angle windows, so a mapping without it would not + # give the same point for the same random numbers + p, _, _, det = map(np.asarray, cut.map_forward([r], [condition])) + momenta[index], accepted[index] = p, det > 0 + + assert accepted[0].any() + # orientation 1 is the pi rotation of orientation 0 + assert momenta[1] == approx(mirror(momenta[0]), abs=1e-7) + # and both are cut on the momenta that come out, not on the other + # orientation (every rambo point of a massless 2 -> 2 is physical) + 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] == passed).all() + # which matters only if the cut can tell the two orientations apart + differ = (accepted[0] != accepted[1]).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) diff --git a/tests/acceptance_tests/test_check_xsec_processes_mg7.py b/tests/acceptance_tests/test_check_xsec_processes_mg7.py index ac64f8781..c604e1d99 100644 --- a/tests/acceptance_tests/test_check_xsec_processes_mg7.py +++ b/tests/acceptance_tests/test_check_xsec_processes_mg7.py @@ -175,7 +175,11 @@ def _edit_run_card(toml_path, events, disable_jet_cuts, cuts=None, 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 diff --git a/tests/unit_tests/interface/test_mg7_launch.py b/tests/unit_tests/interface/test_mg7_launch.py index 9eb5e49aa..9b5180e2e 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') @@ -1141,12 +1161,12 @@ def test_legacy_systematics_forces_lhe_only_without_native(self): def test_scan_card_survives(self): """A scan value is not a valid RunCardMG7 entry: the card is edited as text, and only on its output_format line.""" - self.edit('e_cm = 13000.0', 'e_cm = "scan:[13000, 14000]"') + self.edit('ebeam1 = 6500.0', 'ebeam1 = "scan:[6500, 7000]"') self.launch.force_lhe_output_if_needed({}) self.launch.force_lhe_output_if_needed({'madspin': 'ON'}) self.assertEqual(self.output_format(), 'lhe') with open(self.card) as f: - self.assertIn('e_cm = "scan:[13000, 14000]"', f.read()) + self.assertIn('ebeam1 = "scan:[6500, 7000]"', f.read()) def test_missing_output_format_line(self): with open(self.card, 'w') as f: diff --git a/tests/unit_tests/various/test_banner.py b/tests/unit_tests/various/test_banner.py index 74d358993..65b73fecf 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 @@ -1622,6 +1623,87 @@ 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_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() + 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() @@ -1641,7 +1723,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 @@ -1705,6 +1789,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['beam']['scale_factor'], 0.5) @@ -1730,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() @@ -1882,7 +1981,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_npy') - 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})