diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
index 3188a76..1c82bd7 100644
--- a/FixedOrderGen/src/config.cc
+++ b/FixedOrderGen/src/config.cc
@@ -1,353 +1,357 @@
#include "config.hh"
#include <cctype>
#include "Subleading.hh"
#include "HEJ/config.hh"
#include "HEJ/YAMLreader.hh"
namespace HEJFOG{
using HEJ::set_from_yaml;
using HEJ::set_from_yaml_if_defined;
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 = {
"process", "events", "subleading fraction","subleading channels",
"scales", "scale factors", "max scale ratio", "pdf",
"event output", "analysis", "import scales"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "peak pt", "algorithm", "R", "max rapidity"}){
supported["jets"][jet_opt] = "";
}
for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){
for(auto && particle_opt: {"mass", "width"}){
supported["particle properties"][particle_type][particle_opt] = "";
}
supported["particle properties"][particle_type]["decays"]["into"] = "";
supported["particle properties"][particle_type]["decays"]["branching ratio"] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && beam_opt: {"energy", "particles"}){
supported["beam"][beam_opt] = "";
}
for(auto && unweight_opt: {"sample size", "max deviation"}){
supported["unweight"][unweight_opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
return supported;
}();
return supported;
}
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
){
const auto p = HEJ::get_jet_parameters(node, entry);
JetParameters result;
result.def = p.def;
result.min_pt = p.min_pt;
set_from_yaml(result.max_y, node, entry, "max rapidity");
set_from_yaml_if_defined(result.peak_pt, node, entry, "peak pt");
+ if(result.peak_pt && *result.peak_pt <= result.min_pt)
+ throw std::invalid_argument{
+ "Value of option 'peak pt' has to be larger than 'min pt'."
+ };
return result;
}
Beam get_Beam(
YAML::Node const & node, std::string const & entry
){
Beam beam;
std::vector<HEJ::ParticleID> particles;
set_from_yaml(beam.energy, node, entry, "energy");
set_from_yaml_if_defined(particles, node, entry, "particles");
if(! particles.empty()){
for(HEJ::ParticleID particle: particles){
if(particle != HEJ::pid::p && particle != HEJ::pid::p_bar){
throw std::invalid_argument{
"Unsupported value in option " + entry + ": particles:"
" only proton ('p') and antiproton ('p_bar') beams are supported"
};
}
}
if(particles.size() != 2){
throw std::invalid_argument{"Not exactly two beam particles"};
}
beam.particles.front() = particles.front();
beam.particles.back() = particles.back();
}
return beam;
}
std::vector<std::string> split(
std::string const & str, std::string const & delims
){
std::vector<std::string> result;
for(size_t begin, end = 0; end != str.npos;){
begin = str.find_first_not_of(delims, end);
if(begin == str.npos) break;
end = str.find_first_of(delims, begin + 1);
result.emplace_back(str.substr(begin, end - begin));
}
return result;
}
std::invalid_argument invalid_incoming(std::string const & what){
return std::invalid_argument{
"Incoming particle type " + what + " not supported,"
" incoming particles have to be 'p', 'p_bar' or partons"
};
}
std::invalid_argument invalid_outgoing(std::string const & what){
return std::invalid_argument{
"Outgoing particle type " + what + " not supported,"
" outgoing particles have to be 'j', 'photon', 'W+', 'W-', 'Z', 'H'"
};
}
Process get_process(
YAML::Node const & node, std::string const & entry
){
Process result;
std::string process_string;
set_from_yaml(process_string, node, entry);
assert(! process_string.empty());
const auto particles = split(process_string, " \n\t\v=>");
if(particles.size() < 3){
throw std::invalid_argument{
"Bad format in option process: '" + process_string
+ "', expected format is 'in1 in2 => out1 ...'"
};
}
result.incoming.front() = HEJ::to_ParticleID(particles[0]);
result.incoming.back() = HEJ::to_ParticleID(particles[1]);
for(size_t i = 0; i < result.incoming.size(); ++i){
const HEJ::ParticleID in = result.incoming[i];
if(
in != HEJ::pid::proton && in != HEJ::pid::p_bar
&& !HEJ::is_parton(in)
){
throw invalid_incoming(particles[i]);
}
}
result.njets = 0;
for(size_t i = result.incoming.size(); i < particles.size(); ++i){
assert(! particles[i].empty());
if(particles[i] == "j") ++result.njets;
else if(std::isdigit(particles[i].front())
&& particles[i].back() == 'j')
result.njets += std::stoi(particles[i]);
else{
const auto pid = HEJ::to_ParticleID(particles[i]);
if(!HEJ::is_AWZH_boson(pid)){
throw invalid_outgoing(particles[i]);
}
if(result.boson){
throw std::invalid_argument{
"More than one outgoing boson is not supported"
};
}
result.boson = pid;
}
}
if(result.njets < 2){
throw std::invalid_argument{
"Process has to include at least two jets ('j')"
};
}
return result;
}
HEJFOG::Subleading to_subleading_channel(YAML::Node const & yaml){
std::string name;
using HEJFOG::Subleading;
set_from_yaml(name, yaml);
if(name == "none")
return none;
if(name == "all")
return all;
if(name == "unordered" || name == "uno")
return uno;
throw HEJ::unknown_option("Unknown subleading channel '"+name+"'");
}
unsigned int get_subleading_channels(YAML::Node const & node){
using YAML::NodeType;
using HEJFOG::Subleading;
// all channels allowed by default
if(!node) return all;
switch(node.Type()){
case NodeType::Undefined:
return all;
case NodeType::Null:
return none;
case NodeType::Scalar:
return to_subleading_channel(node);
case NodeType::Map:
throw HEJ::invalid_type{"map is not a valid option for subleading channels"};
case NodeType::Sequence:
unsigned int channels = HEJFOG::Subleading::none;
for(auto && channel_node: node){
channels |= get_subleading_channels(channel_node);
}
return channels;
}
throw std::logic_error{"unreachable"};
}
Decay get_decay(YAML::Node const & node){
Decay decay;
set_from_yaml(decay.products, node, "into");
set_from_yaml(decay.branching_ratio, node, "branching ratio");
return decay;
}
std::vector<Decay> get_decays(YAML::Node const & node){
using YAML::NodeType;
if(!node) return {};
switch(node.Type()){
case NodeType::Null:
case NodeType::Undefined:
return {};
case NodeType::Scalar:
throw HEJ::invalid_type{"value is not a list of decays"};
case NodeType::Map:
return {get_decay(node)};
case NodeType::Sequence:
std::vector<Decay> result;
for(auto && decay_str: node){
result.emplace_back();
set_from_yaml(result.back().products, decay_str, "into");
set_from_yaml(result.back().branching_ratio, decay_str, "branching ratio");
}
return result;
}
throw std::logic_error{"unreachable"};
}
ParticleProperties get_particle_properties(
YAML::Node const & node, std::string const & entry
){
ParticleProperties result;
set_from_yaml(result.mass, node, entry, "mass");
set_from_yaml(result.width, node, entry, "width");
try{
result.decays = get_decays(node[entry]["decays"]);
}
catch(HEJ::missing_option const & ex){
throw HEJ::missing_option{entry + ": decays: " + ex.what()};
}
catch(HEJ::invalid_type const & ex){
throw HEJ::invalid_type{entry + ": decays: " + ex.what()};
}
return result;
}
ParticlesPropMap get_all_particles_properties(YAML::Node const & node){
ParticlesPropMap result;
for(auto const & entry: node) {
const auto name = entry.first.as<std::string>();
const auto id = HEJ::to_ParticleID(name);
result.emplace(id, get_particle_properties(node,name));
}
return result;
}
UnweightSettings get_unweight(
YAML::Node const & node, std::string const & entry
){
UnweightSettings result;
set_from_yaml(result.sample_size, node, entry, "sample size");
if(result.sample_size <= 0){
throw std::invalid_argument{
"negative sample size " + std::to_string(result.sample_size)
};
}
set_from_yaml(result.max_dev, node, entry, "max deviation");
return result;
}
Config to_Config(YAML::Node const & yaml){
try{
HEJ::assert_all_options_known(yaml, get_supported_options());
}
catch(HEJ::unknown_option const & ex){
throw HEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.process = get_process(yaml, "process");
set_from_yaml(config.events, yaml, "events");
config.jets = get_jet_parameters(yaml, "jets");
config.beam = get_Beam(yaml, "beam");
for(size_t i = 0; i < config.process.incoming.size(); ++i){
const auto & in = config.process.incoming[i];
using namespace HEJ::pid;
if( (in == p || in == p_bar) && in != config.beam.particles[i]){
throw std::invalid_argument{
"Particle type of beam " + std::to_string(i+1) + " incompatible"
+ " with type of incoming particle " + std::to_string(i+1)
};
}
}
set_from_yaml(config.pdf_id, yaml, "pdf");
set_from_yaml(config.subleading_fraction, yaml, "subleading fraction");
if(config.subleading_fraction < 0 || config.subleading_fraction > 1){
throw std::invalid_argument{
"subleading fraction has to be between 0 and 1"
};
}
if(config.subleading_fraction == 0)
config.subleading_channels = Subleading::none;
else
config.subleading_channels = get_subleading_channels(yaml["subleading channels"]);
if(!config.process.boson && config.subleading_channels != Subleading::none)
throw HEJ::not_implemented("Subleading processes for pure Jet production not implemented yet");
if(yaml["particle properties"]){
config.particles_properties = get_all_particles_properties(
yaml["particle properties"]);
}
if(config.process.boson
&& config.particles_properties.find(*(config.process.boson))
== config.particles_properties.end())
throw HEJ::missing_option("Process wants to generate boson "
+std::to_string(*(config.process.boson))+", but particle properties are missing");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scales = HEJ::to_ScaleConfig(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = HEJ::to_RNGConfig(yaml, "random generator");
config.Higgs_coupling = HEJ::get_Higgs_coupling(yaml, "Higgs coupling");
if(yaml["unweight"]) config.unweight = get_unweight(yaml, "unweight");
return config;
}
} // namespace anonymous
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;
}
}
}
diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
index 91260ab..569295e 100644
--- a/doc/sphinx/HEJFOG.rst
+++ b/doc/sphinx/HEJFOG.rst
@@ -1,291 +1,292 @@
The HEJ Fixed Order Generator
=============================
For high jet multiplicities event generation with standard fixed-order
generators becomes increasingly cumbersome. For example, the
leading-order production of a Higgs Boson with five or more jets is
computationally prohibitively expensive.
To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
that allows to generate events with high jet multiplicities. To
facilitate the computation the limit of Multi-Regge Kinematics with
large invariant masses between all outgoing particles is assumed in the
matrix elements. The typical use of the ``HEJFOG`` is to supplement
low-multiplicity events from standard generators with high-multiplicity
events before using the HEJ 2 program to add high-energy
resummation.
Installation
------------
The ``HEJFOG`` comes bundled together with HEJ 2 and the
installation is very similar. After downloading HEJ 2 and
installing the prerequisites as described in :ref:`Installation` the
``HEJFOG`` can be installed with::
cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory
make install
where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
subdirectory in the HEJ 2 directory. If HEJ 2 was
installed to a non-standard location, it may be necessary to specify the
directory containing :file:`HEJ-config.cmake`. If the base installation
directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
installing the ``HEJFOG`` would read::
cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory
make install
The installation can be tested with::
make test
provided that the CT10nlo PDF set is installed.
Running the fixed-order generator
---------------------------------
After installing the ``HEJFOG`` you can modify the provided
configuration file :file:`configFO.yml` and run the generator with::
HEJFOG configFO.yml
The resulting event file, by default named :file:`HEJFO.lhe`, can then be
fed into HEJ 2 like any event file generated from a standard
fixed-order generator, see :ref:`Running HEJ 2`.
Settings
--------
Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
<http://yaml.org/>`_ configuration file. The settings are
.. _`process`:
**process**
The scattering process for which events are being generated. The
format is
:code:`in1 in2 => out1 out2 ...`
The incoming particles, :code:`in1`, :code:`in2` can be
- quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
- gluons: :code:`g`
- protons :code:`p` or antiprotons :code:`p_bar`
At most one of the outgoing particles can be a boson, the rest has to be
partonic. At the moment only the Higgs boson :code:`h` is supported. All
other outgoing particles are jets. Multiple jets can be grouped together, so
:code:`p p => h j j` is the same as :code:`p p => h 2j`. There have to be at
least two jets. Further decays of the boson can be added through the
:ref:`particle properties<particle properties: particle: decays>`.
.. _`events`:
**events**
Specifies the number of events to generate.
.. _`jets`:
**jets**
Defines the properties of the generated jets.
.. _`jets: min pt`:
**min pt**
Minimum jet transverse momentum in GeV.
.. _`jets: peak pt`:
**peak pt**
Optional setting to specify the dominant jet transverse momentum
in GeV. If the generated events are used as input for HEJ
resummation, this should be set to the minimum transverse momentum
- of the resummation jets. The effect is that only a small fraction
- of jets will be generated with a transverse momentum below the
- value of this setting.
+ of the resummation jets. In all cases it has to be larger than
+ :code:`min pt`. The effect is that only a small fraction of jets
+ will be generated with a transverse momentum below the value of
+ this setting.
.. _`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.
.. _`jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`jets: max rapidity`:
**max rapidity**
Maximum absolute value of the jet rapidity.
.. _`beam`:
**beam**
Defines various properties of the collider beam.
.. _`beam: energy`:
**energy**
The beam energy in GeV. For example, the 13
TeV LHC corresponds to a value of 6500.
.. _`beam: particles`:
**particles**
A list :code:`[p1, p2]` of two beam particles. The only supported
entries are protons :code:`p` and antiprotons :code:`p_bar`.
.. _`pdf`:
**pdf**
The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.
.. _`subleading fraction`:
**subleading fraction**
This setting is related to the fraction of events that are not FKL
configurations and thus subleading in the high-energy
limit. Currently only unordered emissions are implemented, and only
for Higgs boson plus multijet processes. This value must be positive
and not larger than 1. It should typically be chosen between 0.01 and
0.1. Note that while this parameter influences the chance of
generating subleading configurations, it generally does not
correspond to the actual fraction of subleading events.
.. _`subleading channels`:
**subleading channels**
Optional parameter to select the production of specific channels that
are subleading in the high-energy limit. Only has an effect if
:code:`subleading fraction` is non-zero. Currently three values are
supported:
- :code:`all`: All subleading channels are allowed. This is the default.
- :code:`none`: No subleading contribution, only FKL configurations
are allowed. This is equivalent to :code:`subleading fraction: 0`.
- :code:`unordered`: Unordered emission are allowed.
Unordered emission are any rapidity ordering where exactly one
gluon is emitted outside the FKL rapidity ordering. More
precisely, if at least one of the incoming particles is a quark or
antiquark and there are more than two jets in the final state,
:code:`subleading fraction` states the probability that the
flavours of the outgoing particles are assigned in such a way that
an unordered configuration arises.
.. _`unweight`:
**unweight**
This setting defines the parameters for the partial unweighting of
events. You can disable unweighting by removing this entry from the
configuration file.
.. _`unweight: sample size`:
**sample size**
The number of weighted events used to calibrate the unweighting.
A good default is to set this to the number of target
`events`_. If the number of `events`_ is large this can
lead to significant memory consumption and a lower value should be
chosen. Contrarily, for large multiplicities the unweighting
efficiency becomes worse and the sample size should be increased.
.. _`unweight: max deviation`:
**max deviation**
Controls the range of events to which unweighting is applied. A
larger value means that a larger fraction of events are unweighted.
Typical values are between -1 and 1.
.. _`particle properties`:
**particle properties**
Specifies various properties of the different particles (Higgs, W or Z).
This is only relevant if the chosen `process`_ is the production of the
corresponding particles with jets. E.g. for the `process`_
:code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
:code:`decays` of the :code:`Higgs` boson are required, while all other
particle properties will be ignored.
.. _`particle properties: particle`:
**Higgs, W+, W- or Z**
The particle (Higgs, |W+|, |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.
.. _`particle properties: particle: decays`:
**decays**
Optional setting specifying the decays of the particle. Only the decay
into two particles is implemented. Each decay has the form
:code:`{into: [p1,p2], branching ratio: r}`
where :code:`p1` and :code:`p2` are the particle names of the
decay product (e.g. :code:`photon`) and :code:`r` is the branching
ratio.
Decays of a Higgs boson are treated as the production and subsequent
decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
supported.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. For details, see the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings`. Note that this should usually be a
single value, as the weights resulting from additional scale choices
will simply be ignored when adding high-energy resummation with HEJ 2.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. See the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`RanLux init`:
.. _`random generator`:
**random generator**
Sets parameters for random number generation. See the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`analysis`:
**analysis**
Specifies the name and settings for a custom analysis library. This
can be useful to specify cuts at the fixed-order level. See the
corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. See the corresponding entry in the
HEJ 2 :ref:`HEJ 2 settings` for details.