diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
index 33ba36b..3a17408 100644
--- a/include/HEJ/Config.hh
+++ b/include/HEJ/Config.hh
@@ -1,229 +1,229 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <string>
#include <vector>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "HEJ/Constants.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/ScaleFunction.hh"
namespace HEJ{
//! Jet parameters
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Transverse Momentum */
};
//! Settings for scale variation
struct ScaleConfig{
//! Base scale choices
std::vector<ScaleFunction> base;
//! Factors for multiplicative scale variation
std::vector<double> factors;
//! Maximum ratio between renormalisation and factorisation scale
double max_ratio;
};
//! Settings for random number generator
struct RNGConfig {
//! Random number generator name
std::string name;
//! Optional initial seed
optional<std::string> seed;
};
//! Settings for partial unweighting
struct PartialUnweightConfig {
//! Number of trials for training
int trials;
//! Maximum distance in standard deviations from mean logarithmic weight
double max_dev;
};
/**! Possible treatments for fixed-order input events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Perform resummation */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
//! Container to store the treatments for various event types
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
//! Possible setting for the event weight
enum class WeightType{
weighted, //!< weighted events
unweighted_resum, //!< unweighted only resummation part
partially_unweighted //!< mixed weighted and unweighted
};
/**! Input parameters.
*
* This struct handles stores all configuration parameters
* needed in a HEJ 2 run.
*
* \internal To add a new option:
* 1. Add a member to the Config struct.
* 2. Inside "src/YAMLreader.cc":
* - Add the option name to the "supported" Node in
* get_supported_options.
* - Initialise the new Config member in to_Config.
* The functions set_from_yaml (for mandatory options) and
* set_from_yaml_if_defined (non-mandatory) may be helpful.
* 3. Add a new entry (with short description) to config.yaml
* 4. Update the user documentation in "doc/Sphinx/"
*/
struct Config {
//! %Parameters for scale variation
ScaleConfig scales;
//! Resummation jet properties
JetParameters resummation_jets;
//! Fixed-order jet properties
JetParameters fixed_order_jets;
//! Minimum transverse momentum for extremal partons
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
//! Number of resummation configurations to generate per fixed-order event
int trials;
//! Maximal number of events
optional<size_t> max_events;
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Event output files names and formats
std::vector<OutputFile> output;
//! Parameters for random number generation
RNGConfig rng;
//! Map to decide what to do for different event types
EventTreatMap treat;
//! %Parameters for custom analyses
- YAML::Node analysis_parameters;
+ std::vector<YAML::Node> analyses_parameters;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! Type of event weight e.g. (un)weighted
WeightType weight_type;
//! Settings for partial unweighting
HEJ::optional<PartialUnweightConfig> unweight_config;
};
//! Configuration options for the PhaseSpacePoint class
struct PhaseSpacePointConfig {
//! Properties of resummation jets
JetParameters jet_param;
//! Minimum transverse momentum for extremal partons
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
};
//! Configuration options for the MatrixElement class
struct MatrixElementConfig {
MatrixElementConfig() = default;
MatrixElementConfig(
bool log_correction,
HiggsCouplingSettings Higgs_coupling,
EWConstants ew_parameters,
double regulator_lambda = CLAMBDA
):
log_correction{log_correction},
Higgs_coupling{Higgs_coupling},
ew_parameters{ew_parameters},
regulator_lambda{regulator_lambda}
{}
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
};
//! Configuration options for the EventReweighter class
struct EventReweighterConfig {
//! Settings for phase space point generation
PhaseSpacePointConfig psp_config;
//! Settings for matrix element calculation
MatrixElementConfig ME_config;
//! Access properties of resummation jets
JetParameters & jet_param() {
return psp_config.jet_param;}
//! Access properties of resummation jets (const version)
JetParameters const & jet_param() const {
return psp_config.jet_param;}
//! Treatment of the various event types
EventTreatMap treat;
};
/**! Extract PhaseSpacePointConfig from Config
*
* \internal We do not provide a PhaseSpacePointConfig constructor from Config
* so that PhaseSpacePointConfig remains an aggregate.
* This faciliates writing client code (e.g. the HEJ fixed-order generator)
* that creates a PhaseSpacePointConfig *without* a Config object.
*
* @see to_MatrixElementConfig, to_EventReweighterConfig
*/
inline
PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) {
return {
conf.resummation_jets,
conf.min_extparton_pt,
conf.max_ext_soft_pt_fraction
};
}
/**! Extract MatrixElementConfig from Config
*
* @see to_PhaseSpacePointConfig, to_EventReweighterConfig
*/
inline
MatrixElementConfig to_MatrixElementConfig(Config const & conf) {
return {conf.log_correction, conf.Higgs_coupling,
conf.ew_parameters, conf.regulator_lambda};
}
/**! Extract EventReweighterConfig from Config
*
* @see to_PhaseSpacePointConfig, to_MatrixElementConfig
*/
inline
EventReweighterConfig to_EventReweighterConfig(Config const & conf) {
return {
to_PhaseSpacePointConfig(conf),
to_MatrixElementConfig(conf),
conf.treat
};
}
} // namespace HEJ
diff --git a/include/HEJ/get_analysis.hh b/include/HEJ/get_analysis.hh
index 2dae744..9f36c39 100644
--- a/include/HEJ/get_analysis.hh
+++ b/include/HEJ/get_analysis.hh
@@ -1,33 +1,37 @@
/** \file
* \brief Contains the get_analysis function
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
+#include <vector>
#include "HEJ/Analysis.hh"
namespace YAML {
class Node;
}
namespace HEJ {
//! Load an analysis
/**
* @param parameters Analysis parameters
* @param heprup General run informations
* @returns A pointer to an Analysis instance
*
* If parameters["plugin"] exists, an analysis (deriving from the
* \ref Analysis class) will be loaded from the library parameters["plugin"].
* Otherwise, if parameters["rivet"] exists, the corresponding RivetAnalysis
* will be loaded. If none of these parameters are specified, a pointer to
* the default EmptyAnalysis is returned.
*/
std::unique_ptr<Analysis> get_analysis(
YAML::Node const & parameters, LHEF::HEPRUP const & heprup);
+
+ std::vector<std::unique_ptr<Analysis>> get_analyses(
+ std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup);
}
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index c47dfce..2f79eed 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,555 +1,557 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/YAMLreader.hh"
#include <algorithm>
#include <iostream>
#include <limits>
#include <map>
#include <string>
#include <unordered_map>
#include <vector>
#include <dlfcn.h>
#include "HEJ/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/Constants.hh"
namespace HEJ{
class Event;
namespace{
//! Get YAML tree of supported options
/**
* The configuration file is checked against this tree of options
* in assert_all_options_known.
*/
YAML::Node const & get_supported_options(){
const static YAML::Node supported = [](){
YAML::Node supported;
static const auto opts = {
"trials", "min extparton pt", "max ext soft pt fraction",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "event output", "analysis", "vev",
"regulator parameter", "max events"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "algorithm", "R"}){
supported["resummation jets"][jet_opt] = "";
supported["fixed order jets"][jet_opt] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
for(auto && opt: {"FKL", "unordered", "extremal qqx", "central qqx", "non-resummable"}){
supported["event treatment"][opt] = "";
}
for(auto && particle_type: {"Higgs", "W", "Z"}){
for(auto && particle_opt: {"mass", "width"}){
supported["particle properties"][particle_type][particle_opt] = "";
}
}
for(auto && opt: {"type", "trials", "max deviation"}){
supported["unweight"][opt] = "";
}
return supported;
}();
return supported;
}
fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
using namespace fastjet;
static const std::map<std::string, fastjet::JetAlgorithm> known = {
{"kt", kt_algorithm},
{"cambridge", cambridge_algorithm},
{"antikt", antikt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"plugin", plugin_algorithm}
};
const auto res = known.find(algo);
if(res == known.end()){
throw std::invalid_argument("Unknown jet algorithm \"" + algo + "\"");
}
return res->second;
}
EventTreatment to_EventTreatment(std::string const & name){
static const std::map<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
{"discard", EventTreatment::discard}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment \"" + name + "\"");
}
return res->second;
}
WeightType to_weight_type(std::string const & setting){
if(setting == "weighted")
return WeightType::weighted;
if(setting =="resummation")
return WeightType::unweighted_resum;
if(setting =="partial")
return WeightType::partially_unweighted;
throw std::invalid_argument{"Unknown weight type \"" + setting + "\""};
}
} // namespace anonymous
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
setting = to_JetAlgorithm(yaml.as<std::string>());
}
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
setting = to_EventTreatment(yaml.as<std::string>());
}
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
setting = to_ParticleID(yaml.as<std::string>());
}
void set_from_yaml(WeightType & setting, YAML::Node const & yaml){
setting = to_weight_type(yaml.as<std::string>());
}
} // namespace detail
JetParameters get_jet_parameters(
YAML::Node const & node,
std::string const & entry
){
assert(node);
JetParameters result;
fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
double R;
set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
set_from_yaml(R, node, entry, "R");
result.def = fastjet::JetDefinition{jet_algo, R};
set_from_yaml(result.min_pt, node, entry, "min pt");
return result;
}
RNGConfig to_RNGConfig(
YAML::Node const & node,
std::string const & entry
){
assert(node);
RNGConfig result;
set_from_yaml(result.name, node, entry, "name");
set_from_yaml_if_defined(result.seed, node, entry, "seed");
return result;
}
ParticleProperties get_particle_properties(
YAML::Node const & node, std::string const & entry,
std::string const & boson
){
ParticleProperties result;
set_from_yaml(result.mass, node, entry, boson, "mass");
set_from_yaml(result.width, node, entry, boson, "width");
return result;
}
EWConstants get_ew_parameters(YAML::Node const & node){
EWConstants result;
double vev;
set_from_yaml(vev, node, "vev");
result.set_vevWZH(vev,
get_particle_properties(node, "particle properties", "W"),
get_particle_properties(node, "particle properties", "Z"),
get_particle_properties(node, "particle properties", "Higgs")
);
return result;
}
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node,
std::string const & entry
){
assert(node);
static constexpr double mt_max = 2e4;
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(node[entry]){
throw std::invalid_argument{
"Higgs coupling settings require building HEJ 2 "
"with QCDloop support"
};
}
#endif
HiggsCouplingSettings settings;
set_from_yaml_if_defined(settings.mt, node, entry, "mt");
set_from_yaml_if_defined(settings.mb, node, entry, "mb");
set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
if(settings.use_impact_factors){
if(settings.mt != std::numeric_limits<double>::infinity()){
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
else{
// huge values of the top mass are numerically unstable
settings.mt = std::min(settings.mt, mt_max);
}
return settings;
}
FileFormat to_FileFormat(std::string const & name){
static const std::map<std::string, FileFormat> known = {
{"Les Houches", FileFormat::Les_Houches},
{"HepMC", FileFormat::HepMC},
{"HepMC2", FileFormat::HepMC2},
{"HepMC3", FileFormat::HepMC3},
{"HDF5", FileFormat::HDF5}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown file format \"" + name + "\"");
}
return res->second;
}
std::string extract_suffix(std::string const & filename){
size_t separator = filename.rfind('.');
if(separator == 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;
if(suffix == "hepmc3") return FileFormat::HepMC3;
if(suffix == "hepmc2") return FileFormat::HepMC2;
if(suffix == "hdf5") return FileFormat::HDF5;
throw std::invalid_argument{
"Can't determine format for output file \"" + filename + "\""
};
}
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
){
if(!conf.IsMap()) return;
if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
for(auto const & entry: conf){
const auto name = entry.first.as<std::string>();
if(! supported[name]) throw unknown_option{name};
/* check sub-options, e.g. 'resummation jets: min pt'
* we don't check analysis sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(name != "analysis" && name != "import scales"){
try{
assert_all_options_known(conf[name], supported[name]);
}
catch(unknown_option const & ex){
throw unknown_option{name + ": " + ex.what()};
}
catch(invalid_type const & ex){
throw invalid_type{name + ": " + ex.what()};
}
}
}
}
} // namespace HEJ
namespace YAML {
Node convert<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
}
bool convert<HEJ::OutputFile>::decode(Node const & node, HEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = HEJ::to_FileFormat(it->first.as<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
out.format = HEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
} // namespace YAML
namespace HEJ{
namespace detail{
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
setting = yaml.as<OutputFile>();
}
}
namespace{
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
if(!yaml["fixed order jets"]) return;
set_from_yaml_if_defined(
fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
);
fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
double R = fixed_order_jets.def.R();
set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
fixed_order_jets.def = fastjet::JetDefinition{algo, R};
}
// like std::stod, but throw if not the whole string can be converted
double to_double(std::string const & str){
std::size_t pos;
const double result = std::stod(str, &pos);
if(pos < str.size()){
throw std::invalid_argument(str + " is not a valid double value");
}
return result;
}
using EventScale = double (*)(Event const &);
void import_scale_functions(
std::string const & file,
std::vector<std::string> const & scale_names,
std::unordered_map<std::string, EventScale> & known
) {
auto handle = dlopen(file.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
for(auto const & scale: scale_names) {
void * sym = dlsym(handle, scale.c_str());
error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
known.emplace(scale, reinterpret_cast<EventScale>(sym));
}
}
auto get_scale_map(
YAML::Node const & yaml
) {
std::unordered_map<std::string, EventScale> scale_map;
scale_map.emplace("H_T", H_T);
scale_map.emplace("max jet pperp", max_jet_pt);
scale_map.emplace("jet invariant mass", jet_invariant_mass);
scale_map.emplace("m_j1j2", m_j1j2);
if(yaml["import scales"]) {
if(! yaml["import scales"].IsMap()) {
throw invalid_type{"Entry 'import scales' is not a map"};
}
for(auto const & import: yaml["import scales"]) {
const auto file = import.first.as<std::string>();
const auto scale_names =
import.second.IsSequence()
?import.second.as<std::vector<std::string>>()
:std::vector<std::string>{import.second.as<std::string>()};
import_scale_functions(file, scale_names, scale_map);
}
}
return scale_map;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be H_T,
* H_T/2 is then composite (take H_T and then divide by 2)
*/
ScaleFunction parse_simple_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
) {
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const auto it = known.find(scale_fun);
if(it != end(known)) return {it->first, it->second};
try{
const double scale = to_double(scale_fun);
return {scale_fun, FixedScale{scale}};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: \"" + scale_fun + "\""};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
ScaleFunction parse_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
// parse from right to left => a/b/c gives (a/b)/c
const size_t delim = scale_fun.find_last_of("*/");
if(delim == scale_fun.npos){
return parse_simple_ScaleFunction(scale_fun, known);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
if(scale_fun[delim] == '/'){
return parse_ScaleFunction(first, known)
/ parse_ScaleFunction(second, known);
}
else{
assert(scale_fun[delim] == '*');
return parse_ScaleFunction(first, known)
* parse_ScaleFunction(second, known);
}
}
EventTreatMap get_event_treatment(
YAML::Node const & node, std::string const & entry
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::discard},
{unob, EventTreatment::discard},
{unof, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{non_resummable, EventTreatment::discard}
};
set_from_yaml(treat.at(FKL), node, entry, "FKL");
set_from_yaml(treat.at(unob), node, entry, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(qqxexb), node, entry, "extremal qqx");
treat.at(qqxexf) = treat.at(qqxexb);
set_from_yaml(treat.at(qqxmid), node, entry, "central qqx");
set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable");
if(treat[non_resummable] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-resummable events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
try{
assert_all_options_known(yaml, get_supported_options());
}
catch(unknown_option const & ex){
throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
set_from_yaml_if_defined(config.min_extparton_pt, yaml, "min extparton pt");
if(config.min_extparton_pt!=0)
std::cerr << "WARNING: \"min extparton pt\" is deprecated."
<< " Please use \"max ext soft pt fraction\" instead.\n";
set_from_yaml(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
// Sets the standard value, then changes this if defined
config.regulator_lambda=CLAMBDA;
set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
set_from_yaml_if_defined(config.max_events, yaml, "max events");
set_from_yaml(config.trials, yaml, "trials");
config.weight_type = WeightType::weighted;
set_from_yaml_if_defined(config.weight_type, yaml, "unweight", "type");
if(config.weight_type == WeightType::partially_unweighted) {
config.unweight_config = PartialUnweightConfig{};
set_from_yaml(
config.unweight_config->trials, yaml,
"unweight", "trials"
);
set_from_yaml(
config.unweight_config->max_dev, yaml,
"unweight", "max deviation"
);
}
else if(yaml["unweight"]) {
for(auto && opt: {"trials", "max deviation"}) {
if(yaml["unweight"][opt]) {
throw std::invalid_argument{
"'unweight: " + std::string{opt} + "' "
"is only supported if 'unweight: type' is set to 'partial'"
};
}
}
}
set_from_yaml(config.log_correction, yaml, "log correction");
config.treat = get_event_treatment(yaml, "event treatment");
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = to_RNGConfig(yaml, "random generator");
- set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
+ YAML::Node analysis;
+ set_from_yaml_if_defined(analysis, yaml, "analysis");
+ config.analyses_parameters.push_back(analysis);
config.scales = to_ScaleConfig(yaml);
config.ew_parameters = get_ew_parameters(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace anonymous
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
auto scale_funs = get_scale_map(yaml);
std::vector<std::string> scales;
set_from_yaml(scales, yaml, "scales");
config.base.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.base),
[scale_funs](auto const & entry){
return parse_ScaleFunction(entry, scale_funs);
}
);
set_from_yaml_if_defined(config.factors, yaml, "scale factors");
config.max_ratio = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
return config;
}
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
} // namespace HEJ
diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index 37067c1..8d76fa5 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,376 +1,382 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <array>
#include <chrono>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include "yaml-cpp/yaml.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/CombinedEventWriter.hh"
#include "HEJ/Config.hh"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/BufferedEventReader.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/get_analysis.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/optional.hh"
#include "HEJ/ProgressBar.hh"
#include "HEJ/stream.hh"
#include "HEJ/Unweighter.hh"
#include "HEJ/Version.hh"
#include "HEJ/YAMLreader.hh"
HEJ::Config load_config(char const * filename){
try{
return HEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
-std::unique_ptr<HEJ::Analysis> get_analysis(
- YAML::Node const & parameters, LHEF::HEPRUP const & heprup
+std::vector<std::unique_ptr<HEJ::Analysis>> get_analyses(
+ std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
){
try{
- return HEJ::get_analysis(parameters, heprup);
+ return HEJ::get_analyses(parameters, heprup);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
// unique_ptr is a workaround:
// HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
std::unique_ptr<HEJ::ProgressBar<double>> make_progress_bar(
std::vector<double> const & xs
) {
if(xs.empty()) return {};
const double Born_xs = std::accumulate(begin(xs), end(xs), 0.);
return std::make_unique<HEJ::ProgressBar<double>>(std::cout, Born_xs);
}
std::string time_to_string(const time_t time){
char s[30];
struct tm * p = localtime(&time);
strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
return s;
}
HEJ::Event to_event(
LHEF::HEPEUP const & hepeup,
HEJ::JetParameters const & fixed_order_jets
) {
HEJ::Event::EventData event_data{hepeup};
event_data.reconstruct_intermediate();
return HEJ::Event{
std::move(event_data).cluster(
fixed_order_jets.def, fixed_order_jets.min_pt
)
};
}
void unweight(
HEJ::Unweighter & unweighter,
HEJ::WeightType weight_type,
std::vector<HEJ::Event> & events,
HEJ::RNG & ran
) {
if(weight_type == HEJ::WeightType::unweighted_resum){
unweighter.set_cut_to_maxwt(events);
}
events.erase(
unweighter.unweight(begin(events), end(events), ran),
end(events)
);
}
// peek up to nevents events from reader
std::vector<LHEF::HEPEUP> peek_events(
HEJ::BufferedEventReader & reader,
const int nevents
) {
std::vector<LHEF::HEPEUP> events;
while(
static_cast<int>(events.size()) < nevents
&& reader.read_event()
) {
events.emplace_back(reader.hepeup());
}
// put everything back into the reader
for(auto it = rbegin(events); it != rend(events); ++it) {
reader.emplace(*it);
}
return events;
}
void append_resummed_events(
std::vector<HEJ::Event> & resummation_events,
HEJ::EventReweighter & reweighter,
LHEF::HEPEUP const & hepeup,
const int trials,
HEJ::JetParameters const & fixed_order_jets
) {
const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets);
if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
return;
}
const auto resummed = reweighter.reweight(FO_event, trials);
resummation_events.insert(
end(resummation_events),
begin(resummed), end(resummed)
);
}
void train(
HEJ::Unweighter & unweighter,
HEJ::BufferedEventReader & reader,
HEJ::EventReweighter & reweighter,
const int total_trials,
const double max_dev,
double reweight_factor,
HEJ::JetParameters const & fixed_order_jets
) {
std::cout << "Reading up to " << total_trials << " training events...\n";
auto FO_events = peek_events(reader, total_trials);
if(FO_events.empty()) {
throw std::runtime_error{
"No events generated to calibrate the unweighting weight!"
"Please increase the number \"trials\" or deactivate the unweighting."
};
}
const int trials = total_trials/FO_events.size();
// adjust reweight factor so that the overall normalisation
// is the same as in the full run
reweight_factor *= trials;
for(auto & hepeup: FO_events) {
hepeup.XWGTUP *= reweight_factor;
}
std::cout << "Training unweighter with "
<< trials << '*' << FO_events.size() << " events\n";
auto progress = HEJ::ProgressBar<int>{
std::cout, static_cast<int>(FO_events.size())
};
std::vector<HEJ::Event> resummation_events;
for(auto const & hepeup: FO_events) {
append_resummed_events(
resummation_events,
reweighter, hepeup, trials, fixed_order_jets
);
++progress;
}
unweighter.set_cut_to_peakwt(resummation_events, max_dev);
std::cout << "\nUnweighting events with weight up to "
<< unweighter.get_cut() << '\n';
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn != 3) {
std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
{
std::cout << "Starting " << HEJ::Version::package_name_full()
<< ", revision " << HEJ::Version::revision() << " ("
<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
fastjet::ClusterSequence::print_banner();
// read configuration
const HEJ::Config config = load_config(argv[1]);
auto reader = HEJ::make_reader(argv[2]);
assert(reader);
auto heprup{ reader->heprup() };
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = HEJ::Version::package_name();
heprup.generators.back().version = HEJ::Version::String();
- std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
- config.analysis_parameters, heprup
- );
- assert(analysis != nullptr);
+ auto analyses = get_analyses( config.analyses_parameters, heprup );
+ assert(analyses.empty() || analyses.front() != nullptr);
HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
double global_reweight = 1.;
const auto & max_events = config.max_events;
// if we need the event number:
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){
// try to read from LHE head
auto input_events{reader->number_events()};
if(!input_events) {
// else count manually
auto t_reader = HEJ::make_reader(argv[2]);
input_events = 0;
while(t_reader->read_event()) ++(*input_events);
}
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){
// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
<< "assuming \"cross section = average weight\".\n"
<< "converting to \"cross section = sum of weights\" ";
global_reweight /= *input_events;
}
if(max_events && (*input_events > *max_events)){
// maximal number of events given
global_reweight *= *input_events/static_cast<double>(*max_events);
std::cout << "Processing " << *max_events
<< " out of " << *input_events << " events\n";
}
}
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
std::shared_ptr<HEJ::RNG> ran{
HEJ::make_RNG(config.rng.name, config.rng.seed)};
assert(ran != nullptr);
HEJ::EventReweighter hej{
reader->heprup(),
std::move(scale_gen),
to_EventReweighterConfig(config),
ran
};
HEJ::optional<HEJ::Unweighter> unweighter{};
if(config.weight_type != HEJ::WeightType::weighted) {
unweighter = HEJ::Unweighter{};
}
if(config.weight_type == HEJ::WeightType::partially_unweighted) {
HEJ::BufferedEventReader buffered_reader{std::move(reader)};
assert(config.unweight_config);
train(
*unweighter,
buffered_reader,
hej,
config.unweight_config->trials,
config.unweight_config->max_dev,
global_reweight/config.trials,
config.fixed_order_jets
);
reader = std::make_unique<HEJ::BufferedEventReader>(
std::move(buffered_reader)
);
}
// status infos & eye candy
size_t nevent = 0;
std::array<int, HEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
auto progress = make_progress_bar(reader->heprup().XSECUP);
HEJ::CrossSectionAccumulator xs;
std::map<HEJ::StatusCode, int> status_counter;
size_t total_trials = 0;
size_t total_resum = 0;
// Loop over the events in the input file
while(reader->read_event() && (!max_events || nevent < *max_events) ){
++nevent;
// reweight events so that the total cross section is conserved
auto hepeup = reader->hepeup();
hepeup.XWGTUP *= global_reweight;
const auto FO_event = to_event(hepeup, config.fixed_order_jets);
if(FO_event.central().weight == 0) {
static const bool warned_once = [argv,nevent](){
std::cerr
<< "WARNING: event number " << nevent
<< " in " << argv[2] << " has zero weight. "
"Ignoring this and all further events with vanishing weight.\n";
return true;
}();
(void) warned_once; // shut up compiler warnings
continue;
}
auto resummed_events{ hej.reweight(FO_event, config.trials) };
// some bookkeeping
for(auto const & s: hej.status())
++status_counter[s];
total_trials+=hej.status().size();
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
if(unweighter) {
unweight(*unweighter, config.weight_type, resummed_events, *ran);
}
// analysis
for(auto & 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);
+ bool passed = false;
+ for(auto const & analysis: analyses){
+ if(analysis->pass_cuts(ev, FO_event)){
+ passed = true;
+ analysis->fill(ev, FO_event);
+ };
+ }
+ if(passed){
writer.write(ev);
} else {
ev.parameters()*=0; // do not use discarded events afterwards
}
}
xs.fill_correlated(resummed_events);
total_resum += resummed_events.size();
if(progress) progress->increment(FO_event.central().weight);
} // main event loop
std::cout << '\n';
- analysis->finalise();
+ for(auto const & analysis: analyses){
+ analysis->finalise();
+ }
using namespace HEJ::event_type;
std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
std::cout << '\t' << name(EventType::first_type) << ": "
<< nevent_type[EventType::first_type]
<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
<< '\n';
for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
<< nevent_type[i]
<< ", failed to reconstruct " << nfailed_type[i]
<< '\n';
}
std::cout << '\n' << xs << '\n';
std::cout << "Generation statistic: "
<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
<< " trials successful.\n";
for(auto && entry: status_counter){
const double fraction = static_cast<double>(entry.second)/total_trials;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":")
<< " [";
for(int i = 0; i < percent/2; ++i) std::cout << '#';
for(int i = percent/2; i < 50; ++i) std::cout << ' ';
std::cout << "] " <<std::setw(2)<<std::right<< percent << "%\n";
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nFinished " << HEJ::Version::package_name() << " at "
<< time_to_string(clock::to_time_t(clock::now()))
<< "\n=> Runtime: " << run_time.count() << " sec ("
<< nevent/run_time.count() << " Events/sec).\n";
return EXIT_SUCCESS;
}
diff --git a/src/get_analysis.cc b/src/get_analysis.cc
index 2963fb7..56f56c5 100644
--- a/src/get_analysis.cc
+++ b/src/get_analysis.cc
@@ -1,41 +1,51 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/get_analysis.hh"
#include <dlfcn.h>
#include <string>
#include "yaml-cpp/yaml.h"
#include "HEJ/EmptyAnalysis.hh"
#include "HEJ/RivetAnalysis.hh"
namespace HEJ{
std::unique_ptr<Analysis> get_analysis(
YAML::Node const & parameters, LHEF::HEPRUP const & heprup
){
if(!parameters["plugin"]){
if(parameters["rivet"])
return RivetAnalysis::create(parameters, heprup);
return EmptyAnalysis::create(parameters, heprup);
}
using AnalysisMaker = std::unique_ptr<Analysis> (*)(
YAML::Node const &, LHEF::HEPRUP const &);
const auto plugin_name = parameters["plugin"].as<std::string>();
auto handle = dlopen(plugin_name.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error(error);
void * sym = dlsym(handle, "make_analysis");
error = dlerror();
if(error != nullptr) throw std::runtime_error(error);
auto make_analysis = reinterpret_cast<AnalysisMaker>(sym);
return make_analysis(parameters, heprup);
}
+
+ std::vector<std::unique_ptr<Analysis>> get_analyses(
+ std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
+ ){
+ std::vector<std::unique_ptr<Analysis>> anas;
+ for(auto const & param: parameters){
+ anas.emplace_back(get_analysis(param, heprup));
+ }
+ return anas;
+ }
}