Page MenuHomeHEPForge
Files F7877375

No OneTemporary
Actions

View Options

diff --git a/Changes.md b/Changes.md
index b46b5f3..04fb44a 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,69 +1,72 @@
# Changelog
This is the log for changes to the HEJ program. Further changes to the HEJ API
are documented in `Changes-API.md`. If you are using HEJ as a library, please
also read the changes there.
## Version 2.X
### 2.X.0
* Resummation for W bosons with jets
- New subleading processes `extremal qqx` & `central qqx` 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 mutliple jets together with a (off-shell) W bosons
decaying into lepton & neutrino
* Resummation can now be performed on `unordered` subleading processes
in pure jets.
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Print cross sections at end of run
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
* 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
* 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
- FOG: decays are now specified in `decays` setting (previously under
`particle properties`)
* Allow changing the regulator lambda in input (`regulator parameter`, only for
advanced users)
+* Use relative fraction for soft transverse momentum in extremal jets (`max ext
+ soft pt fraction`). This supersedes `min extparton pt`, which is now
+ deprecated and will be removed in future versions.
* Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`)
* Added support to read `hdf5` event files suggested in
[arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
[HighFive](https://github.com/BlueBrain/HighFive))
* Support input with avarage weight equal to the cross section (`IDWTUP=1 or 4`)
* 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
* Optional setting to specify maximal number of Fixed Order events (`max
events`, default is all)
## 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.
diff --git a/config.yml b/config.yml
index ff01962..8252b87 100644
--- a/config.yml
+++ b/config.yml
@@ -1,116 +1,117 @@
# number of attempted resummation phase space points for each input event
trials: 10
-min extparton pt: 30 # minimum transverse momentum of extremal partons
-
# maximum soft transverse momentum fraction in extremal jets
-#
-# max ext soft pt fraction: 0.1
+max ext soft pt fraction: 0.1
resummation jets: # resummation jet properties
min pt: 35 # minimum jet transverse momentum
algorithm: antikt # jet clustering algorithm
R: 0.4 # jet R parameter
fixed order jets: # properties of input jets
min pt: 30
# 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 qqx: keep
central qqx: keep
non-resummable: keep
# central scale choice or choices
#
# 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]
log correction: false # whether or not to include higher order logs
# event output files
#
# the supported formats are
# - Les Houches (suffix .lhe)
# - HepMC2 (suffix .hepmc)
# - HepMC3 (suffix .hepmc3 or .hepmc)
# TODO: - ROOT ntuples (suffix .root)
#
# 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
# to use a rivet analysis
#
# analysis:
# rivet: MC_XS # rivet analysis name
# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
#
# to use a custom analysis
#
# 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
# 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
## ---------------------------------------------------------------------- ##
## The following settings are only intended for advances users. ##
## Please DO NOT SET them unless you know exactly what you are doing! ##
## ---------------------------------------------------------------------- ##
#
+# minimum transverse momentum of extremal partons
+# deprecated: use "max ext soft pt fraction" instead
+# min extparton pt: 30
+#
# 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 a7bd4ec..b85d499 100644
--- a/doc/sphinx/HEJ.rst
+++ b/doc/sphinx/HEJ.rst
@@ -1,365 +1,365 @@
.. _`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 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
..
- *TODO* Production of a Z boson or photon with jets
where at least two jets are required in each case. For the time being,
only leading-order 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.
-.. _`min extparton pt`:
-
-**min extparton pt**
- Specifies the minimum transverse momentum in GeV of the most forward
- and the most backward parton. This setting is needed to regulate an
- otherwise uncancelled divergence. Its value should be slightly below
- the minimum transverse momentum of jets specified by `resummation
- jets: min pt`_. See also the `max ext soft pt fraction`_ setting.
-
.. _`max ext soft pt fraction`:
**max ext soft pt fraction**
- Specifies the maximum fraction that soft radiation can contribute to
- the transverse momentum of each the most forward and the most backward
- jet. Values between around 0.05 and 0.1 are recommended. See also the
- `min extparton pt`_ setting.
+ Specifies the maximum fraction that soft radiation can contribute to the
+ transverse momentum of each the most forward and the most backward jet. This
+ setting is needed to regulate an otherwise uncancelled divergence. Values
+ between around 0.05 and 0.1 are recommended.
.. _`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`_.
.. _`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.
.. _`FKL`:
**FKL**
Specifies how to treat events respecting FKL rapidity ordering. These
configurations are dominant in the high-energy limit.
.. _`unordered`:
**unordered**
Specifies how to treat events with one 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 qqx`:
**extremal qqx**
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. :code:`reweight` is currently only supported for W boson
plus jets production.
.. _`central qqx`:
**central qqx**
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. :code:`reweight` is currently only supported for W boson
plus jets production.
.. _`non-resummable`:
**non-resummable**
Specifies how to treat events where no resummation is possible. 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`.
.. _`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.
.. _`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`.
.. _`analysis`:
**analysis**
Name and Setting for the event analyses; either a custom
analysis plugin or Rivet. For the first the :code:`plugin` sub-entry
should be set to the analysis file path. All further entries are passed on
to the analysis. To use Rivet a list of Rivet analyses have to be
given in :code:`rivet` and prefix for the yoda file has to be set
through :code:`output`. 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.
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.
+.. _`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 `max ext soft pt fraction`_ 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 ad62d69..79f5089 100644
--- a/include/HEJ/config.hh
+++ b/include/HEJ/config.hh
@@ -1,200 +1,204 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <string>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "HEJ/Constants.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/ScaleFunction.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
optional<std::string> seed;
};
/**! 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<event_type::EventType, EventTreatment>;
/**! 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
+ //! %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
- double min_extparton_pt;
+ //! \deprecated This will be removed in future versions.
+ //! Use \ref max_ext_soft_pt_fraction instead.
+ double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
//! Number of resummation configurations to generate per fixed-order event
int trials;
//! Maximal number of events
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 analyses
+ //! %Parameters for custom analyses
YAML::Node analysis_parameters;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
};
//! Configuration options for the PhaseSpacePoint class
struct PhaseSpacePointConfig {
//! Properties of resummation jets
JetParameters jet_param;
//! Minimum transverse momentum for extremal partons
- double min_extparton_pt;
+ //! \deprecated This will be removed in future versions.
+ //! Use \ref max_ext_soft_pt_fraction instead.
+ double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
};
//! Configuration options for the MatrixElement class
struct MatrixElementConfig {
MatrixElementConfig() = default;
MatrixElementConfig(
bool log_correction,
HiggsCouplingSettings Higgs_coupling,
EWConstants ew_parameters,
double regulator_lambda = CLAMBDA
):
log_correction{log_correction},
Higgs_coupling{Higgs_coupling},
ew_parameters{ew_parameters},
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;
//! 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;
//! Properties of resummation jets
JetParameters 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.min_extparton_pt,
conf.max_ext_soft_pt_fraction
};
}
/**! 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.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.resummation_jets, conf.treat
};
}
} // namespace HEJ
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 9a7922d..e0da9d1 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,510 +1,512 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \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/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/Constants.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",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "event output", "analysis", "vev",
"regulator parameter", "max events"
};
// 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: {"FKL", "unordered", "extremal qqx", "central qqx", "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] = "";
}
}
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}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment \"" + name + "\"");
}
return res->second;
}
} // namespace anonymous
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>());
}
} // 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;
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;
}
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;
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]){
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}
};
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 == filename.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;
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 analysis sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(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;
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
) {
auto 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));
}
}
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"]) {
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); }
);
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 == scale_fun.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);
}
else{
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 {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::discard},
{unob, EventTreatment::discard},
{unof, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, 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(qqxexb), node, entry, "extremal qqx");
treat.at(qqxexf) = treat.at(qqxexb);
set_from_yaml(treat.at(qqxmid), node, entry, "central qqx");
set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable");
if(treat[non_resummable] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight Fixed Order 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(config.min_extparton_pt, yaml, "min extparton pt");
+ 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 use \"max ext soft pt fraction\" instead.\n";
+ set_from_yaml(
+ config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
+ );
// 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");
- config.max_ext_soft_pt_fraction = 0.1;
- set_from_yaml_if_defined(
- config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
- );
set_from_yaml(config.trials, yaml, "trials");
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.analysis_parameters, yaml, "analysis");
config.scales = to_ScaleConfig(yaml);
config.ew_parameters = get_ew_parameters(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace anonymous
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/ME_data/config_mt.yml b/t/ME_data/config_mt.yml
index 25be9fd..5e1be8d 100644
--- a/t/ME_data/config_mt.yml
+++ b/t/ME_data/config_mt.yml
@@ -1,44 +1,44 @@
trials: 1
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.385
width: 2.085
Z:
mass: 91.187
width: 2.495
Higgs coupling:
use impact factors: false
mt: 174
include bottom: false
diff --git a/t/ME_data/config_mtinf.yml b/t/ME_data/config_mtinf.yml
index f1a8f7b..3f87d26 100644
--- a/t/ME_data/config_mtinf.yml
+++ b/t/ME_data/config_mtinf.yml
@@ -1,39 +1,39 @@
trials: 1
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
non-resummable: discard
scales: 125
log correction: false
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.385
width: 2.085
Z:
mass: 91.187
width: 2.495
random generator:
name: mixmax
seed: 1
diff --git a/t/ME_data/config_mtmb.yml b/t/ME_data/config_mtmb.yml
index eb0075a..924f524 100644
--- a/t/ME_data/config_mtmb.yml
+++ b/t/ME_data/config_mtmb.yml
@@ -1,45 +1,45 @@
trials: 1
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.385
width: 2.085
Z:
mass: 91.187
width: 2.495
Higgs coupling:
use impact factors: false
mt: 174
include bottom: true
mb: 4.7
diff --git a/t/ME_data/config_w_ME.yml b/t/ME_data/config_w_ME.yml
index 27d826b..d0da43a 100644
--- a/t/ME_data/config_w_ME.yml
+++ b/t/ME_data/config_w_ME.yml
@@ -1,39 +1,39 @@
trials: 1
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: reweight
extremal qqx: reweight
central qqx: reweight
non-resummable: discard
scales: 125
log correction: false
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.419
width: 2.0476
Z:
mass: 91.187
width: 2.495
random generator:
name: mixmax
seed: 1
diff --git a/t/check_res.cc b/t/check_res.cc
index d33e50b..e3b8fad 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,171 +1,169 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <cmath>
#include <iostream>
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/stream.hh"
#define ASSERT(x) if(!(x)) { \
std::cerr << "Assertion '" #x "' failed.\n"; \
return EXIT_FAILURE; \
}
namespace{
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
const fastjet::JetDefinition Born_jet_def{jet_def};
constexpr double Born_jetptmin = 30;
- constexpr double extpartonptmin = 30;
constexpr double max_ext_soft_pt_fraction = 0.1;
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
using EventTreatment = HEJ::EventTreatment;
using namespace HEJ::event_type;
HEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{non_resummable, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
if(!HEJ::event_type::is_resummable(ev.type()))
return true;
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
};
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "unof"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
if(argn == 5 && std::string(argv[4]) == "unob"){
--argn;
treat[unof] = EventTreatment::discard;
treat[unob] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitf"){
--argn;
treat[qqxexb] = EventTreatment::discard;
treat[qqxexf] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitb"){
--argn;
treat[qqxexb] = EventTreatment::reweight;
treat[qqxexf] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "qqxmid"){
--argn;
treat[qqxmid] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
if(argn != 4){
std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
return EXIT_FAILURE;
}
const double xsec_ref = std::stod(argv[2]);
const double tolerance = std::stod(argv[3]);
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
HEJ::PhaseSpacePointConfig psp_conf;
psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
- psp_conf.min_extparton_pt = extpartonptmin;
psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
HEJ::MatrixElementConfig ME_conf;
ME_conf.log_correction = log_corr;
ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
ME_conf.ew_parameters.set_vevWZH(vev, Wprop, Zprop, Hprop);
HEJ::EventReweighterConfig conf;
conf.psp_config = std::move(psp_conf);
conf.ME_config = std::move(ME_conf);
conf.jet_param = psp_conf.jet_param;
conf.treat = treat;
reader.readEvent();
const bool has_Higgs = std::find(
begin(reader.hepeup.IDUP),
end(reader.hepeup.IDUP),
25
) != end(reader.hepeup.IDUP);
const double mu = has_Higgs?125.:91.188;
HEJ::ScaleGenerator scale_gen{
{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
};
HEJ::Mixmax ran{};
HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};
double xsec = 0.;
double xsec_err = 0.;
do{
auto ev_data = HEJ::Event::EventData{reader.hepeup};
ev_data.reconstruct_intermediate();
HEJ::Event ev{
ev_data.cluster(
Born_jet_def, Born_jetptmin
)
};
auto resummed_events = hej.reweight(ev, 100);
for(auto const & ev: resummed_events) {
ASSERT(correct_colour(ev));
ASSERT(std::isfinite(ev.central().weight));
xsec += ev.central().weight;
xsec_err += ev.central().weight*ev.central().weight;
}
} while(reader.readEvent());
xsec_err = std::sqrt(xsec_err);
const double significance =
std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance );
std::cout << xsec_ref << " +/- " << tolerance << " ~ "
<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";
if(significance > 3.){
std::cerr << "Cross section is off by over 3 sigma!\n";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/t/jet_config.yml b/t/jet_config.yml
index 9d7dc31..79e0363 100644
--- a/t/jet_config.yml
+++ b/t/jet_config.yml
@@ -1,41 +1,41 @@
trials: 10
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: keep
extremal qqx: discard
central qqx: keep
non-resummable: keep
log correction: false
scales: 91.188
random generator:
name: ranlux64
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.385
width: 2.085
Z:
mass: 91.187
width: 2.495
event output:
- tst.lhe
diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml
index 446a58c..027a2e8 100644
--- a/t/jet_config_with_import.yml
+++ b/t/jet_config_with_import.yml
@@ -1,44 +1,44 @@
trials: 10
-min extparton pt: 30
+max ext soft pt fraction: 0.1
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
event treatment:
FKL: reweight
unordered: discard
extremal qqx: discard
central qqx: discard
non-resummable: discard
log correction: false
scales: softest_jet_pt
event output:
- tst.lhe
random generator:
name: ranlux64
vev: 246.2196508
particle properties:
Higgs:
mass: 125
width: 0.004165
W:
mass: 80.385
width: 2.085
Z:
mass: 91.187
width: 2.495
import scales:
./libscales.so: softest_jet_pt
diff --git a/t/test_psp.cc b/t/test_psp.cc
index 92b6a2b..105bc49 100644
--- a/t/test_psp.cc
+++ b/t/test_psp.cc
@@ -1,68 +1,65 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "LHEF/LHEF.h"
#include "HEJ/stream.hh"
#include "HEJ/config.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/PhaseSpacePoint.hh"
#include "HEJ/Ranlux64.hh"
namespace{
constexpr int max_trials = 100;
- constexpr double extpartonptmin = 45.;
- constexpr double max_ext_soft_pt_fraction =
- std::numeric_limits<double>::infinity();
+ constexpr double max_ext_soft_pt_fraction = 0.1;
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
constexpr double min_jet_pt = 50;
};
int main(int argn, char** argv) {
if(argn != 2){
std::cerr << "Usage: " << argv[0] << " eventfile";
return EXIT_FAILURE;
}
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
LHEF::Writer writer{std::cerr};
writer.heprup = reader.heprup;
HEJ::PhaseSpacePointConfig conf;
conf.jet_param = HEJ::JetParameters{jet_def, min_jet_pt};
- conf.min_extparton_pt = extpartonptmin;
conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
HEJ::Ranlux64 ran{};
while(reader.readEvent()){
const HEJ::Event ev{
HEJ::Event::EventData{reader.hepeup}( jet_def, min_jet_pt )
};
for(int trial = 0; trial < max_trials; ++trial){
HEJ::PhaseSpacePoint psp{ev, conf, ran};
if(psp.weight() != 0){
HEJ::Event::EventData tmp_ev;
tmp_ev.incoming = psp.incoming();
tmp_ev.outgoing = psp.outgoing();
tmp_ev.parameters.central = {0,0,0};
HEJ::Event out_ev{ tmp_ev(jet_def, min_jet_pt) };
if(out_ev.type() != ev.type()){
using HEJ::event_type::name;
std::cerr << "Wrong class of phase space point:\n"
"original event of class " << name(ev.type()) << ":\n";
writer.hepeup = reader.hepeup;
writer.writeEvent();
std::cerr << "Phase space point of class " << name(out_ev.type()) << ":\n";
writer.hepeup = to_HEPEUP(out_ev, &writer.heprup);
writer.writeEvent();
return EXIT_FAILURE;
}
}
}
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:30 PM (1 d, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804956
Default Alt Text
(58 KB)

Event Timeline