Page MenuHomeHEPForge
Files F8725491

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
index b600d21..c7278b3 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,343 +1,360 @@
/** \file config.hh
* \brief The file which handles the configuration file parameters
*
* Contains the JetParameters Struct and ScaleConfig Struct. Also
* contains the ScaleGenerator Class and the EventTreatment class.
* Config struct is also defined here. The configuration files
* parameters are read and then stored within this objects.
*/
#pragma once
#include <string>
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/ScaleFunction.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/output_formats.hh"
#include "RHEJ/event_types.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "RHEJ/exceptions.hh"
#include "RHEJ/utility.hh"
#include "yaml-cpp/yaml.h"
namespace RHEJ{
/**! \struct JetParameters config.hh "include/RHEJ/config.hh"
* \brief Contains Jet Definition and Minimum Pt to be a Jet.
*
* Contains the Fastjet definition of a Jet and a threshold (minimum) value for pt
*/
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Pt */
};
/**! \struct ScaleConfig config.hh "include/RHEJ/config.hh"
* \brief Sets up the possible choices for scale variation?
*
* There is a vector which contains the possible scale choices and
* also vector of doubles with the scale factors along with a max
* scale ratio.
*/
struct ScaleConfig{
std::vector<std::unique_ptr<ScaleFunction>> scales; /**< Vector of possible Scale choices */
std::vector<double> scale_factors; /**< Vector of possible Scale Factors */
double max_scale_ratio; /**< Maximum ratio for the scales */
};
/**! \class ScaleGenerator config.hh "include/RHEJ/config.hh"
* \brief A Class used to generate scales
*
* This class has a clustered event class and scale config struct
* within it.
*/
class ScaleGenerator{
public:
ScaleGenerator() = default;
explicit ScaleGenerator(ScaleConfig && settings):
settings_{std::move(settings)}
{}
Event operator()(Event ev) const;
ScaleConfig const & settings() const;
private:
ScaleConfig settings_;
};
/**
* \brief Enumeration which deals with how to treat events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Reweigh the event */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
/**
* \brief An ordered Map with EventType keys and mapped values EventTreatment
*
* If called upon with EventTreatMap.at(EventType) it will return the treatment which has
* been specified for that EventType.
*/
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
- /**! \struct Config config.hh "include/RHEJ/config.hh"
- * \brief Config Struct for user input parameters.
+ /**! \class Config config.hh "include/RHEJ/config.hh"
+ * \brief Config Class for user input parameters.
*
- * The struct which handles the input card given by the user.
+ * The class which handles the input card given by the user.
*
* \internal To add a new option:
- * 1. Add a member to the Config struct.
+ * 1. Add a member to the Config Class.
* 2. Add the option name to the "supported" Node in
* get_supported_options.
* 3. 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.
* 4. Add a new entry (with short description) to config.yaml
*/
- struct Config {
- ScaleGenerator scale_gen; /**< Scale */
- JetParameters resummation_jets; /**< Resummation Jets */
- JetParameters fixed_order_jets; /**< Fixed-Order Jets */
- double min_extparton_pt; /**< Min External Parton Pt*/
- double max_ext_soft_pt_fraction; /**< Min External Parton Pt Fraction */
- int trials; /**< Number of resummation points to generate per FO */
- bool log_correction; /**< Log Correct On or Off */
- bool unweight; /**< Unweight On or Off */
- std::vector<OutputFile> output; /**< Output File Type */
- optional<std::string> RanLux_init; /**< Ranlux File Choice*/
- EventTreatMap treat; /**< Event Treat Map? */
- YAML::Node analysis_parameters; /**< Analysis Parameters */
+ class Config {
+ public:
+ // private:
+ ScaleGenerator _scale_gen; /**< Scale */
+ JetParameters _resummation_jets; /**< Resummation Jets */
+ JetParameters _fixed_order_jets; /**< Fixed-Order Jets */
+ double _min_extparton_pt; /**< Min External Parton Pt*/
+ double _max_ext_soft_pt_fraction; /**< Min External Parton Pt Fraction */
+ int _trials; /**< Number of resummation points to generate per FO */
+ bool _log_correction; /**< Log Correct On or Off */
+ bool _unweight; /**< Unweight On or Off */
+ std::vector<OutputFile> _output; /**< Output File Type */
+ optional<std::string> _RanLux_init; /**< Ranlux File Choice*/
+ EventTreatMap _treat; /**< Event Treat Map? */
+ YAML::Node _analysis_parameters; /**< Analysis Parameters */
+ //! Settings for Effective Higgs-Gluon Coupling
+ HiggsCouplingSettings _Higgs_coupling;
+ public:
+ const ScaleGenerator& scale_gen() const {return _scale_gen;}
+ JetParameters resummation_jets() const {return _resummation_jets;}
+ JetParameters fixed_order_jets() const {return _fixed_order_jets;}
+ double min_extparton_pt() const {return _min_extparton_pt;}
+ double max_ext_soft_pt_fraction() const {return _max_ext_soft_pt_fraction;}
+ int trials() const {return _trials;}
+ bool log_correction() const {return _log_correction;}
+ bool unweight() const {return _unweight;}
+ std::vector<OutputFile> output() const {return _output;}
+ optional<std::string> RanLux_init() const {return _RanLux_init;}
+ EventTreatMap treat() const {return _treat;}
+ YAML::Node analysis_parameters() const {return _analysis_parameters;}
//! Settings for Effective Higgs-Gluon Coupling
- HiggsCouplingSettings Higgs_coupling;
+ HiggsCouplingSettings Higgs_coupling() const {return _Higgs_coupling;}
};
//! Load configuration from file
/**
* @param config_file Name of the YAML configuration file
* @returns The reversed HEJ configuration
*/
Config load_config(std::string const & config_file);
//! Set option using the corresponding YAML entry
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* If the entry does not exist or has the wrong type or format
* an exception is thrown.
*
* For example
* @code
* set_from_yaml(foobar, yaml, "foo", "bar")
* @endcode
* is equivalent to
* @code
* foobar = yaml["foo"]["bar"].as<decltype(foobar)>()
* @endcode
* with improved diagnostics on errors.
*
* @see set_from_yaml_if_defined
*/
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
//! Set option using the corresponding YAML entry, if present
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* This function works similar to set_from_yaml, but does not
* throw any exception if the requested YAML entry does not exist.
*
* @see set_from_yaml
*/
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
);
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node, std::string const & entry
);
ScaleConfig to_ScaleConfig(YAML::Node const & yaml);
ParticleID to_ParticleID(std::string const & particle_name);
//! Check whether all options in configuration are supported
/**
* @param conf Configuration to be checked
* @param supported Tree of supported options
*
* If conf contains an entry that does not appear in supported
* an unknown_option exception is thrown. Sub-entries of "analysis"
* (if present) are not checked.
*
* @see unknown_option
*/
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
);
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml);
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml);
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml);
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml);
inline
void set_from_yaml(YAML::Node & setting, YAML::Node const & yaml){
setting = yaml;
}
template<typename Scalar>
void set_from_yaml(Scalar & setting, YAML::Node const & yaml){
assert(yaml);
if(!yaml.IsScalar()){
throw invalid_type{"value is not a scalar"};
}
try{
setting = yaml.as<Scalar>();
}
catch(...){
throw invalid_type{
"value " + yaml.as<std::string>()
+ " cannot be converted to a " + type_string(setting)
};
}
}
template<typename T>
void set_from_yaml(optional<T> & setting, YAML::Node const & yaml){
T tmp;
set_from_yaml(tmp, yaml);
setting = tmp;
}
template<typename T>
void set_from_yaml(std::vector<T> & setting, YAML::Node const & yaml){
assert(yaml);
// special case: treat a single value like a vector with one element
if(yaml.IsScalar()){
setting.resize(1);
return set_from_yaml(setting.front(), yaml);
}
if(yaml.IsSequence()){
setting.resize(yaml.size());
for(size_t i = 0; i < setting.size(); ++i){
set_from_yaml(setting[i], yaml[i]);
}
return;
}
throw invalid_type{""};
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) throw missing_option{""};
set_from_yaml(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
template<typename T>
void set_from_yaml_if_defined(T & setting, YAML::Node const & yaml){
return set_from_yaml(setting, yaml);
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) return;
set_from_yaml_if_defined(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
}
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": " + ex.what()
};
}
catch(missing_option const &){
throw missing_option{
"No entry for mandatory option " + join(": ", names...)
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml_if_defined(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": " + ex.what()
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
}
namespace YAML {
template<>
struct convert<RHEJ::OutputFile> {
static Node encode(RHEJ::OutputFile const & outfile);
static bool decode(Node const & node, RHEJ::OutputFile & out);
};
}
diff --git a/src/config.cc b/src/config.cc
index 697f6c1..eaeb27b 100644
--- a/src/config.cc
+++ b/src/config.cc
@@ -1,472 +1,472 @@
#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 == "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 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;
}
// 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{}};
}
if(scale_fun == "m12"){
return ret_type{new HardestJetsInvariantMass{}};
}
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();
+ 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"
+ 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");
+ 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;
}
}
diff --git a/src/main.cc b/src/main.cc
index ba3cc85..8766df8 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,146 +1,150 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/stream.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn < 3) {
std::cerr << "\n# Usage:\n./HEJ config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
- config.analysis_parameters
+ config.analysis_parameters()
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
- RHEJ::CombinedEventWriter writer{config.output, reader.heprup};
+ RHEJ::CombinedEventWriter writer{config.output(), reader.heprup};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
global_reweight = input_events/static_cast<double>(max_events);
std::cout << "processing " << max_events
<< " out of " << input_events << " events\n";
}
RHEJ::EventReweighter rhej{
reader.heprup,
- config.resummation_jets.def, config.resummation_jets.min_pt,
- config.treat,
- config.min_extparton_pt, config.max_ext_soft_pt_fraction,
- config.log_correction,
- config.Higgs_coupling
+ config.resummation_jets().def, config.resummation_jets().min_pt,
+ config.treat(),
+ config.min_extparton_pt(), config.max_ext_soft_pt_fraction(),
+ config.log_correction(),
+ config.Higgs_coupling()
};
- if(config.RanLux_init){
- RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
+ if(config.RanLux_init()){
+ RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init());
}
int nevent = 0;
std::array<int, RHEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
// Loop over the events in the inputfile
while(reader.readEvent()){
+ std::cout << "#event " << nevent << std::endl;
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
- if(nevent == max_events) break;
+ if(nevent == max_events) {
+ std::cout << "max #event reached " << nevent << " " << max_events << std::endl;
+ break;
+ }
nevent++;
if (nevent % 10000 == 0)
std::cout << "event number " << nevent << std::endl;
// calculate rHEJ weight
RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
- config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
+ config.fixed_order_jets().def, config.fixed_order_jets().min_pt,
};
auto resummed_events = rhej.reweight(
FO_event,
- config.trials, config.scale_gen
+ config.trials(), config.scale_gen()
);
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
for(auto const & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
if(analysis->pass_cuts(ev, FO_event)){
analysis->fill(ev, FO_event);
writer.write(ev);
}
}
} // main event loop
using namespace RHEJ::event_type;
std::cout<< "Events processed: " << nevent << '\n';
for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
<< ", failed to reconstruct " << nfailed_type[ev_type]
<< '\n';
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 1:46 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243518
Default Alt Text
(36 KB)

Event Timeline