No OneTemporary
Actions

Size

27 KB

Subscribers

None

View Options

	diff --git a/CMakeLists.txt b/CMakeLists.txt
	index 6f5eedf..6456970 100644
	--- a/CMakeLists.txt
	+++ b/CMakeLists.txt
	@@ -1,170 +1,177 @@
	cmake_minimum_required(VERSION 3.1 FATAL_ERROR)

	set(CMAKE_LEGACY_CYGWIN_WIN32 0)
	set(CMAKE_EXPORT_COMPILE_COMMANDS ON)

	project("Reversed HEJ" VERSION 0.0.1 LANGUAGES C CXX)

	## Flags for the compiler. No warning allowed.
	if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Werror")
	elseif (MSVC)
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc")
	endif()
	set(CMAKE_CXX_STANDARD_REQUIRED ON)
	set(CMAKE_CXX_STANDARD 14)

	## Create Version
	set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")


	# Get the latest abbreviated commit hash of the working branch
	execute_process(
	COMMAND git rev-parse HEAD
	WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
	OUTPUT_VARIABLE PROJECT_GIT_REVISION
	OUTPUT_STRIP_TRAILING_WHITESPACE
	)
	# Get the current working branch
	execute_process(
	COMMAND git rev-parse --abbrev-ref HEAD
	WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
	OUTPUT_VARIABLE PROJECT_GIT_BRANCH
	OUTPUT_STRIP_TRAILING_WHITESPACE
	)

	CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
	${PROJECT_BINARY_DIR}/include/RHEJ/Version.hh @ONLY )


	## Add directories and find dependences
	include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)

	find_package(fastjet REQUIRED)
	include_directories(${fastjet_INCLUDE_PATH})
	find_package(clhep REQUIRED)
	include_directories(${clhep_INCLUDE_PATH})
	find_package(lhapdf REQUIRED)
	include_directories(${lhapdf_INCLUDE_PATH})
	find_package(gsl REQUIRED)
	include_directories(${gsl_INCLUDE_PATH})
	find_package(Boost REQUIRED COMPONENTS iostreams)
	include_directories(${Boost_INCLUDE_DIRS})
	# find_package(yaml-cpp REQUIRED)
	# include_directories(${YAML_CPP_INCLUDE_DIR})
	find_package(HepMC 2)
	if(${HepMC_FOUND})
	message (STATUS "HepMC installation found: ${HepMC_INCLUDE_DIRS}")
	set(
	CMAKE_CXX_FLAGS
	"${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_HepMC_VERSION=${HepMC_VERSION_MAJOR}"
	)
	include_directories(${HepMC_INCLUDE_DIRS})
	endif()
	find_package(QCDloop 2)
	if(${QCDloop_FOUND})
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DRHEJ_BUILD_WITH_QCDLOOP")
	include_directories(SYSTEM ${QCDloop_INCLUDE_DIRS})
	endif()

	add_subdirectory(src)

	## define executable
	add_executable(rHEJ src/main.cc)

	## link libraries
	target_link_libraries(rHEJ rhej)

	#plugins
	add_subdirectory("analysis-plugins")

	file(GLOB rhej_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/RHEJ/.hh ${PROJECT_BINARY_DIR}/include/RHEJ/.hh)
	file(GLOB lhef_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/LHEF/*.h)
	install(FILES ${rhej_headers} DESTINATION include/RHEJ/)
	install(FILES ${lhef_headers} DESTINATION include/LHEF/)
	install(TARGETS rHEJ DESTINATION bin)

	enable_testing()
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
	add_executable(test_Matrix ${tst_dir}/test_Matrix.cc)
	add_executable(test_classify ${tst_dir}/test_classify.cc)
	add_executable(test_psp ${tst_dir}/test_psp.cc)
	add_executable(test_ME_h_3j ${tst_dir}/test_ME_h_3j.cc)
	add_executable(test_ME_hjets_mt174 ${tst_dir}/test_ME_hjets_mt174.cc)
	add_executable(check_res ${tst_dir}/check_res.cc)
	add_executable(check_lhe ${tst_dir}/check_lhe.cc)
	+add_library(scales SHARED ${tst_dir}/scales.cc)
	+add_executable(test_scale_import ${tst_dir}/test_scale_import)
	target_link_libraries(test_Matrix rhej)
	target_link_libraries(test_classify rhej)
	target_link_libraries(test_psp rhej)
	target_link_libraries(test_ME_h_3j rhej)
	target_link_libraries(test_ME_hjets_mt174 rhej)
	target_link_libraries(check_res rhej)
	target_link_libraries(check_lhe rhej)
	+target_link_libraries(test_scale_import rhej)

	add_test(
	NAME t_matrix
	COMMAND test_Matrix
	)
	add_test(
	NAME t_classify
	COMMAND test_classify ${tst_dir}/classify.lhe.gz
	)
	add_test(
	NAME t_psp
	COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz
	)
	add_test(
	NAME t_ME
	COMMAND test_ME_h_3j
	)
	if(${QCDloop_FOUND})
	add_test(
	NAME t_ME_mt174
	COMMAND test_ME_hjets_mt174
	)
	endif()
	add_test(
	NAME t_2j
	COMMAND check_res ${tst_dir}/2j.lhe.gz 3.382e7 752159
	)
	add_test(
	NAME t_3j
	COMMAND check_res ${tst_dir}/3j.lhe.gz 2.5019e+06 97075
	)
	add_test(
	NAME t_4j
	COMMAND check_res ${tst_dir}/4j.lhe.gz 479170 19203.9
	)
	add_test(
	NAME t_h_3j
	COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.793107 0.0479054
	)
	add_test(
	NAME t_h_3j_uno
	COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0260428 0.00129855 uno
	)
	if(${HepMC_FOUND})
	file(READ "${tst_dir}/jet_config.yml" config)
	file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc")
	add_test(
	NAME t_main
	COMMAND rHEJ ${tst_dir}/jet_config_withHepMC.yml ${tst_dir}/2j.lhe.gz
	)
	if(${HepMC_VERSION_MAJOR} GREATER 2)
	add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
	target_link_libraries(check_hepmc rhej)
	add_test(
	NAME t_hepmc
	COMMAND check_hepmc tst.hepmc
	)
	endif()
	else()
	add_test(
	NAME t_main
	COMMAND rHEJ ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
	)
	endif()
	add_test(
	NAME t_lhe
	COMMAND check_lhe tst.lhe
	)
	+add_test(
	+ NAME t_scale_import
	+ COMMAND test_scale_import ${tst_dir}/jet_config_with_import.yml
	+ )
	diff --git a/config.yml b/config.yml
	index 1eeaa0f..95e7482 100644
	--- a/config.yml
	+++ b/config.yml
	@@ -1,80 +1,85 @@
	# number of attempted resummation phase space points for each input event
	trials: 10

	min extparton pt: 30 # minimum transverse momentum of extremal partons

	resummation jets: # resummation jet properties
	min pt: 35 # minimum jet transverse momentum
	algorithm: antikt # jet algorithm
	R: 0.4 # jet R parameter

	fixed order jets: # properties of input jets
	min pt: 30
	# by default, algorithm and R are like for resummation jets

	# treatment of he various event classes
	# the supported settings are: reweight, keep, discard
	# non-FKL events cannot be reweighted
	FKL: reweight
	unordered: keep
	non-FKL: keep

	# scale settings similar to original HEJ
	#
	# Use combinations of max jet pperp, input scales, ht/2,
	# and the jet invariant mass and vary all scales by factors
	# of 1, sqrt(2), and 2. Discard combinations where mur and muf
	# differ by a factor of more than two.
	#
	# The weight entries in the final events are ordered as follows:
	# 0-18: max jet pperp
	# 19-37: input scales
	# 38-56: ht/2
	# 57-75: jet invariant mass
	# In each of these groups, the first entry corresponds to the basic
	# scale choice. In the following entries, mur and muf are varied with
	# the above factors. The entries are ordered lexicographically so that
	# mur1 < mur2 or (mur1 == mur2 and muf1 < muf2).
	#
	# Note that in contrast to HEJ, the central choice for the event is always
	# max jet pperp and cannot be configured (yet).
	#
	# scales: [max jet pperp, input, Ht/2, jet invariant mass]
	# scale factors: [0.5, 0.7071, 1, 1.41421, 2]
	# max scale ratio: 2.0001

	scales: 91.188

	+# import scale setting functions
	+#
	+# import scales:
	+# lib_my_scales.so: [scale0,scale1]
	+
	log correction: false # whether or not to include higher order logs
	unweight: false # TODO: whether or not to unweight events

	# event output files
	#
	# the supported formats are
	# - Les Houches (suffix .lhe)
	# - HepMC (suffix .hepmc3)
	# TODO: - ROOT ntuples (suffix .root)
	#
	# An output file's format is deduced either automatically from the suffix
	# or from an explicit specification, e.g.
	# - Les Houches: outfile

	event output:
	- RHEJ.lhe
	# - RHEJ_events.hepmc

	analysis:
	# to use a custom analysis
	# plugin: ./src/analysis-plugins/libVBF.so
	# output: RHEJ.root
	# wtwt cut: # optional cut on (event weight)^2

	#RanLux init: ranlux.0 # file for initialisation of random number engine

	# parameters for Higgs-gluon couplings
	# this requires compilation with looptools
	# Higgs coupling:
	# use impact factors: false
	# mt: 174
	# include bottom: true
	# mb: 4.7
	diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
	index 798092d..ebf9581 100644
	--- a/src/YAMLreader.cc
	+++ b/src/YAMLreader.cc
	@@ -1,427 +1,478 @@
	#include "RHEJ/YAMLreader.hh"

	#include <set>
	#include <string>
	#include <vector>
	#include <iostream>
	#include <stdexcept>

	+#include <dlfcn.h>
	+
	namespace RHEJ{

	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-FKL",
	- "scales", "scale factors", "max scale ratio",
	+ "scales", "scale factors", "max scale ratio", "import scales",
	"log correction", "unweight", "event output", "analysis",
	"RanLux init"
	};
	// 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] = "";
	}
	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 RHEJ::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;
	}

	HiggsCouplingSettings get_Higgs_coupling(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	static constexpr double mt_max = 2e4;
	#ifndef RHEJ_BUILD_WITH_QCDLOOP
	if(node[entry]){
	throw std::invalid_argument{
	"Higgs coupling settings require building Reversed HEJ "
	"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"){
	+ 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 RHEJ

	namespace YAML {

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

	bool convert<RHEJ::OutputFile>::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 YAML

	namespace RHEJ{

	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;
	}

	+ void import_scale_functions(
	+ std::string const & file,
	+ std::vector<std::string> const & scale_names,
	+ std::unordered_map<std::string, ScaleFunction> & known
	+ ) {
	+ using ScaleFunction = double (*)(Event const &);
	+
	+ 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<ScaleFunction>(sym));
	+ }
	+ }
	+
	+ auto get_scale_map(
	+ YAML::Node const & yaml
	+ ) {
	+ std::unordered_map<std::string, ScaleFunction> 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) {
	+ ScaleFunction parse_simple_ScaleFunction(
	+ std::string const & scale_fun,
	+ std::unordered_map<std::string, ScaleFunction> const & known
	+ ) {
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	- if(scale_fun == "H_T") return H_T;
	- if(scale_fun == "max jet pperp") return max_jet_pt;
	- if(scale_fun == "jet invariant mass") return jet_invariant_mass;
	- if(scale_fun == "m_j1j2") return m_j1j2;
	+ const auto it = known.find(scale_fun);
	+ if(it != end(known)) return it->second;
	try{
	const double scale = to_double(scale_fun);
	return 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::string const & scale_fun,
	+ std::unordered_map<std::string, ScaleFunction> 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);
	+ 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});
	double factor;
	ScaleFunction fun;
	if(scale_fun[delim] == '/'){
	factor = 1/to_double(second);
	- fun = parse_simple_ScaleFunction(first);
	+ fun = parse_simple_ScaleFunction(first, known);
	}
	else{
	assert(scale_fun[delim] == '*');
	try{
	factor = to_double(second);
	- fun = parse_simple_ScaleFunction(first);
	+ fun = parse_simple_ScaleFunction(first, known);
	}
	catch(std::invalid_argument const &){
	factor = to_double(first);
	- fun = parse_simple_ScaleFunction(second);
	+ fun = parse_simple_ScaleFunction(second, known);
	}
	}
	assert(fun != nullptr);
	return Product{factor, std::move(fun)};
	}

	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},
	{nonFKL, 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(nonFKL), 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){
	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");
	set_from_yaml_if_defined(config.RanLux_init, yaml, "RanLux init");
	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),
	- parse_ScaleFunction
	+ [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 RHEJ
	diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml
	new file mode 100644
	index 0000000..66a84c0
	--- /dev/null
	+++ b/t/jet_config_with_import.yml
	@@ -0,0 +1,26 @@
	+trials: 10
	+
	+min extparton pt: 30
	+
	+resummation jets:
	+ min pt: 35
	+ algorithm: antikt
	+ R: 0.4
	+
	+fixed order jets:
	+ min pt: 30
	+
	+FKL: reweight
	+unordered: discard
	+non-FKL: discard
	+
	+log correction: false
	+unweight: false
	+
	+scales: softest_jet_pt
	+
	+event output:
	+ - tst.lhe
	+
	+import scales:
	+ ./libscales.so: softest_jet_pt
	diff --git a/t/scales.cc b/t/scales.cc
	new file mode 100644
	index 0000000..4636443
	--- /dev/null
	+++ b/t/scales.cc
	@@ -0,0 +1,8 @@
	+#include "RHEJ/Event.hh"
	+
	+extern "C"
	+__attribute__((visibility("default")))
	+double softest_jet_pt(RHEJ::Event const & ev){
	+ const auto softest_jet = sorted_by_pt(ev.jets()).back();
	+ return softest_jet.perp();
	+}
	diff --git a/t/test_scale_import.cc b/t/test_scale_import.cc
	new file mode 100644
	index 0000000..b16a3b2
	--- /dev/null
	+++ b/t/test_scale_import.cc
	@@ -0,0 +1,31 @@
	+#include <stdexcept>
	+#include <iostream>
	+
	+#include "RHEJ/YAMLreader.hh"
	+#include "RHEJ/Event.hh"
	+
	+int main(int argc, char** argv) {
	+ constexpr double ep = 1e-7;
	+ if (argc != 2) {
	+ throw std::logic_error{"wrong number of args"};
	+ }
	+ const RHEJ::Config config = RHEJ::load_config(argv[1]);
	+
	+ RHEJ::UnclusteredEvent tmp;
	+ tmp.outgoing.push_back(
	+ {RHEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)}
	+ );
	+ tmp.outgoing.push_back(
	+ {RHEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)}
	+ );
	+ RHEJ::Event ev{
	+ std::move(tmp),
	+ fastjet::JetDefinition{fastjet::kt_algorithm, 0.4},
	+ 20.
	+ };
	+
	+ const double softest_pt = config.scales.base[0](ev);
	+ if(std::abs(softest_pt-30.) > ep){
	+ throw std::logic_error{"wrong softest pt"};
	+ }
	+}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:11 PM (21 h, 57 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242403
Default Alt Text: (27 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions