diff --git a/config.yml b/config.yml index 508fcb1..8288cea 100644 --- a/config.yml +++ b/config.yml @@ -1,90 +1,89 @@ # number of attempted resummation phase space points for each input event trials: 10 min extparton pt: 30 # minimum transverse momentum of extremal partons # maximum soft transverse momentum fraction in extremal jets # # max ext soft pt fraction: 0.1 resummation jets: # resummation jet properties min pt: 35 # minimum jet transverse momentum algorithm: antikt # jet clustering algorithm R: 0.4 # jet R parameter fixed order jets: # properties of input jets min pt: 30 # by default, algorithm and R are like for resummation jets # treatment of he various event classes # the supported settings are: reweight, keep, discard # non-HEJ events cannot be reweighted FKL: reweight unordered: keep non-HEJ: keep # central scale choice or choices # # scales: [125, max jet pperp, H_T/2, 2*jet invariant mass, m_j1j2] scales: 91.188 # factors by which the central scales should be multiplied # renormalisation and factorisation scales are varied independently # # scale factors: [0.5, 0.7071, 1, 1.41421, 2] # maximum ratio between renormalisation and factorisation scale # # max scale ratio: 2.0001 # import scale setting functions # # import scales: # lib_my_scales.so: [scale0,scale1] log correction: false # whether or not to include higher order logs -unweight: false # TODO: whether or not to unweight events # event output files # # the supported formats are # - Les Houches (suffix .lhe) # - HepMC (suffix .hepmc3) # TODO: - ROOT ntuples (suffix .root) # # An output file's format is deduced either automatically from the suffix # or from an explicit specification, e.g. # - Les Houches: outfile event output: - HEJ.lhe # - HEJ_events.hepmc # to use a rivet analysis # # analysis: # rivet: MC_XS # rivet analysis name # output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added # # to use a custom analysis # # analysis: # plugin: /path/to/libmyanalysis.so # my analysis parameter: some value # selection of random number generator and seed # the choices are # - mixmax (seed is an integer) # - ranlux64 (seed is a filename containing parameters) random generator: name: mixmax # seed: 1 # parameters for Higgs-gluon couplings # this requires compilation with qcdloop # # Higgs coupling: # use impact factors: false # mt: 174 # include bottom: true # mb: 4.7 diff --git a/doc/sphinx/HEJ.rst b/doc/sphinx/HEJ.rst index b36f709..5beb42e 100644 --- a/doc/sphinx/HEJ.rst +++ b/doc/sphinx/HEJ.rst @@ -1,281 +1,279 @@ .. _`Running HEJ 2`: Running HEJ 2 ============= Quick start ----------- In order to run HEJ 2, you need a configuration file and a file containing fixed-order events. A sample configuration is given by the :file:`config.yml` file distributed together with HEJ 2. Events in the Les Houches Event File format can be generated with standard Monte Carlo generators like `MadGraph5_aMC@NLO <https://launchpad.net/mg5amcnlo>`_ or `Sherpa <https://sherpa.hepforge.org/trac/wiki>`_. HEJ 2 assumes that the cross section is given by the sum of the event weights. Depending on the fixed-order generator it may be necessary to adjust the weights in the Les Houches Event File accordingly. The processes supported by HEJ 2 are - Pure multijet production - Production of a Higgs boson with jets .. - *TODO* Production of a W boson with jets - *TODO* Production of a Z boson or photon with jets where at least two jets are required in each case. For the time being, only leading-order events are supported. After generating an event file :file:`events.lhe` adjust the parameters under the `fixed order jets`_ setting in :file:`config.yml` to the settings in the fixed-order generation. Resummation can then be added by running:: HEJ config.yml events.lhe Using the default settings, this will produce an output event file :file:`HEJ.lhe` with events including high-energy resummation. .. _`HEJ 2 settings`: Settings -------- HEJ 2 configuration files follow the `YAML <http://yaml.org/>`_ format. The following configuration parameters are supported: .. _`trials`: **trials** High-energy resummation is performed by generating a number of resummation phase space configurations corresponding to the input fixed-order event. This parameter specifies how many such configurations HEJ 2 should try to generate for each input event. Typical values vary between 10 and 100. .. _`min extparton pt`: **min extparton pt** Specifies the minimum transverse momentum in GeV of the most forward and the most backward parton. This setting is needed to regulate an otherwise uncancelled divergence. Its value should be slightly below the minimum transverse momentum of jets specified by `resummation jets: min pt`_. See also the `max ext soft pt fraction`_ setting. .. _`max ext soft pt fraction`: **max ext soft pt fraction** Specifies the maximum fraction that soft radiation can contribute to the transverse momentum of each the most forward and the most backward jet. Values between around 0.05 and 0.1 are recommended. See also the `min extparton pt`_ setting. .. _`fixed order jets`: **fixed order jets** This tag collects a number of settings specifying the jet definition in the event input. The settings should correspond to the ones used in the fixed-order Monte Carlo that generated the input events. .. _`fixed order jets: min pt`: **min pt** Minimum transverse momentum in GeV of fixed-order jets. .. _`fixed order jets: algorithm`: **algorithm** The algorithm used to define jets. Allowed settings are :code:`kt`, :code:`cambridge`, :code:`antikt`, :code:`cambridge for passive`. See the `FastJet <http://fastjet.fr/>`_ documentation for a description of these algorithms. .. _`fixed order jets: R`: **R** The R parameter used in the jet algorithm, roughly corresponding to the jet radius in the plane spanned by the rapidity and the azimuthal angle. .. _`resummation jets`: **resummation jets** This tag collects a number of settings specifying the jet definition in the observed, i.e. resummed events. These settings are optional, by default the same values as for the `fixed order jets`_ are assumed. .. _`resummation jets: min pt`: **min pt** Minimum transverse momentum in GeV of resummation jets. This should be between 25% and 50% larger than the minimum transverse momentum of fixed order jets set by `fixed order jets: min pt`_. .. _`resummation jets: algorithm`: **algorithm** The algorithm used to define jets. The HEJ 2 approach to resummation relies on properties of :code:`antikt` jets, so this value is strongly recommended. For a list of possible other values, see the `fixed order jets: algorithm`_ setting. .. _`resummation jets: R`: **R** The R parameter used in the jet algorithm. .. _`FKL`: **FKL** Specifies how to treat FKL events. The possible values are :code:`reweight` to enable resummation, :code:`keep` to keep the events as they are up to a possible change of renormalisation and factorisation scale, and :code:`discard` to discard these events. .. _`unordered`: **unordered** Specifies how to treat events with one emission that does not respect FKL ordering. The possible values are the same as for the `FKL`_ setting, but :code:`reweight` may not be supported for all process types. .. _`non-HEJ`: **non-HEJ** Specifies how to treat events where no resummation is possible. The allowed values are :code:`keep` to keep the events as they are up to a possible change of renormalisation and factorisation scale and :code:`discard` to discard these events. .. _`scales`: **scales** Specifies the renormalisation and factorisation scales for the output events. This can either be a single entry or a list :code:`[scale1, scale2, ...]`. For the case of a list the first entry defines the central scale. Possible values are fixed numbers to set the scale in GeV or the following: - :code:`H_T`: The sum of the scalar transverse momenta of all final-state particles - :code:`max jet pperp`: The maximum transverse momentum of all jets - :code:`jet invariant mass`: Sum of the invariant masses of all jets - :code:`m_j1j2`: Invariant mass between the two hardest jets. Scales can be multiplied or divided by an overall factor, e.g. :code:`H_T/2`. It is also possible to import scales from an external library, see :ref:`Custom scales` .. _`scale factors`: **scale factors** A list of numeric factors by which each of the `scales`_ should be multiplied. Renormalisation and factorisation scales are varied independently. For example, a list with entries :code:`[0.5, 2]` would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`); (0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\ :sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to the order of the final event weights. .. _`max scale ratio`: **max scale ratio** Specifies the maximum factor by which renormalisation and factorisation scales may difer. For a value of :code:`2` and the example given for the `scale factors`_ the scale choices (0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`) will be discarded. .. _`log correction`: **log correction** Whether to include corrections due to the evolution of the strong coupling constant in the virtual corrections. Allowed values are :code:`true` and :code:`false`. -.. TODO: unweight - .. _`event output`: **event output** Specifies the name of a single event output file or a list of such files. The file format is either specified explicitly or derived from the suffix. For example, :code:`events.lhe` or, equivalently :code:`Les Houches: events.lhe` generates an output event file :code:`events.lhe` in the Les Houches format. The supported formats are - :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches event file format. - :code:`file.hepmc` or :code:`HepMC: file`: The HepMC format. .. _`random generator`: **random generator** Sets parameters for random number generation. .. _`random generator: name`: **name** Which random number generator to use. Currently, :code:`mixmax` and :code:`ranlux64` are implemented. Mixmax is recommended. See the `CLHEP documentation <http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for details on the generators. .. _`random generator: seed`: **seed** The seed for random generation. This should be a single number for :code:`mixmax` and the name of a state file for :code:`ranlux64`. .. _`analysis`: **analysis** Name and Setting for the event analyses; either a custom analysis plugin or Rivet. For the first the :code:`plugin` sub-entry should be set to the analysis file path. All further entries are passed on to the analysis. To use Rivet a list of Rivet analyses have to be given in :code:`rivet` and prefix for the yoda file has to be set through :code:`output`. See :ref:`Writing custom analyses` for details. .. _`Higgs coupling`: **Higgs coupling** This collects a number of settings concerning the effective coupling of the Higgs boson to gluons. This is only relevant for the production process of a Higgs boson with jets and only supported if HEJ 2 was compiled with `QCDLoop <https://github.com/scarrazza/qcdloop>`_ support. .. _`Higgs coupling: use impact factors`: **use impact factors** Whether to use impact factors for the coupling to the most forward and most backward partons. Impact factors imply the infinite top-quark mass limit. .. _`Higgs coupling: mt`: **mt** The value of the top-quark mass in GeV. If this is not specified, the limit of an infinite mass is taken. .. _`Higgs coupling: include bottom`: **include bottom** Whether to include the Higgs coupling to bottom quarks. .. _`Higgs coupling: mb`: **mb** The value of the bottom-quark mass in GeV. Only used for the Higgs coupling, external bottom-quarks are always assumed to be massless. diff --git a/include/HEJ/config.hh b/include/HEJ/config.hh index 7c51c05..962bfa6 100644 --- a/include/HEJ/config.hh +++ b/include/HEJ/config.hh @@ -1,173 +1,171 @@ /** \file * \brief HEJ 2 configuration parameters */ #pragma once #include <string> #include <memory> #include "fastjet/JetDefinition.hh" #include "yaml-cpp/yaml.h" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/optional.hh" #include "HEJ/output_formats.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/utility.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; }; /**! 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>; /**! 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 double min_extparton_pt; //! Maximum transverse momentum fraction from soft radiation in extremal jets double max_ext_soft_pt_fraction; //! Number of resummation configurations to generate per fixed-order event int trials; //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction; - //! Whether to perform unweighting - bool unweight; //! 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; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; }; //! Configuration options for the PhaseSpacePoint class struct PhaseSpacePointConfig { //! Properties of resummation jets JetParameters jet_param; //! Minimum transverse momentum for extremal partons double min_extparton_pt; //! Maximum transverse momentum fraction from soft radiation in extremal jets double max_ext_soft_pt_fraction; }; //! Configuration options for the MatrixElement class struct MatrixElementConfig { //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; }; //! 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; //! Properties of resummation jets JetParameters 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}; } /**! Extract EventReweighterConfig from Config * * @see to_PhaseSpacePointConfig, to_MatrixElementConfig */ inline EventReweighterConfig to_EventReweighterConfig(Config const & conf) { return { to_PhaseSpacePointConfig(conf), to_MatrixElementConfig(conf), conf.resummation_jets, conf.treat }; } } // namespace HEJ diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc index 769edaf..608f010 100644 --- a/src/YAMLreader.cc +++ b/src/YAMLreader.cc @@ -1,458 +1,457 @@ #include "HEJ/YAMLreader.hh" #include <set> #include <string> #include <vector> #include <iostream> #include <stdexcept> #include <dlfcn.h> namespace HEJ{ 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-HEJ", "scales", "scale factors", "max scale ratio", "import scales", - "log correction", "unweight", "event output", "analysis" + "log correction", "event output", "analysis" }; // 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] = ""; } 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; } } // 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>()); } } // 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; } 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} }; 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 * 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())) ); const size_t delim = scale_fun.find_first_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_simple_ScaleFunction(first, known)/to_double(second); } else{ assert(scale_fun[delim] == '*'); try{ const double factor = to_double(second); return factor*parse_simple_ScaleFunction(first, known); } catch(std::invalid_argument const &){ const double factor = to_double(first); return factor*parse_simple_ScaleFunction(first, known); } } } 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}, {nonHEJ, 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(nonHEJ), yaml, "non-HEJ"); if(treat[nonHEJ] == EventTreatment::reweight){ throw std::invalid_argument{"Cannot reweight non-HEJ 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"); config.rng = to_RNGConfig(yaml, "random generator"); set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis"); config.scales = to_ScaleConfig(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/t/ME_data/config_mt.yml b/t/ME_data/config_mt.yml index 4656032..5a73374 100644 --- a/t/ME_data/config_mt.yml +++ b/t/ME_data/config_mt.yml @@ -1,29 +1,28 @@ trials: 1 min extparton pt: 30 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 FKL: reweight unordered: reweight non-HEJ: discard scales: 125 log correction: false -unweight: false random generator: name: mixmax seed: 1 Higgs coupling: use impact factors: false mt: 174 include bottom: false diff --git a/t/ME_data/config_mtinf.yml b/t/ME_data/config_mtinf.yml index 0dbc0c1..a7c71b6 100644 --- a/t/ME_data/config_mtinf.yml +++ b/t/ME_data/config_mtinf.yml @@ -1,24 +1,23 @@ trials: 1 min extparton pt: 30 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 FKL: reweight unordered: reweight non-HEJ: discard scales: 125 log correction: false -unweight: false random generator: name: mixmax seed: 1 diff --git a/t/ME_data/config_mtmb.yml b/t/ME_data/config_mtmb.yml index d2c5d7b..c23cb8c 100644 --- a/t/ME_data/config_mtmb.yml +++ b/t/ME_data/config_mtmb.yml @@ -1,30 +1,29 @@ trials: 1 min extparton pt: 30 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 FKL: reweight unordered: reweight non-HEJ: discard scales: 125 log correction: false -unweight: false random generator: name: mixmax seed: 1 Higgs coupling: use impact factors: false mt: 174 include bottom: true mb: 4.7 diff --git a/t/jet_config.yml b/t/jet_config.yml index ced0d0c..1c56c5c 100644 --- a/t/jet_config.yml +++ b/t/jet_config.yml @@ -1,26 +1,25 @@ trials: 10 min extparton pt: 30 resummation jets: min pt: 35 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 FKL: reweight unordered: discard non-HEJ: discard log correction: false -unweight: false scales: 91.188 random generator: name: ranlux64 event output: - tst.lhe diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml index 467f765..4ad014a 100644 --- a/t/jet_config_with_import.yml +++ b/t/jet_config_with_import.yml @@ -1,29 +1,28 @@ trials: 10 min extparton pt: 30 resummation jets: min pt: 35 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 FKL: reweight unordered: discard non-HEJ: discard log correction: false -unweight: false scales: softest_jet_pt event output: - tst.lhe random generator: name: ranlux64 import scales: ./libscales.so: softest_jet_pt