Page MenuHomeHEPForge
Files F8308610

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/ScaleFunction.hh b/include/RHEJ/ScaleFunction.hh
index 0ef96d9..6719e35 100644
--- a/include/RHEJ/ScaleFunction.hh
+++ b/include/RHEJ/ScaleFunction.hh
@@ -1,125 +1,133 @@
/** \file ScaleFunction.hh
* \brief ScaleFunction Header file. Contains the different scale choice structs.
*/
#pragma once
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/Event.hh"
namespace RHEJ{
//! Abstract base class for calculating event scales
/**
* ScaleFunction is an abstract base class for functors
* that calculate the renormalisation and factorisation
* scales for an event.
*
* We use polymorphism here so that the scale can be chosen
* at run time, e.g. in the configuration file
*/
struct ScaleFunction{
virtual EventParameters operator()(Event const & ev) const = 0;
virtual ~ScaleFunction(){};
};
//! Copy scales from input event
/**
* This ScaleFunction functor just returns the scale choices from
* the input event.
*/
struct InputScales: ScaleFunction{
EventParameters operator()(Event const & ev) const override{
return ev.central();
}
};
//! Choose fixed number as scale
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice given by a single fixed number, for
* example \f$\mu_r = \mu_f = 125\f$
*/
class FixedScale: public ScaleFunction{
public:
FixedScale(double scale): scale_{scale} {}
EventParameters operator()(Event const & ev) const override{
return {scale_, scale_, ev.central().weight};
}
private:
double scale_;
};
//! Choose \f$H_T\f$ as scale
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal to \f$H_T\f$, which is defined as the
* sum of the scalar transverse momenta of all outgoing particles.
*/
struct Ht: public ScaleFunction{
public:
EventParameters operator()(Event const & ev) const override;
};
//! Choose maximum jet transverse momentum as scale
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal to the maximum scalar transverse
* momentum of all jets.
*/
struct MaxJetPperp: public ScaleFunction{
EventParameters operator()(Event const & ev) const override;
};
//! Choose sum of jet invariant masses as scale
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal to the sum of all scalar transverse
* momenta of all jets.
*/
struct JetInvariantMass: public ScaleFunction{
EventParameters operator()(Event const & ev) const override;
};
//! Choose invariant mass of the two hardest jets as scale
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal to the invariant mass between the
* two hardest jets.
*/
struct HardestJetsInvariantMass: public ScaleFunction{
EventParameters operator()(Event const & ev) const override;
};
+ /*
+ * This is a scale function that is used for the reversed intro paper
+ * TODO: remove and replace by some plugin functionality (like for analyses)
+ */
+ struct Max125M12: public ScaleFunction{
+ EventParameters operator()(Event const & ev) const override;
+ };
+
//! Choose a scale proportional to exponential of the rapidity separtion
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal \f$m_h \exp(\Delta y_{fb})$\f
*/
struct DeltaRapScale: public ScaleFunction{
EventParameters operator()(Event const & ev) const override;
};
//! Choose scales as a multiple of some base scale.
/**
* This ScaleFunction functor returns a renormalisation and
* factorisation scale choice equal to the given factor times the
* given base scale choice. For example, passing 0.5 and a pointer to
* an Ht object to the constructor, the generated scale will be
* \f$H_T/2\f$.
*/
class Product: public ScaleFunction{
public:
Product(double factor, std::unique_ptr<ScaleFunction> base_scale):
factor_{factor},
base_scale_{std::move(base_scale)}
{}
EventParameters operator()(Event const & ev) const override;
private:
double factor_;
std::unique_ptr<ScaleFunction> base_scale_;
};
}
diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
index a7494d5..4d98a26 100644
--- a/src/ScaleFunction.cc
+++ b/src/ScaleFunction.cc
@@ -1,64 +1,72 @@
#include "RHEJ/ScaleFunction.hh"
#include <numeric>
#include <cassert>
namespace RHEJ{
namespace{
double ht(Event const & ev){
double result = 0.;
for(size_t i = 0; i < ev.outgoing().size(); ++i){
auto const decay_products = ev.decays().find(i);
if(decay_products == end(ev.decays())){
result += ev.outgoing()[i].perp();
}
else{
for(auto const & particle: decay_products->second){
result += particle.perp();
}
}
}
return result;
}
}
EventParameters Ht::operator()(Event const & ev) const{
const double mu = ht(ev);
return {mu, mu, ev.central().weight};
}
EventParameters MaxJetPperp::operator()(Event const & ev) const{
const double mu = sorted_by_pt(ev.jets()).front().pt();
return {mu, mu, ev.central().weight};
}
EventParameters JetInvariantMass::operator()(Event const & ev) const{
fastjet::PseudoJet sum;
for(const auto & jet: ev.jets()) sum+=jet;
const double mu = sum.m();
return {mu, mu, ev.central().weight};
}
EventParameters HardestJetsInvariantMass::operator()(Event const & ev) const{
const auto jets = sorted_by_pt(ev.jets());
assert(jets.size() >= 2);
const double mu = (jets[0] + jets[1]).m();
return {mu, mu, ev.central().weight};
};
+ EventParameters Max125M12::operator()(Event const & ev) const{
+ const auto jets = sorted_by_pt(ev.jets());
+ assert(jets.size() >= 2);
+ const double m12 = (jets[0] + jets[1]).m();
+ const double mu = std::max(125., m12);
+ return {mu, mu, ev.central().weight};
+ };
+
EventParameters DeltaRapScale::operator()(Event const & ev) const{
double mu = 125.;
const auto & out = ev.outgoing();
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
mu*= exp(out.back().rapidity() - out.front().rapidity());
return {mu, mu, ev.central().weight};
};
EventParameters Product::operator()(Event const &ev) const{
EventParameters s = (*base_scale_)(ev);
s.muf *= factor_;
s.mur *= factor_;
return s;
}
}
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 94d78b7..ae729c2 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,445 +1,448 @@
#include "RHEJ/YAMLreader.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", "max order"
};
// 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;
}
} // namespace anonymous
namespace detail{
} // namespace detail
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>());
}
} // 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;
}
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 RHEJ
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 YAML
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{}};
}
+ if(scale_fun == "max_125_m12"){
+ return ret_type{new Max125M12{}};
+ }
if(scale_fun == "Delta rapidity"){
return ret_type{new DeltaRapScale{}};
}
try{
const double scale = to_double(scale_fun);
return ret_type{new FixedScale{scale}};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
std::unique_ptr<ScaleFunction> make_ScaleFunction_ptr(
std::string const & scale_fun
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const size_t delim = scale_fun.find_first_of("*/");
if(delim == scale_fun.npos){
return make_simple_ScaleFunction_ptr(scale_fun);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
double factor;
std::unique_ptr<ScaleFunction> fun;
if(scale_fun[delim] == '/'){
factor = 1/to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
else{
assert(scale_fun[delim] == '*');
try{
factor = to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
catch(std::invalid_argument const &){
factor = to_double(first);
fun = make_simple_ScaleFunction_ptr(second);
}
}
assert(fun != nullptr);
return std::unique_ptr<ScaleFunction>{
new Product{factor, std::move(fun)}
};
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::reweight},
{unob, EventTreatment::keep},
{unof, EventTreatment::keep},
{nonFKL, EventTreatment::keep}
};
set_from_yaml(treat.at(FKL), yaml, "FKL");
set_from_yaml(treat.at(unob), yaml, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(nonFKL), yaml, "non-FKL");
if(treat[nonFKL] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-FKL events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
try{
assert_all_options_known(yaml, get_supported_options());
}
catch(unknown_option const & ex){
throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt");
config.max_ext_soft_pt_fraction = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
set_from_yaml(config.trials, yaml, "trials");
set_from_yaml(config.log_correction, yaml, "log correction");
set_from_yaml(config.unweight, yaml, "unweight");
config.treat = get_event_treatment(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
set_from_yaml_if_defined(config.RanLux_init, yaml, "RanLux init");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scale_gen = ScaleGenerator{to_ScaleConfig(yaml)};
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
set_from_yaml_if_defined(config.max_order, yaml, "max order");
return config;
}
} // namespace anonymous
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;
}
}
} // namespace RHEJ

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:36 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022793
Default Alt Text
(22 KB)

Event Timeline