No OneTemporary
Actions

Size

39 KB

Subscribers

None

View Options

	diff --git a/CMakeLists.txt b/CMakeLists.txt
	index 2a6cf4c..b4d9ef8 100644
	--- a/CMakeLists.txt
	+++ b/CMakeLists.txt
	@@ -1,147 +1,147 @@
	cmake_minimum_required(VERSION 3.1 FATAL_ERROR)

	set(CMAKE_LEGACY_CYGWIN_WIN32 0)
	set(CMAKE_EXPORT_COMPILE_COMMANDS ON)

	project("Reversed HEJ" 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)

	## Add directories and find dependences
	include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
	set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")

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

	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.418e6 67247.8
	)
	add_test(
	NAME t_4j
	COMMAND check_res ${tst_dir}/4j.lhe.gz 551096 30572.5
	)
	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.0246366 0.001 uno
	)
	if(${HepMC_FOUND})
	if(${HepMC_VERSION_MAJOR} GREATER 2)
	file(READ "${tst_dir}/jet_config.yml" config)
	- file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc3")
	+ file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc")
	add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
	target_link_libraries(check_hepmc rhej)
	add_test(
	NAME t_main
	COMMAND rHEJ ${tst_dir}/jet_config_withHepMC.yml ${tst_dir}/2j.lhe.gz
	)
	add_test(
	NAME t_hepmc
	- COMMAND check_hepmc tst.hepmc3
	+ COMMAND check_hepmc tst.hepmc
	)
	else()
	add_test(
	NAME t_main
	COMMAND rHEJ ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
	)
	endif()
	endif()
	add_test(
	NAME t_lhe
	COMMAND check_lhe tst.lhe
	)
	diff --git a/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml
	index ec58822..2112682 100644
	--- a/FixedOrderGen/configFO.yml
	+++ b/FixedOrderGen/configFO.yml
	@@ -1,49 +1,49 @@
	trials : 200

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 230000

	process: p p => h 4j

	# fraction of events with two extremal emissions in one direction
	# that contain an unordered emission
	unordered fraction: 0.05

	scales: max jet pperp

	event output:
	- HEJFO.lhe
	# - RHEJ.lhe
	-# - RHEJ_events.hepmc3
	+# - RHEJ_events.hepmc

	Higgs properties:
	mass: 125
	width: 0.004165
	decays: {into: [photon, photon], branching ratio: 0.0023568762400521404}

	# analysis:
	# # to use a custom analysis
	# plugin: /home/andersen/HEJ/PURE/GitReverse/build/analysis-plugins/src/libVBF.so
	# min dy12: 2.8
	# min m12: 400
	# output: HEJFO.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/config.yml b/config.yml
	index 7edc809..95129b5 100644
	--- a/config.yml
	+++ b/config.yml
	@@ -1,80 +1,80 @@
	# 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: input

	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.hepmc3
	+ - 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/doc/sphinx/analyses.rst b/doc/sphinx/analyses.rst
	index f85a070..eab8ef0 100644
	--- a/doc/sphinx/analyses.rst
	+++ b/doc/sphinx/analyses.rst
	@@ -1,166 +1,166 @@
	.. _`Writing custom analyses`:

	Writing custom analyses
	=======================

	-Both reversed HEJ and HEJ FOG can generate HepMC 3 files, so you can
	+Both reversed HEJ and HEJ FOG can generate HepMC files, so you can
	always run a rivet analysis on these.

	If you want to write your own native analysis, the easiest way is to
	create a new .cc file, for example my_analysis.cc in the
	analysis-plugins/src folder.

	An analysis is a class that derives from the abstract :code:`Analysis`
	base class provided by reversed HEJ. It has to implement two public
	functions:

	* The :code:`pass_cuts` member function return true if and only if the
	given event (first argument) passes the analysis cuts

	* The :code:`fill` member function adds an event to the analysis, which
	for example can be used to fill histograms. Reversed HEJ and HEJ FOG
	will only give you events for which :code:`pass_cuts` has returned
	true.

	Both functions take two arguments: the resummation event generated by
	reversed HEJ and the original fixed-order input event. Usually, the
	second argument can be ignored. It can be used, for example, for
	implementing cuts that depend on the ratio of the weights between the
	fixed-order and the resummation event.

	In addition to the two member functions, there has to be a global
	make_analysis function that takes the analysis parameters in the form of
	a YAML :code:`Node` and returns a :code:`std::unique_ptr` to the
	Analysis.

	The following code creates the simplest conceivable analysis.::

	#include <memory> // for std::unique_ptr

	#include "RHEJ/Analysis.hh"

	class MyAnalysis: public RHEJ::Analysis {
	public:
	MyAnalysis(YAML::Node const & /* config */) {}

	void fill(
	RHEJ::Event const & /* event */,
	RHEJ::Event const & /* FO_event */
	) override{
	};

	bool pass_cuts(
	RHEJ::Event const & /* event */,
	RHEJ::Event const & /* FO_event */
	) override{
	return true;
	};

	};

	extern "C"
	__attribute__((visibility("default")))
	std::unique_ptr<RHEJ::Analysis> make_analysis(
	YAML::Node const & config
	){
	return std::make_unique<MyAnalysis>(config);
	}

	Reversed HEJ and HEJ FOG load analyses dynamically at run time. Since
	the loader cannot deal with the name mangling which is usually done by
	compilers, we have to specify :code:`extern "C"`.

	The standard compiler settings for analyses ensure that symbols (for
	example your class name) are not visible in the compiled analysis. This
	ensures that there can be no name clashes with the remaining
	code. However, since reversed HEJ and HEJ FOG have to be able to call
	:code:`make_analysis`, we use
	:code:`__attribute__((visibility("default")))` to ensure its visibility.

	To build the analysis go to some build directory and run

	.. code-block:: Bash

	cmake some/directory/analysis-plugins
	make my_analysis

	This will create the analysis library libmy_analysis.so in the src
	subdirectory. You can then use it in reversed HEJ or HEJ FOG by adding

	.. code-block:: YAML

	analysis:
	plugin: analysis/build/directory/src/libmy_analysis.so

	to the .yml configuration file.

	As a more interesting example, here is the code for an analysis that
	sums up the total cross section and prints the result to both standard
	output and a file specified in the .yml config with

	.. code-block:: YAML

	analysis:
	plugin: analysis/build/directory/src/libmy_analysis.so
	output: outfile

	To access the configuration at run time, reversed HEJ uses the yaml-cpp
	library; for more details see the `yaml-cpp tutorial
	<https://github.com/jbeder/yaml-cpp/wiki/Tutorial>`_. The analysis code
	itself is::

	#include <memory>
	#include <iostream>
	#include <fstream>
	#include <string>
	#include <cmath>

	#include "RHEJ/Analysis.hh"
	#include "RHEJ/Event.hh"

	#include "yaml-cpp/yaml.h"

	class MyAnalysis: public RHEJ::Analysis {
	public:
	MyAnalysis(YAML::Node const & config):
	xsection_{0.}, xsection_error_{0.},
	outfile_{config["output"].as<std::string>()}
	{}

	void fill(
	RHEJ::Event const & event,
	RHEJ::Event const & /* FO_event */
	) override{
	const double wt = event.central().weight;
	xsection_ += wt;
	xsection_error_ += wt*wt;
	};

	bool pass_cuts(
	RHEJ::Event const & /* event */,
	RHEJ::Event const & /* FO_event */
	) override{
	return true;
	};

	~MyAnalysis(){
	std::cout << "cross section: " << xsection_ << " +- "
	<< std::sqrt(xsection_error_) << "\n";
	std::ofstream fout{outfile_};
	fout << "cross section: " << xsection_ << " +- "
	<< std::sqrt(xsection_error_) << "\n";
	}

	private:
	double xsection_, xsection_error_;
	std::string outfile_;

	};

	extern "C"
	__attribute__((visibility("default")))
	std::unique_ptr<RHEJ::Analysis> make_analysis(
	YAML::Node const & config
	){
	return std::make_unique<MyAnalysis>(config);
	}
	diff --git a/doc/sphinx/rHEJ.rst b/doc/sphinx/rHEJ.rst
	index e7e282d..282f06d 100644
	--- a/doc/sphinx/rHEJ.rst
	+++ b/doc/sphinx/rHEJ.rst
	@@ -1,260 +1,260 @@
	.. _`Running reversed HEJ`:

	Running reversed HEJ
	====================

	Quick start
	-----------

	In order to run reversed HEJ, you need a configuration file and a file
	containing fixed-order events. A sample configuration is given by the
	:file:`config.yml` file distributed together with reversed HEJ. Events
	in the Les Houches Event File format can be generated with standard
	Monte Carlo generators like `MadGraph5_aMC@NLO
	<https://launchpad.net/mg5amcnlo>`_ or `Sherpa
	<https://sherpa.hepforge.org/trac/wiki>`_. The processes supported by
	reversed HEJ are

	- Pure multijet production
	- Production of a Higgs boson with jets

	..
	- TODO Production of a W boson with jets
	- TODO Production of a Z boson or photon with jets

	where at least two jets are required in each case. For the time being,
	only leading-order events are supported.

	After generating an event file :file:`events.lhe` adjust the parameters
	under the `fixed order jets`_ setting in :file:`config.yml` to the
	settings in the fixed-order generation. Resummation can then be added by
	running::

	rHEJ config.yml events.lhe

	Using the default settings, this will produce an output event file
	:file:`RHEJ.lhe` with events including high-energy resummation.

	.. _`reversed HEJ settings`:

	Settings
	--------

	Reversed HEJ configuration files follow the `YAML <http://yaml.org/>`_
	format. The following configuration parameters are supported:

	.. _`trials`:

	trials
	High-energy resummation is performed by generating a number of
	resummation phase space configurations corresponding to the input
	fixed-order event. This parameter specifies how many such
	configurations reversed HEJ should try to generate for each input
	event. Typical values vary between 10 and 100.

	.. _`min extparton pt`:

	min extparton pt
	Specifies the minimum transverse momentum in GeV of the most forward
	and the most backward parton. This setting is needed to regulate an
	otherwise uncancelled divergence. Its value should be slightly below
	the minimum transverse momentum of jets specified by `resummation
	jets: min pt`_. See also the `max ext soft pt fraction`_ setting.

	.. _`max ext soft pt fraction`:

	max ext soft pt fraction
	Specifies the maximum fraction that soft radiation can contribute to
	the transverse momentum of each the most forward and the most backward
	jet. Values between around 0.05 and 0.1 are recommended. See also the
	`min extparton pt`_ setting.

	.. _`fixed order jets`:

	fixed order jets
	This tag collects a number of settings specifying the jet definition
	in the event input. The settings should correspond to the ones used in
	the fixed-order Monte Carlo that generated the input events.

	.. _`fixed order jets: min pt`:

	min pt
	Minimum transverse momentum in GeV of fixed-order jets.

	.. _`fixed order jets: algorithm`:

	algorithm
	The algorithm used to define jets. Allowed settings are
	:code:`kt`, :code:`cambridge`, :code:`antikt`, :code:`genkt`,
	:code:`cambridge for passive`, :code:`genkt for passive`,
	:code:`ee kt`, :code:`ee genkt`. See the `FastJet
	<http://fastjet.fr/>`_ documentation for a description of these
	algorithms.

	.. _`fixed order jets: R`:

	R
	The R parameter used in the jet algorithm, roughly corresponding
	to the jet radius in the plane spanned by the rapidity and the
	azimuthal angle.

	.. _`resummation jets`:

	resummation jets
	This tag collects a number of settings specifying the jet definition
	in the observed, i.e. resummed events. These settings are optional, by
	default the same values as for the `fixed order jets`_ are assumed.

	.. _`resummation jets: min pt`:

	min pt
	Minimum transverse momentum in GeV of resummation jets. This
	should be around 25% larger than the minimum transverse momentum
	of fixed order jets set by `fixed order jets: min pt`_.

	.. _`resummation jets: algorithm`:

	algorithm
	The algorithm used to define jets. The reversed HEJ approach to
	resummation relies on properties of :code:`antikt` jets, so this
	value is strongly recommended. For a list of possible other
	values, see the `fixed order jets: algorithm`_ setting.

	.. _`resummation jets: R`:

	R
	The R parameter used in the jet algorithm.

	.. _`FKL`:

	FKL
	Specifies how to treat FKL events. The possible values are
	:code:`reweight` to enable resummation, :code:`keep` to keep the
	events as they are up to a possible change of renormalisation and
	factorisation scale, and :code:`discard` to discard these events.

	.. _`unordered`:

	unordered
	Specifies how to treat events with one emission that does not respect
	FKL ordering. The possible values are the same as for the `FKL`_
	setting, but :code:`reweight` may not be supported for all process
	types.


	.. TODO: rename to non-HEJ
	.. _`non-FKL`:

	non-FKL
	Specifies how to treat events where no resummation is possible. The
	allowed values are :code:`keep` to keep the events as they are up to
	a possible change of renormalisation and factorisation scale and
	:code:`discard` to discard these events.

	.. _`scales`:

	scales
	Specifies the renormalisation and factorisation scales for the output
	events. This can either be a single entry or a list :code:`[scale1,
	scale2, ...]`. For the case of a list the first entry defines the
	central scale. Possible values are fixed numbers to set the scale in
	GeV or the following:

	- :code:`input`: The same scale as in the input events
	- :code:`Ht`: The sum of the scalar transverse momenta of all
	final-state particles
	- :code:`max jet pperp`: The maximum transverse momentum of all jets
	- :code:`jet invariant mass`: Sum of the invariant masses of all jets

	Scales can be multiplied or divided by an overall factor, e.g. :code:`Ht/2`.

	.. _`scale factors`:

	scale factors
	A list of numeric factors by which each of the `scales`_ should be
	multiplied. Renormalisation and factorisation scales are varied
	independently. For example, a list with entries :code:`[0.5, 2]`
	would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`);
	(0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\
	:sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to
	the order of the final event weights.

	.. _`max scale ratio`:

	max scale ratio
	Specifies the maximum factor by which renormalisation and
	factorisation scales may difer. For a value of :code:`2` and the
	example given for the `scale factors`_ the scale choices
	(0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`)
	will be discarded.

	.. _`log correction`:

	log correction
	Whether to include corrections due to the evolution of the strong
	coupling constant in the virtual corrections. Allowed values are
	:code:`true` and :code:`false`.

	.. TODO: unweight

	.. _`event output`:

	event output
	Specifies the name of a single event output file or a list of such
	files. The file format is either specified explicitly or derived from
	the suffix. For example, :code:`events.lhe` or, equivalently
	:code:`Les Houches: events.lhe` generates an output event file
	:code:`events.lhe` in the Les Houches format. The supported formats
	are

	- :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches
	event file format.
	- - :code:`file.hepmc3` or :code:`HepMC: file`: The HepMC format.
	+ - :code:`file.hepmc` or :code:`HepMC: file`: The HepMC format.

	.. _`RanLux init`:

	RanLux init
	Specifies a file to be loaded for the initialisation of the RanLux64
	random number generation. See the `CLHEP documentation
	<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for
	details.

	.. _`analysis`:

	analysis
	Specifies the name and settings for a custom analysis library. The
	:code:`plugin` sub-entry should be set to the analysis file path. All
	further entries are passed on to the analysis. See
	:ref:`Writing custom analyses` for details.

	.. _`Higgs coupling`:

	Higgs coupling
	This collects a number of settings concerning the effective coupling
	of the Higgs boson to gluons. This is only relevant for the
	production process of a Higgs boson with jets and only supported if
	reversed HEJ was compiled with `QCDloop
	<https://github.com/scarrazza/qcdloop>`_ support.

	.. _`Higgs coupling: use impact factors`:

	use impact factors
	Whether to use impact factors for the coupling to the most forward
	and most backward partons. Impact factors imply the infinite
	top-quark mass limit.

	.. _`Higgs coupling: mt`:

	mt
	The value of the top-quark mass in GeV. If this is not specified,
	the limit of an infinite mass is taken.

	.. _`Higgs coupling: include bottom`:

	include bottom
	Whether to include the Higgs coupling to bottom quarks.

	.. _`Higgs coupling: mb`:

	mb
	The value of the bottom-quark mass in GeV.
	diff --git a/src/config.cc b/src/config.cc
	index 570953b..eef13e2 100644
	--- a/src/config.cc
	+++ b/src/config.cc
	@@ -1,462 +1,462 @@
	#include "RHEJ/config.hh"
	#include "RHEJ/exceptions.hh"

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

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

	}

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

	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 == "hepmc3") return FileFormat::HepMC;
	+ 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
	*/
	if(name != "analysis"){
	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 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 RHEJ{
	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;
	}

	// simple (as in non-composite) scale functions
	/**
	* An example for a simple scale function would be Ht,
	* Ht/2 is then composite (take Ht and then divide by 2)
	*/
	std::unique_ptr<ScaleFunction> make_simple_ScaleFunction_ptr(
	std::string const & scale_fun
	){
	using ret_type = std::unique_ptr<ScaleFunction>;
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	if(scale_fun == "input") return ret_type{new InputScales{}};
	if(scale_fun == "Ht") return ret_type{new Ht{}};
	if(scale_fun == "max jet pperp") return ret_type{new MaxJetPperp{}};
	if(scale_fun == "jet invariant mass"){
	return ret_type{new JetInvariantMass{}};
	}
	try{
	const double scale = to_double(scale_fun);
	return ret_type{new FixedScale{scale}};
	} catch(std::invalid_argument const &){}
	throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
	}

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

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

	std::unique_ptr<ScaleFunction> make_ScaleFunction_ptr(
	std::string const & scale_fun
	){
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	const size_t delim = scale_fun.find_first_of("*/");
	if(delim == scale_fun.npos){
	return make_simple_ScaleFunction_ptr(scale_fun);
	}
	const std::string first = trim_back(std::string{scale_fun, 0, delim});
	const std::string second = trim_front(std::string{scale_fun, delim+1});
	double factor;
	std::unique_ptr<ScaleFunction> fun;
	if(scale_fun[delim] == '/'){
	factor = 1/to_double(second);
	fun = make_simple_ScaleFunction_ptr(first);
	}
	else{
	assert(scale_fun[delim] == '*');
	try{
	factor = to_double(second);
	fun = make_simple_ScaleFunction_ptr(first);
	}
	catch(std::invalid_argument const &){
	factor = to_double(first);
	fun = make_simple_ScaleFunction_ptr(second);
	}
	}
	assert(fun != nullptr);
	return std::unique_ptr<ScaleFunction>{
	new Product{factor, std::move(fun)}
	};
	}

	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.scale_gen = ScaleGenerator{to_ScaleConfig(yaml)};
	config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
	return config;
	}

	} // anonymous namespace

	ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
	ScaleConfig config;
	std::vector<std::string> scales;
	set_from_yaml(scales, yaml, "scales");
	config.scales.reserve(scales.size());
	std::transform(
	begin(scales), end(scales), std::back_inserter(config.scales),
	make_ScaleFunction_ptr
	);
	set_from_yaml_if_defined(config.scale_factors, yaml, "scale factors");
	config.max_scale_ratio = std::numeric_limits<double>::infinity();
	set_from_yaml_if_defined(config.max_scale_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;
	}
	}

	Event ScaleGenerator::operator()(Event ev) const{
	assert(ev.variations().empty());
	assert(! settings_.scales.empty());

	auto const & scales = settings_.scales;
	auto const & scale_factors = settings_.scale_factors;
	const double max_scale_ratio = settings_.max_scale_ratio;

	// check if we are doing scale variation at all
	if(scales.size() == 1 && scale_factors.empty()){
	ev.central() = (*scales.front())(ev);
	return ev;
	}

	for(auto && base_scale: scales){
	const EventParameters cur_base = (*base_scale)(ev);
	ev.variations().emplace_back(cur_base);
	//multiplicative scale variation
	for(double mur_factor: scale_factors){
	const double mur = cur_base.mur*mur_factor;
	for(double muf_factor: scale_factors){
	if(muf_factor == 1. && mur_factor == 1.) continue;
	const double muf = cur_base.muf*muf_factor;
	if(mur/muf < 1/max_scale_ratio \|\| mur/muf > max_scale_ratio) continue;
	ev.variations().emplace_back(
	EventParameters{mur, muf, cur_base.weight}
	);
	}
	}
	};
	ev.central() = (*scales.front())(ev);
	return ev;
	}

	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions