No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	diff --git a/config.yml b/config.yml
	index 410b203..832552d 100644
	--- a/config.yml
	+++ b/config.yml
	@@ -1,134 +1,137 @@
	## Number of attempted resummation phase space points for each input event
	trials: 10

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

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

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

	## Central scale choice or choices
	#
	## multiple scales are allowed, e.g.
	# scales: [125, max jet pperp, H_T/2, 2*jet invariant mass, m_j1j2]
	scales: 91.188

	## Factors by which the central scales should be multiplied
	## renormalisation and factorisation scales are varied independently
	#
	# scale factors: [0.5, 0.7071, 1, 1.41421, 2]

	## Maximum ratio between renormalisation and factorisation scale
	#
	# max scale ratio: 2.0001

	## Import scale setting functions
	#
	# import scales:
	# lib_my_scales.so: [scale0,scale1]

	## Unweighting setting
	## remove to obtain weighted events
	# unweight:
	# # type of unweighting (one of 'weighted', 'resummation', 'partial')
	# type: partial
	# trials: 10000
	# max deviation: 0

	## Event output files
	#
	# the supported formats are
	# - Les Houches (suffix .lhe)
	# - HepMC2 (suffix .hepmc2)
	# - HepMC3 (suffix .hepmc3 or .hepmc)
	# - HDF5 (suffix .hdf5)
	#
	## An output file's format is deduced either automatically from the suffix
	## or from an explicit specification, e.g.
	## - Les Houches: outfile
	#
	event output:
	- HEJ.lhe
	# - HEJ_events.hepmc

	## Analyses
	#
	# analyses:
	## Rivet analysis
	# - rivet: MC_XS # rivet analysis name
	# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
	## Custom analysis
	# - plugin: /path/to/libmyanalysis.so
	# my analysis parameter: some value

	## Selection of random number generator and seed
	## The choices are
	## - mixmax (seed is an integer)
	## - ranlux64 (seed is a filename containing parameters)
	random generator:
	name: mixmax
	# seed: 1

	## Whether or not to include higher order logs
	log correction: false

	## Truncate higher-order corrections at NLO
	NLO truncation:
	enabled: false
	# nlo order: 2

	+## Only keep low pt contributions
	+# lowpt only: false
	+
	## Vacuum expectation value
	vev: 246.2196508

	## Properties of the weak gauge bosons
	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	## Parameters for Higgs-gluon couplings
	## This requires compilation with QCDloop
	#
	# Higgs coupling:
	# use impact factors: false
	# mt: 174
	# include bottom: true
	# mb: 4.7

	## ---------------------------------------------------------------------- ##
	## The following settings are only intended for advanced users. ##
	## Please DO NOT SET them unless you know exactly what you are doing! ##
	## ---------------------------------------------------------------------- ##
	#
	## Maximum soft transverse momentum fraction in any tagging jets, e.g.
	## extremal or qqbar jet
	# soft pt regulator: 0.1
	#
	## Minimum transverse momentum of extremal partons
	## deprecated: use "soft pt regulator" instead
	# min extparton pt: 30
	#
	## deprecated: this cot directly replaced by "soft pt regulator"
	# max ext soft pt fraction: 0.1
	#
	# max events: -1 # Maximal number of fixed order Events to process
	# regulator parameter: 0.2 # The regulator lambda for the subtraction terms
	diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
	index d9b6ad9..67ff219 100644
	--- a/include/HEJ/Config.hh
	+++ b/include/HEJ/Config.hh
	@@ -1,284 +1,289 @@
	/** \file
	* \brief HEJ 2 configuration parameters
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2022
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <map>
	#include <string>
	#include <utility>
	#include <optional>
	#include <vector>

	#include "fastjet/JetDefinition.hh"
	#include "yaml-cpp/yaml.h"

	#include "HEJ/Constants.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/Fraction.hh"
	#include "HEJ/HiggsCouplingSettings.hh"
	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/output_formats.hh"

	namespace HEJ {

	//! Jet parameters
	struct JetParameters{
	fastjet::JetDefinition def; /*< Jet Definition /
	double min_pt{}; /*< Minimum Jet Transverse Momentum /
	};

	//! Settings for scale variation
	struct ScaleConfig{
	//! Base scale choices
	std::vector<ScaleFunction> base;
	//! Factors for multiplicative scale variation
	std::vector<double> factors;
	//! Maximum ratio between renormalisation and factorisation scale
	double max_ratio{};
	};

	//! Settings for random number generator
	struct RNGConfig {
	//! Random number generator name
	std::string name;
	//! Optional initial seed
	std::optional<std::string> seed;
	};

	//! Settings for partial unweighting
	struct PartialUnweightConfig {
	//! Number of trials for training
	size_t trials;
	//! Maximum distance in standard deviations from mean logarithmic weight
	double max_dev;
	};

	//! Settings for HEJ@NLO
	struct NLOConfig {
	//! Settings for HEJ@NLO Truncation
	bool enabled = false;
	//! NLO Born number of jets
	size_t nj = 2;
	};

	/**! Possible treatments for fixed-order input events.
	*
	* The program will decide on how to treat an event based on
	* the value of this enumeration.
	*/
	enum class EventTreatment{
	reweight, /*< Perform resummation /
	keep, /*< Keep the event /
	discard, /*< Discard the event /
	};

	//! Container to store the treatments for various event types
	using EventTreatMap = std::map<event_type::EventType, EventTreatment>;

	//! Possible setting for the event weight
	enum class WeightType{
	weighted, //!< weighted events
	unweighted_resum, //!< unweighted only resummation part
	partially_unweighted //!< mixed weighted and unweighted
	};

	/**! Input parameters.
	*
	* This struct handles stores all configuration parameters
	* needed in a HEJ 2 run.
	*
	* \internal To add a new option:
	* 1. Add a member to the Config struct.
	* 2. Inside "src/YAMLreader.cc":
	* - Add the option name to the "supported" Node in
	* get_supported_options.
	* - Initialise the new Config member in to_Config.
	* The functions set_from_yaml (for mandatory options) and
	* set_from_yaml_if_defined (non-mandatory) may be helpful.
	* 3. Add a new entry (with short description) to config.yaml
	* 4. Update the user documentation in "doc/Sphinx/"
	*/
	struct Config {
	//! %Parameters for scale variation
	ScaleConfig scales;
	//! Resummation jet properties
	JetParameters resummation_jets;
	//! Fixed-order jet properties
	JetParameters fixed_order_jets;
	//! Minimum transverse momentum for extremal partons
	//! \deprecated This will be removed in future versions.
	//! Use \ref soft_pt_regulator instead.
	double min_extparton_pt = 0.;
	//! \deprecated This is equivalent to\ref soft_pt_regulator
	//! and will be removed in future versions.
	std::optional<Fraction<double>> max_ext_soft_pt_fraction{};
	//! @brief Maximum transverse momentum fraction from soft radiation in any
	//! tagging jet (e.g. extremal or qqbar jet)
	Fraction<double> soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR };
	//! The regulator lambda for the subtraction terms
	double regulator_lambda = CLAMBDA;
	//! Number of resummation configurations to generate per fixed-order event
	size_t trials{};
	//! Maximal number of events
	std::optional<size_t> max_events;
	//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
	bool log_correction{};
	//! Event output files names and formats
	std::vector<OutputFile> output;
	//! Parameters for random number generation
	RNGConfig rng;
	//! Map to decide what to do for different event types
	EventTreatMap treat;
	//! %Parameters for custom analysis
	//! @deprecated use analyses_parameters instead
	YAML::Node analysis_parameters;
	//! %Parameters for custom analyses
	std::vector<YAML::Node> analyses_parameters;
	//! Settings for effective Higgs-gluon coupling
	HiggsCouplingSettings Higgs_coupling;
	//! elector weak parameters
	EWConstants ew_parameters;
	//! Type of event weight e.g. (un)weighted
	WeightType weight_type;
	//! Settings for partial unweighting
	std::optional<PartialUnweightConfig> unweight_config;
	//! HEJ@NLO settings
	NLOConfig nlo;
	+ //! LowPT settings
	+ bool lowpt = false;
	};

	//! Configuration options for the PhaseSpacePoint class
	struct PhaseSpacePointConfig {
	PhaseSpacePointConfig() = default;
	PhaseSpacePointConfig(
	JetParameters jet_param,
	NLOConfig nlo,
	double min_extparton_pt = 0.,
	Fraction<double> soft_pt_regulator =
	Fraction<double>{DEFAULT_SOFT_PT_REGULATOR}
	):
	jet_param{std::move(jet_param)},
	nlo{std::move(nlo)},
	min_extparton_pt{min_extparton_pt},
	soft_pt_regulator{std::move(soft_pt_regulator)}
	{}
	//! Properties of resummation jets
	JetParameters jet_param;
	//! HEJ@NLO settings
	NLOConfig nlo;
	//! Minimum transverse momentum for extremal partons
	//! \deprecated This will be removed in future versions.
	//! Use \ref soft_pt_regulator instead.
	double min_extparton_pt = 0.;
	//! \deprecated This is equivalent to\ref soft_pt_regulator
	//! and will be removed in future versions.
	std::optional<Fraction<double>> max_ext_soft_pt_fraction{};
	//! @brief Maximum transverse momentum fraction from soft radiation in any
	//! tagging jet (e.g. extremal or qqbar jet)
	Fraction<double> soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR };
	};

	//! Configuration options for the MatrixElement class
	struct MatrixElementConfig {
	MatrixElementConfig() = default;
	MatrixElementConfig(
	bool log_correction,
	HiggsCouplingSettings Higgs_coupling,
	EWConstants ew_parameters,
	NLOConfig nlo,
	Fraction<double> soft_pt_regulator = Fraction<double>{DEFAULT_SOFT_PT_REGULATOR},
	double regulator_lambda = CLAMBDA
	):
	log_correction{log_correction},
	Higgs_coupling{std::move(Higgs_coupling)},
	ew_parameters{std::move(ew_parameters)},
	nlo{std::move(nlo)},
	soft_pt_regulator{soft_pt_regulator},
	regulator_lambda{regulator_lambda}
	{}

	//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
	bool log_correction{};
	//! Settings for effective Higgs-gluon coupling
	HiggsCouplingSettings Higgs_coupling;
	//! elector weak parameters
	EWConstants ew_parameters;
	//! HEJ@NLO settings
	NLOConfig nlo;
	//! @brief Maximum transverse momentum fraction from soft radiation in any
	//! tagging jet (e.g. extremal or qqbar jet)
	Fraction<double> soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR };
	//! The regulator lambda for the subtraction terms
	double regulator_lambda = CLAMBDA;
	};

	//! Configuration options for the EventReweighter class
	struct EventReweighterConfig {
	//! Settings for phase space point generation
	PhaseSpacePointConfig psp_config;
	//! Settings for matrix element calculation
	MatrixElementConfig ME_config;
	//! Access properties of resummation jets
	JetParameters & jet_param() {
	return psp_config.jet_param;}
	//! Access properties of resummation jets (const version)
	JetParameters const & jet_param() const {
	return psp_config.jet_param;}
	//! Treatment of the various event types
	EventTreatMap treat;
	+ //! Option to only keep lowpt contribution
	+ bool lowpt = false;
	};

	/**! Extract PhaseSpacePointConfig from Config
	*
	* \internal We do not provide a PhaseSpacePointConfig constructor from Config
	* so that PhaseSpacePointConfig remains an aggregate.
	* This faciliates writing client code (e.g. the HEJ fixed-order generator)
	* that creates a PhaseSpacePointConfig without a Config object.
	*
	* @see to_MatrixElementConfig, to_EventReweighterConfig
	*/
	inline
	PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) {
	return {
	conf.resummation_jets,
	conf.nlo,
	conf.min_extparton_pt,
	conf.max_ext_soft_pt_fraction?*conf.max_ext_soft_pt_fraction
	:conf.soft_pt_regulator
	};
	}

	/**! Extract MatrixElementConfig from Config
	*
	* @see to_PhaseSpacePointConfig, to_EventReweighterConfig
	*/
	inline
	MatrixElementConfig to_MatrixElementConfig(Config const & conf) {
	return {
	conf.log_correction,
	conf.Higgs_coupling,
	conf.ew_parameters,
	conf.nlo,
	conf.soft_pt_regulator,
	conf.regulator_lambda
	};
	}

	/**! Extract EventReweighterConfig from Config
	*
	* @see to_PhaseSpacePointConfig, to_MatrixElementConfig
	*/
	inline
	EventReweighterConfig to_EventReweighterConfig(Config const & conf) {
	return {
	to_PhaseSpacePointConfig(conf),
	to_MatrixElementConfig(conf),
	- conf.treat
	+ conf.treat,
	+ conf.lowpt
	};
	}

	} // namespace HEJ
	diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
	index 5c7c8cf..c734d33 100644
	--- a/src/EventReweighter.cc
	+++ b/src/EventReweighter.cc
	@@ -1,289 +1,290 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/EventReweighter.hh"

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <cstddef>
	#include <functional>
	#include <string>
	#include <unordered_map>
	#include <utility>

	#include "fastjet/ClusterSequence.hh"
	#include "fastjet/PseudoJet.hh"

	#include "LHEF/LHEF.h"

	#include "HEJ/Event.hh"
	#include "HEJ/Fraction.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PhaseSpacePoint.hh"
	#include "HEJ/exceptions.hh"

	namespace HEJ {

	EventReweighter::EventReweighter(
	LHEF::HEPRUP const & heprup,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	std::shared_ptr<RNG> ran
	):
	EventReweighter{
	Beam{
	heprup.EBMUP.first,
	{{
	static_cast<ParticleID>(heprup.IDBMUP.first),
	static_cast<ParticleID>(heprup.IDBMUP.second)
	}}
	},
	heprup.PDFSUP.first,
	std::move(scale_gen),
	std::move(conf),
	std::move(ran)
	}
	{
	if(heprup.EBMUP.second != E_beam_){
	throw std::invalid_argument(
	"asymmetric beam: " + std::to_string(E_beam_)
	+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
	);
	}
	if(heprup.PDFSUP.second != pdf_.id()){
	throw std::invalid_argument(
	"conflicting PDF ids: " + std::to_string(pdf_.id())
	+ " vs. " + std::to_string(heprup.PDFSUP.second)
	);
	}
	}

	EventReweighter::EventReweighter(
	Beam const & beam,
	int pdf_id,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	std::shared_ptr<RNG> ran
	):
	param_{std::move(conf)},
	E_beam_{beam.E},
	pdf_{pdf_id, beam.type.front(), beam.type.back()},
	MEt2_{
	[this](double mu){ return pdf_.Halphas(mu); },
	param_.ME_config
	},
	scale_gen_{std::move(scale_gen)},
	ran_{std::move(ran)}
	{
	// legacy code: override new variable with old
	if(param_.psp_config.max_ext_soft_pt_fraction){
	param_.psp_config.soft_pt_regulator = *param_.psp_config.max_ext_soft_pt_fraction;
	param_.psp_config.max_ext_soft_pt_fraction = {};
	}
	assert(ran_);
	}

	PDF const & EventReweighter::pdf() const{
	return pdf_;
	}

	std::vector<Event> EventReweighter::reweight(
	Event const & input_ev, std::size_t num_events
	){
	auto res_events{ gen_res_events(input_ev, num_events) };
	if(res_events.empty()) return {};
	for(auto & event: res_events) event = scale_gen_(std::move(event));
	return rescale(input_ev, std::move(res_events));
	}

	EventTreatment EventReweighter::treatment(EventType type) const {
	return param_.treat.at(type);
	}

	std::vector<Event> EventReweighter::gen_res_events(
	Event const & ev,
	std::size_t phase_space_points
	){
	assert(ev.variations().empty());
	status_.clear();

	switch(treatment(ev.type())){
	case EventTreatment::discard: {
	status_.emplace_back(StatusCode::discard);
	return {};
	}
	case EventTreatment::keep:
	if(! jets_pass_resummation_cuts(ev)) {
	status_.emplace_back(StatusCode::failed_resummation_cuts);
	return {};
	}
	else {
	status_.emplace_back(StatusCode::good);
	return {ev};
	}
	default:;
	}
	const double Born_shat = shat(ev);

	std::vector<Event> resummation_events;
	status_.reserve(phase_space_points);
	for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
	PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
	status_.emplace_back(psp.status());
	assert(psp.status() != StatusCode::unspecified);
	if(psp.status() != StatusCode::good) continue;
	assert(psp.weight() != 0.);
	if(psp.incoming()[0].E() > E_beam_ \|\| psp.incoming()[1].E() > E_beam_) {
	status_.back() = StatusCode::too_much_energy;
	continue;
	}

	resummation_events.emplace_back(
	to_EventData( std::move(psp) ).cluster(
	param_.jet_param().def, param_.jet_param().min_pt
	)
	);
	auto & new_event = resummation_events.back();
	assert( new_event.valid_hej_state(
	param_.psp_config.soft_pt_regulator,
	param_.psp_config.min_extparton_pt ) );
	if( new_event.type() != ev.type() ) {
	throw std::logic_error{
	"Resummation Event does not match Born event: "
	+ name(new_event.type())
	+ " != "
	+ name(ev.type())
	};
	}
	new_event.generate_colours(*ran_);
	assert(new_event.variations().empty());
	new_event.central().mur = ev.central().mur;
	new_event.central().muf = ev.central().muf;
	const double resum_shat = shat(new_event);
	new_event.central().weight = ev.central().weightBorn_shat*Born_shat/
	(phase_space_pointsresum_shatresum_shat);

	- // Option to only keep lowpt events
	- // Keep only events where there is a fixed order event with at least 1
	- // jet below the resummation jet pt but all resummation jets are above
	- // the resummation jet pt
	- bool passes_pt_cut = true;
	- bool jet_below_cut = false;
	+ if(param_.lowpt){
	+ // Keep only events where there is a fixed order event with at least 1
	+ // jet below the resummation jet pt but all resummation jets are above
	+ // the resummation jet pt
	+ bool passes_pt_cut = true;
	+ bool jet_below_cut = false;

	- for (size_t jet_it = 0; jet_it < ev.jets().size(); ++jet_it){
	- jet_below_cut = jet_below_cut \|\| ev.jets()[jet_it].pt()
	- < param_.jet_param().min_pt;
	- passes_pt_cut = passes_pt_cut && new_event.jets()[jet_it].pt()
	- > param_.jet_param().min_pt;
	- }
	+ for (size_t jet_it = 0; jet_it < ev.jets().size(); ++jet_it){
	+ jet_below_cut = jet_below_cut \|\| ev.jets()[jet_it].pt()
	+ < param_.jet_param().min_pt;
	+ passes_pt_cut = passes_pt_cut && new_event.jets()[jet_it].pt()
	+ > param_.jet_param().min_pt;
	+ }

	- if (!passes_pt_cut \|\| !jet_below_cut) {
	- new_event.central().weight *= 0;
	- status_.emplace_back(StatusCode::too_much_energy);
	- return {};
	+ if (!passes_pt_cut \|\| !jet_below_cut) {
	+ new_event.central().weight *= 0;
	+ status_.emplace_back(StatusCode::too_much_energy);
	+ return {};
	+ }
	}

	}
	return resummation_events;
	}

	std::vector<Event> EventReweighter::rescale(
	Event const & Born_ev,
	std::vector<Event> events
	) const{
	const double Born_pdf = pdf_factors(Born_ev).central;
	const double Born_ME = tree_matrix_element(Born_ev);

	for(auto & cur_event: events){
	const auto pdf = pdf_factors(cur_event);
	assert(pdf.variations.size() == cur_event.variations().size());
	const auto ME = matrix_elements(cur_event);
	assert(ME.variations.size() == cur_event.variations().size());
	cur_event.parameters() = pdfME/(Born_pdf*Born_ME);
	}
	return events;
	}

	bool EventReweighter::jets_pass_resummation_cuts(
	Event const & ev
	) const{
	const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
	fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
	return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
	}

	Weights EventReweighter::pdf_factors(Event const & ev) const{
	auto const & a = ev.incoming().front();
	auto const & b = ev.incoming().back();
	const double xa = a.p.e()/E_beam_;
	const double xb = b.p.e()/E_beam_;

	Weights result;
	std::unordered_map<double, double> known_pdf;
	result.central =
	pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
	pdf_.pdfpt(1,xb,ev.central().muf,b.type);
	known_pdf.emplace(ev.central().muf, result.central);

	result.variations.reserve(ev.variations().size());
	for(auto const & ev_param: ev.variations()){
	const double muf = ev_param.muf;
	auto cur_pdf = known_pdf.find(muf);
	if(cur_pdf == known_pdf.end()){
	cur_pdf = known_pdf.emplace(
	muf,
	pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
	).first;
	}
	result.variations.emplace_back(cur_pdf->second);
	}
	assert(result.variations.size() == ev.variations().size());
	return result;
	}

	Weights
	EventReweighter::matrix_elements(Event const & ev) const{
	assert(param_.treat.count(ev.type()) > 0);
	if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev);
	}

	return MEt2_(ev);
	}

	double EventReweighter::tree_matrix_element(Event const & ev) const{
	assert(ev.variations().empty());
	assert(param_.treat.count(ev.type()) > 0);
	if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev).central;
	}
	return MEt2_.tree(ev).central;
	}

	Weights
	EventReweighter::fixed_order_scale_ME(Event const & ev) const{
	int alpha_s_power = 0;
	for(auto const & part: ev.outgoing()){
	if(is_parton(part))
	++alpha_s_power;
	else if(part.type == pid::Higgs) {
	alpha_s_power += 2;
	}
	// nothing to do for other uncoloured particles
	}

	Weights result;
	result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
	for(auto const & var: ev.variations()){
	result.variations.emplace_back(
	std::pow(pdf_.Halphas(var.mur), alpha_s_power)
	);
	}
	return result;
	}

	} // namespace HEJ
	diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
	index bf14a98..a9229cf 100644
	--- a/src/YAMLreader.cc
	+++ b/src/YAMLreader.cc
	@@ -1,592 +1,593 @@
	/**
	* \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"
	+ "regulator parameter", "max events", "lowpt only"
	};
	// 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"}){
	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}
	};
	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 {
	{not_enough_jets, EventTreatment::discard},
	{bad_final_state, EventTreatment::discard},
	{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}
	};
	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"};
	}
	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, "lowpt only");

	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");
	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, 2:29 PM (1 d, 11 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804741
Default Alt Text: (45 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions