diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 85b06c5..d9cc041 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,375 +1,374 @@
 /** \file
  *  \brief Declares the Event class and helpers
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <array>
 #include <cstddef>
 #include <iosfwd>
 #include <iterator>
 #include <unordered_map>
 #include <utility>
 #include <vector>
-#include <numeric>
 
 #include "boost/iterator/filter_iterator.hpp"
 
 #include "fastjet/ClusterSequence.hh"
 #include "fastjet/PseudoJet.hh"
 
 #include "HEJ/Constants.hh"
 #include "HEJ/event_types.hh"
 #include "HEJ/Parameters.hh"
 #include "HEJ/Particle.hh"
 
 namespace LHEF {
   class HEPEUP;
   class HEPRUP;
 }
 
 namespace fastjet {
   class JetDefinition;
 }
 
 namespace HEJ {
   struct RNG;
   struct UnclusteredEvent;
 
   /** @brief An event with clustered jets
     *
     * This is the main HEJ 2 event class.
     * It contains kinematic information including jet clustering,
     * parameter (e.g. scale) settings and the event weight.
     */
   class Event {
   public:
     class EventData;
 
     //! Iterator over partons
     using ConstPartonIterator = boost::filter_iterator<
       bool (*)(Particle const &),
       std::vector<Particle>::const_iterator
       >;
     //! Reverse Iterator over partons
     using ConstReversePartonIterator = std::reverse_iterator<
                                                 ConstPartonIterator>;
     //! No default Constructor
     Event() = delete;
     //! Event Constructor adding jet clustering to an unclustered event
     //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
     [[deprecated("UnclusteredEvent will be replaced by EventData")]]
     Event(
       UnclusteredEvent const & ev,
       fastjet::JetDefinition const & jet_def, double min_jet_pt
     );
 
     //! @name Particle Access
     //! @{
 
     //! Incoming particles
     std::array<Particle, 2> const &  incoming() const{
       return incoming_;
     }
     //! Outgoing particles
     std::vector<Particle> const &  outgoing() const{
       return outgoing_;
     }
     //! Iterator to the first outgoing parton
     ConstPartonIterator begin_partons() const;
     //! Iterator to the first outgoing parton
     ConstPartonIterator cbegin_partons() const;
 
     //! Iterator to the end of the outgoing partons
     ConstPartonIterator end_partons() const;
     //! Iterator to the end of the outgoing partons
     ConstPartonIterator cend_partons() const;
 
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator rbegin_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator crbegin_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator rend_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator crend_partons() const;
 
     //! Particle decays
     /**
      *  The key in the returned map corresponds to the index in the
      *  vector returned by outgoing()
      */
     std::unordered_map<std::size_t, std::vector<Particle>> const & decays()
     const {
       return decays_;
     }
     //! The jets formed by the outgoing partons, sorted in rapidity
     std::vector<fastjet::PseudoJet> const & jets() const{
       return jets_;
     }
     //! @}
 
     //! @name Weight variations
     //! @{
 
     //! All chosen parameter, i.e. scale choices (const version)
     Parameters<EventParameters> const & parameters() const{
       return parameters_;
     }
     //! All chosen parameter, i.e. scale choices
     Parameters<EventParameters> & parameters(){
       return parameters_;
     }
 
     //! Central parameter choice (const version)
     EventParameters const & central() const{
       return parameters_.central;
     }
     //! Central parameter choice
     EventParameters & central(){
       return parameters_.central;
     }
 
     //! Parameter (scale) variations (const version)
     std::vector<EventParameters> const & variations() const{
       return parameters_.variations;
     }
     //! Parameter (scale) variations
     std::vector<EventParameters> & variations(){
       return parameters_.variations;
     }
 
     //! Parameter (scale) variation (const version)
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters const & variations(std::size_t i) const{
       return parameters_.variations.at(i);
     }
     //! Parameter (scale) variation
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters & variations(std::size_t i){
       return parameters_.variations.at(i);
     }
     //! @}
 
     //! Indices of the jets the outgoing partons belong to
     /**
      *  @param jets   Jets to be tested
      *  @returns      A vector containing, for each outgoing parton,
      *                the index in the vector of jets the considered parton
      *                belongs to. If the parton is not inside any of the
      *                passed jets, the corresponding index is set to -1.
      */
     std::vector<int> particle_jet_indices(
         std::vector<fastjet::PseudoJet> const & jets
     ) const {
       return cs_.particle_jet_indices(jets);
     }
     //! particle_jet_indices() of the Event jets()
     std::vector<int> particle_jet_indices() const {
       return particle_jet_indices(jets());
     }
 
     //! Jet definition used for clustering
     fastjet::JetDefinition const & jet_def() const{
       return cs_.jet_def();
     }
 
     //! Minimum jet transverse momentum
     double min_jet_pt() const{
       return min_jet_pt_;
     }
 
     //! Event type
     event_type::EventType type() const{
       return type_;
     }
 
     //! Give colours to each particle
     /**
      * @returns true if new colours are generated, i.e. same as is_resummable()
      * @details Colour ordering is done according to leading colour in the MRK
      *          limit, see \cite Andersen:2011zd. This only affects \ref
      *          is_resummable() "HEJ" configurations, all other \ref event_type
      *          "EventTypes" will be ignored.
      * @note    This overwrites all previously set colours.
      */
     bool generate_colours(RNG & /*ran*/);
 
     //! Check that current colours are leading in the high energy limit
     /**
      * @details Checks that the colour configuration can be split up in
      *          multiple, rapidity ordered, non-overlapping ladders. Such
      *          configurations are leading in the MRK limit, see
      *          \cite Andersen:2011zd
      *
      * @note This is _not_ to be confused with \ref is_resummable(), however
      *       for all resummable states it is possible to create a leading colour
      *       configuration, see generate_colours()
      */
     bool is_leading_colour() const;
 
     /**
      * @brief Check if given event could have been produced by HEJ
      * @details A HEJ state has to fulfil:
      *          1. type() has to be \ref is_resummable() "resummable"
      *          2. Soft radiation in the tagging jets contributes at most to
      *             `soft_pt_regulator` of the total jet \f$ p_\perp \f$
      *
      * @note This is true for any resummed stated produced by the
      *       EventReweighter or any \ref is_resummable() "resummable" Leading
      *       Order state.
      *
      * @param soft_pt_regulator  Maximum transverse momentum fraction from soft
      *                           radiation in tagging jets
      * @param min_pt             Absolute minimal \f$ p_\perp \f$,
      *                           \b deprecated use soft_pt_regulator instead
      * @return True if this state could have been produced by HEJ
      */
     bool valid_hej_state(
       double soft_pt_regulator = DEFAULT_SOFT_PT_REGULATOR,
       double min_pt = 0.
     ) const;
 
   private:
     //! \internal
     //! @brief Construct Event explicitly from input.
     /** This is only intended to be called from EventData.
      *
      * \warning The input is taken _as is_, sorting and classification has to be
      *          done externally, i.e. by EventData
      */
     Event(
       std::array<Particle, 2> && incoming,
       std::vector<Particle> && outgoing,
       std::unordered_map<std::size_t, std::vector<Particle>> && decays,
       Parameters<EventParameters> && parameters,
       fastjet::JetDefinition const & jet_def,
       double min_jet_pt
     );
 
     //! Iterator over partons (non-const)
     using PartonIterator = boost::filter_iterator<
       bool (*)(Particle const &),
       std::vector<Particle>::iterator
       >;
     //! Reverse Iterator over partons (non-const)
     using ReversePartonIterator = std::reverse_iterator<PartonIterator>;
 
     //! Iterator to the first outgoing parton (non-const)
     PartonIterator begin_partons();
     //! Iterator to the end of the outgoing partons (non-const)
     PartonIterator end_partons();
 
     //! Reverse Iterator to the first outgoing parton (non-const)
     ReversePartonIterator rbegin_partons();
     //! Reverse Iterator to the first outgoing parton (non-const)
     ReversePartonIterator rend_partons();
 
     std::array<Particle, 2> incoming_;
     std::vector<Particle> outgoing_;
     std::unordered_map<std::size_t, std::vector<Particle>> decays_;
     std::vector<fastjet::PseudoJet> jets_;
     Parameters<EventParameters> parameters_;
     fastjet::ClusterSequence cs_;
     double min_jet_pt_;
     event_type::EventType type_;
   }; // end class Event
 
   //! Class to store general Event setup, i.e. Phase space and weights
   class Event::EventData {
   public:
     //! Default Constructor
     EventData() = default;
     //! Constructor from LesHouches event information
     EventData(LHEF::HEPEUP const & hepeup);
     //! Constructor with all values given
     EventData(
       std::array<Particle, 2> incoming,
       std::vector<Particle> outgoing,
       std::unordered_map<std::size_t, std::vector<Particle>> decays,
       Parameters<EventParameters> parameters
     ):
       incoming(std::move(incoming)), outgoing(std::move(outgoing)),
       decays(std::move(decays)), parameters(std::move(parameters))
     {}
 
     //! Generate an Event from the stored EventData.
     /**
      * @details          Do jet clustering and classification.
      *                   Use this to generate an Event.
      *
      * @note             Calling this function destroys EventData
      *
      * @param jet_def    Jet definition
      * @param min_jet_pt minimal \f$p_T\f$ for each jet
      *
      * @returns          Full clustered and classified event.
      */
     Event cluster(
       fastjet::JetDefinition const & jet_def, double min_jet_pt);
 
     //! Alias for cluster()
     Event operator()(
       fastjet::JetDefinition const & jet_def, double const min_jet_pt){
       return cluster(jet_def, min_jet_pt);
     }
 
     //! Sort particles in rapidity
     void sort();
 
     //! Reconstruct intermediate particles from final-state leptons
     /**
      *  Final-state leptons are created from virtual photons, W, or Z bosons.
      *  This function tries to reconstruct such intermediate bosons if they
      *  are not part of the event record.
      */
     void reconstruct_intermediate();
 
     //! Incoming particles
     std::array<Particle, 2> incoming;
     //! Outcoing particles
     std::vector<Particle> outgoing;
     //! Particle decays in the format {outgoing index, decay products}
     std::unordered_map<std::size_t, std::vector<Particle>> decays;
     //! Parameters, e.g. scale or inital weight
     Parameters<EventParameters> parameters;
   }; // end class EventData
 
   //! Print Event
   std::ostream& operator<<(std::ostream & os, Event const & ev);
 
   //! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
   double shat(Event const & ev);
 
   //! Convert an event to a LHEF::HEPEUP
   LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * /*heprup*/);
 
   // put deprecated warning at the end, so don't get the warning inside Event.hh,
   // additionally doxygen can not identify [[deprecated]] correctly
   struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
     UnclusteredEvent;
   //! An event before jet clustering
   //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
   struct UnclusteredEvent{
     //! Default Constructor
     UnclusteredEvent() = default;
     //! Constructor from LesHouches event information
     UnclusteredEvent(LHEF::HEPEUP const & hepeup);
 
     std::array<Particle, 2> incoming;          /**< Incoming Particles */
     std::vector<Particle> outgoing;            /**< Outgoing Particles */
     //! Particle decays in the format {outgoing index, decay products}
     std::unordered_map<std::size_t, std::vector<Particle>> decays;
     //! Central parameter (e.g. scale) choice
     EventParameters central;
     std::vector<EventParameters> variations;    /**< For parameter variation */
   };
 
 } // namespace HEJ