No OneTemporary
Actions

Size

21 KB

Subscribers

None

View Options

	diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
	index 39a24c4..2c4c899 100644
	--- a/src/YAMLreader.cc
	+++ b/src/YAMLreader.cc
	@@ -1,606 +1,608 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/YAMLreader.hh"

	#include <algorithm>
	#include <iostream>
	#include <limits>
	#include <map>
	#include <string>
	#include <unordered_map>
	#include <vector>

	#include <dlfcn.h>

	#include "HEJ/ConfigFlags.hh"
	#include "HEJ/Constants.hh"
	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/output_formats.hh"

	namespace HEJ {
	class Event;

	namespace {
	//! Get YAML tree of supported options
	/**
	* The configuration file is checked against this tree of options
	* in assert_all_options_known.
	*/
	YAML::Node const & get_supported_options(){
	const static YAML::Node supported = [](){
	YAML::Node supported;
	static const auto opts = {
	"trials", "min extparton pt", "max ext soft pt fraction",
	"soft pt regulator",
	"scales", "scale factors", "max scale ratio", "import scales",
	"log correction", "event output", "analysis", "analyses", "vev",
	"regulator parameter", "max events", "off-shell tolerance",
	"require low pt jet"
	};
	// add subnodes to "supported" - the assigned value is irrelevant
	for(auto && opt: opts) supported[opt] = "";
	for(auto && jet_opt: {"min pt", "algorithm", "R"}){
	supported["resummation jets"][jet_opt] = "";
	supported["fixed order jets"][jet_opt] = "";
	}
	for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
	supported["Higgs coupling"][opt] = "";
	}
	for(auto && opt: {"name", "seed"}){
	supported["random generator"][opt] = "";
	}
	for(auto && opt: {"enabled", "nlo order"}){
	supported["NLO truncation"][opt] = "";
	}
	for(
	auto && opt: {
	"FKL", "unordered", "extremal qqbar", "central qqbar", "non-resummable",
	"unknown", "invalid"
	}){
	supported["event treatment"][opt] = "";
	}
	for(auto && particle_type: {"Higgs", "W", "Z"}){
	for(auto && particle_opt: {"mass", "width"}){
	supported["particle properties"][particle_type][particle_opt] = "";
	}
	}
	for(auto && opt: {"type", "trials", "max deviation"}){
	supported["unweight"][opt] = "";
	}
	return supported;
	}();
	return supported;
	}

	fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
	using namespace fastjet;
	static const std::map<std::string, fastjet::JetAlgorithm> known = {
	{"kt", kt_algorithm},
	{"cambridge", cambridge_algorithm},
	{"antikt", antikt_algorithm},
	{"cambridge for passive", cambridge_for_passive_algorithm},
	{"plugin", plugin_algorithm}
	};
	const auto res = known.find(algo);
	if(res == known.end()){
	throw std::invalid_argument("Unknown jet algorithm \"" + algo + "\"");
	}
	return res->second;
	}

	EventTreatment to_EventTreatment(std::string const & name){
	static const std::map<std::string, EventTreatment> known = {
	{"reweight", EventTreatment::reweight},
	{"keep", EventTreatment::keep},
	{"discard", EventTreatment::discard},
	{"abort", EventTreatment::abort}
	};
	const auto res = known.find(name);
	if(res == known.end()){
	throw std::invalid_argument("Unknown event treatment \"" + name + "\"");
	}
	return res->second;
	}

	WeightType to_weight_type(std::string const & setting){
	if(setting == "weighted")
	return WeightType::weighted;
	if(setting =="resummation")
	return WeightType::unweighted_resum;
	if(setting =="partial")
	return WeightType::partially_unweighted;
	throw std::invalid_argument{"Unknown weight type \"" + setting + "\""};
	}

	} // namespace

	namespace detail{
	void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
	setting = to_JetAlgorithm(yaml.as<std::string>());
	}

	void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
	setting = to_EventTreatment(yaml.as<std::string>());
	}

	void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
	setting = to_ParticleID(yaml.as<std::string>());
	}

	void set_from_yaml(WeightType & setting, YAML::Node const & yaml){
	setting = to_weight_type(yaml.as<std::string>());
	}

	} // namespace detail

	JetParameters get_jet_parameters(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	JetParameters result;
	fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
	double R = NAN;
	set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
	set_from_yaml(R, node, entry, "R");
	result.def = fastjet::JetDefinition{jet_algo, R};
	set_from_yaml(result.min_pt, node, entry, "min pt");
	return result;
	}

	RNGConfig to_RNGConfig(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	RNGConfig result;
	set_from_yaml(result.name, node, entry, "name");
	set_from_yaml_if_defined(result.seed, node, entry, "seed");
	return result;
	}

	NLOConfig to_NLOConfig(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	NLOConfig result;
	set_from_yaml_if_defined(result.enabled, node, entry, "enabled");
	set_from_yaml_if_defined(result.nj, node, entry, "nlo order");
	return result;
	}

	ParticleProperties get_particle_properties(
	YAML::Node const & node, std::string const & entry,
	std::string const & boson
	){
	ParticleProperties result{};
	set_from_yaml(result.mass, node, entry, boson, "mass");
	set_from_yaml(result.width, node, entry, boson, "width");
	return result;
	}

	EWConstants get_ew_parameters(YAML::Node const & node){
	EWConstants result;
	double vev = NAN;
	set_from_yaml(vev, node, "vev");
	result.set_vevWZH(vev,
	get_particle_properties(node, "particle properties", "W"),
	get_particle_properties(node, "particle properties", "Z"),
	get_particle_properties(node, "particle properties", "Higgs")
	);
	return result;
	}

	HiggsCouplingSettings get_Higgs_coupling(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	static constexpr double mt_max = 2e4;
	#ifndef HEJ_BUILD_WITH_QCDLOOP
	if(node[entry].IsDefined()){
	throw std::invalid_argument{
	"Higgs coupling settings require building HEJ 2 "
	"with QCDloop support"
	};
	}
	#endif
	HiggsCouplingSettings settings;
	set_from_yaml_if_defined(settings.mt, node, entry, "mt");
	set_from_yaml_if_defined(settings.mb, node, entry, "mb");
	set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
	set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
	if(settings.use_impact_factors){
	if(settings.mt != std::numeric_limits<double>::infinity()){
	throw std::invalid_argument{
	"Conflicting settings: "
	"impact factors may only be used in the infinite top mass limit"
	};
	}
	}
	else{
	// huge values of the top mass are numerically unstable
	settings.mt = std::min(settings.mt, mt_max);
	}
	return settings;
	}

	FileFormat to_FileFormat(std::string const & name){
	static const std::map<std::string, FileFormat> known = {
	{"Les Houches", FileFormat::Les_Houches},
	{"HepMC", FileFormat::HepMC},
	{"HepMC2", FileFormat::HepMC2},
	{"HepMC3", FileFormat::HepMC3},
	{"HDF5", FileFormat::HDF5}
	};
	const auto res = known.find(name);
	if(res == known.end()){
	throw std::invalid_argument("Unknown file format \"" + name + "\"");
	}
	return res->second;
	}

	std::string extract_suffix(std::string const & filename){
	size_t separator = filename.rfind('.');
	if(separator == std::string::npos) return {};
	return filename.substr(separator + 1);
	}

	FileFormat format_from_suffix(std::string const & filename){
	const std::string suffix = extract_suffix(filename);
	if(suffix == "lhe") return FileFormat::Les_Houches;
	if(suffix == "hepmc") return FileFormat::HepMC;
	if(suffix == "hepmc3") return FileFormat::HepMC3;
	if(suffix == "hepmc2") return FileFormat::HepMC2;
	if(suffix == "hdf5") return FileFormat::HDF5;
	throw std::invalid_argument{
	"Can't determine format for output file \"" + filename + "\""
	};
	}

	void assert_all_options_known(
	YAML::Node const & conf, YAML::Node const & supported
	){
	if(!conf.IsMap()) return;
	if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
	for(auto const & entry: conf){
	const auto name = entry.first.as<std::string>();
	if(! supported[name]) throw unknown_option{name};
	/* check sub-options, e.g. 'resummation jets: min pt'
	* we don't check analyses sub-options
	* those depend on the analysis being used and should be checked there
	* similar for "import scales"
	*/
	if(name != "analyses" && name != "analysis" && name != "import scales"){
	try{
	assert_all_options_known(conf[name], supported[name]);
	}
	catch(unknown_option const & ex){
	throw unknown_option{name + ": " + ex.what()};
	}
	catch(invalid_type const & ex){
	throw invalid_type{name + ": " + ex.what()};
	}
	}
	}
	}

	} // namespace HEJ

	namespace YAML {

	Node convert<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
	Node node;
	node[to_string(outfile.format)] = outfile.name;
	return node;
	}

	bool convert<HEJ::OutputFile>::decode(Node const & node, HEJ::OutputFile & out) {
	switch(node.Type()){
	case NodeType::Map: {
	YAML::const_iterator it = node.begin();
	out.format = HEJ::to_FileFormat(it->first.as<std::string>());
	out.name = it->second.as<std::string>();
	return true;
	}
	case NodeType::Scalar:
	out.name = node.as<std::string>();
	out.format = HEJ::format_from_suffix(out.name);
	return true;
	default:
	return false;
	}
	}
	} // namespace YAML

	namespace HEJ {

	namespace detail{
	void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
	setting = yaml.as<OutputFile>();
	}
	}

	namespace {
	void update_fixed_order_jet_parameters(
	JetParameters & fixed_order_jets, YAML::Node const & yaml
	){
	if(!yaml["fixed order jets"]) return;
	set_from_yaml_if_defined(
	fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
	);
	fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
	set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
	double R = fixed_order_jets.def.R();
	set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
	fixed_order_jets.def = fastjet::JetDefinition{algo, R};
	}

	// like std::stod, but throw if not the whole string can be converted
	double to_double(std::string const & str){
	std::size_t pos = 0;
	const double result = std::stod(str, &pos);
	if(pos < str.size()){
	throw std::invalid_argument(str + " is not a valid double value");
	}
	return result;
	}

	using EventScale = double (*)(Event const &);

	void import_scale_functions(
	std::string const & file,
	std::vector<std::string> const & scale_names,
	std::unordered_map<std::string, EventScale> & known
	) {
	void * handle = dlopen(file.c_str(), RTLD_NOW);
	char * error = dlerror();
	if(error != nullptr) throw std::runtime_error{error};

	for(auto const & scale: scale_names) {
	void * sym = dlsym(handle, scale.c_str());
	error = dlerror();
	if(error != nullptr) throw std::runtime_error{error};
	known.emplace(scale, reinterpret_cast<EventScale>(sym)); // NOLINT
	}
	}

	auto get_scale_map(
	YAML::Node const & yaml
	) {
	std::unordered_map<std::string, EventScale> scale_map;
	scale_map.emplace("H_T", H_T);
	scale_map.emplace("max jet pperp", max_jet_pt);
	scale_map.emplace("jet invariant mass", jet_invariant_mass);
	scale_map.emplace("m_j1j2", m_j1j2);
	if(yaml["import scales"].IsDefined()) {
	if(! yaml["import scales"].IsMap()) {
	throw invalid_type{"Entry 'import scales' is not a map"};
	}
	for(auto const & import: yaml["import scales"]) {
	const auto file = import.first.as<std::string>();
	const auto scale_names =
	import.second.IsSequence()
	?import.second.as<std::vector<std::string>>()
	:std::vector<std::string>{import.second.as<std::string>()};
	import_scale_functions(file, scale_names, scale_map);
	}
	}
	return scale_map;
	}

	// simple (as in non-composite) scale functions
	/**
	* An example for a simple scale function would be H_T,
	* H_T/2 is then composite (take H_T and then divide by 2)
	*/
	ScaleFunction parse_simple_ScaleFunction(
	std::string const & scale_fun,
	std::unordered_map<std::string, EventScale> const & known
	) {
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	const auto it = known.find(scale_fun);
	if(it != end(known)) return {it->first, it->second};
	try{
	const double scale = to_double(scale_fun);
	return {scale_fun, FixedScale{scale}};
	} catch(std::invalid_argument const &){}
	throw std::invalid_argument{"Unknown scale choice: \"" + scale_fun + "\""};
	}

	std::string trim_front(std::string const & str){
	const auto new_begin = std::find_if(
	begin(str), end(str), [](char c){ return std::isspace(c) == 0; }
	);
	return std::string(new_begin, end(str));
	}

	std::string trim_back(std::string str){
	size_t pos = str.size() - 1;
	// use guaranteed wrap-around behaviour to check whether we have
	// traversed the whole string
	for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
	str.resize(pos + 1); // note that pos + 1 can be 0
	return str;
	}

	ScaleFunction parse_ScaleFunction(
	std::string const & scale_fun,
	std::unordered_map<std::string, EventScale> const & known
	){
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	// parse from right to left => a/b/c gives (a/b)/c
	const size_t delim = scale_fun.find_last_of("*/");
	if(delim == std::string::npos){
	return parse_simple_ScaleFunction(scale_fun, known);
	}
	const std::string first = trim_back(std::string{scale_fun, 0, delim});
	const std::string second = trim_front(std::string{scale_fun, delim+1});
	if(scale_fun[delim] == '/'){
	return parse_ScaleFunction(first, known)
	/ parse_ScaleFunction(second, known);
	}
	assert(scale_fun[delim] == '*');
	return parse_ScaleFunction(first, known)
	* parse_ScaleFunction(second, known);

	}

	EventTreatMap get_event_treatment(
	YAML::Node const & node, std::string const & entry
	){
	using namespace event_type;
	EventTreatMap treat {
	{FKL, EventTreatment::discard},
	{unob, EventTreatment::discard},
	{unof, EventTreatment::discard},
	{qqbar_exb, EventTreatment::discard},
	{qqbar_exf, EventTreatment::discard},
	{qqbar_mid, EventTreatment::discard},
	{non_resummable, EventTreatment::discard},
	{unknown, EventTreatment::abort},
	{invalid, EventTreatment::abort}
	};
	set_from_yaml(treat.at(FKL), node, entry, "FKL");
	set_from_yaml(treat.at(unob), node, entry, "unordered");
	treat.at(unof) = treat.at(unob);
	set_from_yaml(treat.at(qqbar_exb), node, entry, "extremal qqbar");
	treat.at(qqbar_exf) = treat.at(qqbar_exb);
	set_from_yaml(treat.at(qqbar_mid), node, entry, "central qqbar");
	set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable");
	- if(treat[non_resummable] == EventTreatment::reweight){
	- throw std::invalid_argument{"Cannot reweight non-resummable events"};
	- }
	set_from_yaml_if_defined(treat.at(unknown), node, entry, "unknown");
	set_from_yaml_if_defined(treat.at(invalid), node, entry, "invalid");
	+ for(auto type: {non_resummable, unknown, invalid}) {
	+ if(treat[type] == EventTreatment::reweight){
	+ throw std::invalid_argument{"Cannot reweight " + name(type) + " events"};
	+ }
	+ }
	return treat;
	}

	Config to_Config(YAML::Node const & yaml){
	try{
	assert_all_options_known(yaml, get_supported_options());
	}
	catch(unknown_option const & ex){
	throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
	}

	Config config;
	config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
	config.fixed_order_jets = config.resummation_jets;
	update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
	set_from_yaml_if_defined(config.min_extparton_pt, yaml, "min extparton pt");
	if(config.min_extparton_pt!=0)
	std::cerr << "WARNING: \"min extparton pt\" is deprecated."
	<< " Please remove this entry or set \"soft pt regulator\" instead.\n";
	set_from_yaml_if_defined(
	config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
	);
	if(config.max_ext_soft_pt_fraction){
	std::cerr << "WARNING: \"max ext soft pt fraction\" is deprecated."
	<< " Please remove this entry or set \"soft pt regulator\" instead.\n";
	config.soft_pt_regulator = *config.max_ext_soft_pt_fraction;
	} else {
	set_from_yaml_if_defined(
	config.soft_pt_regulator, yaml, "soft pt regulator"
	);
	}
	// Sets the standard value, then changes this if defined
	config.regulator_lambda=CLAMBDA;
	set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
	set_from_yaml_if_defined(config.max_events, yaml, "max events");
	set_from_yaml(config.trials, yaml, "trials");
	config.weight_type = WeightType::weighted;
	set_from_yaml_if_defined(config.weight_type, yaml, "unweight", "type");
	if(config.weight_type == WeightType::partially_unweighted) {
	config.unweight_config = PartialUnweightConfig{};
	set_from_yaml(
	config.unweight_config->trials, yaml,
	"unweight", "trials"
	);
	set_from_yaml(
	config.unweight_config->max_dev, yaml,
	"unweight", "max deviation"
	);
	}
	else if(yaml["unweight"].IsDefined()) {
	for(auto && opt: {"trials", "max deviation"}) {
	if(yaml["unweight"][opt].IsDefined()) {
	throw std::invalid_argument{
	"'unweight: " + std::string{opt} + "' "
	"is only supported if 'unweight: type' is set to 'partial'"
	};
	}
	}
	}
	set_from_yaml(config.log_correction, yaml, "log correction");
	config.treat = get_event_treatment(yaml, "event treatment");
	set_from_yaml_if_defined(config.output, yaml, "event output");
	config.rng = to_RNGConfig(yaml, "random generator");
	set_from_yaml_if_defined(config.lowpt, yaml, "require low pt jet");

	set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses");
	if(yaml["analysis"].IsDefined()){
	std::cerr <<
	"WARNING: Configuration entry 'analysis' is deprecated. "
	" Use 'analyses' instead.\n";
	set_from_yaml(config.analysis_parameters, yaml, "analysis");
	if(!config.analysis_parameters.IsNull()){
	config.analyses_parameters.push_back(config.analysis_parameters);
	}
	}
	config.scales = to_ScaleConfig(yaml);
	config.ew_parameters = get_ew_parameters(yaml);
	config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");

	//HEJ@NLO Truncation
	config.nlo = to_NLOConfig(yaml, "NLO truncation");
	set_from_yaml_if_defined(
	config.off_shell_tolerance,
	yaml,
	"off-shell tolerance"
	);
	return config;
	}

	} // namespace

	ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
	ScaleConfig config;
	auto scale_funs = get_scale_map(yaml);
	std::vector<std::string> scales;
	set_from_yaml(scales, yaml, "scales");
	config.base.reserve(scales.size());
	std::transform(
	begin(scales), end(scales), std::back_inserter(config.base),
	[scale_funs](auto const & entry){
	return parse_ScaleFunction(entry, scale_funs);
	}
	);
	set_from_yaml_if_defined(config.factors, yaml, "scale factors");
	config.max_ratio = std::numeric_limits<double>::infinity();
	set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
	return config;
	}

	Config load_config(std::string const & config_file){
	try{
	return to_Config(YAML::LoadFile(config_file));
	}
	catch(...){
	std::cerr << "Error reading " << config_file << ":\n ";
	throw;
	}
	}

	} // namespace HEJ

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:55 PM (1 d, 5 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805938
Default Alt Text: (21 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions