diff --git a/Changes.md b/Changes.md index 0c3d748..fd9edeb 100644 --- a/Changes.md +++ b/Changes.md @@ -1,191 +1,197 @@ # Changelog This is the log for changes to the HEJ program. Further changes to the HEJ API are documented in [`Changes-API.md`](Changes-API.md). If you are using HEJ as a library, please also read the changes there. ## In Development (To be moved) * Added HEJ@NLO truncation option to yaml config. ## Version 2.2 This release adds support for the new processes: * Wp+Wp with jets. * Wm+Wm with jets. This release includes minor changes which affect users of HEJ as a library (see [`Changes-API.md`](Changes-API.md)). Compilation now requires a compiler supporting C++17. ### 2.2.0 #### New Processes * Resummation for WpWp/WmWm with jets, including interference between configurations. * Significantly improved description of Higgs boson plus jet production: - Processes with extremal Higgs boson emission are now treated as leading-log. - Resummation is now enabled for Higgs boson production with a single jet. * `HEJFOG` can generate multiple jets together with a charged lepton-antilepton pair via a virtual Z boson or photon. +#### Updates to configuration file + +* With the new option `off-shell tolerance` HEJ repairs the momenta + of massless particles that are slightly off-shell and incoming + particles with small but non-vanishing transverse momenta. + ## Version 2.1 This release adds support for two new processes: * W boson with jets. * Jets with a charged lepton-antilepton pair via a virtual Z boson or photon. In addition, the complete set of first subleading processes (unordered gluon, central and extremal quark-antiquark pair) is implemented for pure jets and W + jets, see [arXiv:2012.10310](https://arxiv.org/abs/2012.10310). Unordered gluon emission is also supported for Higgs boson + jets and Z boson/photon + jets. This release include many changes to the code, which affect users of HEJ as a library (see [`Changes-API.md`](Changes-API.md)). ### 2.1.3 * Updated documentation. ### 2.1.2 * Updated `cxxopts.hpp` dependency. ### 2.1.1 * Fixed invalid iterator accesses. * Reorganised automated tests. * Updated documentation. ### 2.1.0 #### New Processes * Resummation for W bosons with jets - New subleading processes `extremal qqbar` & `central qqbar` for a quark and anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other subleading processes also work with W's) - `HEJFOG` can generate multiple jets together with a (off-shell) W bosons decaying into lepton & neutrino * Resummation for jets with a charged lepton-antilepton pair via a virtual Z boson or photon. Includes the `unordered` subleading process. * Resummation can now be performed on all subleading processes within pure jets also. This includes `unordered`, `extremal qqbar` and `central qqbar` processes. #### More Physics implementation * Partons now have a Colour charge - Colours are read from and written to LHE files - For reweighted events the colours are created according to leading colour in the FKL limit * Use relative fraction for soft transverse momentum in tagging jets (`soft pt regulator`) as new (optional) parameter. - This supersedes `min extparton pt`, which is marked **deprecated** and will be removed in version 2.2.0 - This is a direct replacement for the old `max ext soft pt fraction`, which is also **deprecated**. * Changed the redistribution of the momenta of soft emissions. Now also bosons can take part of the recoil (previously only jets). #### Updates to Runcard * Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2` * Grouped `event treatment` for subleading channels together in runcard - Rename `non-HEJ` processes to `non-resummable` * Read electro-weak constants from input - new mandatory setting `vev` to change vacuum expectation value - new mandatory settings `particle properties` to specify mass & width of bosons - `HEJFOG`: decays are now specified in `decays` setting (previously under `particle properties`) * Allow loading multiple analyses with `analyses`. The old `analysis` (with "i") is marked **deprecated**. * Optional setting to specify maximal number of Fixed Order events (`max events`, default is all) * Allow changing the regulator lambda in input (`regulator parameter`, only for advanced users) #### Changes to Input/Output * Added support to read & write `hdf5` event files suggested in [arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs [HighFive](https://github.com/BlueBrain/HighFive)) * Support input with average weight equal to the cross section (`IDWTUP=1 or 4`) * Support unmodified Les Houches Event Files written by Sherpa with `cross section = sum(weights)/sum(trials)` * Analyses now get general run information (`LHEF::HEPRUP`) in the constructor. **This might break previously written, external analyses!** - external analyses should now be created with `make_analysis(YAML::Node const & config, LHEF::HEPRUP const & heprup)` * Support `rivet` version 3 with both `HepMC` version 2 and 3 - Multiple weights with `rivet 3` will only create one `.yoda` file (instead of one per weight/scale) * Added option to unweight only resummation events (`unweight: {type: resummation}`) * Added option for partially unweighting resummation events, similar to the fixed-order generator. * Improved unweighting algorithm. * Follow HepMC convention for particle Status codes: incoming = 11, decaying = 2, outgoing = 1 (unchanged) #### Miscellaneous * Print cross sections at end of run * Added example analysis & scale to `examples/`. Everything in `examples/` will be build when the flag `-DBUILD_EXAMPLES=TRUE` is set in `cmake`. * Dropped support for HepMC 3.0.0, either HepMC version 2 or >3.1 is required - It is now possible to write out both HepMC 2 and HepMC 3 events at the same time * Require LHADPF version 6. Dropped support for all other versions. * Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`) * Currents are now generated with [`FORM`](https://github.com/vermaseren/form) - `FORM` is included as a `git submodule`, use `git submodule update --init` to download `FORM` * Create [Sphinx](http://sphinx-doc.org/) and [Doxygen](http://doxygen.org/) documentation by `make sphinx` or `make doxygen` in your `build/` folder ## Version 2.0 First release of HEJ 2. Complete code rewrite compared to HEJ 1. Improved matching to Fixed Order ([arXiv:1805.04446](https://arxiv.org/abs/1805.04446)). Implemented processes: Higgs boson with jets (FKL and unordered gluon emission, with finite quark mass loop, [arXiv:1812.08072](https://arxiv.org/abs/1812.08072)), and pure jets (only FKL). See [arXiv:1902.08430](https://arxiv.org/abs/1902.08430) ## 2.0.7 * Added missing `#include` directives. ### 2.0.6 * Fixed compiling rivet when YODA headers are _outside_ of rivet directory. ### 2.0.5 * Fixed event classification for input not ordered in rapidity. ### 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`. ### 2.0.3 * Fixed parsing of (numerical factor) * (base scale) in configuration. * Don't change scale names, but sanitise Rivet output file names instead. ### 2.0.2 * Changed scale names to `"_over_"` and `"_times_"` for proper file names (was `"/"` and `"*"` before). ### 2.0.1 * Fixed name of fixed-order generator in error message. ### 2.0.0 * First release. diff --git a/config.yml b/config.yml index 410b203..83cf2e3 100644 --- a/config.yml +++ b/config.yml @@ -1,134 +1,137 @@ ## Number of attempted resummation phase space points for each input event trials: 10 resummation jets: # resummation jet properties min pt: 30 # minimum jet transverse momentum algorithm: antikt # jet clustering algorithm R: 0.4 # jet R parameter fixed order jets: # properties of input jets min pt: 20 # by default, algorithm and R are like for resummation jets ## Treatment of he various event classes ## the supported settings are: reweight, keep, discard ## non-resummable events cannot be reweighted event treatment: FKL: reweight unordered: keep extremal qqbar: keep central qqbar: keep non-resummable: keep ## Central scale choice or choices # ## multiple scales are allowed, e.g. # scales: [125, max jet pperp, H_T/2, 2*jet invariant mass, m_j1j2] scales: 91.188 ## Factors by which the central scales should be multiplied ## renormalisation and factorisation scales are varied independently # # scale factors: [0.5, 0.7071, 1, 1.41421, 2] ## Maximum ratio between renormalisation and factorisation scale # # max scale ratio: 2.0001 ## Import scale setting functions # # import scales: # lib_my_scales.so: [scale0,scale1] ## Unweighting setting ## remove to obtain weighted events # unweight: # # type of unweighting (one of 'weighted', 'resummation', 'partial') # type: partial # trials: 10000 # max deviation: 0 ## Event output files # # the supported formats are # - Les Houches (suffix .lhe) # - HepMC2 (suffix .hepmc2) # - HepMC3 (suffix .hepmc3 or .hepmc) # - HDF5 (suffix .hdf5) # ## An output file's format is deduced either automatically from the suffix ## or from an explicit specification, e.g. ## - Les Houches: outfile # event output: - HEJ.lhe # - HEJ_events.hepmc ## Analyses # # analyses: ## Rivet analysis # - rivet: MC_XS # rivet analysis name # output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added ## Custom analysis # - plugin: /path/to/libmyanalysis.so # my analysis parameter: some value ## Selection of random number generator and seed ## The choices are ## - mixmax (seed is an integer) ## - ranlux64 (seed is a filename containing parameters) random generator: name: mixmax # seed: 1 ## Whether or not to include higher order logs log correction: false ## Truncate higher-order corrections at NLO NLO truncation: enabled: false # nlo order: 2 ## Vacuum expectation value vev: 246.2196508 ## Properties of the weak gauge bosons particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 ## Parameters for Higgs-gluon couplings ## This requires compilation with QCDloop # # Higgs coupling: # use impact factors: false # mt: 174 # include bottom: true # mb: 4.7 +# # Tolerance towards numerical inaccuracies in input momenta +# off-shell tolerance: 0.01 + ## ---------------------------------------------------------------------- ## ## The following settings are only intended for advanced users. ## ## Please DO NOT SET them unless you know exactly what you are doing! ## ## ---------------------------------------------------------------------- ## # ## Maximum soft transverse momentum fraction in any tagging jets, e.g. ## extremal or qqbar jet # soft pt regulator: 0.1 # ## Minimum transverse momentum of extremal partons ## deprecated: use "soft pt regulator" instead # min extparton pt: 30 # ## deprecated: this cot directly replaced by "soft pt regulator" # max ext soft pt fraction: 0.1 # # max events: -1 # Maximal number of fixed order Events to process # regulator parameter: 0.2 # The regulator lambda for the subtraction terms diff --git a/doc/sphinx/HEJ.rst b/doc/sphinx/HEJ.rst index 3198fe7..7dbd804 100644 --- a/doc/sphinx/HEJ.rst +++ b/doc/sphinx/HEJ.rst @@ -1,463 +1,472 @@ .. _`Running HEJ 2`: Running HEJ 2 ============= Quick start ----------- In order to run HEJ 2, you need a configuration file and a file containing fixed-order events. A sample configuration is given by the :file:`config.yml` file distributed together with HEJ 2. Events in the Les Houches Event File format can be generated with standard Monte Carlo generators like `MadGraph5_aMC@NLO `_ or `Sherpa `_. If HEJ 2 was compiled with `HDF5 `_ support, it can also read and write event files in the format suggested in `arXiv:1905.05120 `_. HEJ 2 assumes that the cross section is given by the sum of the event weights. Depending on the fixed-order generator it may be necessary to adjust the weights in the Les Houches Event File accordingly. The processes supported by HEJ 2 are - Pure multijet production - Production of a Higgs boson with jets - Production of a W boson with jets - Production of jets with a charged lepton-antilepton pair, via a virtual Z boson and/or photon - Production of two same-sign W bosons with jets where at least two jets are required in each case. For the time being, only leading-order input events are supported. After generating an event file :file:`events.lhe` adjust the parameters under the `fixed order jets`_ setting in :file:`config.yml` to the settings in the fixed-order generation. Resummation can then be added by running:: HEJ config.yml events.lhe Using the default settings, this will produce an output event file :file:`HEJ.lhe` with events including high-energy resummation. When using the `Docker image `_, HEJ can be run with .. code-block:: bash docker run -v $PWD:$PWD -w $PWD hejdock/hej HEJ config.yml events.lhe .. _`HEJ 2 settings`: Settings -------- HEJ 2 configuration files follow the `YAML `_ format. The following configuration parameters are supported: .. _`trials`: **trials** High-energy resummation is performed by generating a number of resummation phase space configurations corresponding to an input fixed-order event. This parameter specifies how many such configurations HEJ 2 should try to generate for each input event. Typical values vary between 10 and 100. .. _`fixed order jets`: **fixed order jets** This tag collects a number of settings specifying the jet definition in the event input. The settings should correspond to the ones used in the fixed-order Monte Carlo that generated the input events. .. _`fixed order jets: min pt`: **min pt** Minimum transverse momentum in GeV of fixed-order jets. .. _`fixed order jets: algorithm`: **algorithm** The algorithm used to define jets. Allowed settings are :code:`kt`, :code:`cambridge`, :code:`antikt`, :code:`cambridge for passive`. See the `FastJet `_ documentation for a description of these algorithms. .. _`fixed order jets: R`: **R** The R parameter used in the jet algorithm, roughly corresponding to the jet radius in the plane spanned by the rapidity and the azimuthal angle. .. _`resummation jets`: **resummation jets** This tag collects a number of settings specifying the jet definition in the observed, i.e. resummed events. These settings are optional, by default the same values as for the `fixed order jets`_ are assumed. .. _`resummation jets: min pt`: **min pt** Minimum transverse momentum in GeV of resummation jets. This should be between 25% and 50% larger than the minimum transverse momentum of fixed order jets set by `fixed order jets: min pt`_. .. _`resummation jets: algorithm`: **algorithm** The algorithm used to define jets. The HEJ 2 approach to resummation relies on properties of :code:`antikt` jets, so this value is strongly recommended. For a list of possible other values, see the `fixed order jets: algorithm`_ setting. .. _`resummation jets: R`: **R** The R parameter used in the jet algorithm. .. _`event treatment`: **event treatment** Specify how to treat different event types. The different event types contribute to different orders in the high-energy limit. The possible values are :code:`reweight` to enable resummation, :code:`keep` to keep the events as they are up to a possible change of renormalisation and factorisation scale, and :code:`discard` to discard these events. The following types are implemented for the different bosons: .. csv-table:: :header: , "FKL", "unordered", "extremal qqbar", "central qqbar" :widths: auto :align: center :stub-columns: 1 "pure jets", "Yes", "Yes", "Yes", "Yes" "Higgs + jets", "Yes", "Yes", "No", "No" "W + jets", "Yes", "Yes", "Yes", "Yes" "Z/γ + jets", "Yes", "Yes", "No", "No" "same-sign W + jets", "Yes", "No", "No", "No" Non-implemented process will always be classified as :code:`non-resummable`. The different types are: .. _`FKL`: **FKL** Specifies how to treat events respecting FKL rapidity ordering, where all but the two partons extremal in rapidity have to be gluons, e.g. :code:`u d => u g d`. These configurations are dominant in the high-energy limit. .. _`unordered`: **unordered** Specifies how to treat events with one gluon emission that does not respect FKL ordering, e.g. :code:`u d => g u d`. In the high-energy limit, such configurations are logarithmically suppressed compared to FKL configurations. .. _`extremal qqbar`: **extremal qqbar** Specifies how to treat events with a quark-antiquark pair as extremal partons in rapidity, e.g. :code:`g d => u u_bar d`. In the high-energy limit, such configurations are logarithmically suppressed compared to FKL configurations. .. _`central qqbar`: **central qqbar** Specifies how to treat events with a quark-antiquark pair central in rapidity, e.g. :code:`g g => g u u_bar g`. In the high-energy limit, such configurations are logarithmically suppressed compared to FKL configurations. .. _`non-resummable`: **non-resummable** Specifies how to treat events that do not fall into any of the above categories or that are not yet implemented. Only :code:`keep` or :code:`discard` are valid options, *not* :code:`reweight` for obvious reasons. .. _`scales`: **scales** Specifies the renormalisation and factorisation scales for the output events. This can either be a single entry or a list :code:`[scale1, scale2, ...]`. For the case of a list the first entry defines the central scale. Possible values are fixed numbers to set the scale in GeV or the following: - :code:`H_T`: The sum of the scalar transverse momenta of all final-state particles - :code:`max jet pperp`: The maximum transverse momentum of all jets - :code:`jet invariant mass`: Sum of the invariant masses of all jets - :code:`m_j1j2`: Invariant mass between the two hardest jets. Scales can be multiplied or divided by overall factors, e.g. :code:`H_T/2`. It is also possible to import scales from an external library, see :ref:`Custom scales` .. _`scale factors`: **scale factors** A list of numeric factors by which each of the `scales`_ should be multiplied. Renormalisation and factorisation scales are varied independently. For example, a list with entries :code:`[0.5, 2]` would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`); (0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\ :sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to the order of the final event weights. .. _`max scale ratio`: **max scale ratio** Specifies the maximum factor by which renormalisation and factorisation scales may difer. For a value of :code:`2` and the example given for the `scale factors`_ the scale choices (0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`) will be discarded. .. _`log correction`: **log correction** Whether to include corrections due to the evolution of the strong coupling constant in the virtual corrections. Allowed values are :code:`true` and :code:`false`. .. _`NLO truncation`: **NLO truncation** Options to truncate the HEJ resummation at next-to-leading order. Used for bin-by-bin NLO reweighting. .. _`NLO truncation: enabled`: **enabled** Enable truncation. Allowed values are :code:`true` and :code:`false` (default). .. _`NLO truncation: nlo order`: **nlo order** Set the (base) number of jets in the NLO sample. Allowed values are integers (default: 2). .. _`unweight`: **unweight** Settings for unweighting events. Unweighting can greatly reduce the number of resummation events, speeding up analyses and shrinking event file sizes. .. _`type`: **type** How to unweight events. The supported settings are - :code:`weighted`: Generate weighted events. Default, if nothing else specified. - :code:`resummation`: Unweight only resummation events. Each set of resummation events coming from *a single fixed order event* are unweighted separately according to the largest weight in the current chunk of events. - :code:`partial`: Unweight only resummation events with weights below a certain threshold. The weight threshold is determined automatically in a calibration run prior to the usual event generation. .. _`unweight: trials`: **trials** Maximum number of trial resummation events generated in the calibration run for partial unweighting. This option should only be set for partial unweighting. If possible, each trial is generated from a different input fixed-order event. If there are not sufficiently many input events, more than one trial event may be generated for each of them and the actual number of trial events may be smaller than requested. Increasing the number of trials generally leads to better unweighting calibration but increases the run time. Between 1000 and 10000 trials are usually sufficient. .. _`unweight: max deviation`: **max deviation** Controls the range of events to which unweighting is applied. This option should only be set for partial unweighting. A larger value means that a larger fraction of events are unweighted. Typical values are between -1 and 1. .. _`event output`: **event output** Specifies the name of a single event output file or a list of such files. The file format is either specified explicitly or derived from the suffix. For example, :code:`events.lhe` or, equivalently :code:`Les Houches: events.lhe` generates an output event file :code:`events.lhe` in the Les Houches format. The supported formats are - :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches event file format. - :code:`file.hepmc2` or :code:`HepMC2: file`: HepMC format version 2. - :code:`file.hepmc3` or :code:`HepMC3: file`: HepMC format version 3. - :code:`file.hepmc` or :code:`HepMC: file`: The latest supported version of the HepMC format, currently version 3. - :code:`file.hdf5` or :code:`HDF5: file`: The HDF5-based format of `arXiv:1905.05120 `_. .. _`random generator`: **random generator** Sets parameters for random number generation. .. _`random generator: name`: **name** Which random number generator to use. Currently, :code:`mixmax` and :code:`ranlux64` are supported. Mixmax is recommended. See the `CLHEP documentation `_ for details on the generators. .. _`random generator: seed`: **seed** The seed for random generation. This should be a single number for :code:`mixmax` and the name of a state file for :code:`ranlux64`. .. _`analyses`: **analyses** Names and settings for one or more custom and Rivet event analyses. Entries containing the :code:`rivet` key are interpreted as Rivet analyses; the values corresponding to this key should be the analyses names. In addition, there is a mandatory :code:`output` key which determines the prefix for the yoda output file. For a custom analysis the :code:`plugin` sub-entry should be set to the analysis file path. All further entries are passed on to the analysis. See :ref:`Writing custom analyses` for details. .. _`vev`: **vev** Higgs vacuum expectation value in GeV. All electro-weak constants are derived from this together with the `particle properties`_. .. _`particle properties`: **particle properties** Specifies various properties of the different particles (Higgs, W or Z). All electro-weak constants are derived from these together with the :ref:`vacuum expectation value`. .. _`particle properties: particle`: **Higgs, W or Z** The particle (Higgs, |W+| or |W-|, Z) for which the following properties are defined. .. |W+| replace:: W\ :sup:`+` .. |W-| replace:: W\ :sup:`-` .. _`particle properties: particle: mass`: **mass** The mass of the particle in GeV. .. _`particle properties: particle: width`: **width** The total decay width of the particle in GeV. .. _`Higgs coupling`: **Higgs coupling** This collects a number of settings concerning the effective coupling of the Higgs boson to gluons. This is only relevant for the production process of a Higgs boson with jets and only supported if HEJ 2 was compiled with `QCDLoop `_ support. .. _`Higgs coupling: use impact factors`: **use impact factors** Whether to use impact factors for the coupling to the most forward and most backward partons. Impact factors imply the infinite top-quark mass limit. .. _`Higgs coupling: mt`: **mt** The value of the top-quark mass in GeV. If this is not specified, the limit of an infinite mass is taken. .. _`Higgs coupling: include bottom`: **include bottom** Whether to include the Higgs coupling to bottom quarks. .. _`Higgs coupling: mb`: **mb** The value of the bottom-quark mass in GeV. Only used for the Higgs coupling, external bottom-quarks are always assumed to be massless. +.. _`off-shell tolerance`: + +**off-shell tolerance** + Tolerance for numerical inaccuracies in input momenta. Momenta of + massless particles with an invariant mass below the given value + are rescaled to be on-shell. Transverse momentum components of + incoming particles that are smaller than the given tolerance are + set to zero. The default value is 0, leaving input momenta unchanged. + Advanced Settings ~~~~~~~~~~~~~~~~~ All of the following settings are optional. Please **do not set** any of the following options, unless you know exactly what you are doing. The default behaviour gives the most reliable results for a wide range of observables. .. _`soft pt regulator`: **soft pt regulator** Specifies the maximum fraction that soft radiation can contribute to the transverse momentum of each the tagging jets, i.e. any jet that affects the event classification, like the most forward and most backward jet or the jets of the central qqbar pair. This setting is needed to regulate an otherwise cancelled divergence. Default is 0.1 .. _`max ext soft pt fraction`: **max ext soft pt fraction** This is the same as `soft pt regulator`_ and will be removed in future versions. .. _`min extparton pt`: **min extparton pt** Specifies the minimum transverse momentum in GeV of the most forward and the most backward parton. Its value should be slightly below the minimum transverse momentum of jets specified by `resummation jets: min pt`_. This setting got superseded by `soft pt regulator`_ and will be removed in future versions. .. _`max events`: **max events** Maximal number of (input) Fixed Order events. HEJ will stop after processing `max events` many events. Default considers all events. .. _`regulator parameter`: **regulator parameter** Slicing parameter to regularise the subtraction term, called :math:`\lambda` in `arxiv:1706.01002 `_. Default is 0.2 diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh index d9b6ad9..73e1bf9 100644 --- a/include/HEJ/Config.hh +++ b/include/HEJ/Config.hh @@ -1,284 +1,286 @@ /** \file * \brief HEJ 2 configuration parameters * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include "fastjet/JetDefinition.hh" #include "yaml-cpp/yaml.h" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/output_formats.hh" namespace HEJ { //! Jet parameters struct JetParameters{ fastjet::JetDefinition def; /**< Jet Definition */ double min_pt{}; /**< Minimum Jet Transverse Momentum */ }; //! Settings for scale variation struct ScaleConfig{ //! Base scale choices std::vector base; //! Factors for multiplicative scale variation std::vector factors; //! Maximum ratio between renormalisation and factorisation scale double max_ratio{}; }; //! Settings for random number generator struct RNGConfig { //! Random number generator name std::string name; //! Optional initial seed std::optional seed; }; //! Settings for partial unweighting struct PartialUnweightConfig { //! Number of trials for training size_t trials; //! Maximum distance in standard deviations from mean logarithmic weight double max_dev; }; //! Settings for HEJ@NLO struct NLOConfig { //! Settings for HEJ@NLO Truncation bool enabled = false; //! NLO Born number of jets size_t nj = 2; }; /**! Possible treatments for fixed-order input events. * * The program will decide on how to treat an event based on * the value of this enumeration. */ enum class EventTreatment{ reweight, /**< Perform resummation */ keep, /**< Keep the event */ discard, /**< Discard the event */ }; //! Container to store the treatments for various event types using EventTreatMap = std::map; //! Possible setting for the event weight enum class WeightType{ weighted, //!< weighted events unweighted_resum, //!< unweighted only resummation part partially_unweighted //!< mixed weighted and unweighted }; /**! Input parameters. * * This struct handles stores all configuration parameters * needed in a HEJ 2 run. * * \internal To add a new option: * 1. Add a member to the Config struct. * 2. Inside "src/YAMLreader.cc": * - Add the option name to the "supported" Node in * get_supported_options. * - Initialise the new Config member in to_Config. * The functions set_from_yaml (for mandatory options) and * set_from_yaml_if_defined (non-mandatory) may be helpful. * 3. Add a new entry (with short description) to config.yaml * 4. Update the user documentation in "doc/Sphinx/" */ struct Config { //! %Parameters for scale variation ScaleConfig scales; //! Resummation jet properties JetParameters resummation_jets; //! Fixed-order jet properties JetParameters fixed_order_jets; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. //! Use \ref soft_pt_regulator instead. double min_extparton_pt = 0.; //! \deprecated This is equivalent to\ref soft_pt_regulator //! and will be removed in future versions. std::optional> max_ext_soft_pt_fraction{}; //! @brief Maximum transverse momentum fraction from soft radiation in any //! tagging jet (e.g. extremal or qqbar jet) Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; //! Number of resummation configurations to generate per fixed-order event size_t trials{}; //! Maximal number of events std::optional max_events; //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Event output files names and formats std::vector output; //! Parameters for random number generation RNGConfig rng; //! Map to decide what to do for different event types EventTreatMap treat; //! %Parameters for custom analysis //! @deprecated use analyses_parameters instead YAML::Node analysis_parameters; //! %Parameters for custom analyses std::vector analyses_parameters; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; //! Type of event weight e.g. (un)weighted WeightType weight_type; //! Settings for partial unweighting std::optional unweight_config; //! HEJ@NLO settings NLOConfig nlo; + //! Tolerance towards numerical inaccuracies in input momenta + double off_shell_tolerance = 0.; }; //! Configuration options for the PhaseSpacePoint class struct PhaseSpacePointConfig { PhaseSpacePointConfig() = default; PhaseSpacePointConfig( JetParameters jet_param, NLOConfig nlo, double min_extparton_pt = 0., Fraction soft_pt_regulator = Fraction{DEFAULT_SOFT_PT_REGULATOR} ): jet_param{std::move(jet_param)}, nlo{std::move(nlo)}, min_extparton_pt{min_extparton_pt}, soft_pt_regulator{std::move(soft_pt_regulator)} {} //! Properties of resummation jets JetParameters jet_param; //! HEJ@NLO settings NLOConfig nlo; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. //! Use \ref soft_pt_regulator instead. double min_extparton_pt = 0.; //! \deprecated This is equivalent to\ref soft_pt_regulator //! and will be removed in future versions. std::optional> max_ext_soft_pt_fraction{}; //! @brief Maximum transverse momentum fraction from soft radiation in any //! tagging jet (e.g. extremal or qqbar jet) Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; }; //! Configuration options for the MatrixElement class struct MatrixElementConfig { MatrixElementConfig() = default; MatrixElementConfig( bool log_correction, HiggsCouplingSettings Higgs_coupling, EWConstants ew_parameters, NLOConfig nlo, Fraction soft_pt_regulator = Fraction{DEFAULT_SOFT_PT_REGULATOR}, double regulator_lambda = CLAMBDA ): log_correction{log_correction}, Higgs_coupling{std::move(Higgs_coupling)}, ew_parameters{std::move(ew_parameters)}, nlo{std::move(nlo)}, soft_pt_regulator{soft_pt_regulator}, regulator_lambda{regulator_lambda} {} //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; //! HEJ@NLO settings NLOConfig nlo; //! @brief Maximum transverse momentum fraction from soft radiation in any //! tagging jet (e.g. extremal or qqbar jet) Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; }; //! Configuration options for the EventReweighter class struct EventReweighterConfig { //! Settings for phase space point generation PhaseSpacePointConfig psp_config; //! Settings for matrix element calculation MatrixElementConfig ME_config; //! Access properties of resummation jets JetParameters & jet_param() { return psp_config.jet_param;} //! Access properties of resummation jets (const version) JetParameters const & jet_param() const { return psp_config.jet_param;} //! Treatment of the various event types EventTreatMap treat; }; /**! Extract PhaseSpacePointConfig from Config * * \internal We do not provide a PhaseSpacePointConfig constructor from Config * so that PhaseSpacePointConfig remains an aggregate. * This faciliates writing client code (e.g. the HEJ fixed-order generator) * that creates a PhaseSpacePointConfig *without* a Config object. * * @see to_MatrixElementConfig, to_EventReweighterConfig */ inline PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) { return { conf.resummation_jets, conf.nlo, conf.min_extparton_pt, conf.max_ext_soft_pt_fraction?*conf.max_ext_soft_pt_fraction :conf.soft_pt_regulator }; } /**! Extract MatrixElementConfig from Config * * @see to_PhaseSpacePointConfig, to_EventReweighterConfig */ inline MatrixElementConfig to_MatrixElementConfig(Config const & conf) { return { conf.log_correction, conf.Higgs_coupling, conf.ew_parameters, conf.nlo, conf.soft_pt_regulator, conf.regulator_lambda }; } /**! Extract EventReweighterConfig from Config * * @see to_PhaseSpacePointConfig, to_MatrixElementConfig */ inline EventReweighterConfig to_EventReweighterConfig(Config const & conf) { return { to_PhaseSpacePointConfig(conf), to_MatrixElementConfig(conf), conf.treat }; } } // namespace HEJ diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh index c9cd6dd..dae110a 100644 --- a/include/HEJ/Event.hh +++ b/include/HEJ/Event.hh @@ -1,403 +1,413 @@ /** \file * \brief Declares the Event class and helpers * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include #include #include #include "boost/iterator/filter_iterator.hpp" #include "fastjet/ClusterSequence.hh" #include "fastjet/PseudoJet.hh" #include "HEJ/Constants.hh" #include "HEJ/Parameters.hh" #include "HEJ/Particle.hh" #include "HEJ/event_types.hh" namespace LHEF { class HEPEUP; class HEPRUP; } namespace fastjet { class JetDefinition; } namespace HEJ { class EWConstants; struct RNG; struct UnclusteredEvent; /** @brief An event with clustered jets * * This is the main HEJ 2 event class. * It contains kinematic information including jet clustering, * parameter (e.g. scale) settings and the event weight. */ class Event { public: class EventData; //! Iterator over partons using ConstPartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::const_iterator >; //! Reverse Iterator over partons using ConstReversePartonIterator = std::reverse_iterator< ConstPartonIterator>; //! No default Constructor Event() = delete; //! Event Constructor adding jet clustering to an unclustered event //! @deprecated UnclusteredEvent got superseded by EventData //! and will be removed in HEJ 2.2.0 [[deprecated("UnclusteredEvent got superseded by EventData")]] Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! @name Particle Access //! @{ //! Incoming particles std::array const & incoming() const{ return incoming_; } //! Outgoing particles std::vector const & outgoing() const{ return outgoing_; } //! Iterator to the first outgoing parton ConstPartonIterator begin_partons() const; //! Iterator to the first outgoing parton ConstPartonIterator cbegin_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator end_partons() const; //! Iterator to the end of the outgoing partons ConstPartonIterator cend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crbegin_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator rend_partons() const; //! Reverse Iterator to the first outgoing parton ConstReversePartonIterator crend_partons() const; //! Particle decays /** * The key in the returned map corresponds to the index in the * vector returned by outgoing() */ std::unordered_map> const & decays() const { return decays_; } //! The jets formed by the outgoing partons, sorted in rapidity std::vector const & jets() const{ return jets_; } //! @} //! @name Weight variations //! @{ //! All chosen parameter, i.e. scale choices (const version) Parameters const & parameters() const{ return parameters_; } //! All chosen parameter, i.e. scale choices Parameters & parameters(){ return parameters_; } //! Central parameter choice (const version) EventParameters const & central() const{ return parameters_.central; } //! Central parameter choice EventParameters & central(){ return parameters_.central; } //! Parameter (scale) variations (const version) std::vector const & variations() const{ return parameters_.variations; } //! Parameter (scale) variations std::vector & variations(){ return parameters_.variations; } //! Parameter (scale) variation (const version) /** * @param i Index of the requested variation */ EventParameters const & variations(std::size_t i) const{ return parameters_.variations.at(i); } //! Parameter (scale) variation /** * @param i Index of the requested variation */ EventParameters & variations(std::size_t i){ return parameters_.variations.at(i); } //! @} //! Indices of the jets the outgoing partons belong to /** * @param jets Jets to be tested * @returns A vector containing, for each outgoing parton, * the index in the vector of jets the considered parton * belongs to. If the parton is not inside any of the * passed jets, the corresponding index is set to -1. */ std::vector particle_jet_indices( std::vector const & jets ) const { return cs_.particle_jet_indices(jets); } //! particle_jet_indices() of the Event jets() std::vector particle_jet_indices() const { return particle_jet_indices(jets()); } //! Jet definition used for clustering fastjet::JetDefinition const & jet_def() const{ return cs_.jet_def(); } //! Minimum jet transverse momentum double min_jet_pt() const{ return min_jet_pt_; } //! Event type event_type::EventType type() const{ return type_; } //! Give colours to each particle /** * @returns true if new colours are generated, i.e. same as is_resummable() * @details Colour ordering is done according to leading colour in the MRK * limit, see \cite Andersen:2011zd. This only affects \ref * is_resummable() "HEJ" configurations, all other \ref event_type * "EventTypes" will be ignored. * @note This overwrites all previously set colours. */ bool generate_colours(RNG & /*ran*/); //! Check that current colours are leading in the high energy limit /** * @details Checks that the colour configuration can be split up in * multiple, rapidity ordered, non-overlapping ladders. Such * configurations are leading in the MRK limit, see * \cite Andersen:2011zd * * @note This is _not_ to be confused with \ref is_resummable(), however * for all resummable states it is possible to create a leading colour * configuration, see generate_colours() */ bool is_leading_colour() const; /** * @brief Check if given event could have been produced by HEJ * @details A HEJ state has to fulfil: * 1. type() has to be \ref is_resummable() "resummable" * 2. Soft radiation in the tagging jets contributes at most to * `soft_pt_regulator` of the total jet \f$ p_\perp \f$ * * @note This is true for any resummed stated produced by the * EventReweighter or any \ref is_resummable() "resummable" Leading * Order state. * * @param soft_pt_regulator Maximum transverse momentum fraction from soft * radiation in tagging jets * @param min_pt Absolute minimal \f$ p_\perp \f$, * \b deprecated use soft_pt_regulator instead * @return True if this state could have been produced by HEJ */ bool valid_hej_state( double soft_pt_regulator = DEFAULT_SOFT_PT_REGULATOR, double min_pt = 0. ) const; //! Check that the incoming momenta are valid /** * @details Checks that the incoming parton momenta are on-shell and have * vanishing transverse components. * */ bool valid_incoming() const; private: //! \internal //! @brief Construct Event explicitly from input. /** This is only intended to be called from EventData. * * \warning The input is taken _as is_, sorting and classification has to be * done externally, i.e. by EventData */ Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double min_jet_pt ); //! Iterator over partons (non-const) using PartonIterator = boost::filter_iterator< bool (*)(Particle const &), std::vector::iterator >; //! Reverse Iterator over partons (non-const) using ReversePartonIterator = std::reverse_iterator; //! Iterator to the first outgoing parton (non-const) PartonIterator begin_partons(); //! Iterator to the end of the outgoing partons (non-const) PartonIterator end_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rbegin_partons(); //! Reverse Iterator to the first outgoing parton (non-const) ReversePartonIterator rend_partons(); std::array incoming_; std::vector outgoing_; std::unordered_map> decays_; std::vector jets_; Parameters parameters_; fastjet::ClusterSequence cs_; double min_jet_pt_; event_type::EventType type_; }; // end class Event //! Detect if a backward incoming gluon turns into a backward outgoing Higgs boson inline bool is_backward_g_to_h(Event const & ev) { return ev.outgoing().front().type == pid::higgs && ev.incoming().front().type == pid::gluon; } //! Detect if a forward incoming gluon turns into a forward outgoing Higgs boson inline bool is_forward_g_to_h(Event const & ev) { return ev.outgoing().back().type == pid::higgs && ev.incoming().back().type == pid::gluon; } //! Class to store general Event setup, i.e. Phase space and weights class Event::EventData { public: //! Default Constructor EventData() = default; //! Constructor from LesHouches event information EventData(LHEF::HEPEUP const & hepeup); //! Constructor with all values given EventData( std::array incoming, std::vector outgoing, std::unordered_map> decays, Parameters parameters ): incoming(std::move(incoming)), outgoing(std::move(outgoing)), decays(std::move(decays)), parameters(std::move(parameters)) {} //! Generate an Event from the stored EventData. /** * @details Do jet clustering and classification. * Use this to generate an Event. * * @note Calling this function destroys EventData * * @param jet_def Jet definition * @param min_jet_pt minimal \f$p_T\f$ for each jet * * @returns Full clustered and classified event. */ Event cluster( fastjet::JetDefinition const & jet_def, double min_jet_pt); //! Alias for cluster() Event operator()( fastjet::JetDefinition const & jet_def, double const min_jet_pt){ return cluster(jet_def, min_jet_pt); } //! Sort particles in rapidity void sort(); //! Reconstruct intermediate particles from final-state leptons /** * Final-state leptons are created from virtual photons, W, or Z bosons. * This function tries to reconstruct such intermediate bosons if they * are not part of the event record. */ void reconstruct_intermediate(EWConstants const & /*ew_parameters*/); + //! Repair momenta of massless particles + /** + * This function changes the momenta of massless particles as follows: + * - Close-to-zero incoming transverse momenta are set to zero. + * - Nearly on-shell momenta are made lightlike by rescaling energy + * and spatial components. This rescaling is applied to both incoming + * and outgoing particles, including decay products. + */ + void repair_momenta(double tolerance); + //! Incoming particles std::array incoming; //! Outcoing particles std::vector outgoing; //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Parameters, e.g. scale or inital weight Parameters parameters; }; // end class EventData //! Print Event std::ostream& operator<<(std::ostream & os, Event const & ev); //! Square of the partonic centre-of-mass energy \f$\hat{s}\f$ double shat(Event const & ev); //! Tolerance parameter for validity check on incoming momenta static constexpr double TOL = 1e-6; //! Convert an event to a LHEF::HEPEUP LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * /*heprup*/); // put deprecated warning at the end, so don't get the warning inside Event.hh, // additionally doxygen can not identify [[deprecated]] correctly struct [[deprecated("UnclusteredEvent will be replaced by EventData")]] UnclusteredEvent; //! An event before jet clustering //! @deprecated UnclusteredEvent got superseded by EventData //! and will be removed in HEJ 2.2.0 struct UnclusteredEvent{ //! Default Constructor UnclusteredEvent() = default; //! Constructor from LesHouches event information UnclusteredEvent(LHEF::HEPEUP const & hepeup); std::array incoming; /**< Incoming Particles */ std::vector outgoing; /**< Outgoing Particles */ //! Particle decays in the format {outgoing index, decay products} std::unordered_map> decays; //! Central parameter (e.g. scale) choice EventParameters central; std::vector variations; /**< For parameter variation */ }; } // namespace HEJ diff --git a/include/HEJ/PDG_codes.hh b/include/HEJ/PDG_codes.hh index 23694de..f520144 100644 --- a/include/HEJ/PDG_codes.hh +++ b/include/HEJ/PDG_codes.hh @@ -1,241 +1,268 @@ /** \file PDG_codes.hh * \brief Contains the Particle IDs of all relevant SM particles. * * Large enumeration included which has multiple entries for potential * alternative names of different particles. There are also functions * which can be used to determine if a particle is a parton or if * it is a non-gluon boson. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include namespace HEJ { //! particle ids according to PDG namespace pid { //! The possible particle identities. We use PDG IDs as standard. enum ParticleID: int{ //! Unspecified type, should never be used!, debug only unspecified = 0, d = 1, /*!< Down Quark */ down = d, /*!< Down Quark */ u = 2, /*!< Up Quark */ up = u, /*!< Up Quark */ s = 3, /*!< Strange Quark */ strange = s, /*!< Strange Quark */ c = 4, /*!< Charm Quark */ charm = c, /*!< Charm Quark */ b = 5, /*!< Bottom Quark */ bottom = b, /*!< Bottom Quark */ t = 6, /*!< Top Quark */ top = t, /*!< Top Quark */ e = 11, /*!< Electron */ electron = e, /*!< Electron */ nu_e = 12, /*!< Electron Neutrino */ electron_neutrino = nu_e, /*!< Electron neutrino */ mu = 13, /*!< Muon */ muon = mu, /*!< Muon */ nu_mu = 14, /*!< Muon Neutrino */ muon_neutrino = nu_mu, /*!< Muon Neutrino */ tau = 15, /*!< Tau */ nu_tau = 16, /*!< Tau Neutrino */ tau_neutrino = nu_tau, /*!< Tau Neutrino */ d_bar = -d, /*!< Anti-Down Quark */ antidown = d_bar, /*!< Anti-Down Quark */ u_bar = -u, /*!< Anti-Up quark */ antiup = -u, /*!< Anti-Up quark */ s_bar = -s, /*!< Anti-Strange Quark */ antistrange = -s, /*!< Anti-Strange Quark */ c_bar = -c, /*!< Anti-Charm Quark */ anticharm = -c, /*!< Anti-Charm Quark */ b_bar = -b, /*!< Anti-Bottom Quark */ antibottom = -b, /*!< Anti-Bottom Quark */ t_bar = -t, /*!< Anti-Top Quark */ antitop = -t, /*!< Anti-Top Quark */ e_bar = -e, /*!< Positron */ positron = e_bar, /*!< Positron */ antielectron = positron, /*!< Positron */ nu_e_bar = -nu_e, /*!< Electron Anti-Neutrino */ electron_antineutrino = nu_e_bar,/*!< Electron Anti-Neutrino */ mu_bar = -mu, /*!< Anti-Muon */ antimuon = -mu, /*!< Anti-Muon */ nu_mu_bar = -nu_mu, /*!< Muon Anti-Neutrino */ muon_antineutrino = nu_mu_bar, /*!< Muon Anti-Neutrino */ tau_bar = -tau, /*!< Anti-Tau */ antitau = tau_bar, /*!< Anti-Tau */ nu_tau_bar = -nu_tau, /*!< Tau Anti-Neutrino */ tau_antineutrino = nu_tau_bar, /*!< Tau Anti-Neutrino */ gluon = 21, /*!< Gluon */ g = gluon, /*!< Gluon */ photon = 22, /*!< Photon */ gamma = photon, /*!< Photon */ Z = 23, /*!< Z Boson */ Z_photon_mix = 81, /*!< Z/photon superposition */ Z_gamma_mix = Z_photon_mix, /*!< Z/photon superposition */ Wp = 24, /*!< W- Boson */ Wm = -Wp, /*!< W+ Boson */ h = 25, /*!< Higgs Boson */ Higgs = h, /*!< Higgs Boson */ higgs = h, /*!< Higgs Boson */ p = 2212, /*!< Proton */ proton = p, /*!< Proton */ p_bar = -p, /*!< Anti-Proton */ antiproton = p_bar, /*!< Anti-Proton */ }; //! Get the of the particle with the given PDG ID std::string name(ParticleID id); //! return the negative flavour of the given PDG ID ParticleID anti(ParticleID id); } // namespace pid using ParticleID = pid::ParticleID; //! Convert a particle name to the corresponding PDG particle ID ParticleID to_ParticleID(std::string const & name); /** * \brief Function to determine if particle is a quark * @param id PDG ID of particle * @returns true if the particle is a quark, false otherwise */ inline constexpr bool is_quark(ParticleID id){ return (id >= pid::down && id <= pid::top); } /** * \brief Function to determine if particle is an antiquark * @param id PDG ID of particle * @returns true if the particle is an antiquark, false otherwise */ inline constexpr bool is_antiquark(ParticleID id){ return (id <= pid::d_bar && id >= pid::t_bar); } /** * \brief Function to determine if particle is an (anti-)quark * @param id PDG ID of particle * @returns true if the particle is a quark or antiquark, false otherwise */ inline constexpr bool is_anyquark(ParticleID id){ return is_quark(id) || is_antiquark(id); } /** * \brief Function to determine if particle is a gluon * @param id PDG ID of particle * @returns true if the particle is a gluon, false otherwise */ inline constexpr bool is_gluon(ParticleID id){ return id == pid::gluon; } /** * \brief Function to determine if particle is a parton * @param id PDG ID of particle * @returns true if the particle is a parton, false otherwise */ inline constexpr bool is_parton(ParticleID id){ return is_gluon(id) || (is_anyquark(id) && std::abs(id) != pid::top); } /** * \brief function to determine if the particle is a photon, W or Z * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZ_boson(ParticleID id){ return id == pid::Wm || (id >= pid::photon && id <= pid::Wp) || id == pid::Z_photon_mix; } /** * \brief function to determine if the particle is a photon, W, Z, or Higgs * boson * @param id PDG ID of particle * @returns true if the partice is an A,W,Z, or H, false otherwise */ inline constexpr bool is_AWZH_boson(ParticleID id){ return is_AWZ_boson(id) || (id == pid::Higgs); } /** * \brief Function to determine if particle is a lepton * @param id PDG ID of particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(ParticleID id){ return (id >= pid::electron && id <= pid::tau_neutrino); } /** * \brief Function to determine if particle is an antilepton * @param id PDG ID of particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(ParticleID id){ return (id <= pid::positron && id >= pid::nu_tau_bar); } /** * \brief Function to determine if particle is an (anti-)lepton * @param id PDG ID of particle * @returns true if the particle is a lepton or antilepton, * false otherwise */ inline constexpr bool is_anylepton(ParticleID id){ return ( is_lepton(id) || is_antilepton(id)); } /** * \brief Function to determine if particle is a neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(ParticleID id){ return (id == pid::nu_e || id == pid::tau_neutrino || id == pid::muon_neutrino); } /** * \brief Function to determine if particle is an antineutrino * @param id PDG ID of particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(ParticleID id){ return (id == pid::nu_e_bar || id == pid::nu_tau_bar || id == pid::nu_mu_bar); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param id PDG ID of particle * @returns true if the particle is a neutrino or antineutrino, * false otherwise */ inline constexpr bool is_anyneutrino(ParticleID id){ return ( is_neutrino(id) || is_antineutrino(id)); } + + //! Check if a particle is massless + inline + constexpr bool is_massless(ParticleID id){ + // cannot use `std::abs` because it's not `constexpr` + const int abs_id = (id >= 0)?id:-id; + switch(abs_id){ + case pid::bottom: + case pid::top: + case pid::tau: + case pid::Z: + case pid::Z_photon_mix: + case pid::Wp: + case pid::Higgs: + case pid::proton: + return false; + default: + return true; + } + } + + //! Check if a particle is massive + inline + constexpr bool is_massive(ParticleID id){ + return !is_massless(id); + } + } // namespace HEJ diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh index 94f0b4e..62e70da 100644 --- a/include/HEJ/Particle.hh +++ b/include/HEJ/Particle.hh @@ -1,239 +1,251 @@ /** * \file Particle.hh * \brief Contains the particle struct * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/PDG_codes.hh" namespace HEJ { using Colour = std::pair; //! Class representing a particle struct Particle { //! particle type ParticleID type = pid::unspecified; //! particle momentum fastjet::PseudoJet p; //! (optional) colour & anti-colour std::optional colour; //! get rapidity double rapidity() const{ return p.rapidity(); } //! get transverse momentum double perp() const{ return p.perp(); } //! get transverse momentum double pt() const{ return perp(); } //! get momentum in x direction double px() const{ return p.px(); } //! get momentum in y direction double py() const{ return p.py(); } //! get momentum in z direction double pz() const{ return p.pz(); } //! get energy double E() const{ return p.E(); } //! get mass double m() const{ return p.m(); } }; //! Functor to compare rapidities /** * This can be used whenever a rapidity comparison function is needed, * for example in many standard library functions. * * @see pz_less */ struct rapidity_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.rapidity() < p2.rapidity(); } }; //! Functor to compare momenta in z direction /** * This can be used whenever a pz comparison function is needed, * for example in many standard library functions. * * @see rapidity_less */ struct pz_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.pz() < p2.pz(); } }; //! Convert a vector of Particles to a vector of particle momenta inline std::vector to_PseudoJet( std::vector const & v ){ std::vector result; result.reserve(v.size()); for(auto && sp: v) result.emplace_back(sp.p); return result; } //! Functor to compare particle type (PDG) /** * This can be used whenever a particle-type comparison is needed * for example in many standard library functions. */ struct type_less{ template bool operator()(Particle const & p1, Particle const & p2){ return p1.type < p2.type; } }; //! Check if a particle is a parton, i.e. quark, antiquark, or gluon inline constexpr bool is_parton(Particle const & p){ return is_parton(p.type); } //! Check if a particle is a quark inline constexpr bool is_quark(Particle const & p){ return is_quark(p.type); } //! Check if a particle is an anti-quark inline constexpr bool is_antiquark(Particle const & p){ return is_antiquark(p.type); } //! Check if a particle is a quark or anit-quark inline constexpr bool is_anyquark(Particle const & p){ return is_anyquark(p.type); } /** * \brief Function to determine if particle is a lepton * @param p the particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(Particle const & p){ return is_lepton(p.type); } /** * \brief Function to determine if particle is an antilepton * @param p the particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(Particle const & p){ return is_antilepton(p.type); } /** * \brief Function to determine if particle is an (anti-)lepton * @param p the particle * @returns true if the particle is a lepton or antilepton, false otherwise */ inline constexpr bool is_anylepton(Particle const & p){ return is_anylepton(p.type); } /** * \brief Function to determine if particle is a neutrino * @param p the particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(Particle const & p){ return is_neutrino(p.type); } /** * \brief Function to determine if particle is an antineutrino * @param p the particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(Particle const & p){ return is_antineutrino(p.type); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param p the particle * @returns true if the particle is a neutrino or antineutrino, false otherwise */ inline constexpr bool is_anyneutrino(Particle const & p){ return is_anyneutrino(p.type); } + //! Check if a particle is massless + inline + constexpr bool is_massless(Particle const & p){ + return is_massless(p.type); + } + + //! Check if a particle is massive + inline + constexpr bool is_massive(Particle const & p){ + return is_massive(p.type); + } + //! Check if a particle is a photon, W or Z boson inline constexpr bool is_AWZ_boson(Particle const & particle){ return is_AWZ_boson(particle.type); } //! Check if a particle is a photon, W, Z, or Higgs boson inline constexpr bool is_AWZH_boson(Particle const & particle){ return is_AWZH_boson(particle.type); } //! Extract all partons from a vector of particles inline std::vector filter_partons( std::vector const & v ){ std::vector result; result.reserve(v.size()); std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_parton(p); } ); return result; } //! Extract all AWZH bosons from a vector of particles inline std::vector filter_AWZH_bosons( std::vector const & v ){ std::vector result; std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_AWZH_boson(p); } ); return result; } } // namespace HEJ diff --git a/src/Event.cc b/src/Event.cc index cf62654..7e935ac 100644 --- a/src/Event.cc +++ b/src/Event.cc @@ -1,1361 +1,1394 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/Event.hh" #include #include #include #include #include #include #include #include #include #include #include #include "fastjet/ClusterSequence.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" #include "LHEF/LHEF.h" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/RNG.hh" #include "HEJ/exceptions.hh" +#include "HEJ/LorentzVector.hh" #include "HEJ/utility.hh" namespace HEJ { namespace { using std::size_t; //! LHE status codes namespace lhe_status { enum Status: int { in = -1, decay = 2, out = 1, }; } using LHE_Status = lhe_status::Status; //! true if leptonic W decay bool valid_W_decay( int const w_type, // sign of W std::vector const & decays ){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) if( w_type == 1 ) // W+ return is_antilepton(decays[0]) || is_neutrino(decays[0]); // W- return is_lepton(decays[0]) || is_antineutrino(decays[0]); } //! true for Z decay to charged leptons bool valid_Z_decay(std::vector const & decays){ if(decays.size() != 2) // no 1->2 decay return false; const int pidsum = decays[0].type + decays[1].type; if( std::abs(pidsum) != 0 ) // correct charge return false; // leptonic decay (only check first, second follows from pidsum) return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]); } //! true if supported decay bool valid_decay(std::vector const & decays){ return valid_W_decay(+1, decays) || // Wp valid_W_decay(-1, decays) || // Wm valid_Z_decay( decays) // Z/gamma ; } /// @name helper functions to determine event type //@{ /** * \brief check if final state valid for HEJ * * check final state has the expected number of valid decays for bosons * and all the rest are quarks or gluons */ bool final_state_ok(Event const & ev){ size_t invalid_decays = ev.decays().size(); std::vector const & outgoing = ev.outgoing(); for( size_t i=0; i second; fastjet::PseudoJet res = std::accumulate( progeny.cbegin(), progeny.cend(), fastjet::PseudoJet(), [](fastjet::PseudoJet & sum, Particle const & p) -> fastjet::PseudoJet { return std::move(sum) + p.p; } ); if(!nearby(out.p, res, out.E())){ return false; } } // W decays (required) if(std::abs(out.type) == ParticleID::Wp){ if( decay != ev.decays().cend() && valid_W_decay(out.type>0?+1:-1, decay->second) ){ --invalid_decays; } else return false; } // Higgs decays (optional) else if(out.type == ParticleID::h){ if(decay != ev.decays().cend()) --invalid_decays; } // Z decays (required) else if(out.type == ParticleID::Z_photon_mix){ if( decay != ev.decays().cend() && valid_Z_decay(decay->second) ){ --invalid_decays; } else return false; } } else if(! is_parton(out.type)) return false; } // any invalid decays? return invalid_decays == 0; } /** * returns all EventTypes implemented in HEJ */ size_t implemented_types(std::vector const & bosons){ using namespace event_type; // no bosons if(bosons.empty()) return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid; // 1 boson if(bosons.size()== 1) { switch (bosons[0].type) { case ParticleID::Wp: case ParticleID::Wm: return FKL | unob | unof | qqbar_exb | qqbar_exf | qqbar_mid; case ParticleID::Z_photon_mix: return FKL | unob | unof; case ParticleID::h: return FKL | unob | unof; default: return non_resummable; } } // 2 bosons if(bosons.size() == 2) { // Resum only samesign W events if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) { return FKL; } else if(bosons[0].type == ParticleID::Wm && bosons[1].type == ParticleID::Wm) { return FKL; } } return non_resummable; } /** * \brief function which determines if type change is consistent with Wp emission. * @param in incoming Particle id * @param out outgoing Particle id * @param is_qqbar Current both incoming/both outgoing? * * \see is_Wm_Change */ bool is_Wp_Change(ParticleID in, ParticleID out, bool is_qqbar){ using namespace pid; if(!is_qqbar && (in==d_bar || in==u || in==s_bar || in==c)) return out == (in-1); if( is_qqbar && (in==d || in==u_bar || in==s || in==c_bar)) return out == -(in+1); return false; } /** * \brief function which determines if type change is consistent with Wm emission. * @param in incoming Particle id * @param out outgoing Particle id * @param is_qqbar Current both incoming/both outgoing? * * Ensures that change type of quark line is possible by a flavour changing * Wm emission. Allows checking of is_qqbar currents also. */ bool is_Wm_Change(ParticleID in, ParticleID out, bool is_qqbar){ using namespace pid; if(!is_qqbar && (in==d || in==u_bar || in==s || in==c_bar)) return out == (in+1); if( is_qqbar && (in==d_bar || in==u || in==s_bar || in==c)) return out == -(in-1); return false; } /** * \brief checks if particle type remains same from incoming to outgoing * @param in incoming Particle * @param out outgoing Particle * @param is_qqbar Current both incoming/outgoing? */ bool no_flavour_change(ParticleID in, ParticleID out, bool is_qqbar){ const int qqbarCurrent = is_qqbar?-1:1; if(std::abs(in)<=pid::top || in==pid::gluon) return (in==out*qqbarCurrent); return false; } bool has_enough_jets(Event const & event){ if(event.jets().size() >= 2) return true; if(event.jets().empty()) return false; // check for h+jet const auto the_higgs = std::find_if( begin(event.outgoing()), end(event.outgoing()), [](const auto & particle) { return particle.type == pid::higgs; } ); return the_higgs != end(event.outgoing()); } bool is_gluon_to_Higgs(const ParticleID in, const ParticleID out) { return in == pid::gluon && out == pid::Higgs; } /** * \brief check if we have a valid Impact factor * @param in incoming Particle * @param out outgoing Particle * @param is_qqbar Current both incoming/outgoing? * @param W_change returns +1 if Wp, -1 if Wm, else 0 */ bool is_valid_impact_factor( ParticleID in, ParticleID out, bool is_qqbar, int & W_change ){ if( no_flavour_change(in, out, is_qqbar) || is_gluon_to_Higgs(in, out)) { return true; } if( is_Wp_Change(in, out, is_qqbar) ) { W_change+=1; return true; } if( is_Wm_Change(in, out, is_qqbar) ) { W_change-=1; return true; } return false; } bool is_extremal_higgs_off_quark( const ParticleID in, const ParticleID extremal_out, const ParticleID out ) { return in == out && extremal_out == pid::higgs && is_anyquark(in); } //! Returns all possible classifications from the impact factors // the beginning points are changed s.t. after the the classification they // point to the beginning of the (potential) FKL chain // sets W_change: + if Wp change // 0 if no change // - if Wm change // This function can be used with forward & backwards iterators template size_t possible_impact_factors( ParticleID incoming_id, // incoming OutIterator & begin_out, OutIterator const & end_out, // outgoing int & W_change, std::vector const & boson, bool const backward // backward? ){ using namespace event_type; if(begin_out == end_out) return non_resummable; // keep track of all states that we don't test size_t not_tested = qqbar_mid; if(backward) not_tested |= unof | qqbar_exf; else not_tested |= unob | qqbar_exb; // Is this LL current? if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){ ++begin_out; return not_tested | FKL; } // q -> H q and qbar -> H qbar are technically not LL, // but we treat them as such anyway const auto next = std::next(begin_out); if( // first ensure that the next particle is not part of the *other* impact factor next != end_out && is_extremal_higgs_off_quark(incoming_id, begin_out->type, next->type) ) { std::advance(begin_out, 2); return not_tested | FKL; } // or NLL current? // -> needs two partons in two different jets if( std::distance(begin_out, end_out)>=2 ){ auto next = std::next(begin_out); // Is this unordered emisson? if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){ if( is_valid_impact_factor( incoming_id, next->type, false, W_change ) ){ // veto Higgs inside uno assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?unob:unof); } } // Is this QQbar? else if( incoming_id==pid::gluon ){ if( is_valid_impact_factor( begin_out->type, next->type, true, W_change ) ){ // veto Higgs inside qqbar assert(next!=end_out); if( !boson.empty() && boson.front().type == ParticleID::h ){ if( (backward && boson.front().rapidity() < next->rapidity()) ||(!backward && boson.front().rapidity() > next->rapidity())) return non_resummable; } begin_out = std::next(next); return not_tested | (backward?qqbar_exb:qqbar_exf); } } } return non_resummable; } //! Returns all possible classifications from central emissions // the beginning points are changed s.t. after the the classification they // point to the end of the emission chain // sets W_change: + if Wp change // 0 if no change // - if Wm change template size_t possible_central( OutIterator & begin_out, OutIterator const & end_out, int & W_change, std::vector const & boson ){ using namespace event_type; // if we already passed the central chain, // then it is not a valid all-order state if(std::distance(begin_out, end_out) < 0) return non_resummable; // keep track of all states that we don't test size_t possible = unob | unof | qqbar_exb | qqbar_exf; // Find the first quark or antiquark emission begin_out = std::find_if( begin_out, end_out, [](Particle const & p) { return is_anyquark(p); } ); // end of chain -> FKL if( begin_out==end_out ){ return possible | FKL; } // is this a qqbar-pair? // needs two partons in two separate jets auto next = std::next(begin_out); if( next != end_out && is_valid_impact_factor(begin_out->type, next->type, true, W_change) ){ // veto Higgs inside qqbar if( !boson.empty() && boson.front().type == ParticleID::h && boson.front().rapidity() > begin_out->rapidity() && boson.front().rapidity() < next->rapidity() ){ return non_resummable; } begin_out = std::next(next); // remaining chain should be pure FKL (gluon or higgs) if(std::any_of( begin_out, end_out, [](Particle const & p) { return is_anyquark(p); } )) { return non_resummable; } return possible | qqbar_mid; } return non_resummable; } namespace { bool is_parton_or_higgs(Particle const & p) { return is_parton(p) || p.type == pid::higgs; } } /** * \brief Checks for all event types * @param ev Event * @returns Event Type * */ event_type::EventType classify(Event const & ev){ using namespace event_type; if(! has_enough_jets(ev)) return not_enough_jets; // currently we can't handle multiple boson states in the ME. So they are // considered "bad_final_state" even though the "classify" could work with // them. if(! final_state_ok(ev)) return bad_final_state; // initialise variables auto const & in = ev.incoming(); // range for current checks auto begin_out = boost::make_filter_iterator( is_parton_or_higgs, cbegin(ev.outgoing()), cend(ev.outgoing()) ); auto rbegin_out = std::make_reverse_iterator( boost::make_filter_iterator( is_parton_or_higgs, cend(ev.outgoing()), cend(ev.outgoing()) ) ); assert(std::distance(begin(in), end(in)) == 2); assert(std::distance(begin_out, rbegin_out.base()) >= 2); assert(std::is_sorted(begin_out, rbegin_out.base(), rapidity_less{})); auto const bosons{ filter_AWZH_bosons(ev.outgoing()) }; // keep track of potential W couplings, at the end the sum should be 0 int remaining_Wp = 0; int remaining_Wm = 0; for(auto const & boson : bosons){ if(boson.type == ParticleID::Wp) ++remaining_Wp; else if(boson.type == ParticleID::Wm) ++remaining_Wm; } size_t final_type = ~(not_enough_jets | bad_final_state); // check forward impact factor int W_change = 0; final_type &= possible_impact_factors( in.front().type, begin_out, rbegin_out.base(), W_change, bosons, true ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // check backward impact factor W_change = 0; final_type &= possible_impact_factors( in.back().type, rbegin_out, std::make_reverse_iterator(begin_out), W_change, bosons, false ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // check central emissions W_change = 0; final_type &= possible_central( begin_out, rbegin_out.base(), W_change, bosons ); if( final_type == non_resummable ) return non_resummable; if(W_change>0) remaining_Wp-=W_change; else if(W_change<0) remaining_Wm+=W_change; // Check whether the right number of Ws are present if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable; // result has to be unique if( (final_type & (final_type-1)) != 0) return non_resummable; // check that each sub processes is implemented // (has to be done at the end) if( (final_type & ~implemented_types(bosons)) != 0 ) return non_resummable; return static_cast(final_type); } //@} Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){ auto id = static_cast(hepeup.IDUP[i]); auto colour = is_parton(id)?hepeup.ICOLUP[i]:std::optional(); return { id, { hepeup.PUP[i][0], hepeup.PUP[i][1], hepeup.PUP[i][2], hepeup.PUP[i][3] }, colour }; } bool is_decay_product(std::pair const & mothers){ if(mothers.first == 0) return false; return mothers.second == 0 || mothers.first == mothers.second; } } // namespace Event::EventData::EventData(LHEF::HEPEUP const & hepeup){ parameters.central = EventParameters{ hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP }; size_t in_idx = 0; for (int i = 0; i < hepeup.NUP; ++i) { // skip decay products // we will add them later on, but we have to ensure that // the decayed particle is added before if(is_decay_product(hepeup.MOTHUP[i])) continue; auto particle = extract_particle(hepeup, i); // needed to identify mother particles for decay products particle.p.set_user_index(i+1); if(hepeup.ISTUP[i] == LHE_Status::in){ if(in_idx > incoming.size()) { throw std::invalid_argument{ "Event has too many incoming particles" }; } incoming[in_idx++] = std::move(particle); } else outgoing.emplace_back(std::move(particle)); } // add decay products for (int i = 0; i < hepeup.NUP; ++i) { if(!is_decay_product(hepeup.MOTHUP[i])) continue; const int mother_id = hepeup.MOTHUP[i].first; const auto mother = std::find_if( begin(outgoing), end(outgoing), [mother_id](Particle const & particle){ return particle.p.user_index() == mother_id; } ); if(mother == end(outgoing)){ throw std::invalid_argument{"invalid decay product parent"}; } const int mother_idx = std::distance(begin(outgoing), mother); assert(mother_idx >= 0); decays[mother_idx].emplace_back(extract_particle(hepeup, i)); } } Event::Event( UnclusteredEvent const & ev, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): Event( Event::EventData{ ev.incoming, ev.outgoing, ev.decays, Parameters{ev.central, ev.variations} }.cluster(jet_def, min_jet_pt) ) {} //! @TODO remove in HEJ 2.2.0 UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){ Event::EventData const evData{hepeup}; incoming = evData.incoming; outgoing = evData.outgoing; decays = evData.decays; central = evData.parameters.central; variations = evData.parameters.variations; } void Event::EventData::sort(){ // sort particles std::sort( begin(incoming), end(incoming), [](Particle const & o1, Particle const & o2){return o1.p.pz() idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){ return old_outgoing[i].rapidity() < old_outgoing[j].rapidity(); }); outgoing.clear(); outgoing.reserve(old_outgoing.size()); for(size_t i: idx) { outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!decays.empty()){ auto old_decays = std::move(decays); decays.clear(); for(size_t i=0; isecond)); } assert(old_decays.size() == decays.size()); } } namespace { // use valid_X_decay to determine boson type ParticleID reconstruct_type(std::vector const & progeny) { if(valid_W_decay(+1, progeny)) { return ParticleID::Wp; } if(valid_W_decay(-1, progeny)) { return ParticleID::Wm; } if(valid_Z_decay(progeny)) { return ParticleID::Z_photon_mix; } throw not_implemented{ "final state with decay X -> " + name(progeny[0].type) + " + " + name(progeny[1].type) }; } // reconstruct particle with explicit ParticleID Particle reconstruct_boson( std::vector const & progeny, ParticleID const & type ) { Particle progenitor; progenitor.p = progeny[0].p + progeny[1].p; progenitor.type = type; return progenitor; } // reconstruct via call to reconstruct_type Particle reconstruct_boson(std::vector const & progeny) { Particle progenitor {reconstruct_boson(progeny, reconstruct_type(progeny))}; assert(is_AWZH_boson(progenitor)); return progenitor; } using GroupedParticles = std::vector >; using Decay = std::pair >; using Decays = std::vector; // return groups of reconstructable progeny std::vector group_progeny(std::vector & leptons) { /** Warning: The partition in to charged/neutral leptons is valid ONLY for WW. **/ assert(leptons.size() == 4); auto const begin_neutrino = std::partition( begin(leptons), end(leptons), [](Particle const & p) {return !is_anyneutrino(p);} ); std::vector neutrinos (begin_neutrino, end(leptons)); leptons.erase(begin_neutrino, end(leptons)); if(leptons.size() != 2 || neutrinos.size() != 2) { return {}; } assert(leptons.size() == 2 && neutrinos.size() == 2); std::vector valid_groupings; GroupedParticles candidate_grouping{{leptons[0], neutrinos[0]}, {leptons[1], neutrinos[1]}}; if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) { valid_groupings.emplace_back(std::move(candidate_grouping)); } candidate_grouping = {{leptons[1], neutrinos[0]}, {leptons[0], neutrinos[1]}}; if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) { valid_groupings.emplace_back(std::move(candidate_grouping)); } return valid_groupings; } // 'best' decay ordering measure double decay_measure(const Particle& reconstructed, EWConstants const & params) { ParticleProperties ref = params.prop(reconstructed.type); return std::abs(reconstructed.p.m() - ref.mass); } // decay_measure accumulated over decays double decay_measure(Decays const & decays, EWConstants const & params) { return std::accumulate( cbegin(decays), cend(decays), 0., [¶ms] (double dm, Decay const & decay) -> double { return dm + decay_measure(decay.first, params); } ); } // select best combination of decays for the event Decays select_decays ( std::vector & leptons, EWConstants const & ew_parameters ) { std::vector groupings = group_progeny(leptons); std::vector valid_decays; valid_decays.reserve(groupings.size()); // Reconstruct all groupings for(GroupedParticles const & group : groupings) { Decays decays; for(auto const & progeny : group) { decays.emplace_back(make_pair(reconstruct_boson(progeny), progeny)); } valid_decays.emplace_back(decays); } if (valid_decays.empty()) { throw not_implemented{"No supported intermediate reconstruction available"}; } if (valid_decays.size() == 1) { return valid_decays[0]; } // select decay with smallest decay_measure auto selected = std::min_element(cbegin(valid_decays), cend(valid_decays), [&ew_parameters] (auto const & d1, auto const & d2) -> bool { return decay_measure(d1, ew_parameters) < decay_measure(d2, ew_parameters); } ); return *selected; } } // namespace void Event::EventData::reconstruct_intermediate(EWConstants const & ew_parameters) { auto const begin_leptons = std::partition( begin(outgoing), end(outgoing), [](Particle const & p) {return !is_anylepton(p);} ); std::vector leptons(begin_leptons, end(outgoing)); outgoing.erase(begin_leptons, end(outgoing)); if(leptons.empty()) { return; } // nothing to do if(leptons.size() == 2) { outgoing.emplace_back(reconstruct_boson(leptons)); std::sort(begin(leptons), end(leptons), type_less{}); decays.emplace(outgoing.size()-1, std::move(leptons)); } else if(leptons.size() == 4) { Decays valid_decays = select_decays(leptons, ew_parameters); for(auto &decay : valid_decays) { outgoing.emplace_back(decay.first); std::sort(begin(decay.second), end(decay.second), type_less{}); decays.emplace(outgoing.size()-1, std::move(decay.second)); } } else { throw not_implemented { std::to_string(leptons.size()) + " leptons in the final state" }; } } + namespace { + void repair_momentum(fastjet::PseudoJet & p, const double tolerance) { + if(p.e() > 0. && p.m2() != 0. && (p.m2() < tolerance * tolerance)) { + const double rescale = std::sqrt(p.modp() / p.e()); + const double e = p.e() * rescale; + const double px = p.px() / rescale; + const double py = p.py() / rescale; + const double pz = p.pz() / rescale; + p.reset(px, py, pz, e); + } + } + } + + void Event::EventData::repair_momenta(const double tolerance) { + for(auto & in: incoming) { + if(is_massless(in)) { + const double px = (std::abs(in.px()) < tolerance)?0.:in.px(); + const double py = (std::abs(in.py()) < tolerance)?0.:in.py(); + in.p.reset(px, py, in.p.pz(), in.p.e()); + repair_momentum(in.p, tolerance); + } + } + for(auto & out: outgoing) { + if(is_massless(out)) repair_momentum(out.p, tolerance); + } + for(auto & decay: decays) { + for(auto & out: decay.second) { + if(is_massless(out)) repair_momentum(out.p, tolerance); + } + } + } + Event Event::EventData::cluster( fastjet::JetDefinition const & jet_def, double const min_jet_pt ){ sort(); return Event{ std::move(incoming), std::move(outgoing), std::move(decays), std::move(parameters), jet_def, min_jet_pt }; } Event::Event( std::array && incoming, std::vector && outgoing, std::unordered_map> && decays, Parameters && parameters, fastjet::JetDefinition const & jet_def, double const min_jet_pt ): incoming_{std::move(incoming)}, outgoing_{std::move(outgoing)}, decays_{std::move(decays)}, parameters_{std::move(parameters)}, cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def }, min_jet_pt_{min_jet_pt} { jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_)); assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{})); type_ = classify(*this); } namespace { //! check that Particles have a reasonable colour bool correct_colour(Particle const & part){ ParticleID id{ part.type }; if(!is_parton(id)) return !part.colour; if(!part.colour) return false; Colour const & col{ *part.colour }; if(is_quark(id)) return col.first != 0 && col.second == 0; if(is_antiquark(id)) return col.first == 0 && col.second != 0; assert(id==ParticleID::gluon); return col.first != 0 && col.second != 0 && col.first != col.second; } //! Connect parton to t-channel colour line & update the line //! returns false if connection not possible template bool try_connect_t(OutIterator const & it_part, Colour & line_colour){ if( line_colour.first == it_part->colour->second ){ line_colour.first = it_part->colour->first; return true; } if( line_colour.second == it_part->colour->first ){ line_colour.second = it_part->colour->second; return true; } return false; } //! Connect parton to u-channel colour line & update the line //! returns false if connection not possible template bool try_connect_u(OutIterator & it_part, Colour & line_colour){ auto it_next = std::next(it_part); if( try_connect_t(it_next, line_colour) && try_connect_t(it_part, line_colour) ){ it_part=it_next; return true; } return false; } } // namespace bool Event::is_leading_colour() const { if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) ) return false; Colour line_colour = *incoming()[0].colour; std::swap(line_colour.first, line_colour.second); // reasonable colour if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour)) return false; for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){ switch (type()) { case event_type::FKL: if( !try_connect_t(it_part, line_colour) ) return false; break; case event_type::unob: case event_type::qqbar_exb: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(cbegin_partons(), it_part)!=0 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::unof: case event_type::qqbar_exf: { if( !try_connect_t(it_part, line_colour) // u-channel only allowed at impact factor && (std::distance(it_part, cend_partons())!=2 || !try_connect_u(it_part, line_colour))) return false; break; } case event_type::qqbar_mid:{ auto it_next = std::next(it_part); if( !try_connect_t(it_part, line_colour) // u-channel only allowed at q-qbar/qbar-q pair && ( ( !(is_quark(*it_part) && is_antiquark(*it_next)) && !(is_antiquark(*it_part) && is_quark(*it_next))) || !try_connect_u(it_part, line_colour)) ) return false; break; } default: throw std::logic_error{"unreachable"}; } // no colour singlet exchange/disconnected diagram if(line_colour.first == line_colour.second) return false; } return (incoming()[1].colour->first == line_colour.first) && (incoming()[1].colour->second == line_colour.second); } namespace { //! connect incoming Particle to colour flow void connect_incoming(Particle & in, int & colour, int & anti_colour){ in.colour = std::make_pair(anti_colour, colour); // gluon if(in.type == pid::gluon) return; if(in.type > 0){ // quark assert(is_quark(in)); in.colour->second = 0; colour*=-1; return; } // anti-quark assert(is_antiquark(in)); in.colour->first = 0; anti_colour*=-1; } //! connect outgoing Particle to t-channel colour flow template void connect_tchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ assert(colour>0 || anti_colour>0); if(it_part->type == ParticleID::gluon){ // gluon if(colour>0 && anti_colour>0){ // on g line => connect to colour OR anti-colour (random) if(ran.flat() < 0.5){ it_part->colour = std::make_pair(colour+2,colour); colour+=2; } else { it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(colour > 0){ // on q line => connect to available colour it_part->colour = std::make_pair(colour+2, colour); colour+=2; } else { assert(colour<0 && anti_colour>0); // on qbar line => connect to available anti-colour it_part->colour = std::make_pair(anti_colour, anti_colour+2); anti_colour+=2; } } else if(is_quark(*it_part)) { // quark assert(anti_colour>0); if(colour>0){ // on g line => connect and remove anti-colour it_part->colour = std::make_pair(anti_colour, 0); anti_colour+=2; anti_colour*=-1; } else { // on qbar line => new colour colour*=-1; it_part->colour = std::make_pair(colour, 0); } } else if(is_antiquark(*it_part)) { // anti-quark assert(colour>0); if(anti_colour>0){ // on g line => connect and remove colour it_part->colour = std::make_pair(0, colour); colour+=2; colour*=-1; } else { // on q line => new anti-colour anti_colour*=-1; it_part->colour = std::make_pair(0, anti_colour); } } else { // not a parton assert(!is_parton(*it_part)); it_part->colour = {}; } } //! connect to t- or u-channel colour flow template void connect_utchannel( OutIterator & it_part, int & colour, int & anti_colour, RNG & ran ){ OutIterator it_first = it_part++; if(ran.flat()<.5) {// t-channel connect_tchannel(it_first, colour, anti_colour, ran); connect_tchannel(it_part, colour, anti_colour, ran); } else { // u-channel connect_tchannel(it_part, colour, anti_colour, ran); connect_tchannel(it_first, colour, anti_colour, ran); } } } // namespace bool Event::generate_colours(RNG & ran){ // generate only for HEJ events if(!event_type::is_resummable(type())) return false; assert(std::is_sorted( begin(outgoing()), end(outgoing()), rapidity_less{})); assert(incoming()[0].pz() < incoming()[1].pz()); // positive (anti-)colour -> can connect // negative (anti-)colour -> not available/used up by (anti-)quark int colour = COLOUR_OFFSET; int anti_colour = colour+1; // initialise first connect_incoming(incoming_[0], colour, anti_colour); // reset outgoing colours std::for_each(outgoing_.begin(), outgoing_.end(), [](Particle & part){ part.colour = {};}); for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){ switch (type()) { // subleading can connect to t- or u-channel case event_type::unob: case event_type::qqbar_exb: { if( std::distance(begin_partons(), it_part)==0) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::unof: case event_type::qqbar_exf: { if( std::distance(it_part, end_partons())==2) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } case event_type::qqbar_mid:{ auto it_next = std::next(it_part); if( std::distance(begin_partons(), it_part)>0 && std::distance(it_part, end_partons())>2 && ( (is_quark(*it_part) && is_antiquark(*it_next)) || (is_antiquark(*it_part) && is_quark(*it_next)) ) ) connect_utchannel(it_part, colour, anti_colour, ran); else connect_tchannel(it_part, colour, anti_colour, ran); break; } default: // rest has to be t-channel connect_tchannel(it_part, colour, anti_colour, ran); } } // Connect last connect_incoming(incoming_[1], anti_colour, colour); assert(is_leading_colour()); return true; } // generate_colours namespace { bool valid_parton( std::vector const & jets, Particle const & parton, int const idx, double const soft_pt_regulator, double const min_extparton_pt ){ // TODO code overlap with PhaseSpacePoint::pass_extremal_cuts if(min_extparton_pt > parton.pt()) return false; if(idx<0) return false; assert(static_cast(jets.size())>=idx); auto const & jet{ jets[idx] }; return (parton.p - jet).pt()/jet.pt() <= soft_pt_regulator; } } // namespace // this should work with multiple types bool Event::valid_hej_state(double const soft_pt_regulator, double const min_pt ) const { using namespace event_type; if(!is_resummable(type())) return false; auto const & jet_idx{ particle_jet_indices() }; auto idx_begin{ jet_idx.cbegin() }; auto idx_end{ jet_idx.crbegin() }; auto part_begin{ cbegin_partons() }; auto part_end{ crbegin_partons() }; // always seperate extremal jets if(!is_backward_g_to_h(*this)) { if(! valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt)) { return false; } ++part_begin; ++idx_begin; // unob -> second parton in own jet if( type() & (unob | qqbar_exb) ){ if( !valid_parton(jets(), *part_begin, *idx_begin, soft_pt_regulator, min_pt) ) return false; ++part_begin; ++idx_begin; } } if(!is_forward_g_to_h(*this)) { if(!valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt)) { return false; } ++part_end; ++idx_end; if( type() & (unof | qqbar_exf) ){ if( !valid_parton(jets(), *part_end, *idx_end, soft_pt_regulator, min_pt) ) return false; ++part_end; // ++idx_end; // last check, we don't need idx_end afterwards } } if( type() & qqbar_mid ){ // find qqbar pair auto begin_qqbar{ std::find_if( part_begin, part_end.base(), [](Particle const & part) -> bool { return part.type != ParticleID::gluon; } )}; assert(begin_qqbar != part_end.base()); long int qqbar_pos{ std::distance(part_begin, begin_qqbar) }; assert(qqbar_pos >= 0); idx_begin+=qqbar_pos; if( !( valid_parton(jets(), *begin_qqbar, *idx_begin, soft_pt_regulator, min_pt) && valid_parton(jets(), *std::next(begin_qqbar), *std::next(idx_begin), soft_pt_regulator, min_pt) )) return false; } return true; } bool Event::valid_incoming() const{ for(std::size_t i=0; i < incoming_.size(); ++i){ if(!(HEJ::nearby_ep(std::abs(incoming_[i].pz()), incoming_[i].E(), TOL*incoming_[i].E()) && (incoming_[i].pt()==0.))) return false; } return true; } Event::ConstPartonIterator Event::begin_partons() const { return cbegin_partons(); } Event::ConstPartonIterator Event::cbegin_partons() const { return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())}; } Event::ConstPartonIterator Event::end_partons() const { return cend_partons(); } Event::ConstPartonIterator Event::cend_partons() const { return {HEJ::is_parton, cend(outgoing()), cend(outgoing())}; } Event::ConstReversePartonIterator Event::rbegin_partons() const { return crbegin_partons(); } Event::ConstReversePartonIterator Event::crbegin_partons() const { return std::reverse_iterator( cend_partons() ); } Event::ConstReversePartonIterator Event::rend_partons() const { return crend_partons(); } Event::ConstReversePartonIterator Event::crend_partons() const { return std::reverse_iterator( cbegin_partons() ); } Event::PartonIterator Event::begin_partons() { return {HEJ::is_parton, begin(outgoing_), end(outgoing_)}; } Event::PartonIterator Event::end_partons() { return {HEJ::is_parton, end(outgoing_), end(outgoing_)}; } Event::ReversePartonIterator Event::rbegin_partons() { return std::reverse_iterator( end_partons() ); } Event::ReversePartonIterator Event::rend_partons() { return std::reverse_iterator( begin_partons() ); } namespace { void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){ constexpr int prec = 6; const std::streamsize orig_prec = os.precision(); os < const & col){ constexpr int width = 3; if(!col) os << "(no color)"; // American spelling for better alignment else os << "(" <first << ", " <second << ")"; } } // namespace std::ostream& operator<<(std::ostream & os, Event const & ev){ constexpr int prec = 4; constexpr int wtype = 3; // width for types const std::streamsize orig_prec = os.precision(); os < rapidity=" < rapidity=" < rapidity=" < incoming{ event.incoming() }; // First incoming should be positive pz according to LHE standard // (or at least most (everyone?) do it this way, and Pythia assumes it) if(incoming[0].pz() < incoming[1].pz()) std::swap(incoming[0], incoming[1]); for(Particle const & in: incoming){ result.IDUP.emplace_back(in.type); result.ISTUP.emplace_back(LHE_Status::in); result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()}); result.MOTHUP.emplace_back(0, 0); assert(in.colour); result.ICOLUP.emplace_back(*in.colour); } // outgoing for(size_t i = 0; i < event.outgoing().size(); ++i){ Particle const & out = event.outgoing()[i]; result.IDUP.emplace_back(out.type); const int status = event.decays().count(i) != 0u ?LHE_Status::decay :LHE_Status::out; result.ISTUP.emplace_back(status); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); result.MOTHUP.emplace_back(1, 2); if(out.colour) result.ICOLUP.emplace_back(*out.colour); else{ result.ICOLUP.emplace_back(std::make_pair(0,0)); } } // decays for(auto const & decay: event.decays()){ for(auto const & out: decay.second){ result.IDUP.emplace_back(out.type); result.ISTUP.emplace_back(LHE_Status::out); result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()}); const size_t mother_idx = 1 + event.incoming().size() + decay.first; result.MOTHUP.emplace_back(mother_idx, mother_idx); result.ICOLUP.emplace_back(0,0); } } assert(result.ICOLUP.size() == num_particles); static constexpr double unknown_spin = 9.; //per Les Houches accord result.VTIMUP = std::vector(num_particles, unknown_spin); result.SPINUP = result.VTIMUP; return result; } } // namespace HEJ diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc index bf14a98..94557d2 100644 --- a/src/YAMLreader.cc +++ b/src/YAMLreader.cc @@ -1,592 +1,597 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/YAMLreader.hh" #include #include #include #include #include #include #include #include #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/output_formats.hh" namespace HEJ { class Event; namespace { //! Get YAML tree of supported options /** * The configuration file is checked against this tree of options * in assert_all_options_known. */ YAML::Node const & get_supported_options(){ const static YAML::Node supported = [](){ YAML::Node supported; static const auto opts = { "trials", "min extparton pt", "max ext soft pt fraction", "soft pt regulator", "scales", "scale factors", "max scale ratio", "import scales", "log correction", "event output", "analysis", "analyses", "vev", - "regulator parameter", "max events" + "regulator parameter", "max events", "off-shell tolerance" }; // add subnodes to "supported" - the assigned value is irrelevant for(auto && opt: opts) supported[opt] = ""; for(auto && jet_opt: {"min pt", "algorithm", "R"}){ supported["resummation jets"][jet_opt] = ""; supported["fixed order jets"][jet_opt] = ""; } for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){ supported["Higgs coupling"][opt] = ""; } for(auto && opt: {"name", "seed"}){ supported["random generator"][opt] = ""; } for(auto && opt: {"enabled", "nlo order"}){ supported["NLO truncation"][opt] = ""; } for(auto && opt: {"FKL", "unordered", "extremal qqbar", "central qqbar", "non-resummable"}){ supported["event treatment"][opt] = ""; } for(auto && particle_type: {"Higgs", "W", "Z"}){ for(auto && particle_opt: {"mass", "width"}){ supported["particle properties"][particle_type][particle_opt] = ""; } } for(auto && opt: {"type", "trials", "max deviation"}){ supported["unweight"][opt] = ""; } return supported; }(); return supported; } fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){ using namespace fastjet; static const std::map known = { {"kt", kt_algorithm}, {"cambridge", cambridge_algorithm}, {"antikt", antikt_algorithm}, {"cambridge for passive", cambridge_for_passive_algorithm}, {"plugin", plugin_algorithm} }; const auto res = known.find(algo); if(res == known.end()){ throw std::invalid_argument("Unknown jet algorithm \"" + algo + "\""); } return res->second; } EventTreatment to_EventTreatment(std::string const & name){ static const std::map known = { {"reweight", EventTreatment::reweight}, {"keep", EventTreatment::keep}, {"discard", EventTreatment::discard} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown event treatment \"" + name + "\""); } return res->second; } WeightType to_weight_type(std::string const & setting){ if(setting == "weighted") return WeightType::weighted; if(setting =="resummation") return WeightType::unweighted_resum; if(setting =="partial") return WeightType::partially_unweighted; throw std::invalid_argument{"Unknown weight type \"" + setting + "\""}; } } // namespace namespace detail{ void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){ setting = to_JetAlgorithm(yaml.as()); } void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){ setting = to_EventTreatment(yaml.as()); } void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){ setting = to_ParticleID(yaml.as()); } void set_from_yaml(WeightType & setting, YAML::Node const & yaml){ setting = to_weight_type(yaml.as()); } } // namespace detail JetParameters get_jet_parameters( YAML::Node const & node, std::string const & entry ){ assert(node); JetParameters result; fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm; double R = NAN; set_from_yaml_if_defined(jet_algo, node, entry, "algorithm"); set_from_yaml(R, node, entry, "R"); result.def = fastjet::JetDefinition{jet_algo, R}; set_from_yaml(result.min_pt, node, entry, "min pt"); return result; } RNGConfig to_RNGConfig( YAML::Node const & node, std::string const & entry ){ assert(node); RNGConfig result; set_from_yaml(result.name, node, entry, "name"); set_from_yaml_if_defined(result.seed, node, entry, "seed"); return result; } NLOConfig to_NLOConfig( YAML::Node const & node, std::string const & entry ){ assert(node); NLOConfig result; set_from_yaml_if_defined(result.enabled, node, entry, "enabled"); set_from_yaml_if_defined(result.nj, node, entry, "nlo order"); return result; } ParticleProperties get_particle_properties( YAML::Node const & node, std::string const & entry, std::string const & boson ){ ParticleProperties result{}; set_from_yaml(result.mass, node, entry, boson, "mass"); set_from_yaml(result.width, node, entry, boson, "width"); return result; } EWConstants get_ew_parameters(YAML::Node const & node){ EWConstants result; double vev = NAN; set_from_yaml(vev, node, "vev"); result.set_vevWZH(vev, get_particle_properties(node, "particle properties", "W"), get_particle_properties(node, "particle properties", "Z"), get_particle_properties(node, "particle properties", "Higgs") ); return result; } HiggsCouplingSettings get_Higgs_coupling( YAML::Node const & node, std::string const & entry ){ assert(node); static constexpr double mt_max = 2e4; #ifndef HEJ_BUILD_WITH_QCDLOOP if(node[entry].IsDefined()){ throw std::invalid_argument{ "Higgs coupling settings require building HEJ 2 " "with QCDloop support" }; } #endif HiggsCouplingSettings settings; set_from_yaml_if_defined(settings.mt, node, entry, "mt"); set_from_yaml_if_defined(settings.mb, node, entry, "mb"); set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom"); set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors"); if(settings.use_impact_factors){ if(settings.mt != std::numeric_limits::infinity()){ throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } else{ // huge values of the top mass are numerically unstable settings.mt = std::min(settings.mt, mt_max); } return settings; } FileFormat to_FileFormat(std::string const & name){ static const std::map known = { {"Les Houches", FileFormat::Les_Houches}, {"HepMC", FileFormat::HepMC}, {"HepMC2", FileFormat::HepMC2}, {"HepMC3", FileFormat::HepMC3}, {"HDF5", FileFormat::HDF5} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown file format \"" + name + "\""); } return res->second; } std::string extract_suffix(std::string const & filename){ size_t separator = filename.rfind('.'); if(separator == std::string::npos) return {}; return filename.substr(separator + 1); } FileFormat format_from_suffix(std::string const & filename){ const std::string suffix = extract_suffix(filename); if(suffix == "lhe") return FileFormat::Les_Houches; if(suffix == "hepmc") return FileFormat::HepMC; if(suffix == "hepmc3") return FileFormat::HepMC3; if(suffix == "hepmc2") return FileFormat::HepMC2; if(suffix == "hdf5") return FileFormat::HDF5; throw std::invalid_argument{ "Can't determine format for output file \"" + filename + "\"" }; } void assert_all_options_known( YAML::Node const & conf, YAML::Node const & supported ){ if(!conf.IsMap()) return; if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"}; for(auto const & entry: conf){ const auto name = entry.first.as(); if(! supported[name]) throw unknown_option{name}; /* check sub-options, e.g. 'resummation jets: min pt' * we don't check analyses sub-options * those depend on the analysis being used and should be checked there * similar for "import scales" */ if(name != "analyses" && name != "analysis" && name != "import scales"){ try{ assert_all_options_known(conf[name], supported[name]); } catch(unknown_option const & ex){ throw unknown_option{name + ": " + ex.what()}; } catch(invalid_type const & ex){ throw invalid_type{name + ": " + ex.what()}; } } } } } // namespace HEJ namespace YAML { Node convert::encode(HEJ::OutputFile const & outfile) { Node node; node[to_string(outfile.format)] = outfile.name; return node; } bool convert::decode(Node const & node, HEJ::OutputFile & out) { switch(node.Type()){ case NodeType::Map: { YAML::const_iterator it = node.begin(); out.format = HEJ::to_FileFormat(it->first.as()); out.name = it->second.as(); return true; } case NodeType::Scalar: out.name = node.as(); out.format = HEJ::format_from_suffix(out.name); return true; default: return false; } } } // namespace YAML namespace HEJ { namespace detail{ void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){ setting = yaml.as(); } } namespace { void update_fixed_order_jet_parameters( JetParameters & fixed_order_jets, YAML::Node const & yaml ){ if(!yaml["fixed order jets"]) return; set_from_yaml_if_defined( fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt" ); fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm(); set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm"); double R = fixed_order_jets.def.R(); set_from_yaml_if_defined(R, yaml, "fixed order jets", "R"); fixed_order_jets.def = fastjet::JetDefinition{algo, R}; } // like std::stod, but throw if not the whole string can be converted double to_double(std::string const & str){ std::size_t pos = 0; const double result = std::stod(str, &pos); if(pos < str.size()){ throw std::invalid_argument(str + " is not a valid double value"); } return result; } using EventScale = double (*)(Event const &); void import_scale_functions( std::string const & file, std::vector const & scale_names, std::unordered_map & known ) { void * handle = dlopen(file.c_str(), RTLD_NOW); char * error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; for(auto const & scale: scale_names) { void * sym = dlsym(handle, scale.c_str()); error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; known.emplace(scale, reinterpret_cast(sym)); // NOLINT } } auto get_scale_map( YAML::Node const & yaml ) { std::unordered_map scale_map; scale_map.emplace("H_T", H_T); scale_map.emplace("max jet pperp", max_jet_pt); scale_map.emplace("jet invariant mass", jet_invariant_mass); scale_map.emplace("m_j1j2", m_j1j2); if(yaml["import scales"].IsDefined()) { if(! yaml["import scales"].IsMap()) { throw invalid_type{"Entry 'import scales' is not a map"}; } for(auto const & import: yaml["import scales"]) { const auto file = import.first.as(); const auto scale_names = import.second.IsSequence() ?import.second.as>() :std::vector{import.second.as()}; import_scale_functions(file, scale_names, scale_map); } } return scale_map; } // simple (as in non-composite) scale functions /** * An example for a simple scale function would be H_T, * H_T/2 is then composite (take H_T and then divide by 2) */ ScaleFunction parse_simple_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ) { assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); const auto it = known.find(scale_fun); if(it != end(known)) return {it->first, it->second}; try{ const double scale = to_double(scale_fun); return {scale_fun, FixedScale{scale}}; } catch(std::invalid_argument const &){} throw std::invalid_argument{"Unknown scale choice: \"" + scale_fun + "\""}; } std::string trim_front(std::string const & str){ const auto new_begin = std::find_if( begin(str), end(str), [](char c){ return std::isspace(c) == 0; } ); return std::string(new_begin, end(str)); } std::string trim_back(std::string str){ size_t pos = str.size() - 1; // use guaranteed wrap-around behaviour to check whether we have // traversed the whole string for(; pos < str.size() && std::isspace(str[pos]); --pos) {} str.resize(pos + 1); // note that pos + 1 can be 0 return str; } ScaleFunction parse_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ){ assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); // parse from right to left => a/b/c gives (a/b)/c const size_t delim = scale_fun.find_last_of("*/"); if(delim == std::string::npos){ return parse_simple_ScaleFunction(scale_fun, known); } const std::string first = trim_back(std::string{scale_fun, 0, delim}); const std::string second = trim_front(std::string{scale_fun, delim+1}); if(scale_fun[delim] == '/'){ return parse_ScaleFunction(first, known) / parse_ScaleFunction(second, known); } assert(scale_fun[delim] == '*'); return parse_ScaleFunction(first, known) * parse_ScaleFunction(second, known); } EventTreatMap get_event_treatment( YAML::Node const & node, std::string const & entry ){ using namespace event_type; EventTreatMap treat { {not_enough_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {FKL, EventTreatment::discard}, {unob, EventTreatment::discard}, {unof, EventTreatment::discard}, {qqbar_exb, EventTreatment::discard}, {qqbar_exf, EventTreatment::discard}, {qqbar_mid, EventTreatment::discard}, {non_resummable, EventTreatment::discard} }; set_from_yaml(treat.at(FKL), node, entry, "FKL"); set_from_yaml(treat.at(unob), node, entry, "unordered"); treat.at(unof) = treat.at(unob); set_from_yaml(treat.at(qqbar_exb), node, entry, "extremal qqbar"); treat.at(qqbar_exf) = treat.at(qqbar_exb); set_from_yaml(treat.at(qqbar_mid), node, entry, "central qqbar"); set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable"); if(treat[non_resummable] == EventTreatment::reweight){ throw std::invalid_argument{"Cannot reweight non-resummable events"}; } return treat; } Config to_Config(YAML::Node const & yaml){ try{ assert_all_options_known(yaml, get_supported_options()); } catch(unknown_option const & ex){ throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"}; } Config config; config.resummation_jets = get_jet_parameters(yaml, "resummation jets"); config.fixed_order_jets = config.resummation_jets; update_fixed_order_jet_parameters(config.fixed_order_jets, yaml); set_from_yaml_if_defined(config.min_extparton_pt, yaml, "min extparton pt"); if(config.min_extparton_pt!=0) std::cerr << "WARNING: \"min extparton pt\" is deprecated." << " Please remove this entry or set \"soft pt regulator\" instead.\n"; set_from_yaml_if_defined( config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction" ); if(config.max_ext_soft_pt_fraction){ std::cerr << "WARNING: \"max ext soft pt fraction\" is deprecated." << " Please remove this entry or set \"soft pt regulator\" instead.\n"; config.soft_pt_regulator = *config.max_ext_soft_pt_fraction; } else { set_from_yaml_if_defined( config.soft_pt_regulator, yaml, "soft pt regulator" ); } // Sets the standard value, then changes this if defined config.regulator_lambda=CLAMBDA; set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter"); set_from_yaml_if_defined(config.max_events, yaml, "max events"); set_from_yaml(config.trials, yaml, "trials"); config.weight_type = WeightType::weighted; set_from_yaml_if_defined(config.weight_type, yaml, "unweight", "type"); if(config.weight_type == WeightType::partially_unweighted) { config.unweight_config = PartialUnweightConfig{}; set_from_yaml( config.unweight_config->trials, yaml, "unweight", "trials" ); set_from_yaml( config.unweight_config->max_dev, yaml, "unweight", "max deviation" ); } else if(yaml["unweight"].IsDefined()) { for(auto && opt: {"trials", "max deviation"}) { if(yaml["unweight"][opt].IsDefined()) { throw std::invalid_argument{ "'unweight: " + std::string{opt} + "' " "is only supported if 'unweight: type' is set to 'partial'" }; } } } set_from_yaml(config.log_correction, yaml, "log correction"); config.treat = get_event_treatment(yaml, "event treatment"); set_from_yaml_if_defined(config.output, yaml, "event output"); config.rng = to_RNGConfig(yaml, "random generator"); set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses"); if(yaml["analysis"].IsDefined()){ std::cerr << "WARNING: Configuration entry 'analysis' is deprecated. " " Use 'analyses' instead.\n"; set_from_yaml(config.analysis_parameters, yaml, "analysis"); if(!config.analysis_parameters.IsNull()){ config.analyses_parameters.push_back(config.analysis_parameters); } } config.scales = to_ScaleConfig(yaml); config.ew_parameters = get_ew_parameters(yaml); config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling"); //HEJ@NLO Truncation config.nlo = to_NLOConfig(yaml, "NLO truncation"); + set_from_yaml_if_defined( + config.off_shell_tolerance, + yaml, + "off-shell tolerance" + ); return config; } } // namespace ScaleConfig to_ScaleConfig(YAML::Node const & yaml){ ScaleConfig config; auto scale_funs = get_scale_map(yaml); std::vector scales; set_from_yaml(scales, yaml, "scales"); config.base.reserve(scales.size()); std::transform( begin(scales), end(scales), std::back_inserter(config.base), [scale_funs](auto const & entry){ return parse_ScaleFunction(entry, scale_funs); } ); set_from_yaml_if_defined(config.factors, yaml, "scale factors"); config.max_ratio = std::numeric_limits::infinity(); set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio"); return config; } Config load_config(std::string const & config_file){ try{ return to_Config(YAML::LoadFile(config_file)); } catch(...){ std::cerr << "Error reading " << config_file << ":\n "; throw; } } } // namespace HEJ diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc index 94478e2..42f09dc 100644 --- a/src/bin/HEJ.cc +++ b/src/bin/HEJ.cc @@ -1,410 +1,425 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include #include "yaml-cpp/yaml.h" #include "fastjet/ClusterSequence.hh" #include "HEJ/CombinedEventWriter.hh" #include "HEJ/Config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/EventReader.hh" #include "HEJ/BufferedEventReader.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/get_analysis.hh" #include "HEJ/make_RNG.hh" #include "HEJ/ProgressBar.hh" #include "HEJ/stream.hh" #include "HEJ/Unweighter.hh" #include "HEJ/Version.hh" #include "HEJ/YAMLreader.hh" HEJ::Config load_config(char const * filename){ try{ return HEJ::load_config(filename); } catch(std::exception const & exc){ std::cerr << "Error: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } std::vector> get_analyses( std::vector const & parameters, LHEF::HEPRUP const & heprup ){ try{ return HEJ::get_analyses(parameters, heprup); } catch(std::exception const & exc){ std::cerr << "Failed to load analysis: " << exc.what() << '\n'; std::exit(EXIT_FAILURE); } } // unique_ptr is a workaround: // std::optional is a better fit, but gives spurious errors with g++ 7.3.0 std::unique_ptr> make_progress_bar( std::vector const & xs ) { if(xs.empty()) return {}; const double Born_xs = std::accumulate(begin(xs), end(xs), 0.); return std::make_unique>(std::cout, Born_xs); } std::string time_to_string(const time_t time){ char s[30]; struct tm * p = localtime(&time); strftime(s, 30, "%a %b %d %Y %H:%M:%S", p); return s; } void check_one_parton_per_jet(HEJ::Event const & ev) { const std::size_t npartons = std::count_if( ev.outgoing().begin(), ev.outgoing().end(), [](HEJ::Particle const & p) { return is_parton(p); } ); if(ev.jets().size() < npartons) { throw std::invalid_argument{ "Number of partons (" + std::to_string(npartons) + ") in input event" " does not match number of jets (" + std::to_string(ev.jets().size())+ ")" }; } } HEJ::Event to_event( LHEF::HEPEUP const & hepeup, HEJ::JetParameters const & fixed_order_jets, - HEJ::EWConstants const & ew_parameters + HEJ::EWConstants const & ew_parameters, + const double off_shell_tolerance ) { HEJ::Event::EventData event_data{hepeup}; event_data.reconstruct_intermediate(ew_parameters); - + event_data.repair_momenta(off_shell_tolerance); HEJ::Event ev{ std::move(event_data).cluster( fixed_order_jets.def, fixed_order_jets.min_pt ) }; check_one_parton_per_jet(ev); return ev; } void unweight( HEJ::Unweighter & unweighter, HEJ::WeightType weight_type, std::vector & events, HEJ::RNG & ran ) { if(weight_type == HEJ::WeightType::unweighted_resum){ unweighter.set_cut_to_maxwt(events); } events.erase( unweighter.unweight(begin(events), end(events), ran), end(events) ); } // peek up to nevents events from reader std::vector peek_events( HEJ::BufferedEventReader & reader, const int nevents ) { std::vector events; while( static_cast(events.size()) < nevents && reader.read_event() ) { events.emplace_back(reader.hepeup()); } // put everything back into the reader for(auto it = rbegin(events); it != rend(events); ++it) { reader.emplace(*it); } return events; } void append_resummed_events( std::vector & resummation_events, HEJ::EventReweighter & reweighter, LHEF::HEPEUP const & hepeup, const size_t trials, HEJ::JetParameters const & fixed_order_jets, - HEJ::EWConstants const & ew_parameters + HEJ::EWConstants const & ew_parameters, + const double off_shell_tolerance ) { - const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets, ew_parameters); + const HEJ::Event FO_event = to_event( + hepeup, + fixed_order_jets, + ew_parameters, + off_shell_tolerance + ); if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) { return; } const auto resummed = reweighter.reweight(FO_event, trials); resummation_events.insert( end(resummation_events), begin(resummed), end(resummed) ); } void train( HEJ::Unweighter & unweighter, HEJ::BufferedEventReader & reader, HEJ::EventReweighter & reweighter, const size_t total_trials, const double max_dev, double reweight_factor, HEJ::JetParameters const & fixed_order_jets, - HEJ::EWConstants const & ew_parameters + HEJ::EWConstants const & ew_parameters, + const double off_shell_tolerance ) { std::cout << "Reading up to " << total_trials << " training events...\n"; auto FO_events = peek_events(reader, total_trials); if(FO_events.empty()) { throw std::runtime_error{ "No events generated to calibrate the unweighting weight!" "Please increase the number \"trials\" or deactivate the unweighting." }; } const size_t trials = total_trials/FO_events.size(); // adjust reweight factor so that the overall normalisation // is the same as in the full run reweight_factor *= trials; for(auto & hepeup: FO_events) { hepeup.XWGTUP *= reweight_factor; } std::cout << "Training unweighter with " << trials << '*' << FO_events.size() << " events\n"; auto progress = HEJ::ProgressBar{ std::cout, static_cast(FO_events.size()) }; std::vector resummation_events; for(auto const & hepeup: FO_events) { append_resummed_events( resummation_events, reweighter, hepeup, trials, fixed_order_jets, - ew_parameters + ew_parameters, + off_shell_tolerance ); ++progress; } unweighter.set_cut_to_peakwt(resummation_events, max_dev); std::cout << "\nUnweighting events with weight up to " << unweighter.get_cut() << '\n'; } int main(int argn, char** argv) { using clock = std::chrono::system_clock; if (argn != 3) { std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n"; return EXIT_FAILURE; } const auto start_time = clock::now(); { std::cout << "Starting " << HEJ::Version::package_name_full() << ", revision " << HEJ::Version::revision() << " (" << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl; } fastjet::ClusterSequence::print_banner(); // read configuration const HEJ::Config config = load_config(argv[1]); auto reader = HEJ::make_reader(argv[2]); assert(reader); auto heprup{ reader->heprup() }; heprup.generators.emplace_back(LHEF::XMLTag{}); heprup.generators.back().name = HEJ::Version::package_name(); heprup.generators.back().version = HEJ::Version::String(); auto analyses = get_analyses( config.analyses_parameters, heprup ); assert(analyses.empty() || analyses.front() != nullptr); HEJ::CombinedEventWriter writer{config.output, std::move(heprup)}; double global_reweight = 1.; auto const & max_events = config.max_events; // if we need the event number: if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){ // try to read from LHE head auto input_events{reader->number_events()}; if(!input_events) { // else count manually auto t_reader = HEJ::make_reader(argv[2]); input_events = 0; while(t_reader->read_event()) ++(*input_events); } if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){ // IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs std::cout << "Found IDWTUP " << heprup.IDWTUP << ": " << "assuming \"cross section = average weight\".\n" << "converting to \"cross section = sum of weights\" "; global_reweight /= *input_events; } if(max_events && (*input_events > *max_events)){ // maximal number of events given global_reweight *= *input_events/static_cast(*max_events); std::cout << "Processing " << *max_events << " out of " << *input_events << " events\n"; } } HEJ::ScaleGenerator scale_gen{ config.scales.base, config.scales.factors, config.scales.max_ratio }; std::shared_ptr ran{ HEJ::make_RNG(config.rng.name, config.rng.seed)}; assert(ran != nullptr); HEJ::EventReweighter hej{ reader->heprup(), std::move(scale_gen), to_EventReweighterConfig(config), ran }; std::optional unweighter{}; if(config.weight_type != HEJ::WeightType::weighted) { unweighter = HEJ::Unweighter{}; } if(config.weight_type == HEJ::WeightType::partially_unweighted) { HEJ::BufferedEventReader buffered_reader{std::move(reader)}; assert(config.unweight_config); train( *unweighter, buffered_reader, hej, config.unweight_config->trials, config.unweight_config->max_dev, global_reweight/config.trials, config.fixed_order_jets, - config.ew_parameters + config.ew_parameters, + config.off_shell_tolerance ); reader = std::make_unique( std::move(buffered_reader) ); } // status infos & eye candy if (config.nlo.enabled){ std::cout << "HEJ@NLO Truncation Enabled for NLO order: " << config.nlo.nj << std::endl; } size_t nevent = 0; std::array nevent_type{0}, nfailed_type{0}; auto progress = make_progress_bar(reader->heprup().XSECUP); HEJ::CrossSectionAccumulator xs; std::map status_counter; size_t total_trials = 0; size_t total_resum = 0; // Loop over the events in the input file while(reader->read_event() && (!max_events || nevent < *max_events) ){ ++nevent; // reweight events so that the total cross section is conserved auto hepeup = reader->hepeup(); hepeup.XWGTUP *= global_reweight; - const auto FO_event = to_event(hepeup, config.fixed_order_jets, config.ew_parameters); + const auto FO_event = to_event( + hepeup, + config.fixed_order_jets, + config.ew_parameters, + config.off_shell_tolerance + ); if(FO_event.central().weight == 0) { static const bool warned_once = [argv,nevent](){ std::cerr << "WARNING: event number " << nevent << " in " << argv[2] << " has zero weight. " "Ignoring this and all further events with vanishing weight.\n"; return true; }(); (void) warned_once; // shut up compiler warnings continue; } auto resummed_events{ hej.reweight(FO_event, config.trials) }; // some bookkeeping for(auto const & s: hej.status()) ++status_counter[s]; total_trials+=hej.status().size(); ++nevent_type[FO_event.type()]; if(resummed_events.empty()) ++nfailed_type[FO_event.type()]; if(unweighter) { unweight(*unweighter, config.weight_type, resummed_events, *ran); } // analysis for(auto & ev: resummed_events){ //TODO: move pass_cuts to after phase space point generation bool passed = analyses.empty(); for(auto const & analysis: analyses){ if(analysis->pass_cuts(ev, FO_event)){ passed = true; analysis->fill(ev, FO_event); }; } if(passed){ writer.write(ev); } else { ev.parameters()*=0; // do not use discarded events afterwards } } xs.fill_correlated(resummed_events); total_resum += resummed_events.size(); if(progress) progress->increment(FO_event.central().weight); } // main event loop std::cout << '\n'; for(auto const & analysis: analyses){ analysis->finalise(); } writer.finish(); using namespace HEJ::event_type; std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n'; std::cout << '\t' << name(EventType::first_type) << ": " << nevent_type[EventType::first_type] << ", failed to reconstruct " << nfailed_type[EventType::first_type] << '\n'; for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){ std::cout << '\t' << name(static_cast(i)) << ": " << nevent_type[i] << ", failed to reconstruct " << nfailed_type[i] << '\n'; } std::cout << '\n' << xs << '\n'; std::cout << "Generation statistic: " << status_counter[HEJ::StatusCode::good] << "/" << total_trials << " trials successful.\n"; for(auto && entry: status_counter){ const double fraction = static_cast(entry.second)/total_trials; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":") << " ["; for(int i = 0; i < percent/2; ++i) std::cout << '#'; for(int i = percent/2; i < 50; ++i) std::cout << ' '; std::cout << "] " < run_time = (clock::now() - start_time); std::cout << "\nFinished " << HEJ::Version::package_name() << " at " << time_to_string(clock::to_time_t(clock::now())) << "\n=> Runtime: " << run_time.count() << " sec (" << nevent/run_time.count() << " Events/sec).\n"; return EXIT_SUCCESS; } diff --git a/t/CMakeLists.txt b/t/CMakeLists.txt index 1ee1ec2..a02ffc2 100644 --- a/t/CMakeLists.txt +++ b/t/CMakeLists.txt @@ -1,624 +1,632 @@ set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}") set(tst_ME_data_dir "${tst_dir}/ME_data") # small library for common test functions add_library(hej_test SHARED hej_test.cc) target_include_directories(hej_test PUBLIC ${tst_dir}) target_link_libraries(hej_test HEJ) # test event classification # test explicit configurations add_executable(test_classify ${tst_dir}/test_classify.cc) target_compile_options(test_classify PRIVATE "-O0") # avoid compiler optimisation target_link_libraries(test_classify HEJ hej_test) add_test( NAME classify COMMAND test_classify ) # test against reference data add_executable(test_classify_ref ${tst_dir}/test_classify_ref.cc) target_link_libraries(test_classify_ref HEJ hej_test) add_test( NAME classify_ref COMMAND test_classify_ref ${tst_dir}/classify_ref ${tst_dir}/classify.lhe.gz ) add_test( NAME classify_ref_4j COMMAND test_classify_ref ${tst_dir}/classify_ref_4j ${tst_dir}/4j.lhe.gz ) add_test( NAME classify_ref_W4j COMMAND test_classify_ref ${tst_dir}/classify_ref_W4j ${tst_dir}/MGemnubar4j_unweighted.lhe.gz ) add_test( NAME t_classify_ref_Z4j COMMAND test_classify_ref ${tst_dir}/classify_ref_Z4j ${tst_dir}/Z_4j.lhe.gz ) # test for valid W decays add_executable(test_decay ${tst_dir}/test_decay.cc) target_link_libraries(test_decay HEJ hej_test) add_test( NAME valid_decay COMMAND test_decay ) # test valid jet cuts on tagging jets add_executable(test_jet_cuts ${tst_dir}/test_jet_cuts.cc) target_link_libraries(test_jet_cuts HEJ hej_test) add_test( NAME jet_cuts COMMAND test_jet_cuts ) # test phase space point add_executable(test_psp ${tst_dir}/test_psp.cc) target_link_libraries(test_psp HEJ hej_test) add_test( NAME PhaseSpace COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz ) # test importing analyses file(COPY "${tst_dir}/analysis_config.yml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") set(test_config "${CMAKE_CURRENT_BINARY_DIR}/analysis_config.yml") get_target_property(ANALYSIS_PATH AnalysisTemplate_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisTemplate_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS "") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_simple.yml @ONLY ) add_test( NAME analysis_simple COMMAND $ analysis_config_simple.yml ${tst_dir}/2j.lhe.gz ) get_target_property(ANALYSIS_PATH AnalysisPrint_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisPrint_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS " output: ana_output") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_print.yml @ONLY ) add_test( NAME analysis_print COMMAND $ analysis_config_print.yml ${tst_dir}/2j.lhe.gz ) if(ROOT_FOUND) get_target_property(ANALYSIS_PATH AnalysisROOT_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisROOT_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS "") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_root.yml @ONLY ) add_test( NAME analysis_root COMMAND $ analysis_config_root.yml ${tst_dir}/2j.lhe.gz ) endif() if(RIVET_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - rivet: MC_XS\n output: ana_rivet\n") endif() add_test( NAME analysis_all COMMAND $ ${test_config} ${tst_dir}/2j.lhe.gz ) # test importing scales (from examples/softestptScale) add_executable(test_scale_import ${tst_dir}/test_scale_import) target_link_libraries(test_scale_import HEJ) get_target_property(SCALE_PATH softestptScale_lib BINARY_DIR) get_target_property(SCALE_LIB softestptScale_lib OUTPUT_NAME) set(SCALE_NAME "softest_jet_pt") configure_file( ${tst_dir}/jet_config_with_import.yml.in jet_config_with_import.yml @ONLY ) add_test( NAME scale_import COMMAND test_scale_import jet_config_with_import.yml ) # test scale arithmetic (e.g. 2*H_T/4) add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics) target_link_libraries(test_scale_arithmetics HEJ hej_test) add_test( NAME scale_arithmetics COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz ) # test "ParameterDescription" add_executable(test_descriptions ${tst_dir}/test_descriptions) target_link_libraries(test_descriptions HEJ hej_test) add_test( NAME descriptions COMMAND test_descriptions ) # test "EventParameters*Weight" add_executable(test_parameters ${tst_dir}/test_parameters) target_link_libraries(test_parameters HEJ hej_test) add_test( NAME parameters COMMAND test_parameters ) # test Particle IDs add_executable(test_pid ${tst_dir}/test_pid) target_link_libraries(test_pid HEJ hej_test) add_test( NAME pid COMMAND test_pid ) # test EW constants add_executable(test_EWConstants ${tst_dir}/test_EWConstants) target_link_libraries(test_EWConstants HEJ hej_test) add_test( NAME EWConstants COMMAND test_EWConstants ) # test unweighting add_executable(test_unweighter ${tst_dir}/test_unweighter) target_link_libraries(test_unweighter HEJ hej_test) add_test( NAME unweighter COMMAND test_unweighter ${tst_dir}/4j.lhe.gz ) # test colour generation add_executable(test_colours ${tst_dir}/test_colours) target_link_libraries(test_colours HEJ hej_test) add_test( NAME colour_flow COMMAND test_colours ) add_executable(test_colours2 ${tst_dir}/test_colours2) target_link_libraries(test_colours2 HEJ hej_test) add_test( NAME colour_flow2 COMMAND test_colours2 ) # test matrix elements add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc) target_link_libraries(test_ME_generic HEJ hej_test) add_test( NAME ME_j COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_j_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree_virt.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_h COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_h_subl COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_subl_mtinf_tree.dat ${tst_ME_data_dir}/PSP_h_subl.lhe.gz ) add_test( NAME ME_hj COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_hj_mtinf_tree.dat ${tst_ME_data_dir}/PSP_hj.lhe.gz ) add_test( NAME ME_h_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_virt.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_hj_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_hj_mtinf_virt.dat ${tst_ME_data_dir}/PSP_hj.lhe.gz ) if(QCDloop_FOUND) add_test( NAME ME_h_mt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_h_mtmb COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_h_subl_mtmb COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_subl_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h_subl.lhe.gz ) add_test( NAME ME_hj_mt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_hj_mt_tree.dat ${tst_ME_data_dir}/PSP_hj.lhe.gz ) add_test( NAME ME_hj_mtmb COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_hj_mtmb_tree.dat ${tst_ME_data_dir}/PSP_hj.lhe.gz ) endif() add_test( NAME ME_j_subl COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree.dat ${tst_ME_data_dir}/PSP_jets_subl.lhe.gz ) add_test( NAME ME_j_subl_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_virt.dat ${tst_ME_data_dir}/PSP_jets_subl.lhe.gz ) add_test( NAME ME_j_subl_4j COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_new.dat ${tst_dir}/4j.lhe.gz ) add_test( NAME ME_j_subl_4j_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_new_virt.dat ${tst_dir}/4j.lhe.gz ) add_test( NAME ME_w_FKL COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_tree.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz ) add_test( NAME ME_w_FKL_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_virt.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz ) add_test( NAME ME_Wp COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz ) add_test( NAME ME_Wp_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp_virt.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz ) add_test( NAME ME_Wm COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz ) add_test( NAME ME_Wm_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm_virt.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz ) add_test( NAME ME_Z COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z.dat ${tst_ME_data_dir}/PSP_Z.lhe.gz ) add_test( NAME ME_Z_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_virt.dat ${tst_ME_data_dir}/PSP_Z.lhe.gz ) add_test( NAME ME_Z_uno COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_uno_tree.dat ${tst_ME_data_dir}/PSP_Z_uno.lhe.gz ) add_test( NAME ME_Z_uno_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_uno_virt.dat ${tst_ME_data_dir}/PSP_Z_uno.lhe.gz ) add_test( NAME ME_WpWp COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME.yml ${tst_ME_data_dir}/ME_WpWp.dat ${tst_ME_data_dir}/PSP_WpWp.lhe.gz ) add_test( NAME ME_WpWp_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME.yml ${tst_ME_data_dir}/ME_WpWp_virt.dat ${tst_ME_data_dir}/PSP_WpWp.lhe.gz ) add_test( NAME ME_WmWm COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME.yml ${tst_ME_data_dir}/ME_WmWm.dat ${tst_ME_data_dir}/PSP_WmWm.lhe.gz ) add_test( NAME ME_WmWm_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME.yml ${tst_ME_data_dir}/ME_WmWm_virt.dat ${tst_ME_data_dir}/PSP_WmWm.lhe.gz ) # HEJ@NLO Virtuals add_test( NAME ME_WmWm_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME_hejnlo.yml ${tst_ME_data_dir}/ME_WmWm_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_WmWm.lhe.gz ) add_test( NAME ME_WpWp_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_ww_ME_hejnlo.yml ${tst_ME_data_dir}/ME_WpWp_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_WpWp.lhe.gz ) add_test( NAME ME_Wp_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_wzhpure_ME_hejnlo.yml ${tst_ME_data_dir}/ME_Wp_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz ) add_test( NAME ME_Wm_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_wzhpure_ME_hejnlo.yml ${tst_ME_data_dir}/ME_Wm_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz ) add_test( NAME ME_Z_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_wzhpure_ME_hejnlo.yml ${tst_ME_data_dir}/ME_Z_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_Z.lhe.gz ) add_test( NAME ME_jets_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_wzhpure_ME_hejnlo.yml ${tst_ME_data_dir}/ME_jets_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_h_hejnlo COMMAND test_ME_generic ${tst_ME_data_dir}/config_wzhpure_ME_hejnlo.yml ${tst_ME_data_dir}/ME_h_virt_hejnlo.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_j_log_corr COMMAND test_ME_generic ${tst_ME_data_dir}/config_log_corr.yml ${tst_ME_data_dir}/ME_jets_log_corr.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_j_lambda COMMAND test_ME_generic ${tst_ME_data_dir}/config_lambda.yml ${tst_ME_data_dir}/ME_jets_lambda.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) # test main executable file(COPY "${tst_dir}/jet_config.yml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") set(test_config "${CMAKE_CURRENT_BINARY_DIR}/jet_config.yml") if(HighFive_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hdf5\n") endif() if(HepMC3_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hepmc\n") endif() if(HepMC_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hepmc2\n") endif() if(rivet_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config}\nanalysis:\n rivet: MC_XS\n output: tst\n") endif() set(test_cmd_main "$\\\;${test_config}\\\;${tst_dir}/2j.lhe.gz") # check that HepMC3 output is correct if(HepMC3_FOUND) add_executable(check_hepmc ${tst_dir}/check_hepmc.cc) target_link_libraries(check_hepmc HEJ ${HEPMC3_LIBRARIES}) target_include_directories(check_hepmc PRIVATE ${HEPMC3_INCLUDE_DIR}) set(test_cmd_hepmc "$\\\;tst.hepmc") else() set(test_cmd_hepmc "") endif() # check that LHEF output is correct add_executable(check_lhe ${tst_dir}/check_lhe.cc) target_link_libraries(check_lhe HEJ hej_test) set(test_cmd_lhe "$\\\;tst.lhe") # check that rivet interface is consistent with naive rivet if(rivet_FOUND) # this assumes "rivet" and "yodadiff" are found in PATH if(rivet_USE_HEPMC3) set(hepmc_file "tst.hepmc") else() set(hepmc_file "tst.hepmc2") endif() if(rivet_USE_HEPMC3 OR (rivet_VERSION VERSION_LESS 3)) set(histo_exclude "") else() # rivet 3 with HepMC 2 is inconsistent in order of weights # -> interface != direct call (by permutation) # REQUIRES Yoda 1.7.5 set(histo_exclude "-M\\\;\\\\d") endif() set(test_cmd_rivet "rivet\\\;-a\\\;MC_XS\\\;${hepmc_file}\\\;-o\\\;tst_direct.yoda\ \;yodadiff\\\;${histo_exclude}\\\;tst.yoda\\\;tst_direct.yoda") else() set(test_cmd_rivet "") endif() # Run dependent tests in one command to ensure correct execution order # Note: The commands are concatenated with "\;" to escape CMake lists. # Thus arguments have to be escaped twice "\\\;". # e.g. "cmd1\;cmd2\\\;arg1\\\;arg2" is executed like "cmd1 && cmd2 arg1 arg2" add_test( NAME main COMMAND ${CMAKE_COMMAND} -DCMDS=${test_cmd_main}\;${test_cmd_hepmc}\;${test_cmd_lhe}\;${test_cmd_rivet} -P ${PROJECT_SOURCE_DIR}/cmake/run_multiple_tests.cmake ) add_test( NAME main_example COMMAND $ ${PROJECT_SOURCE_DIR}/config.yml ${tst_dir}/h_3j.lhe.gz ) # check that Sherpas LHE input can be read add_executable(check_lhe_sherpa ${tst_dir}/check_lhe_sherpa.cc) target_link_libraries(check_lhe_sherpa HEJ hej_test) add_test( NAME sherpa_reader COMMAND check_lhe_sherpa ${tst_dir}/SherpaLHE.lhe 1.62624e+08 ) # check HDF5 reader & writer if(HighFive_FOUND) add_executable(test_hdf5 ${tst_dir}/test_hdf5.cc) target_link_libraries(test_hdf5 HEJ) add_test( NAME hdf5_read COMMAND test_hdf5 ${tst_dir}/Wm9-g4-repack.hdf5 ) add_executable(test_hdf5_write ${tst_dir}/test_hdf5_write.cc) target_link_libraries(test_hdf5_write HEJ hej_test) add_test( NAME hdf5_write COMMAND test_hdf5_write ${tst_dir}/Wm9-g4-repack.hdf5 output.hdf5 ) endif() # check rivet interface if(RIVET_FOUND) add_executable(check_rivet ${tst_dir}/check_rivet.cc) target_link_libraries(check_rivet HEJ rivet::rivet) add_test( NAME rivet COMMAND check_rivet ) endif() # test boson reconstruction add_executable(cmp_events ${tst_dir}/cmp_events.cc) target_link_libraries(cmp_events HEJ) add_test( NAME reconstruct_W COMMAND cmp_events ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz ${tst_dir}/epnu2jLOFKL_unweight_noW.lhe.gz ) +# test repairing slightly off-shell momenta +add_executable(repair_off-shell ${tst_dir}/repair_off-shell.cc) +target_link_libraries(repair_off-shell HEJ) +add_test( + NAME repair_off-shell + COMMAND repair_off-shell ${tst_dir}/off-shell.lhe.gz + ) + # test resummed result add_executable(check_res ${tst_dir}/check_res.cc) target_link_libraries(check_res HEJ hej_test) add_test( NAME xs_h_1j COMMAND check_res ${tst_dir}/h_1j.lhe.gz 4.97087 0.0302519 ) add_test( NAME xs_2j COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684 ) add_test( NAME xs_3j COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6 ) add_test( NAME xs_3j_unof COMMAND check_res ${tst_dir}/3j.lhe.gz 133399 4688.83 unof ) add_test( NAME xs_3j_unob COMMAND check_res ${tst_dir}/3j.lhe.gz 105247 3449.45 unob ) add_test( NAME xs_3j_splitf COMMAND check_res ${tst_dir}/3j.lhe.gz 97659.9 2748.65 splitf ) add_test( NAME xs_3j_splitb COMMAND check_res ${tst_dir}/3j.lhe.gz 107150 2799.8 splitb ) add_test( NAME xs_4j COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6 ) add_test( NAME xs_4j_qqbar_mid COMMAND check_res ${tst_dir}/4j.lhe.gz 21866.7 1716.96 qqbar_mid ) add_test( NAME xs_h_3j COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.667678 0.0126859 ) add_test( NAME xs_h_3j_unof COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0110041 0.000274214 unof ) add_test( NAME xs_h_3j_unob COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0116516 0.000376737 unob ) add_test( NAME xs_epnu_3j COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 84394.7 3770.35 ) add_test( NAME xs_emnubar_3j COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted.lhe.gz 53489.1 977.234 ) add_test( NAME xs_epnu_4j COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.lhe.gz 19616.2 774.336 ) add_test( NAME xs_emnubar_4j COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted.lhe.gz 10169.6 248.741 ) if(TEST_ALL) # deactivate long tests by default add_test( NAME xs_h_1j_long COMMAND check_res ${tst_dir}/h_1j_large.lhe.gz 4.70112 0.00415353 ) add_test( NAME xs_epnu_2j COMMAND check_res ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz 262.7 3 ) add_test( NAME xs_epnu_3j_long COMMAND check_res ${tst_dir}/MGepnu3j_unweighted_large.lhe.gz 38.9512 1 ) add_test( NAME xs_emnubar_3j_long COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted_large.lhe.gz 24.1575 1 ) add_test( NAME xs_epnu_3j_unof COMMAND check_res ${tst_dir}/MGepnu3j_unweighted_large.lhe.gz 9.63702 0.128355 unof ) add_test( NAME xs_epnu_3j_unob COMMAND check_res ${tst_dir}/MGepnu3j_unweighted_large.lhe.gz 9.70119 0.108436 unob ) add_test( NAME xs_epnu_3j_splitf COMMAND check_res ${tst_dir}/MGepnu3j_unweighted_large.lhe.gz 2.91995 0.0463182 splitf ) add_test( NAME xs_epnu_3j_splitb COMMAND check_res ${tst_dir}/MGepnu3j_unweighted_large.lhe.gz 3.40708 0.0550975 splitb ) add_test( NAME xs_epnu_4j_long COMMAND check_res ${tst_dir}/MGepnu4j_unweighted_large.lhe.gz 10.2542 0.135106 ) add_test( NAME xs_emnubar_4j_long COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted_large.lhe.gz 5.73282 0.0583738 ) add_test( NAME xs_epnu_4j_qqbar_mid COMMAND check_res ${tst_dir}/MGepnu4j_unweighted_large.lhe.gz 0.802603 0.0126908 qqbar_mid ) endif() diff --git a/t/off-shell.lhe.gz b/t/off-shell.lhe.gz new file mode 100644 index 0000000..6a544b0 --- /dev/null +++ b/t/off-shell.lhe.gz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:87990f06c3c57c26e0b72a9105e258d09158cdae39eefa7bd484e090e5543c40 +size 1152 diff --git a/t/repair_off-shell.cc b/t/repair_off-shell.cc new file mode 100644 index 0000000..3182a83 --- /dev/null +++ b/t/repair_off-shell.cc @@ -0,0 +1,52 @@ +/** + * \authors The HEJ collaboration (see AUTHORS for details) + * \date 2022 + * \copyright GPLv2 or later + */ +#include +#include +#include + +#include "hej_test.hh" + +#include "HEJ/Event.hh" +#include "HEJ/Particle.hh" +#include "HEJ/EventReader.hh" + +namespace { + constexpr double TOLERANCE = 0.1; + constexpr double EP = 1e-5; + + void test_repair(LHEF::HEPEUP const & hepeup, const double tolerance) { + HEJ::Event::EventData ev{hepeup}; + ev.repair_momenta(tolerance); + for(auto const & in: ev.incoming) { + ASSERT(in.px() == 0. || std::abs(in.px()) > tolerance); + ASSERT(in.py() == 0. || std::abs(in.py()) > tolerance); + ASSERT(in.m() < EP || in.m() > tolerance); + } + for(auto const & out: ev.outgoing) { + ASSERT(is_massive(out) || out.m() < EP || out.m() > tolerance); + } + for(auto const & decay: ev.decays) { + for(auto const & out: decay.second) { + ASSERT(is_massive(out) || out.m() < EP || out.m() > tolerance); + } + } + } +} + +int main(int argc, char** argv) { + if(argc != 2) { + std::cerr << "Usage: repair_off-shell EVENTS.LHE\n"; + return EXIT_FAILURE; + } + + auto reader = HEJ::make_reader(argv[1]); + while(reader->read_event()) { + for(auto tolerance: {0., TOLERANCE, std::numeric_limits::max()}) { + test_repair(reader->hepeup(), tolerance); + } + } + +}