Page MenuHomeHEPForge
Files F8725662

No OneTemporary
Actions

View Options

diff --git a/src/config.cc b/src/config.cc
index 694d75c..0128c0f 100644
--- a/src/config.cc
+++ b/src/config.cc
@@ -1,599 +1,622 @@
#include "RHEJ/config.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
namespace RHEJ{
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 == "hepmc3") return FileFormat::HepMC;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
}
namespace YAML {
+
+ std::string to_string(NodeType::value type){
+ switch (type) {
+ case NodeType::Null: return "Null";
+ case NodeType::Scalar: return "Scalar";
+ case NodeType::Sequence: return "Sequence";
+ case NodeType::Map: return "Map";
+ case NodeType::Undefined: return "Undefined";
+ default:
+ throw std::logic_error{"Missing case in switch statement"};
+ }
+ }
+
template<>
struct convert<RHEJ::OutputFile> {
static Node encode(RHEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
}
static bool decode(Node const & node, RHEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = RHEJ::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 = RHEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
};
}
namespace RHEJ{
namespace{
struct invalid_type: std::invalid_argument {
explicit invalid_type(std::string const & what):
std::invalid_argument{what} {};
explicit invalid_type(char const * what):
std::invalid_argument{what} {};
};
struct missing_option: std::logic_error {
explicit missing_option(std::string const & what):
std::logic_error{what} {};
explicit missing_option(char const * what):
std::logic_error{what} {};
};
static const std::set<std::string> known = {
{"trials"},
{"min extparton pt"}, {"max ext soft pt fraction"},
{"resummation jets"}, {"fixed order jets"},
{"FKL"}, {"non-FKL"}, {"unordered"},
{"log correction"},
{"event output"},
{"analysis"},
{"unweight"}, {"RanLux init"},
{"scales"}, {"scale factors"}, {"max scale ratio"},
{"Higgs coupling"}
};
static const std::set<std::string> known_jet = {
{"min pt"}, {"algorithm"}, {"R"}
};
template<typename T>
std::string type_string();
template<>
std::string type_string<std::vector<double>>(){
return "array of doubles";
}
template<>
std::string type_string<std::vector<std::string>>(){
return "array of strings";
}
template<>
std::string type_string<int>(){
return "integer";
}
template<>
std::string type_string<double>(){
return "double";
}
template<>
std::string type_string<bool>(){
return "bool";
}
template<>
std::string type_string<std::string>(){
return "string";
}
template<>
std::string type_string<std::vector<OutputFile>>(){
return "array of output files";
}
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;
}
// helper struct to allow partial template specialisation
template<typename T>
struct existing_as_helper{
static T as(YAML::Node const & config, std::string const & entry){
assert(config[entry]);
try{
return config[entry].as<T>();
}
catch(std::exception const &){
throw invalid_type{
"Entry " + entry + " is not of type " + type_string<T>()
};
}
}
};
template<>
struct existing_as_helper<fastjet::JetAlgorithm>{
static fastjet::JetAlgorithm as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
const std::string algo_name =
existing_as_helper<std::string>::as(config, entry);
return to_JetAlgorithm(algo_name);
}
};
template<>
struct existing_as_helper<EventTreatment>{
static EventTreatment as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
const std::string name =
existing_as_helper<std::string>::as(config, entry);
return to_EventTreatment(name);
}
};
template<>
struct existing_as_helper<OutputFile>{
static OutputFile as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
YAML::convert<RHEJ::OutputFile> converter{};
OutputFile out;
if(converter.decode(config[entry], out)) return out;
throw std::invalid_argument{
"Bad output file setting: " + config[entry].as<std::string>()
};
}
};
template<typename T>
struct existing_as_helper<std::vector<T>>{
static std::vector<T> as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
// special case: treat a single value like a vector with one element
if(config[entry].IsScalar()){
return {existing_as_helper<T>::as(config, entry)};
}
try{
return config[entry].as<std::vector<T>>();
}
catch(std::exception const & p){
std::cerr << "Error: " << p.what() << '\n';
throw invalid_type{
"Entry " + entry + " is not of type " + type_string<std::vector<T>>()
};
}
}
};
template<typename T>
T existing_as(YAML::Node const & config, std::string const & entry){
return existing_as_helper<T>::as(config, entry);
}
template<typename T>
T as(YAML::Node const & config, std::string const & entry);
template<>
JetParameters existing_as<JetParameters>(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
YAML::Node const & jet_config = config[entry];
JetParameters result;
for(auto const & opt: jet_config){
auto const & opt_name = opt.first.as<std::string>();
if(! known_jet.count(opt_name)){
std::cerr << "In entry " + entry + ": "
+ "Unknown option " + opt_name + "\n";
}
}
try{
result.def = fastjet::JetDefinition{
as<fastjet::JetAlgorithm>(jet_config, "algorithm"),
as<double>(jet_config, "R")
};
result.min_pt = as<double>(jet_config, "min pt");
return result;
}
catch(missing_option const & exc){
throw missing_option{"In entry " + entry + ":\n" + exc.what()};
}
catch(invalid_type const & exc){
throw invalid_type{"In entry " + entry + ":\n" + exc.what()};
}
}
template<typename T>
T as(YAML::Node const & config, std::string const & entry){
if(!config[entry]){
throw missing_option{"No entry for option " + entry};
}
return existing_as<T>(config, entry);
}
template<typename T>
optional<T> as_optional(YAML::Node const & config, std::string const & entry){
if(!config[entry]) return {};
return {existing_as<T>(config, entry)};
}
template<>
HiggsCouplingSettings as<HiggsCouplingSettings>(
YAML::Node const & config, std::string const & entry
){
static constexpr double mt_max = 2e4;
HiggsCouplingSettings settings;
if(!config[entry]) return settings;
#ifndef RHEJ_BUILD_WITH_LT
throw std::invalid_argument{
"Higgs coupling settings require building Reversed HEJ "
"with looptools support"
};
#endif
- if(config[entry]["mt"]){
- settings.mt = config[entry]["mt"].as<double>();
- }
- if(config[entry]["mb"]){
- settings.mb = config[entry]["mb"].as<double>();
- }
- if(config[entry]["include bottom"]){
- settings.include_bottom = config[entry]["include bottom"].as<bool>();
+ YAML::Node const & node = config[entry];
+ if(!node.IsMap()){
+ throw std::invalid_argument{
+ "entry '" + entry + "' is a " + to_string(node.Type())
+ + ", not a Map of properties"
+ };
}
- if(config[entry]["use impact factors"]){
- settings.use_impact_factors =
- config[entry]["use impact factors"].as<bool>();
+ for(auto const & property: node){
+ const auto key = property.first.as<std::string>();
+ auto const & value = property.second;
+ if(key == "mt") settings.mt = value.as<double>();
+ else if(key == "mb") settings.mb = value.as<double>();
+ else if(key == "include bottom"){
+ settings.include_bottom = value.as<bool>();
+ }
+ else if(key == "use impact factors"){
+ settings.use_impact_factors = value.as<bool>();
+ }
+ else{
+ std::cerr << "WARNING: unknown option '" << key
+ << "' in '" << entry << "'\n";
+ }
}
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;
}
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
YAML::Node const & yaml_jets = yaml["fixed order jets"];
assert(yaml_jets);
for(auto const & opt: yaml_jets){
auto const & opt_name = opt.first.as<std::string>();
if(! known_jet.count(opt_name)){
std::cerr << "In entry fixed order jets: "
"Unknown option " + opt_name + "\n";
}
}
try{
auto algo = as_optional<fastjet::JetAlgorithm>(yaml_jets, "algorithm");
if(algo){
fixed_order_jets.def = fastjet::JetDefinition{
*algo, fixed_order_jets.def.R()
};
}
auto R = as_optional<double>(yaml_jets, "R");
if(R){
fixed_order_jets.def = fastjet::JetDefinition{
fixed_order_jets.def.jet_algorithm(), *R
};
}
auto min_pt = as_optional<double>(yaml_jets, "min pt");
if(min_pt) fixed_order_jets.min_pt = *min_pt;
}
catch(missing_option const & exc){
throw missing_option{
std::string{"In entry fixed order jets:\n"} + exc.what()
};
}
catch(invalid_type const & exc){
throw invalid_type{
std::string{"In entry fixed order jets:\n"} + exc.what()
};
}
}
// 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;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be Ht,
* Ht/2 is then composite (take Ht and then divide by 2)
*/
std::unique_ptr<ScaleFunction> make_simple_ScaleFunction_ptr(
std::string const & scale_fun
){
using ret_type = std::unique_ptr<ScaleFunction>;
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
if(scale_fun == "input") return ret_type{new InputScales{}};
if(scale_fun == "Ht") return ret_type{new Ht{}};
if(scale_fun == "max jet pperp") return ret_type{new MaxJetPperp{}};
if(scale_fun == "jet invariant mass"){
return ret_type{new JetInvariantMass{}};
}
try{
const double scale = to_double(scale_fun);
return ret_type{new 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;
}
std::unique_ptr<ScaleFunction> make_ScaleFunction_ptr(
std::string const & scale_fun
){
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 make_simple_ScaleFunction_ptr(scale_fun);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
double factor;
std::unique_ptr<ScaleFunction> fun;
if(scale_fun[delim] == '/'){
factor = 1/to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
else{
assert(scale_fun[delim] == '*');
try{
factor = to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
catch(std::invalid_argument const &){
factor = to_double(first);
fun = make_simple_ScaleFunction_ptr(second);
}
}
assert(fun != nullptr);
return std::unique_ptr<ScaleFunction>{
new Product{factor, std::move(fun)}
};
}
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
const auto scales = as<std::vector<std::string>>(yaml, "scales");
config.scales.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.scales),
make_ScaleFunction_ptr
);
auto scale_factors = as_optional<std::vector<double>>(
yaml, "scale factors"
);
if(scale_factors) config.scale_factors = std::move(*scale_factors);
const auto max_scale_ratio = as_optional<double>(yaml, "max scale ratio");
static_assert(
std::numeric_limits<double>::has_infinity, "infinity not supported"
);
config.max_scale_ratio = max_scale_ratio?
*max_scale_ratio:
std::numeric_limits<double>::infinity()
;
return config;
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_class;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard}
};
treat.emplace(FKL, as<EventTreatment>(yaml, "FKL"));
const auto unordered_treatment = as<EventTreatment>(yaml, "unordered");
treat.emplace(unordered_forward, unordered_treatment);
treat.emplace(unordered_backward, unordered_treatment);
treat.emplace(nonFKL, as<EventTreatment>(yaml, "non-FKL"));
if(treat[nonFKL] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-FKL events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
for(auto const & opt: yaml){
auto const & opt_name = opt.first.as<std::string>();
if(! known.count(opt_name)){
std::cerr << "WARNING: unknown option " + opt_name + "\n";
}
}
Config config;
config.resummation_jets = as<JetParameters>(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
if(yaml["fixed order jets"]){
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
}
config.min_extparton_pt = as<double>(yaml, "min extparton pt");
config.max_ext_soft_pt_fraction = yaml["max ext soft pt fraction"]?
as<double>(yaml, "max ext soft pt fraction"):
std::numeric_limits<double>::infinity()
;
config.trials = as<int>(yaml, "trials");
config.log_correction = as<bool>(yaml, "log correction");
config.unweight = as<bool>(yaml, "unweight");
config.treat = get_event_treatment(yaml);
if(yaml["event output"]){
config.output = as<std::vector<OutputFile>>(yaml, "event output");
}
config.RanLux_init = as_optional<std::string>(yaml, "RanLux init");
config.analysis_parameters = yaml["analysis"]?yaml["analysis"]:YAML::Node{};
config.scale_gen = ScaleGenerator{to_ScaleConfig(yaml)};
config.Higgs_coupling = as<HiggsCouplingSettings>(yaml, "Higgs coupling");
return config;
}
} // anonymous namespace
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;
}
}
Event ScaleGenerator::operator()(Event ev) const{
assert(ev.variations().empty());
assert(! settings_.scales.empty());
auto const & scales = settings_.scales;
auto const & scale_factors = settings_.scale_factors;
const double max_scale_ratio = settings_.max_scale_ratio;
// check if we are doing scale variation at all
if(scales.size() == 1 && scale_factors.empty()){
ev.central() = (*scales.front())(ev);
return ev;
}
for(auto && base_scale: scales){
const EventParameters cur_base = (*base_scale)(ev);
ev.variations().emplace_back(cur_base);
//multiplicative scale variation
for(double mur_factor: scale_factors){
const double mur = cur_base.mur*mur_factor;
for(double muf_factor: scale_factors){
if(muf_factor == 1. && mur_factor == 1.) continue;
const double muf = cur_base.muf*muf_factor;
if(mur/muf < 1/max_scale_ratio || mur/muf > max_scale_ratio) continue;
ev.variations().emplace_back(
EventParameters{mur, muf, cur_base.weight}
);
}
}
};
ev.central() = (*scales.front())(ev);
return ev;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:16 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243624
Default Alt Text
(20 KB)

Event Timeline