No OneTemporary
Actions

Size

39 KB

Subscribers

None

View Options

	diff --git a/CMakeLists.txt b/CMakeLists.txt
	index 9f87541..11fc0b5 100644
	--- a/CMakeLists.txt
	+++ b/CMakeLists.txt
	@@ -1,249 +1,255 @@
	cmake_minimum_required(VERSION 3.1 FATAL_ERROR)

	set(CMAKE_LEGACY_CYGWIN_WIN32 0)
	set(CMAKE_EXPORT_COMPILE_COMMANDS ON)

	project("HEJ" VERSION 2.0.3 LANGUAGES C CXX)

	# Set a default build type if none was specified
	set(default_build_type "Release")

	if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
	message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
	set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
	STRING "Choose the type of build." FORCE)
	# Set the possible values of build type for cmake-gui
	set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
	"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
	endif()

	## 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")
	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
	)

	## target directories for install
	set(INSTALL_INCLUDE_DIR_BASE include)
	set(INSTALL_INCLUDE_DIR ${INSTALL_INCLUDE_DIR_BASE}/HEJ)
	set(INSTALL_BIN_DIR bin)
	set(INSTALL_LIB_DIR lib)
	set(INSTALL_CONFIG_DIR lib/cmake/HEJ)

	## Template files
	configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
	${PROJECT_BINARY_DIR}/include/HEJ/Version.hh @ONLY )

	configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/HEJ-config.cc.in
	${PROJECT_BINARY_DIR}/src/bin/HEJ-config.cc @ONLY )

	# Generate CMake config file
	include(CMakePackageConfigHelpers)
	configure_package_config_file(
	cmake/Templates/hej-config.cmake.in
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
	INSTALL_DESTINATION ${INSTALL_CONFIG_DIR}
	PATH_VARS INSTALL_INCLUDE_DIR_BASE INSTALL_LIB_DIR
	)
	write_basic_package_version_file(
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
	COMPATIBILITY SameMajorVersion
	)
	install(
	FILES ${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
	DESTINATION ${INSTALL_CONFIG_DIR})

	## 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 2.3 REQUIRED)
	include_directories(${clhep_INCLUDE_PATH})
	find_package(lhapdf REQUIRED)
	include_directories(${lhapdf_INCLUDE_PATH})
	## Amend unintuitive behaviour of FindBoost.cmake
	## Priority of BOOST_ROOT over BOOSTROOT matches FindBoost.cmake
	## at least for cmake 3.12
	if(DEFINED BOOST_ROOT)
	message("BOOST_ROOT set - only looking for Boost in ${BOOST_ROOT}")
	set(Boost_NO_BOOST_CMAKE ON)
	elseif(DEFINED BOOSTROOT)
	message("BOOSTROOT set - only looking for Boost in ${BOOSTROOT}")
	set(Boost_NO_BOOST_CMAKE ON)
	endif()
	find_package(Boost REQUIRED COMPONENTS iostreams)
	include_directories(${Boost_INCLUDE_DIRS})
	find_package(yaml-cpp)
	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} -DHEJ_BUILD_WITH_HepMC_VERSION=${HepMC_VERSION_MAJOR}"
	)
	include_directories(${HepMC_INCLUDE_DIRS})
	find_package(rivet)
	if(${rivet_FOUND})
	include_directories(${rivet_INCLUDE_PATH})
	set(
	CMAKE_CXX_FLAGS
	"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_RIVET"
	)
	endif()
	endif()
	find_package(QCDloop 2)
	if(${QCDloop_FOUND})
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_QCDLOOP")
	include_directories(SYSTEM ${QCDloop_INCLUDE_DIRS})
	endif()

	add_subdirectory(src)

	## define executable
	add_executable(HEJ src/bin/HEJ.cc)

	## link libraries
	target_link_libraries(HEJ hejlib)

	add_executable(HEJ-config src/bin/HEJ-config.cc)

	file(GLOB hej_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/HEJ/.hh ${PROJECT_BINARY_DIR}/include/HEJ/.hh)
	file(GLOB lhef_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/LHEF/*.h)
	install(FILES ${hej_headers} DESTINATION ${INSTALL_INCLUDE_DIR})
	install(FILES ${lhef_headers} DESTINATION include/LHEF/)
	install(TARGETS HEJ HEJ-config DESTINATION ${INSTALL_BIN_DIR})

	## tests
	enable_testing()
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
	add_executable(test_classify ${tst_dir}/test_classify.cc)
	add_executable(test_psp ${tst_dir}/test_psp.cc)
	add_executable(test_ME_generic ${tst_dir}/test_ME_generic.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)
	add_executable(test_descriptions ${tst_dir}/test_descriptions)
	+add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics)
	target_link_libraries(test_classify hejlib)
	target_link_libraries(test_psp hejlib)
	target_link_libraries(test_ME_generic hejlib)
	target_link_libraries(check_res hejlib)
	target_link_libraries(check_lhe hejlib)
	target_link_libraries(test_scale_import hejlib)
	target_link_libraries(test_descriptions hejlib)
	+target_link_libraries(test_scale_arithmetics hejlib)

	## add tests
	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
	)
	set(tst_ME_data_dir "${tst_dir}/ME_data")
	add_test(
	NAME t_ME_j
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
	)
	add_test(
	NAME t_ME_h
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	if(${QCDloop_FOUND})
	add_test(
	NAME t_ME_h_mt
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	add_test(
	NAME t_ME_h_mtmb
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	endif()
	add_test(
	NAME t_2j
	COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684
	)
	add_test(
	NAME t_3j
	COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
	)
	add_test(
	NAME t_4j
	COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
	)
	add_test(
	NAME t_h_3j
	COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334
	)
	add_test(
	NAME t_h_3j_uno
	COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0261968 0.000341549 uno
	)
	if(${HepMC_FOUND})
	file(READ "${tst_dir}/jet_config.yml" config)
	file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc")
	if(${rivet_FOUND})
	file(READ "${tst_dir}/jet_config_withHepMC.yml" config)
	file(WRITE "${tst_dir}/jet_config_withRivet.yml" "${config}\n\nanalysis:\n rivet: MC_XS\n output: tst")
	add_test(
	NAME t_main
	COMMAND HEJ ${tst_dir}/jet_config_withRivet.yml ${tst_dir}/2j.lhe.gz
	)
	else()
	add_test(
	NAME t_main
	COMMAND HEJ ${tst_dir}/jet_config_withHepMC.yml ${tst_dir}/2j.lhe.gz
	)
	endif()
	if(${HepMC_VERSION_MAJOR} GREATER 2)
	add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
	target_link_libraries(check_hepmc hejlib)
	add_test(
	NAME t_hepmc
	COMMAND check_hepmc tst.hepmc
	)
	endif()
	else()
	add_test(
	NAME t_main
	COMMAND HEJ ${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
	)
	add_test(
	NAME t_descriptions
	COMMAND test_descriptions
	)
	+add_test(
	+ NAME t_scale_arithmetics
	+ COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
	+ )
	diff --git a/Changes.md b/Changes.md
	index d6ccf10..ebdf244 100644
	--- a/Changes.md
	+++ b/Changes.md
	@@ -1,15 +1,20 @@
	# Version 2.0

	+## 2.X.0
	+
	+* Made `MatrixElement::tree_kin` and `MatrixElement::tree_param` public
	+* Allow multiplication and division of multiple scale functions e.g. `H_T/2*m_j1j2`
	+
	## 2.0.3

	* Fixed parsing of (numerical factor) * (base scale) in configuration
	* Don't change scale names, but sanitise Rivet output file names instead

	## 2.0.2

	* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
	`"/"` and `"*"` before)

	## 2.0.1

	* Fixed name of fixed-order generator in error message.
	diff --git a/include/HEJ/ScaleFunction.hh b/include/HEJ/ScaleFunction.hh
	index 5132a4e..42c2c99 100644
	--- a/include/HEJ/ScaleFunction.hh
	+++ b/include/HEJ/ScaleFunction.hh
	@@ -1,160 +1,174 @@
	/** \file
	* \brief Functions to calculate the (renormalisation and factorisation) scales for an event
	*
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <functional>
	#include <memory>
	#include <string>
	#include <utility>
	#include <vector>

	namespace HEJ{
	class Event;

	//! Class to calculate the scale associated with an event
	class ScaleFunction {
	public:
	//! Constructor
	/**
	* @param name Name of the scale choice (e.g. H_T)
	* @param fun Function used to calculate the scale
	*/
	ScaleFunction(std::string name, std::function<double(Event const &)> fun):
	name_{std::move(name)},
	fun_{std::move(fun)}
	{}

	//! Name of the scale choice
	std::string const & name() const {
	return name_;
	}

	//! Calculate the scale associated with an event
	double operator()(Event const & ev) const {
	return fun_(ev);
	}

	private:
	friend ScaleFunction operator*(double factor, ScaleFunction base_scale);
	friend ScaleFunction operator/(ScaleFunction base_scale, double denom);
	+ friend ScaleFunction operator*(ScaleFunction func1, ScaleFunction func2);
	+ friend ScaleFunction operator/(ScaleFunction func1, ScaleFunction func2);

	std::string name_;
	std::function<double(Event const &)> fun_;
	};

	//! Multiply a scale choice by a constant factor
	/**
	* For example, multiplying 0.5 and a scale function for H_T
	* will result in a scale function for H_T/2.
	*/
	ScaleFunction operator*(double factor, ScaleFunction base_scale);
	+ //! Multiply a scale choice by a second one
	+ /**
	+ * For example, multiplying H_T and m_j1j2
	+ * will result in a scale function for H_T*m_j1j2.
	+ */
	+ ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale);
	//! Divide a scale choice by a constant factor
	/**
	* For example, dividing a scale function for H_T by 2
	* will result in a scale function for H_T/2.
	*/
	ScaleFunction operator/(ScaleFunction base_scale, double denom);
	+ //! Divide a scale choice by a second one
	+ /**
	+ * For example, dividing a scale function for H_T by m_j1j2
	+ * will result in a scale function for H_T/m_j1j2.
	+ */
	+ ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom);

	//! Calculate \f$H_T\f$ for the input event
	/**
	* \f$H_T\f$ is the sum of the (scalar) transverse momenta of all
	* final-state particles
	*/
	double H_T(Event const &);
	//! The maximum of all (scalar) jet transverse momentum
	double max_jet_pt(Event const &);
	//! The invariant mass of the sum of all jet momenta
	double jet_invariant_mass(Event const &);
	//! Invariant mass of the two hardest jets
	double m_j1j2(Event const &);

	//! Functor that returns a fixed scale regardless of the input event
	class FixedScale {
	public:
	explicit FixedScale(double mu): mu_{mu} {}
	double operator()(Event const &) const {
	return mu_;
	}
	private:
	double mu_;
	};

	class ParameterDescription;

	//! Generate combinations of renormalisation and factorisation scales
	class ScaleGenerator{
	public:
	ScaleGenerator() = default;

	/** \brief Constructor
	* @param scale_functions_begin Iterator to first base scale
	* @param scale_functions_end Iterator past last base scale
	* @param scale_factors_begin Iterator to first scale factor
	* @param scale_factors_end Iterator past last scale factor
	* @param max_scale_ratio Maximum ratio between renormalisation
	* and factorisation scale
	*/
	template<class ScaleFunIterator, class FactorIterator>
	ScaleGenerator(
	ScaleFunIterator scale_functions_begin,
	ScaleFunIterator scale_functions_end,
	FactorIterator scale_factors_begin,
	FactorIterator scale_factors_end,
	double max_scale_ratio
	):
	scales_(scale_functions_begin, scale_functions_end),
	scale_factors_(scale_factors_begin, scale_factors_end),
	max_scale_ratio_{max_scale_ratio}
	{
	gen_descriptions();
	}

	/** \brief Constructor
	* @param scales Base scales
	* @param scale_factors Factors to multiply the base scales
	* @param max_scale_ratio Maximum ratio between renormalisation
	* and factorisation scale
	*/
	ScaleGenerator(
	std::vector<ScaleFunction> scales,
	std::vector<double> scale_factors,
	double max_scale_ratio
	):
	scales_(std::move(scales)),
	scale_factors_(std::move(scale_factors)),
	max_scale_ratio_{max_scale_ratio}
	{
	gen_descriptions();
	}

	/** \brief Adjust event parameters, adding scale variation
	*
	* The central renormalisation and factorisation scale of the returned
	* event is given be the first base scale passed to the constructor.
	* The scale variation (stored in event.variation()) is constructed as
	* follows: For each base scale according to the arguments of the
	* constructor we add one variation where both renormalisation and
	* factorisation scale are set according to the current base scale.
	* Then, all combinations where the base renormalisation and factorisation
	* scales are multiplied by one of the scale factors are added.
	* The case were both scales are multiplied by one is skipped.
	* Scale combinations where the ratio is larger than the maximum scale ratio
	* set in the constructor are likewise discarded.
	*/
	Event operator()(Event event) const;

	private:
	void gen_descriptions();

	std::vector<ScaleFunction> scales_;
	std::vector<double> scale_factors_;
	std::vector<std::shared_ptr<ParameterDescription>> descriptions_;
	double max_scale_ratio_;
	};

	}
	diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
	index 67a8e99..5fc7d7e 100644
	--- a/src/ScaleFunction.cc
	+++ b/src/ScaleFunction.cc
	@@ -1,149 +1,171 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/ScaleFunction.hh"

	#include <cassert>

	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"

	namespace HEJ{

	double H_T(Event const & ev){
	double result = 0.;
	for(size_t i = 0; i < ev.outgoing().size(); ++i){
	auto const decay_products = ev.decays().find(i);
	if(decay_products == end(ev.decays())){
	result += ev.outgoing()[i].perp();
	}
	else{
	for(auto const & particle: decay_products->second){
	result += particle.perp();
	}
	}
	}
	return result;
	}

	double max_jet_pt(Event const & ev) {
	return sorted_by_pt(ev.jets()).front().pt();
	}

	double jet_invariant_mass(Event const & ev) {
	fastjet::PseudoJet sum;
	for(const auto & jet: ev.jets()) sum+=jet;
	return sum.m();
	}

	double m_j1j2(Event const & ev) {
	const auto jets = sorted_by_pt(ev.jets());
	assert(jets.size() >= 2);
	return (jets[0] + jets[1]).m();
	}

	ScaleFunction operator*(double factor, ScaleFunction base_scale) {
	base_scale.name_.insert(0, std::to_string(factor) + '*');
	auto new_fun =
	[factor,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
	return factor*fun(ev);
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	+ ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale) {
	+ base_scale.name_.insert(0, factor.name_ + '*');
	+ auto new_fun =
	+ [fun1{std::move(factor.fun_)}, fun2{std::move(base_scale.fun_)}, name{base_scale.name_}]
	+ (HEJ::Event const & ev) {
	+ return fun1(ev)*fun2(ev);
	+ };
	+ base_scale.fun_ = std::move(new_fun);
	+ return base_scale;
	+ }
	+
	ScaleFunction operator/(ScaleFunction base_scale, double denom) {
	base_scale.name_.append('/' + std::to_string(denom));
	auto new_fun =
	[denom,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
	return fun(ev)/denom;
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	+ ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom) {
	+ base_scale.name_.append('/' + denom.name_);
	+ auto new_fun =
	+ [fun2{std::move(denom.fun_)}, fun1{std::move(base_scale.fun_)}]
	+ (HEJ::Event const & ev) {
	+ return fun1(ev)/fun2(ev);
	+ };
	+ base_scale.fun_ = std::move(new_fun);
	+ return base_scale;
	+ }
	+
	// TODO: significant logic duplication with operator()
	void ScaleGenerator::gen_descriptions() {
	if(scales_.empty()) {
	throw std::logic_error{"Need at least one scale"};
	}
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(scales_.front().name(), 1., 1.)
	);

	for(auto & base_scale: scales_){
	const auto base_name = base_scale.name();
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(base_name, 1., 1.)
	);
	//multiplicative scale variation
	for(double mur_factor: scale_factors_){
	for(double muf_factor: scale_factors_){
	if(muf_factor == 1. && mur_factor == 1.) continue;
	if(
	mur_factor/muf_factor < 1/max_scale_ratio_
	\|\| mur_factor/muf_factor > max_scale_ratio_
	) continue;
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(
	base_name, mur_factor, muf_factor
	)
	);
	}
	}
	}
	}

	Event ScaleGenerator::operator()(Event ev) const {
	if(! ev.variations().empty()) {
	throw std::invalid_argument{"Input event already has scale variation"};
	}
	assert(!scales_.empty());
	assert(!descriptions_.empty());

	size_t descr_idx = 0;
	const double mu_central = (scales_.front())(ev);
	ev.central().mur = mu_central;
	ev.central().muf = mu_central;
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.central().description = descriptions_[descr_idx++];

	// check if we are doing scale variation at all
	if(scales_.size() == 1 && scale_factors_.empty()) return ev;

	for(auto & base_scale: scales_){
	const double mu_base = base_scale(ev);
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.variations().emplace_back(
	EventParameters{
	mu_base, mu_base, ev.central().weight, descriptions_[descr_idx++]
	}
	);
	//multiplicative scale variation
	for(double mur_factor: scale_factors_){
	const double mur = mu_base*mur_factor;
	for(double muf_factor: scale_factors_){
	if(muf_factor == 1. && mur_factor == 1.) continue;
	const double muf = mu_base*muf_factor;
	if(
	mur/muf < 1/max_scale_ratio_
	\|\| mur/muf > max_scale_ratio_
	) continue;
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.variations().emplace_back(
	EventParameters{
	mur, muf, ev.central().weight, descriptions_[descr_idx++]
	}
	);
	}
	}
	};
	return ev;
	}

	}
	diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
	index 5336e8d..dc2e36c 100644
	--- a/src/YAMLreader.cc
	+++ b/src/YAMLreader.cc
	@@ -1,469 +1,465 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \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/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",
	"FKL", "unordered", "non-HEJ",
	"scales", "scale factors", "max scale ratio", "import scales",
	"log correction", "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},
	{"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}
	};
	const auto res = known.find(name);
	if(res == known.end()){
	throw std::invalid_argument("Unknown event treatment " + name);
	}
	return res->second;
	}

	} // namespace anonymous

	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("*/");
	+ // parse from right to left => a/b/c gives (a/b)/c
	+ const size_t delim = scale_fun.find_last_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);
	+ return parse_ScaleFunction(first, known)
	+ / parse_ScaleFunction(second, known);
	}
	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(second, known);
	- }
	+ return parse_ScaleFunction(first, known)
	+ * parse_ScaleFunction(second, 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");
	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
	diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
	new file mode 100644
	index 0000000..015dd6a
	--- /dev/null
	+++ b/t/test_scale_arithmetics.cc
	@@ -0,0 +1,87 @@
	+
	+// Generic tester for the ME for a given set of PSP
	+// reference weights and PSP (as LHE file) have to be given as _individual_ files
	+
	+#include <fstream>
	+
	+#include "LHEF/LHEF.h"
	+
	+#include "HEJ/EventReweighter.hh"
	+#include "HEJ/make_RNG.hh"
	+#include "HEJ/Event.hh"
	+#include "HEJ/YAMLreader.hh"
	+#include "HEJ/stream.hh"
	+
	+constexpr double alpha_s = 0.118;
	+constexpr double ep = 1e-13;
	+
	+void dump(HEJ::Event const & ev){
	+ {
	+ LHEF::Writer writer{std::cout};
	+ std::cout << std::setprecision(6);
	+ writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
	+ writer.writeEvent();
	+ }
	+ std::cout << "Rapidity ordering:\n";
	+ for(const auto & part: ev.outgoing()){
	+ std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
	+ }
	+}
	+
	+int main(int argn, char** argv){
	+ if(argn != 3){
	+ std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
	+ return EXIT_FAILURE;
	+ }
	+ HEJ::Config config = HEJ::load_config(argv[1]);
	+ config.scales = HEJ::to_ScaleConfig(
	+ YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2H_T/2, H_T/22, H_T/2/24, H_TH_T/H_T]")
	+ );
	+
	+ HEJ::istream in{argv[2]};
	+ LHEF::Reader reader{in};
	+
	+ auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
	+
	+ HEJ::ScaleGenerator scale_gen{
	+ config.scales.base,
	+ config.scales.factors,
	+ config.scales.max_ratio
	+ };
	+
	+ HEJ::EventReweighter resum{
	+ reader.heprup,
	+ std::move(scale_gen),
	+ to_EventReweighterConfig(config),
	+ *ran
	+ };
	+
	+ size_t i = 0;
	+ while(reader.readEvent()){
	+ ++i;
	+
	+ HEJ::Event event{
	+ HEJ::UnclusteredEvent{reader.hepeup},
	+ config.resummation_jets.def,
	+ config.resummation_jets.min_pt
	+ };
	+
	+ auto resummed = resum.reweight(event, config.trials);
	+ for(auto && ev: resummed) {
	+ for(auto &&var: ev.variations()) {
	+ if(std::abs(var.muf - ev.central().muf) > ep) {
	+ std::cerr
	+ << std::setprecision(15)
	+ << "unequal scales: " << var.muf
	+ << " != " << ev.central().muf << '\n'
	+ << "in resummed event:\n";
	+ dump(ev);
	+ std::cerr << "\noriginal event:\n";
	+ dump(event);
	+ return EXIT_FAILURE;
	+ }
	+ }
	+ }
	+ }
	+
	+}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 11:00 PM (1 d, 20 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242767
Default Alt Text: (39 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions