Page MenuHomeHEPForge
Files F8725689

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
index bab25fb..e61f466 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,114 +1,260 @@
/** \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_classes.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 EventClass keys and mapped values EventTreatment
*
* If called upon with EventTreatMap.at(EventClass) it will return the treatment which has
* been specified for that EventClass.
*/
using EventTreatMap = std::map<event_class::EventClass, EventTreatment>;
/**! \struct Config config.hh "include/RHEJ/config.hh"
* \brief Config Struct for user input parameters.
*
* The struct which handles the input card given by the user.
*/
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 */
//! Settings for Effective Higgs-Gluon Coupling
HiggsCouplingSettings Higgs_coupling;
};
Config load_config(std::string const & config_file);
+
+ template<typename T, typename... YamlNames>
+ void set_from_yaml(
+ T & setting,
+ YAML::Node const & yaml, YamlNames const & ... names
+ );
+
+ template<typename T, typename... YamlNames>
+ void set_from_yaml_if_defined(
+ T & setting,
+ YAML::Node const & yaml, YamlNames const & ... names
+ );
+
+ namespace detail{
+ 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{""};
+ }
+
+ void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml);
+ void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml);
+
+ inline
+ void set_from_yaml(YAML::Node & setting, YAML::Node const & yaml){
+ setting = yaml;
+ }
+
+ 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...) + ": value " + ex.what()
+ + " cannot be converted to a " + type_string(setting)
+ };
+ }
+ 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()
+ };
+ }
+ }
+
}
diff --git a/include/RHEJ/utility.hh b/include/RHEJ/utility.hh
index 542e710..9caab71 100644
--- a/include/RHEJ/utility.hh
+++ b/include/RHEJ/utility.hh
@@ -1,118 +1,141 @@
/**
* \file utility.hh
* Author: Tuomas Hapola
* \brief Contains the struct Sparticle which contains particle information and parameters.
*
* Also contains some functions which are useful for comparison on particles parameters.
*/
#pragma once
#include <algorithm>
+#include <boost/core/demangle.hpp>
+
// FastJet Includes
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "RHEJ/PDG_codes.hh"
namespace RHEJ{
/** \struct Sparticle utility.hh "include/RHEJ/utility.hh"
* \brief The struct which contains all the parameters of a particle.
*
* Constituted by PiD and PseudoJet Objects.
*/
struct Sparticle {
ParticleID type;
fastjet::PseudoJet p;
//! get rapidity function
double rapidity() const{
return p.rapidity();
}
//! get perp function
double perp() const{
return p.perp();
}
//! get Px function
double px() const{
return p.px();
}
//! get Py function
double py() const{
return p.py();
}
//! Get Pz function
double pz() const{
return p.pz();
}
//! Get P0 function
double E() const{
return p.E();
}
};
/** \struct rapidity_less utility.hh "include/RHEJ/utility.hh
- * \brief A struct which allows quick comparison of two different jets rapidities.
+ * \brief A struct which allows quick comparison of two different jets rapidities.
*/
struct rapidity_less{
template<class FourVector>
bool operator()(FourVector const & p1, FourVector const & p2){
return p1.rapidity() < p2.rapidity();
}
};
-
+
/** \struct pz_less utility.hh "include/RHEJ/utility.hh
- * \brief A struct which allows quick comparison of Pz between two jets.
+ * \brief A struct which allows quick comparison of Pz between two jets.
*/
struct pz_less{
template<class FourVector>
bool operator()(FourVector const & p1, FourVector const & p2){
return p1.pz() < p2.pz();
}
};
-
+
inline
std::vector<fastjet::PseudoJet> to_PseudoJet(
std::vector<Sparticle> const & v
){
std::vector<fastjet::PseudoJet> result;
for(auto && sp: v) result.emplace_back(sp.p);
return result;
}
inline
bool is_parton(Sparticle const & p){
return is_parton(p.type);
}
inline bool is_AWZH_boson(Sparticle const & particle){
return is_AWZH_boson(particle.type);
}
inline
std::vector<Sparticle> filter_partons(
std::vector<Sparticle> const & v
){
std::vector<Sparticle> result;
result.reserve(v.size());
std::copy_if(
begin(v), end(v), std::back_inserter(result),
[](Sparticle const & p){ return is_parton(p); }
);
return result;
}
template<class T, class... Args>
std::unique_ptr<T> make_unique(Args&&... a){
return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
}
template<typename T, typename... U>
constexpr
std::array<T, 1 + sizeof...(U)> make_array(T t, U&&... rest){
return {{t, std::forward<U>(rest)...}};
}
+ inline
+ std::string join(
+ std::string const & /* delim */, std::string const & str
+ ){
+ return str;
+ }
+
+ template<typename... Strings>
+ std::string join(
+ std::string const & delim,
+ std::string const & first, std::string const & second,
+ Strings&&... rest
+ ){
+ return join(delim, first + delim + second, std::forward<Strings>(rest)...);
+ }
+
+ template<typename T>
+ std::string type_string(T&&){
+ return boost::core::demangle(typeid(T).name());
+ }
+
}
diff --git a/src/config.cc b/src/config.cc
index e7e9f5e..a5851ce 100644
--- a/src/config.cc
+++ b/src/config.cc
@@ -1,606 +1,412 @@
#include "RHEJ/config.hh"
#include "RHEJ/exceptions.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
namespace RHEJ{
+
+ namespace{
+ fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
+ using namespace fastjet;
+ static const std::map<std::string, fastjet::JetAlgorithm> known = {
+ {"kt", kt_algorithm},
+ {"cambridge", cambridge_algorithm},
+ {"antikt", antikt_algorithm},
+ {"genkt", genkt_algorithm},
+ {"cambridge for passive", cambridge_for_passive_algorithm},
+ {"genkt for passive", genkt_for_passive_algorithm},
+ {"ee kt", ee_kt_algorithm},
+ {"ee genkt", ee_genkt_algorithm},
+ {"plugin", plugin_algorithm}
+ };
+ const auto res = known.find(algo);
+ if(res == known.end()){
+ throw std::invalid_argument("Unknown jet algorithm " + algo);
+ }
+ return res->second;
+ }
+
+ EventTreatment to_EventTreatment(std::string const & name){
+ static const std::map<std::string, EventTreatment> known = {
+ {"reweight", EventTreatment::reweight},
+ {"keep", EventTreatment::keep},
+ {"discard", EventTreatment::discard}
+ };
+ const auto res = known.find(name);
+ if(res == known.end()){
+ throw std::invalid_argument("Unknown event treatment " + name);
+ }
+ return res->second;
+ }
+
+ }
+
+ namespace 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>());
+ }
+ }
+
+ namespace{
+ 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;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
}
namespace YAML {
std::string to_string(NodeType::value type){
switch (type) {
case NodeType::Null: return "Null";
case NodeType::Scalar: return "Scalar";
case NodeType::Sequence: return "Sequence";
case NodeType::Map: return "Map";
case NodeType::Undefined: return "Undefined";
default:
throw std::logic_error{"Missing case in switch statement"};
}
}
template<>
struct convert<RHEJ::OutputFile> {
static Node encode(RHEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
}
static bool 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{
static const std::set<std::string> known = {
{"trials"},
{"min extparton pt"}, {"max ext soft pt fraction"},
{"resummation jets"}, {"fixed order jets"},
{"FKL"}, {"non-FKL"}, {"unordered"},
{"log correction"},
{"event output"},
{"analysis"},
{"unweight"}, {"RanLux init"},
{"scales"}, {"scale factors"}, {"max scale ratio"},
{"Higgs coupling"}
};
static const std::set<std::string> known_jet = {
{"min pt"}, {"algorithm"}, {"R"}
};
- template<typename T>
- std::string type_string();
-
- template<>
- std::string type_string<std::vector<double>>(){
- return "array of doubles";
- }
-
- template<>
- std::string type_string<std::vector<std::string>>(){
- return "array of strings";
- }
-
- template<>
- std::string type_string<int>(){
- return "integer";
- }
-
- template<>
- std::string type_string<double>(){
- return "double";
- }
-
- template<>
- std::string type_string<bool>(){
- return "bool";
- }
-
- template<>
- std::string type_string<std::string>(){
- return "string";
- }
-
- template<>
- std::string type_string<std::vector<OutputFile>>(){
- return "array of output files";
- }
-
- 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;
- }
-
- // helper struct to allow partial template specialisation
- template<typename T>
- struct existing_as_helper{
-
- static T as(YAML::Node const & config, std::string const & entry){
- assert(config[entry]);
- try{
- return config[entry].as<T>();
- }
- catch(std::exception const &){
- throw invalid_type{
- "Entry " + entry + " is not of type " + type_string<T>()
- };
- }
- }
-
- };
-
- template<>
- struct existing_as_helper<fastjet::JetAlgorithm>{
-
- static fastjet::JetAlgorithm as(
- YAML::Node const & config, std::string const & entry
- ){
- assert(config[entry]);
- const std::string algo_name =
- existing_as_helper<std::string>::as(config, entry);
- return to_JetAlgorithm(algo_name);
- }
-
- };
-
- template<>
- struct existing_as_helper<EventTreatment>{
-
- static EventTreatment as(
- YAML::Node const & config, std::string const & entry
- ){
- assert(config[entry]);
- const std::string name =
- existing_as_helper<std::string>::as(config, entry);
- return to_EventTreatment(name);
- }
-
- };
-
- template<>
- struct existing_as_helper<OutputFile>{
-
- static OutputFile as(
- YAML::Node const & config, std::string const & entry
- ){
- assert(config[entry]);
- YAML::convert<RHEJ::OutputFile> converter{};
- OutputFile out;
- if(converter.decode(config[entry], out)) return out;
- throw std::invalid_argument{
- "Bad output file setting: " + config[entry].as<std::string>()
- };
- }
- };
-
- template<typename T>
- struct existing_as_helper<std::vector<T>>{
-
- static std::vector<T> as(
- YAML::Node const & config, std::string const & entry
- ){
- assert(config[entry]);
- // special case: treat a single value like a vector with one element
- if(config[entry].IsScalar()){
- return {existing_as_helper<T>::as(config, entry)};
- }
- try{
- return config[entry].as<std::vector<T>>();
- }
- catch(std::exception const & p){
- std::cerr << "Error: " << p.what() << '\n';
- throw invalid_type{
- "Entry " + entry + " is not of type " + type_string<std::vector<T>>()
- };
- }
- }
-
- };
-
- template<typename T>
- T existing_as(YAML::Node const & config, std::string const & entry){
- return existing_as_helper<T>::as(config, entry);
- }
-
- template<typename T>
- T as(YAML::Node const & config, std::string const & entry);
-
- template<>
- JetParameters existing_as<JetParameters>(
- YAML::Node const & config, std::string const & entry
- ){
- assert(config[entry]);
- YAML::Node const & jet_config = config[entry];
- JetParameters result;
- for(auto const & opt: jet_config){
- auto const & opt_name = opt.first.as<std::string>();
- if(! known_jet.count(opt_name)){
- throw unknown_option{"Unknown option " + opt_name};
- }
- }
- try{
- result.def = fastjet::JetDefinition{
- as<fastjet::JetAlgorithm>(jet_config, "algorithm"),
- as<double>(jet_config, "R")
- };
- result.min_pt = as<double>(jet_config, "min pt");
- return result;
- }
- catch(missing_option const & exc){
- throw missing_option{"In entry " + entry + ":\n" + exc.what()};
- }
- catch(invalid_type const & exc){
- throw invalid_type{"In entry " + entry + ":\n" + exc.what()};
- }
- }
-
- template<typename T>
- T as(YAML::Node const & config, std::string const & entry){
- if(!config[entry]){
- throw missing_option{"No entry for option " + entry};
- }
- return existing_as<T>(config, entry);
- }
-
- template<typename T>
- optional<T> as_optional(YAML::Node const & config, std::string const & entry){
- if(!config[entry]) return {};
- return {existing_as<T>(config, entry)};
- }
-
- template<>
- HiggsCouplingSettings as<HiggsCouplingSettings>(
- YAML::Node const & config, std::string const & entry
- ){
- static constexpr double mt_max = 2e4;
- HiggsCouplingSettings settings;
- if(!config[entry]) return settings;
-#ifndef RHEJ_BUILD_WITH_QCDLOOP
- throw std::invalid_argument{
- "Higgs coupling settings require building Reversed HEJ "
- "with QCDloop support"
- };
-#endif
- YAML::Node const & node = config[entry];
- if(!node.IsMap()){
- throw std::invalid_argument{
- "entry '" + entry + "' is a " + to_string(node.Type())
- + ", not a Map of properties"
- };
- }
- for(auto const & property: node){
- const auto key = property.first.as<std::string>();
- auto const & value = property.second;
- if(key == "mt") settings.mt = value.as<double>();
- else if(key == "mb") settings.mb = value.as<double>();
- else if(key == "include bottom"){
- settings.include_bottom = value.as<bool>();
- }
- else if(key == "use impact factors"){
- settings.use_impact_factors = value.as<bool>();
- }
- else{
- throw unknown_option{"Unknown option '" + key + "'"};
- }
- }
- 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;
- }
-
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
- YAML::Node const & yaml_jets = yaml["fixed order jets"];
- assert(yaml_jets);
- for(auto const & opt: yaml_jets){
- auto const & opt_name = opt.first.as<std::string>();
- if(! known_jet.count(opt_name)){
- throw unknown_option{"Unknown option '" + opt_name + "'"};
- }
- }
- try{
- auto algo = as_optional<fastjet::JetAlgorithm>(yaml_jets, "algorithm");
- if(algo){
- fixed_order_jets.def = fastjet::JetDefinition{
- *algo, fixed_order_jets.def.R()
- };
- }
- auto R = as_optional<double>(yaml_jets, "R");
- if(R){
- fixed_order_jets.def = fastjet::JetDefinition{
- fixed_order_jets.def.jet_algorithm(), *R
- };
- }
- auto min_pt = as_optional<double>(yaml_jets, "min pt");
- if(min_pt) fixed_order_jets.min_pt = *min_pt;
- }
- catch(missing_option const & exc){
- throw missing_option{
- std::string{"In entry fixed order jets:\n"} + exc.what()
- };
- }
- catch(invalid_type const & exc){
- throw invalid_type{
- std::string{"In entry fixed order jets:\n"} + exc.what()
- };
- }
+ 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)}
};
}
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
- const auto scales = as<std::vector<std::string>>(yaml, "scales");
+ 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
);
- auto scale_factors = as_optional<std::vector<double>>(
- yaml, "scale factors"
- );
- if(scale_factors) config.scale_factors = std::move(*scale_factors);
- const auto max_scale_ratio = as_optional<double>(yaml, "max scale ratio");
- static_assert(
- std::numeric_limits<double>::has_infinity, "infinity not supported"
- );
- config.max_scale_ratio = max_scale_ratio?
- *max_scale_ratio:
- std::numeric_limits<double>::infinity()
- ;
+ 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;
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_class;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
- {bad_final_state, EventTreatment::discard}
+ {bad_final_state, EventTreatment::discard},
+ {FKL, EventTreatment::reweight},
+ {unob, EventTreatment::keep},
+ {unof, EventTreatment::keep},
+ {nonFKL, EventTreatment::keep}
};
- treat.emplace(FKL, as<EventTreatment>(yaml, "FKL"));
- const auto unordered_treatment = as<EventTreatment>(yaml, "unordered");
- treat.emplace(unordered_forward, unordered_treatment);
- treat.emplace(unordered_backward, unordered_treatment);
- treat.emplace(nonFKL, as<EventTreatment>(yaml, "non-FKL"));
+ 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){
for(auto const & opt: yaml){
auto const & opt_name = opt.first.as<std::string>();
if(! known.count(opt_name)){
throw unknown_option{"Unknown option '" + opt_name + "'"};
}
}
Config config;
- config.resummation_jets = as<JetParameters>(yaml, "resummation jets");
+ config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
- if(yaml["fixed order jets"]){
- update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
- }
- config.min_extparton_pt = as<double>(yaml, "min extparton pt");
- config.max_ext_soft_pt_fraction = yaml["max ext soft pt fraction"]?
- as<double>(yaml, "max ext soft pt fraction"):
- std::numeric_limits<double>::infinity()
- ;
- config.trials = as<int>(yaml, "trials");
- config.log_correction = as<bool>(yaml, "log correction");
- config.unweight = as<bool>(yaml, "unweight");
+ 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);
- if(yaml["event output"]){
- config.output = as<std::vector<OutputFile>>(yaml, "event output");
- }
- config.RanLux_init = as_optional<std::string>(yaml, "RanLux init");
- config.analysis_parameters = yaml["analysis"]?yaml["analysis"]:YAML::Node{};
+ 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 = as<HiggsCouplingSettings>(yaml, "Higgs coupling");
+ config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // anonymous namespace
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, Jan 21, 2:18 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243641
Default Alt Text
(37 KB)

Event Timeline