diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index b09dabf..149ff0c 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,257 +1,263 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/EventReweighter.hh"
 
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <functional>
 #include <string>
 #include <unordered_map>
 #include <utility>
 
 #include "fastjet/ClusterSequence.hh"
 #include "fastjet/PseudoJet.hh"
 
 #include "LHEF/LHEF.h"
 
 #include "HEJ/Event.hh"
 #include "HEJ/Fraction.hh"
 #include "HEJ/PDG_codes.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PhaseSpacePoint.hh"
 #include "HEJ/exceptions.hh"
 
 namespace HEJ {
 
   EventReweighter::EventReweighter(
       LHEF::HEPRUP const & heprup,
       ScaleGenerator scale_gen,
       EventReweighterConfig conf,
       std::shared_ptr<RNG> ran
   ):
     EventReweighter{
       Beam{
         heprup.EBMUP.first,
         {{
           static_cast<ParticleID>(heprup.IDBMUP.first),
           static_cast<ParticleID>(heprup.IDBMUP.second)
         }}
       },
       heprup.PDFSUP.first,
       std::move(scale_gen),
       std::move(conf),
       std::move(ran)
     }
   {
     if(heprup.EBMUP.second != E_beam_){
       throw std::invalid_argument(
           "asymmetric beam: " + std::to_string(E_beam_)
           + " ---> <--- " + std::to_string(heprup.EBMUP.second)
       );
     }
     if(heprup.PDFSUP.second != pdf_.id()){
       throw std::invalid_argument(
           "conflicting PDF ids: " + std::to_string(pdf_.id())
           + " vs. " + std::to_string(heprup.PDFSUP.second)
       );
     }
   }
 
   EventReweighter::EventReweighter(
       Beam const & beam,
       int pdf_id,
       ScaleGenerator scale_gen,
       EventReweighterConfig conf,
       std::shared_ptr<RNG> ran
   ):
     param_{std::move(conf)},
     E_beam_{beam.E},
     pdf_{pdf_id, beam.type.front(), beam.type.back()},
     MEt2_{
       [this](double mu){ return pdf_.Halphas(mu); },
       param_.ME_config
     },
     scale_gen_{std::move(scale_gen)},
     ran_{std::move(ran)}
   {
     assert(ran_);
   }
 
   PDF const & EventReweighter::pdf() const{
     return pdf_;
   }
 
   std::vector<Event> EventReweighter::reweight(
       Event const & input_ev, std::size_t num_events
   ){
     auto res_events{ gen_res_events(input_ev, num_events) };
     if(res_events.empty()) return {};
     for(auto & event: res_events) event = scale_gen_(std::move(event));
     return rescale(input_ev, std::move(res_events));
   }
 
   EventTreatment EventReweighter::treatment(EventType type) const {
     return param_.treat.at(type);
   }
 
   std::vector<Event> EventReweighter::gen_res_events(
       Event const & ev,
       std::size_t phase_space_points
   ){
     assert(ev.variations().empty());
     status_.clear();
 
     switch(treatment(ev.type())){
     case EventTreatment::discard: {
       status_.emplace_back(StatusCode::discard);
       return {};
     }
     case EventTreatment::keep:
       if(! jets_pass_resummation_cuts(ev)) {
         status_.emplace_back(StatusCode::failed_resummation_cuts);
         return {};
       }
       else {
         status_.emplace_back(StatusCode::good);
         return {ev};
       }
     default:;
     }
     const double Born_shat = shat(ev);
 
     std::vector<Event> resummation_events;
     status_.reserve(phase_space_points);
     for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
       PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
       status_.emplace_back(psp.status());
       assert(psp.status() != StatusCode::unspecified);
       if(psp.status() != StatusCode::good) continue;
       assert(psp.weight() != 0.);
       if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
         status_.back() = StatusCode::too_much_energy;
         continue;
       }
 
       resummation_events.emplace_back(
         to_EventData( std::move(psp) ).cluster(
           param_.jet_param().def, param_.jet_param().min_pt
         )
       );
       auto & new_event = resummation_events.back();
       assert( new_event.valid_hej_state(
         param_.psp_config.max_ext_soft_pt_fraction,
         param_.psp_config.min_extparton_pt ) );
-      if( new_event.type() != ev.type() )
-        throw std::logic_error{"Resummation Event does not match Born event"};
+      if( new_event.type() != ev.type() ) {
+        throw std::logic_error{
+          "Resummation Event does not match Born event: "
+          + name(new_event.type())
+          + " != "
+          + name(ev.type())
+        };
+      }
       new_event.generate_colours(*ran_);
       assert(new_event.variations().empty());
       new_event.central().mur = ev.central().mur;
       new_event.central().muf = ev.central().muf;
       const double resum_shat = shat(new_event);
       new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
         (phase_space_points*resum_shat*resum_shat);
     }
     return resummation_events;
   }
 
   std::vector<Event> EventReweighter::rescale(
       Event const & Born_ev,
       std::vector<Event> events
   ) const{
     const double Born_pdf = pdf_factors(Born_ev).central;
     const double Born_ME = tree_matrix_element(Born_ev);
 
     for(auto & cur_event: events){
       const auto pdf = pdf_factors(cur_event);
       assert(pdf.variations.size() == cur_event.variations().size());
       const auto ME = matrix_elements(cur_event);
       assert(ME.variations.size() == cur_event.variations().size());
       cur_event.parameters() *= pdf*ME/(Born_pdf*Born_ME);
     }
     return events;
   }
 
   bool EventReweighter::jets_pass_resummation_cuts(
       Event const & ev
   ) const{
     const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
     fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
     return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
   }
 
   Weights EventReweighter::pdf_factors(Event const & ev) const{
     auto const & a = ev.incoming().front();
     auto const & b = ev.incoming().back();
     const double xa = a.p.e()/E_beam_;
     const double xb = b.p.e()/E_beam_;
 
     Weights result;
     std::unordered_map<double, double> known_pdf;
     result.central =
       pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
       pdf_.pdfpt(1,xb,ev.central().muf,b.type);
     known_pdf.emplace(ev.central().muf, result.central);
 
     result.variations.reserve(ev.variations().size());
     for(auto const & ev_param: ev.variations()){
       const double muf = ev_param.muf;
       auto cur_pdf = known_pdf.find(muf);
       if(cur_pdf == known_pdf.end()){
         cur_pdf = known_pdf.emplace(
             muf,
             pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
         ).first;
       }
       result.variations.emplace_back(cur_pdf->second);
     }
     assert(result.variations.size() == ev.variations().size());
     return result;
   }
 
   Weights
   EventReweighter::matrix_elements(Event const & ev) const{
     assert(param_.treat.count(ev.type()) > 0);
     if(param_.treat.find(ev.type())->second == EventTreatment::keep){
       return fixed_order_scale_ME(ev);
     }
 
     return MEt2_(ev);
   }
 
   double EventReweighter::tree_matrix_element(Event const & ev) const{
     assert(ev.variations().empty());
     assert(param_.treat.count(ev.type()) > 0);
     if(param_.treat.find(ev.type())->second == EventTreatment::keep){
       return fixed_order_scale_ME(ev).central;
     }
     return MEt2_.tree(ev).central;
   }
 
   Weights
   EventReweighter::fixed_order_scale_ME(Event const & ev) const{
     int alpha_s_power = 0;
     for(auto const & part: ev.outgoing()){
       if(is_parton(part))
         ++alpha_s_power;
       else if(part.type == pid::Higgs) {
         alpha_s_power += 2;
       }
       // nothing to do for other uncoloured particles
     }
 
     Weights result;
     result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
     for(auto const & var: ev.variations()){
       result.variations.emplace_back(
           std::pow(pdf_.Halphas(var.mur), alpha_s_power)
       );
     }
     return result;
   }
 
 } // namespace HEJ