Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F10664257
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
42 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
index fe4b9dc..38a4c88 100644
--- a/FixedOrderGen/src/config.cc
+++ b/FixedOrderGen/src/config.cc
@@ -1,305 +1,302 @@
#include "config.hh"
#include <cctype>
#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", "unordered fraction", "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", "Wm", "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");
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())){
- if(particles[i].back() != 'j'){
- throw invalid_outgoing(particles[i]);
- }
+ 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;
}
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;
using namespace HEJ;
// @TODO allow more synonyms
if(node["Higgs"])
result[pid::Higgs] = get_particle_properties(node,"Higgs");
if(node["Wp"])
result[pid::Wp] = get_particle_properties(node,"Wp");
if(node["Wm"])
result[pid::Wm] = get_particle_properties(node,"Wm");
if(node["Z"])
result[pid::Z] = get_particle_properties(node,"Z");
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.unordered_fraction, yaml, "unordered fraction");
if(config.unordered_fraction < 0 || config.unordered_fraction > 1){
throw std::invalid_argument{
"unordered fraction has to be between 0 and 1"
};
}
config.particles_properties = get_all_particles_properties(
yaml["Particle properties"]);
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/include/HEJ/PDG_codes.hh b/include/HEJ/PDG_codes.hh
index 8b93ca3..01dc521 100644
--- a/include/HEJ/PDG_codes.hh
+++ b/include/HEJ/PDG_codes.hh
@@ -1,92 +1,97 @@
/** \file PDG_codes.hh
* \brief Contains the Particle IDs of all relevant SM particles.
*
* Large enumeration included which has multiple entries for potential
* alternative names of different particles. There are also functions
* which can be used to determine if a particle is a parton or if
* it is a non-gluon boson.
*/
#pragma once
+#include <string>
+
namespace HEJ {
//! particle ids according to PDG
namespace pid {
//! The possible particle identities. We use PDG IDs as standard.
enum ParticleID{
d = 1, /*!< Down Quark */
down = d, /*!< Down Quark */
u = 2, /*!< Up Quark */
up = u, /*!< Up Quark */
s = 3, /*!< Strange Quark */
strange = s, /*!< Strange Quark */
c = 4, /*!< Charm Quark */
charm = c, /*!< Charm Quark */
b = 5, /*!< Bottom Quark */
bottom = b, /*!< Bottom Quark */
t = 6, /*!< Top Quark */
top = t, /*!< Top Quark */
e = 11, /*!< Electron */
electron = e, /*!< Electron */
nu_e = 12, /*!< Electron Neutrino */
electron_neutrino = nu_e, /*!< Electron neutrino */
mu = 13, /*!< Muon */
muon = mu, /*!< Muon */
nu_mu = 14, /*!< Muon Neutrino */
muon_neutrino = nu_mu, /*!< Muon Neutrino */
tau = 15, /*!< Tau */
nu_tau = 16, /*!< Tau Neutrino */
tau_neutrino = nu_tau, /*!< Tau Neutrino */
d_bar = -d, /*!< Anti-Down Quark */
u_bar = -u, /*!< Anti-Up quark */
s_bar = -s, /*!< Anti-Strange Quark */
c_bar = -c, /*!< Anti-Charm Quark */
b_bar = -b, /*!< Anti-Bottom Quark */
t_bar = -t, /*!< Anti-Top Quark */
e_bar = -e, /*!< Positron */
nu_e_bar = -nu_e, /*!< Anti-Electron Neutrino */
mu_bar = -mu, /*!< Anti-Muon */
nu_mu_bar = -nu_mu, /*!< Anti-Muon Neutrino */
tau_bar = -tau, /*!< Anti-Tau */
nu_tau_bar = -nu_tau, /*!< Anti-Tau Neutrino */
gluon = 21, /*!< Gluon */
g = gluon, /*!< Gluon */
photon = 22, /*!< Photon */
gamma = photon, /*!< Photon */
Z = 23, /*!< Z Boson */
Wp = 24, /*!< W- Boson */
Wm = -Wp, /*!< W+ Boson */
h = 25, /*!< Higgs Boson */
Higgs = h, /*!< Higgs Boson */
higgs = h, /*!< Higgs Boson */
p = 2212, /*!< Proton */
proton = p, /*!< Proton */
p_bar = -p, /*!< Anti-Proton */
};
}
using ParticleID = pid::ParticleID;
+ //! Convert a particle name to the corresponding PDG particle ID
+ ParticleID to_ParticleID(std::string const & name);
+
/**
* \brief Function to determine if particle is a parton
* @param p PDG ID of particle
* @returns true if the particle is a parton, false otherwise
*/
inline
constexpr bool is_parton(ParticleID p){
return p == pid::gluon || std::abs(p) <= pid::top;
}
/**
* \brief function to determine if the particle is a photon, W, Z, or Higgs boson
* @param id PDG ID of particle
* @returns true if the partice is a A,W,Z, or H, false otherwise
*/
inline
constexpr bool is_AWZH_boson(ParticleID id){
return id == pid::Wm || (id >= pid::photon && id <= pid::Higgs);
}
}
diff --git a/include/HEJ/YAMLreader.hh b/include/HEJ/YAMLreader.hh
index d2ed33e..8740f14 100644
--- a/include/HEJ/YAMLreader.hh
+++ b/include/HEJ/YAMLreader.hh
@@ -1,244 +1,241 @@
/** \file
* \brief The file which handles the configuration file parameters
*
* The configuration files parameters are read and then stored
* within this objects.
*/
#pragma once
#include "yaml-cpp/yaml.h"
#include "HEJ/config.hh"
#include "HEJ/exceptions.hh"
namespace HEJ{
//! Load configuration from file
/**
* @param config_file Name of the YAML configuration file
* @returns The HEJ 2 configuration
*/
Config load_config(std::string const & config_file);
//! Set option using the corresponding YAML entry
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* If the entry does not exist or has the wrong type or format
* an exception is thrown.
*
* For example
* @code
* set_from_yaml(foobar, yaml, "foo", "bar")
* @endcode
* is equivalent to
* @code
* foobar = yaml["foo"]["bar"].as<decltype(foobar)>()
* @endcode
* with improved diagnostics on errors.
*
* @see set_from_yaml_if_defined
*/
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
//! Set option using the corresponding YAML entry, if present
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* This function works similar to set_from_yaml, but does not
* throw any exception if the requested YAML entry does not exist.
*
* @see set_from_yaml
*/
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
//! Extract jet parameters from YAML configuration
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
);
//! Extract Higgs coupling settings from YAML configuration
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node, std::string const & entry
);
//! Extract scale setting parameters from YAML configuration
ScaleConfig to_ScaleConfig(YAML::Node const & yaml);
//! Extract random number generator settings from YAML configuration
RNGConfig to_RNGConfig(YAML::Node const & node, std::string const & entry);
- //! Convert a particle name to the corresponding PDG particle ID
- ParticleID to_ParticleID(std::string const & particle_name);
-
//! Check whether all options in configuration are supported
/**
* @param conf Configuration to be checked
* @param supported Tree of supported options
*
* If conf contains an entry that does not appear in supported
* an unknown_option exception is thrown. Sub-entries of "analysis"
* (if present) are not checked.
*
* @see unknown_option
*/
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
);
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml);
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml);
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml);
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml);
inline
void set_from_yaml(YAML::Node & setting, YAML::Node const & yaml){
setting = yaml;
}
template<typename Scalar>
void set_from_yaml(Scalar & setting, YAML::Node const & yaml){
assert(yaml);
if(!yaml.IsScalar()){
throw invalid_type{"value is not a scalar"};
}
try{
setting = yaml.as<Scalar>();
}
catch(...){
throw invalid_type{
"value " + yaml.as<std::string>()
+ " cannot be converted to a " + type_string(setting)
};
}
}
template<typename T>
void set_from_yaml(optional<T> & setting, YAML::Node const & yaml){
T tmp;
set_from_yaml(tmp, yaml);
setting = tmp;
}
template<typename T>
void set_from_yaml(std::vector<T> & setting, YAML::Node const & yaml){
assert(yaml);
// special case: treat a single value like a vector with one element
if(yaml.IsScalar()){
setting.resize(1);
return set_from_yaml(setting.front(), yaml);
}
if(yaml.IsSequence()){
setting.resize(yaml.size());
for(size_t i = 0; i < setting.size(); ++i){
set_from_yaml(setting[i], yaml[i]);
}
return;
}
throw invalid_type{""};
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) throw missing_option{""};
set_from_yaml(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
template<typename T>
void set_from_yaml_if_defined(T & setting, YAML::Node const & yaml){
return set_from_yaml(setting, yaml);
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) return;
set_from_yaml_if_defined(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
}
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": " + ex.what()
};
}
catch(missing_option const &){
throw missing_option{
"No entry for mandatory option " + join(": ", names...)
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml_if_defined(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": " + ex.what()
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
}
namespace YAML {
template<>
struct convert<HEJ::OutputFile> {
static Node encode(HEJ::OutputFile const & outfile);
static bool decode(Node const & node, HEJ::OutputFile & out);
};
}
diff --git a/src/PDG_codes.cc b/src/PDG_codes.cc
new file mode 100644
index 0000000..3dfc576
--- /dev/null
+++ b/src/PDG_codes.cc
@@ -0,0 +1,48 @@
+#include "HEJ/PDG_codes.hh"
+
+#include <map>
+#include <iostream>
+
+namespace HEJ{
+ ParticleID to_ParticleID(std::string const & name){
+ using namespace HEJ::pid;
+ static const std::map<std::string, ParticleID> known = {
+ {"d", d}, {"down", down}, {"1",static_cast<ParticleID>(1)},
+ {"u", u}, {"up", up}, {"2",static_cast<ParticleID>(2)},
+ {"s", s}, {"strange", strange}, {"3",static_cast<ParticleID>(3)},
+ {"c", c}, {"charm", charm}, {"4",static_cast<ParticleID>(4)},
+ {"b", b}, {"bottom", bottom}, {"5",static_cast<ParticleID>(5)},
+ {"t", t}, {"top", top}, {"6",static_cast<ParticleID>(6)},
+ {"e", e}, {"electron", electron}, {"11",static_cast<ParticleID>(11)},
+ {"nu_e", nu_e}, {"electron_neutrino", electron_neutrino}, {"12",static_cast<ParticleID>(12)},
+ {"mu", mu}, {"muon", muon}, {"13",static_cast<ParticleID>(13)},
+ {"nu_mu", nu_mu}, {"muon_neutrino", muon_neutrino}, {"14",static_cast<ParticleID>(14)},
+ {"tau", tau}, {"15",static_cast<ParticleID>(15)},
+ {"nu_tau", nu_tau}, {"tau_neutrino", tau_neutrino}, {"16",static_cast<ParticleID>(16)},
+ {"d_bar", d_bar}, {"-1",static_cast<ParticleID>(-1)},
+ {"u_bar", u_bar}, {"-2",static_cast<ParticleID>(-2)},
+ {"s_bar", s_bar}, {"-3",static_cast<ParticleID>(-3)},
+ {"c_bar", c_bar}, {"-4",static_cast<ParticleID>(-4)},
+ {"b_bar", b_bar}, {"-5",static_cast<ParticleID>(-5)},
+ {"t_bar", t_bar}, {"-6",static_cast<ParticleID>(-6)},
+ {"e_bar", e_bar}, {"-11",static_cast<ParticleID>(-11)},
+ {"nu_e_bar", nu_e_bar}, {"-12",static_cast<ParticleID>(-12)},
+ {"mu_bar", mu_bar}, {"-13",static_cast<ParticleID>(-13)},
+ {"nu_mu_bar", nu_mu_bar}, {"-14",static_cast<ParticleID>(-14)},
+ {"tau_bar", tau_bar}, {"-15",static_cast<ParticleID>(-15)},
+ {"nu_tau_bar", nu_tau_bar}, {"-16",static_cast<ParticleID>(-16)},
+ {"gluon", gluon}, {"g", g}, {"21",static_cast<ParticleID>(21)},
+ {"photon", photon}, {"gamma", gamma}, {"22",static_cast<ParticleID>(22)},
+ {"Z", Z}, {"23",static_cast<ParticleID>(23)},
+ {"Wp", Wp}, {"W+", Wp}, {"24",static_cast<ParticleID>(24)},
+ {"Wm", Wm}, {"W-", Wm}, {"-24",static_cast<ParticleID>(-24)},
+ {"h", h}, {"H", h}, {"Higgs", Higgs}, {"higgs", higgs}, {"25",static_cast<ParticleID>(25)},
+ {"p", p}, {"proton", proton}, {"p_bar", p_bar}, {"2212",static_cast<ParticleID>(2212)}
+ };
+ const auto res = known.find(name);
+ if(res == known.end()){
+ throw std::invalid_argument("Unknown particle " + name);
+ }
+ return res->second;
+ }
+}
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index ef503de..13f3888 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,486 +1,462 @@
#include "HEJ/YAMLreader.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
#include <dlfcn.h>
namespace HEJ{
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",
"FKL", "unordered", "non-HEJ",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "unweight", "event output", "analysis"
};
// 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] = "";
}
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},
{"genkt", genkt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"genkt for passive", genkt_for_passive_algorithm},
{"ee kt", ee_kt_algorithm},
{"ee genkt", ee_genkt_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
- ParticleID to_ParticleID(std::string const & name){
- using namespace HEJ::pid;
- static const std::map<std::string, ParticleID> known = {
- {"d", d}, {"down", down}, {"u", u}, {"up", up}, {"s", s}, {"strange", strange},
- {"c", c}, {"charm", charm}, {"b", b}, {"bottom", bottom}, {"t", t}, {"top", top},
- {"e", e}, {"electron", electron}, {"nu_e", nu_e}, {"electron_neutrino", electron_neutrino},
- {"mu", mu}, {"muon", muon}, {"nu_mu", nu_mu}, {"muon_neutrino", muon_neutrino},
- {"tau", tau}, {"nu_tau", nu_tau}, {"tau_neutrino", tau_neutrino},
- {"d_bar", d_bar}, {"u_bar", u_bar}, {"s_bar", s_bar}, {"c_bar", c_bar},
- {"b_bar", b_bar}, {"t_bar", t_bar}, {"e_bar", e_bar},
- {"nu_e_bar", nu_e_bar}, {"mu_bar", mu_bar}, {"nu_mu_bar", nu_mu_bar},
- {"tau_bar", tau_bar}, {"nu_tau_bar", nu_tau_bar},
- {"gluon", gluon}, {"g", g}, {"photon", photon}, {"gamma", gamma},
- {"Z", Z}, {"Wp", Wp}, {"Wm", Wm}, {"W+", Wp}, {"W-", Wm},
- {"h", h}, {"H", h}, {"Higgs", Higgs}, {"higgs", higgs},
- {"p", p}, {"proton", proton}, {"p_bar", p_bar}
- };
- const auto res = known.find(name);
- if(res == known.end()){
- throw std::invalid_argument("Unknown particle " + name);
- }
- return res->second;
- }
-
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;
}
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}
};
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;
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()))
);
const size_t delim = scale_fun.find_first_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_simple_ScaleFunction(first, known)/to_double(second);
}
else{
assert(scale_fun[delim] == '*');
try{
const double factor = to_double(second);
return factor*parse_simple_ScaleFunction(first, known);
}
catch(std::invalid_argument const &){
const double factor = to_double(first);
return factor*parse_simple_ScaleFunction(first, known);
}
}
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::reweight},
{unob, EventTreatment::keep},
{unof, EventTreatment::keep},
{nonHEJ, EventTreatment::keep}
};
set_from_yaml(treat.at(FKL), yaml, "FKL");
set_from_yaml(treat.at(unob), yaml, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(nonHEJ), yaml, "non-HEJ");
if(treat[nonHEJ] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-HEJ 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");
config.max_ext_soft_pt_fraction = std::numeric_limits<double>::infinity();
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");
set_from_yaml(config.unweight, yaml, "unweight");
config.treat = get_event_treatment(yaml);
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.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
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:34 AM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887755
Default Alt Text
(42 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment