Page MenuHomeHEPForge
Files F7878687

No OneTemporary
Actions

View Options

diff --git a/config.yml b/config.yml
index 0bfdcce..d26d29c 100644
--- a/config.yml
+++ b/config.yml
@@ -1,140 +1,142 @@
## 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
+## the supported settings are: reweight, keep, discard, abort
+## non-resummable, unknown, and invalid events cannot be reweighted
event treatment:
FKL: reweight
unordered: keep
extremal qqbar: keep
central qqbar: keep
non-resummable: keep
+ unknown: abort
+ invalid: abort
## 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
## Only keep low pt contributions
# require low pt jet: false
## 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 72247e2..0c18246 100644
--- a/doc/sphinx/HEJ.rst
+++ b/doc/sphinx/HEJ.rst
@@ -1,488 +1,513 @@
.. _`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 <https://launchpad.net/mg5amcnlo>`_
or `Sherpa <https://sherpa.hepforge.org/trac/wiki>`_. If HEJ 2 was
compiled with `HDF5 <https://www.hdfgroup.org/>`_ support, it can also
read and write event files in the format suggested in
`arXiv:1905.05120 <https://arxiv.org/abs/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 <https://hub.docker.com/r/hejdock/hej>`_,
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 <http://yaml.org/>`_
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
<http://fastjet.fr/>`_ 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`_. The
contribution from the lower min pt fixed order jets can also be
calculated separately using the `require low pt jet`_ option.
.. _`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:
+ contribute to different orders in the high-energy limit. The 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.
+ Specifies how to treat events where resummation is currently not
+ supported for the given rapidity ordering and/or assignment of
+ parton flavours.
+
+ .. _`unknown`:
+
+ **unknown**
+ Specifies how to treat events for processes that are not
+ implemented, e.g. di-photon production.
+
+ .. _`invalid`:
+
+ **invalid**
+ Specifies how to treat events that are considered
+ unphysical. Events in this category violate one or more of the
+ following requirements:
+
+ - Charge conservation.
+ - Four-momentum conservation.
+ - Momenta of massless particles must be lightlike.
+ - Incoming particles must not have transverse momentum.
+
+ The possible treatments 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
+ - :code:`discard` to discard these events
+ - :code:`abort` to abort the program if one of these events is encountered
+
+ The settings :code:`keep`, :code:`discard`, and :code:`abort` are
+ supported for all event types. The :code:`reweight` treatment is
+ implemented for different event types depending on the considered
+ process:
+
+ .. 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"
+
+
.. _`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 <https://arxiv.org/abs/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
<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ 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<vev>`.
.. _`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
<https://github.com/scarrazza/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.
.. _`require low pt jet`:
**require low pt jet**
Restrict output to only include events which have fixed order jets with pT
less than the `resummation jets: min pt`_. Used to calculate the
contribution from a lower min pt in the fixed order jets. This option allows
this calculation to be done separately from the calculation where the fixed
order min pt is equal to the resummation jets min pt. This is usually a small
correction.
.. _`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 <https://arxiv.org/abs/1706.01002>`_. Default is 0.2.
diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
index 3b399de..e1f0b26 100644
--- a/include/HEJ/Config.hh
+++ b/include/HEJ/Config.hh
@@ -1,291 +1,292 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <string>
#include <utility>
#include <optional>
#include <vector>
#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<ScaleFunction> base;
//! Factors for multiplicative scale variation
std::vector<double> 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<std::string> 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 */
+ abort, /**< Throw an exception */
};
//! Container to store the treatments for various event types
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
//! 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<Fraction<double>> max_ext_soft_pt_fraction{};
//! @brief Maximum transverse momentum fraction from soft radiation in any
//! tagging jet (e.g. extremal or qqbar jet)
Fraction<double> 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<size_t> 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<OutputFile> 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<YAML::Node> 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<PartialUnweightConfig> unweight_config;
//! HEJ@NLO settings
NLOConfig nlo;
//! LowPT settings
bool lowpt = false;
//! 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<double> soft_pt_regulator =
Fraction<double>{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<Fraction<double>> max_ext_soft_pt_fraction{};
//! @brief Maximum transverse momentum fraction from soft radiation in any
//! tagging jet (e.g. extremal or qqbar jet)
Fraction<double> 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<double> soft_pt_regulator = Fraction<double>{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<double> 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;
//! Option to only keep lowpt contribution
bool lowpt = false;
};
/**! 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,
conf.lowpt
};
}
} // namespace HEJ
diff --git a/include/HEJ/EventReweighter.hh b/include/HEJ/EventReweighter.hh
index d65e545..2a4d213 100644
--- a/include/HEJ/EventReweighter.hh
+++ b/include/HEJ/EventReweighter.hh
@@ -1,196 +1,206 @@
/** \file
* \brief Declares the EventReweighter class
*
* EventReweighter is the main class used within HEJ 2. It reweights the
* resummation events.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include <memory>
#include <utility>
#include <vector>
#include "HEJ/Config.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/ScaleFunction.hh"
#include "HEJ/StatusCode.hh"
#include "HEJ/event_types.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
struct RNG;
//! Beam parameters
/**
* Currently, only symmetric beams are supported,
* so there is a single beam energy.
*/
struct Beam{
double E; /**< Beam energy */
std::array<ParticleID, 2> type; /**< Beam particles */
};
//! Main class for reweighting events in HEJ.
class EventReweighter{
using EventType = event_type::EventType;
public:
EventReweighter(
Beam const & beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
std::shared_ptr<RNG> ran /**< Random number generator */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< LHEF event header */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
std::shared_ptr<RNG> ran /**< Random number generator */
);
//! Get the used pdf
PDF const & pdf() const;
//! Check the lowpt only restriction passes for lowpt runs
bool pass_low_pt(
HEJ::Event const & input_ev
);
//! Get event treatment
EventTreatment treatment(EventType type) const;
//! Generate resummation events for a given fixed-order event
/**
* @param ev Fixed-order event corresponding
* to the resummation events
* @param num_events Number of trial resummation configurations.
* @returns A vector of resummation events.
*
* The result vector depends on the type of the input event and the
* \ref EventTreatment of different types as specified in the constructor:
*
* - EventTreatment::reweight: The result vector contains between 0 and
* num_events resummation events.
* - EventTreatment::keep: If the input event passes the resummation
* jet cuts the result vector contains one
* event. Otherwise it is empty.
* - EventTreatment::discard: The result vector is empty
*/
std::vector<Event> reweight(
Event const & ev,
std::size_t num_events
);
//! Gives all StatusCodes of the last reweight()
/**
* Each StatusCode corresponds to one tried generation. Only good
* StatusCodes generated an event.
*/
std::vector<StatusCode> const & status() const {
return status_;
}
private:
/** \internal
* \brief main generation/reweighting function:
* generate phase space points and divide out Born factors
*/
std::vector<Event> gen_res_events(
Event const & ev, std::size_t phase_space_points
);
std::vector<Event> rescale(
Event const & Born_ev, std::vector<Event> events
) const;
/** \internal
* \brief Do the Jets pass the resummation Cuts?
*
* @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
bool jets_pass_resummation_cuts(Event const & ev) const;
/** \internal
* \brief pdf_factors Function
*
* @param ev Event in Question
* @returns EventFactor due to PDFs
*
* Calculates the Central value and the variation due
* to the PDF choice made.
*/
Weights pdf_factors(Event const & ev) const;
/** \internal
* \brief matrix_elements Function
*
* @param ev Event in question
* @returns EventFactor due to MatrixElements
*
* Calculates the Central value and the variation due
* to the Matrix Element.
*/
Weights matrix_elements(Event const & ev) const;
/** \internal
* \brief Scale-dependent part of fixed-order matrix element
*
* @param ev Event in question
* @returns EventFactor scale variation due to FO-ME.
*
* This is only called to compute the scale variation for events where
* we don't do resummation (e.g. non-FKL).
* Since at tree level the scale dependence is just due to alpha_s,
* it is enough to return the alpha_s(mur) factors in the matrix element.
* The rest drops out in the ratio of (output event ME)/(input event ME),
* so we never have to compute it.
*/
Weights fixed_order_scale_ME(Event const & ev) const;
/** \internal
* \brief Computes the tree level matrix element
*
* @param ev Event in Question
* @returns HEJ approximation to Tree level Matrix Element
*
* This computes the HEJ approximation to the tree level FO
* Matrix element which is used within the LO weighting process.
*/
double tree_matrix_element(Event const & ev) const;
//! \internal General parameters
EventReweighterConfig param_;
//! \internal Beam energy
double E_beam_;
//! \internal PDF
PDF pdf_;
//! \internal Object to calculate the square of the matrix element
MatrixElement MEt2_;
//! \internal Object to calculate event renormalisation and factorisation scales
ScaleGenerator scale_gen_;
//! \internal random number generator
std::shared_ptr<RNG> ran_;
//! \internal StatusCode of each attempt
std::vector<StatusCode> status_;
};
+ //! Exception indicating that an event is not accepted by the reweighter
+ /**
+ * When constructing an EventReweighter one can specify the \ref abort treatment
+ * for a certain EventType. If such an Event is passed to EventReweighter::reweight(),
+ * this exception is thrown.
+ */
+ struct abort_event: std::invalid_argument {
+ explicit abort_event(Event const & ev);
+ };
+
} // namespace HEJ
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 14ecdb1..7ee97df 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,293 +1,301 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <string>
#include <unordered_map>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/Fraction.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PhaseSpacePoint.hh"
#include "HEJ/exceptions.hh"
namespace HEJ {
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
EventReweighter{
Beam{
heprup.EBMUP.first,
{{
static_cast<ParticleID>(heprup.IDBMUP.first),
static_cast<ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
std::move(ran)
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
}
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam const & beam,
int pdf_id,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
param_{std::move(conf)},
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
param_.ME_config
},
scale_gen_{std::move(scale_gen)},
ran_{std::move(ran)}
{
// legacy code: override new variable with old
if(param_.psp_config.max_ext_soft_pt_fraction){
param_.psp_config.soft_pt_regulator = *param_.psp_config.max_ext_soft_pt_fraction;
param_.psp_config.max_ext_soft_pt_fraction = {};
}
assert(ran_);
}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
bool EventReweighter::pass_low_pt(
HEJ::Event const & input_ev
){
// Keep only events where there is a fixed order event with at least 1
// jet below the resummation jet pt but all resummation jets are above
// the resummation jet pt
if(param_.treat.at(EventType::non_resummable)
!= EventTreatment::discard){
throw std::logic_error{
"Non-resummable events should be discarded for lowpt runs"
};
}
return std::any_of(begin(input_ev.jets()),
end(input_ev.jets()),
[&](fastjet::PseudoJet jet)
{return jet.pt() < param_.jet_param().min_pt;});
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, std::size_t num_events
){
if(param_.lowpt && !EventReweighter::pass_low_pt(input_ev)){
return {};
}
auto res_events{ gen_res_events(input_ev, num_events) };
if(res_events.empty()) return {};
for(auto & event: res_events) event = scale_gen_(std::move(event));
return rescale(input_ev, std::move(res_events));
}
EventTreatment EventReweighter::treatment(EventType type) const {
return param_.treat.at(type);
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
std::size_t phase_space_points
){
assert(ev.variations().empty());
status_.clear();
switch(treatment(ev.type())){
case EventTreatment::discard: {
status_.emplace_back(StatusCode::discard);
return {};
}
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) {
status_.emplace_back(StatusCode::failed_resummation_cuts);
return {};
}
else {
status_.emplace_back(StatusCode::good);
return {ev};
}
+ case EventTreatment::abort:
+ throw abort_event{ev};
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
status_.reserve(phase_space_points);
for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
status_.emplace_back(psp.status());
assert(psp.status() != StatusCode::unspecified);
if(psp.status() != StatusCode::good) continue;
assert(psp.weight() != 0.);
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
status_.back() = StatusCode::too_much_energy;
continue;
}
resummation_events.emplace_back(
to_EventData( std::move(psp) ).cluster(
param_.jet_param().def, param_.jet_param().min_pt
)
);
auto & new_event = resummation_events.back();
assert( new_event.valid_hej_state(
param_.psp_config.soft_pt_regulator,
param_.psp_config.min_extparton_pt ) );
if( new_event.type() != ev.type() ) {
throw std::logic_error{
"Resummation Event does not match Born event: "
+ name(new_event.type())
+ " != "
+ name(ev.type())
};
}
new_event.generate_colours(*ran_);
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.parameters() *= pdf*ME/(Born_pdf*Born_ME);
}
return events;
}
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
}
Weights EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
Weights result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & ev_param: ev.variations()){
const double muf = ev_param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
Weights
EventReweighter::matrix_elements(Event const & ev) const{
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
return MEt2_(ev);
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(ev).central;
}
Weights
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
int alpha_s_power = 0;
for(auto const & part: ev.outgoing()){
if(is_parton(part))
++alpha_s_power;
else if(part.type == pid::Higgs) {
alpha_s_power += 2;
}
// nothing to do for other uncoloured particles
}
Weights result;
result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
std::pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
+ abort_event::abort_event(Event const & ev):
+ std::invalid_argument{
+ "Encountered `" + name(ev.type()) + "` event:\n"
+ + to_string(ev)
+ } {}
+
} // namespace HEJ
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 6cb3d17..5d301af 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,603 +1,608 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/YAMLreader.hh"
#include <algorithm>
#include <iostream>
#include <limits>
#include <map>
#include <string>
#include <unordered_map>
#include <vector>
#include <dlfcn.h>
#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", "off-shell tolerance",
"require low pt jet"
};
// 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"}){
+ for(
+ auto && opt: {
+ "FKL", "unordered", "extremal qqbar", "central qqbar", "non-resummable",
+ "unknown", "invalid"
+ }){
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<std::string, fastjet::JetAlgorithm> 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<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
- {"discard", EventTreatment::discard}
+ {"discard", EventTreatment::discard},
+ {"abort", EventTreatment::abort}
};
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<std::string>());
}
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
setting = to_EventTreatment(yaml.as<std::string>());
}
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
setting = to_ParticleID(yaml.as<std::string>());
}
void set_from_yaml(WeightType & setting, YAML::Node const & yaml){
setting = to_weight_type(yaml.as<std::string>());
}
} // 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<double>::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<std::string, FileFormat> 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<std::string>();
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<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
}
bool convert<HEJ::OutputFile>::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<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
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<OutputFile>();
}
}
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<std::string> const & scale_names,
std::unordered_map<std::string, EventScale> & 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<EventScale>(sym)); // NOLINT
}
}
auto get_scale_map(
YAML::Node const & yaml
) {
std::unordered_map<std::string, EventScale> 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<std::string>();
const auto scale_names =
import.second.IsSequence()
?import.second.as<std::vector<std::string>>()
:std::vector<std::string>{import.second.as<std::string>()};
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<std::string, EventScale> 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<std::string, EventScale> 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 {
- {invalid, EventTreatment::discard},
- {unknown, EventTreatment::discard},
{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}
+ {non_resummable, EventTreatment::discard},
+ {unknown, EventTreatment::abort},
+ {invalid, EventTreatment::abort}
};
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");
set_from_yaml(treat.at(unknown), node, entry, "unknown");
set_from_yaml(treat.at(invalid), node, entry, "invalid");
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.lowpt, yaml, "require low pt jet");
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<std::string> 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<double>::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/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
index f94e4e4..aec8ab8 100644
--- a/t/test_scale_arithmetics.cc
+++ b/t/test_scale_arithmetics.cc
@@ -1,104 +1,105 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <memory>
#include <cstdlib>
#include <utility>
#include "HEJ/Config.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/RNG.hh"
#include "HEJ/stream.hh"
#include "HEJ/YAMLreader.hh"
#include "LHEF/LHEF.h"
namespace {
constexpr double EP = 1e-13;
void dump(HEJ::Event const & ev){
{
LHEF::Writer writer{std::cout};
std::cout << std::setprecision(6);
writer.hepeup = to_HEPEUP(ev, nullptr);
writer.writeEvent();
}
std::cout << "Rapidity ordering:\n";
for(auto const & part: ev.outgoing()){
std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
}
}
}
int main(int argn, char** argv){
if(argn != 3){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
return EXIT_FAILURE;
}
HEJ::Config config = HEJ::load_config(argv[1]);
config.scales = HEJ::to_ScaleConfig(
YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]")
);
+ config.treat[HEJ::event_type::unknown] = HEJ::EventTreatment::keep;
HEJ::istream in{argv[2]};
LHEF::Reader reader{in};
std::shared_ptr<HEJ::RNG> ran{
HEJ::make_RNG(config.rng.name, config.rng.seed)};
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
HEJ::EventReweighter resum{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
ran
};
std::size_t i = 0;
while(reader.readEvent()){
++i;
HEJ::Event::EventData data{reader.hepeup};
shuffle_particles(data);
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
auto resummed = resum.reweight(event, config.trials);
for(auto && ev: resummed) {
for(auto &&var: ev.variations()) {
if(std::abs(var.muf - ev.central().muf) > EP) {
std::cerr
<< std::setprecision(15)
<< "unequal scales: " << var.muf
<< " != " << ev.central().muf << '\n'
<< "in resummed event:\n";
dump(ev);
std::cerr << "\noriginal event:\n";
dump(event);
return EXIT_FAILURE;
}
}
}
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:40 PM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805642
Default Alt Text
(74 KB)

Event Timeline