diff --git a/include/HEJ/ScaleFunction.hh b/include/HEJ/ScaleFunction.hh
index 5132a4e..42c2c99 100644
--- a/include/HEJ/ScaleFunction.hh
+++ b/include/HEJ/ScaleFunction.hh
@@ -1,160 +1,174 @@
 /** \file
  *  \brief Functions to calculate the (renormalisation and factorisation) scales for an event
  *
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <functional>
 #include <memory>
 #include <string>
 #include <utility>
 #include <vector>
 
 namespace HEJ{
   class Event;
 
   //! Class to calculate the scale associated with an event
   class ScaleFunction {
   public:
     //! Constructor
     /**
      *  @param name     Name of the scale choice (e.g. H_T)
      *  @param fun      Function used to calculate the scale
      */
     ScaleFunction(std::string name, std::function<double(Event const &)> fun):
       name_{std::move(name)},
       fun_{std::move(fun)}
     {}
 
     //! Name of the scale choice
     std::string const & name() const {
       return name_;
     }
 
     //! Calculate the scale associated with an event
     double operator()(Event const & ev) const {
       return fun_(ev);
     }
 
   private:
     friend ScaleFunction operator*(double factor, ScaleFunction base_scale);
     friend ScaleFunction operator/(ScaleFunction base_scale, double denom);
+    friend ScaleFunction operator*(ScaleFunction func1, ScaleFunction func2);
+    friend ScaleFunction operator/(ScaleFunction func1, ScaleFunction func2);
 
     std::string name_;
     std::function<double(Event const &)> fun_;
   };
 
   //! Multiply a scale choice by a constant factor
   /**
    * For example, multiplying 0.5 and a scale function for H_T
    * will result in a scale function for H_T/2.
    */
   ScaleFunction operator*(double factor, ScaleFunction base_scale);
+  //! Multiply a scale choice by a second one
+  /**
+   * For example, multiplying H_T and m_j1j2
+   * will result in a scale function for H_T*m_j1j2.
+   */
+  ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale);
   //! Divide a scale choice by a constant factor
   /**
    * For example, dividing a scale function for H_T by 2
    * will result in a scale function for H_T/2.
    */
   ScaleFunction operator/(ScaleFunction base_scale, double denom);
+  //! Divide a scale choice by a second one
+  /**
+   * For example, dividing a scale function for H_T by m_j1j2
+   * will result in a scale function for H_T/m_j1j2.
+   */
+  ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom);
 
   //! Calculate \f$H_T\f$ for the input event
   /**
    * \f$H_T\f$ is the sum of the (scalar) transverse momenta of all
    * final-state particles
    */
   double H_T(Event const &);
   //! The maximum of all (scalar) jet transverse momentum
   double max_jet_pt(Event const &);
   //! The invariant mass of the sum of all jet momenta
   double jet_invariant_mass(Event const &);
   //! Invariant mass of the two hardest jets
   double m_j1j2(Event const &);
 
   //! Functor that returns a fixed scale regardless of the input event
   class FixedScale {
   public:
     explicit FixedScale(double mu): mu_{mu} {}
     double operator()(Event const &) const {
       return mu_;
     }
   private:
     double mu_;
   };
 
   class ParameterDescription;
 
   //! Generate combinations of renormalisation and factorisation scales
   class ScaleGenerator{
   public:
     ScaleGenerator() = default;
 
     /** \brief Constructor
      * @param scale_functions_begin   Iterator to first base scale
      * @param scale_functions_end     Iterator past last base scale
      * @param scale_factors_begin     Iterator to first scale factor
      * @param scale_factors_end       Iterator past last scale factor
      * @param max_scale_ratio         Maximum ratio between renormalisation
      *                                and factorisation scale
      */
     template<class ScaleFunIterator, class FactorIterator>
     ScaleGenerator(
         ScaleFunIterator scale_functions_begin,
         ScaleFunIterator scale_functions_end,
         FactorIterator scale_factors_begin,
         FactorIterator scale_factors_end,
         double max_scale_ratio
     ):
       scales_(scale_functions_begin, scale_functions_end),
       scale_factors_(scale_factors_begin, scale_factors_end),
       max_scale_ratio_{max_scale_ratio}
     {
       gen_descriptions();
     }
 
     /** \brief Constructor
      * @param scales            Base scales
      * @param scale_factors     Factors to multiply the base scales
      * @param max_scale_ratio   Maximum ratio between renormalisation
      *                          and factorisation scale
      */
     ScaleGenerator(
         std::vector<ScaleFunction> scales,
         std::vector<double> scale_factors,
         double max_scale_ratio
     ):
       scales_(std::move(scales)),
       scale_factors_(std::move(scale_factors)),
       max_scale_ratio_{max_scale_ratio}
     {
       gen_descriptions();
     }
 
     /** \brief Adjust event parameters, adding scale variation
      *
      * The central renormalisation and factorisation scale of the returned
      * event is given be the first base scale passed to the constructor.
      * The scale variation (stored in event.variation()) is constructed as
      * follows: For each base scale according to the arguments of the
      * constructor we add one variation where both renormalisation and
      * factorisation scale are set according to the current base scale.
      * Then, all combinations where the base renormalisation and factorisation
      * scales are multiplied by one of the scale factors are added.
      * The case were both scales are multiplied by one is skipped.
      * Scale combinations where the ratio is larger than the maximum scale ratio
      * set in the constructor are likewise discarded.
      */
     Event operator()(Event event) const;
 
   private:
     void gen_descriptions();
 
     std::vector<ScaleFunction> scales_;
     std::vector<double> scale_factors_;
     std::vector<std::shared_ptr<ParameterDescription>> descriptions_;
     double max_scale_ratio_;
   };
 
 }
diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
index 67a8e99..a258a05 100644
--- a/src/ScaleFunction.cc
+++ b/src/ScaleFunction.cc
@@ -1,149 +1,172 @@
 /**
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 #include "HEJ/ScaleFunction.hh"
 
 #include <cassert>
 
 #include "HEJ/Event.hh"
 #include "HEJ/exceptions.hh"
 
 namespace HEJ{
 
   double H_T(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;
   }
 
   double max_jet_pt(Event const & ev) {
     return sorted_by_pt(ev.jets()).front().pt();
   }
 
   double jet_invariant_mass(Event const & ev) {
     fastjet::PseudoJet sum;
     for(const auto & jet: ev.jets()) sum+=jet;
     return sum.m();
   }
 
   double m_j1j2(Event const & ev) {
     const auto jets = sorted_by_pt(ev.jets());
     assert(jets.size() >= 2);
     return (jets[0] + jets[1]).m();
   }
 
   ScaleFunction operator*(double factor, ScaleFunction base_scale) {
     base_scale.name_.insert(0, std::to_string(factor) + '*');
     auto new_fun =
       [factor,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
       return factor*fun(ev);
     };
     base_scale.fun_ = std::move(new_fun);
     return base_scale;
   }
 
+  ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale) {
+    std::cout << factor.name() << " * " << base_scale.name() << std::endl;
+    base_scale.name_.insert(0, factor.name_ + '*');
+    auto new_fun =
+      [fun1{std::move(factor.fun_)}, fun2{std::move(base_scale.fun_)}, name{base_scale.name_}]
+      (HEJ::Event const & ev) {
+        return fun1(ev)*fun2(ev);
+    };
+    base_scale.fun_ = std::move(new_fun);
+    return base_scale;
+  }
+
   ScaleFunction operator/(ScaleFunction base_scale, double denom) {
     base_scale.name_.append('/' + std::to_string(denom));
     auto new_fun =
       [denom,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
       return fun(ev)/denom;
     };
     base_scale.fun_ = std::move(new_fun);
     return base_scale;
   }
 
+  ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom) {
+    base_scale.name_.append('/' + denom.name_);
+    auto new_fun =
+      [fun2{std::move(denom.fun_)}, fun1{std::move(base_scale.fun_)}]
+      (HEJ::Event const & ev) {
+        return fun1(ev)/fun2(ev);
+    };
+    base_scale.fun_ = std::move(new_fun);
+    return base_scale;
+  }
+
   // TODO: significant logic duplication with operator()
   void ScaleGenerator::gen_descriptions() {
     if(scales_.empty()) {
       throw std::logic_error{"Need at least one scale"};
     }
     descriptions_.emplace_back(
         std::make_shared<ParameterDescription>(scales_.front().name(), 1., 1.)
     );
 
     for(auto & base_scale: scales_){
       const auto base_name = base_scale.name();
       descriptions_.emplace_back(
           std::make_shared<ParameterDescription>(base_name, 1., 1.)
       );
       //multiplicative scale variation
       for(double mur_factor: scale_factors_){
         for(double muf_factor: scale_factors_){
           if(muf_factor == 1. && mur_factor == 1.) continue;
           if(
               mur_factor/muf_factor < 1/max_scale_ratio_
               || mur_factor/muf_factor > max_scale_ratio_
           ) continue;
           descriptions_.emplace_back(
               std::make_shared<ParameterDescription>(
                   base_name, mur_factor, muf_factor
               )
           );
         }
       }
     }
   }
 
   Event ScaleGenerator::operator()(Event ev) const {
     if(! ev.variations().empty()) {
       throw std::invalid_argument{"Input event already has scale variation"};
     }
     assert(!scales_.empty());
     assert(!descriptions_.empty());
 
     size_t descr_idx = 0;
     const double mu_central = (scales_.front())(ev);
     ev.central().mur = mu_central;
     ev.central().muf = mu_central;
     assert(descr_idx < descriptions_.size());
     assert(descriptions_[descr_idx]);
     ev.central().description = descriptions_[descr_idx++];
 
     // check if we are doing scale variation at all
     if(scales_.size() == 1 && scale_factors_.empty()) return ev;
 
     for(auto & base_scale: scales_){
       const double mu_base = base_scale(ev);
       assert(descr_idx < descriptions_.size());
       assert(descriptions_[descr_idx]);
       ev.variations().emplace_back(
           EventParameters{
             mu_base, mu_base, ev.central().weight, descriptions_[descr_idx++]
           }
       );
       //multiplicative scale variation
       for(double mur_factor: scale_factors_){
         const double mur = mu_base*mur_factor;
         for(double muf_factor: scale_factors_){
           if(muf_factor == 1. && mur_factor == 1.) continue;
           const double muf = mu_base*muf_factor;
           if(
               mur/muf < 1/max_scale_ratio_
               || mur/muf > max_scale_ratio_
           ) continue;
           assert(descr_idx < descriptions_.size());
           assert(descriptions_[descr_idx]);
           ev.variations().emplace_back(
               EventParameters{
                 mur, muf, ev.central().weight, descriptions_[descr_idx++]
               }
           );
         }
       }
     };
     return ev;
   }
 
 }
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 5336e8d..e5ce822 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,469 +1,464 @@
 /**
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \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"
 
 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",
           "FKL", "unordered", "non-HEJ",
           "scales", "scale factors", "max scale ratio", "import scales",
           "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);
+        return parse_ScaleFunction(first, known)
+          / parse_ScaleFunction(second, known);
       }
       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(second, known);
-        }
+        return parse_ScaleFunction(first, known)
+          * parse_ScaleFunction(second, 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");
       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