diff --git a/Changes.md b/Changes.md
index 1e17c7e..edf7b10 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,50 +1,51 @@
 # Version 2.0
 
 ## 2.X.0
 
 * Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
 * Allow multiplication and division of multiple scale functions e.g.
   `H_T/2*m_j1j2`
 * Follow HepMC convention for particle Status codes: incoming = 11,
   decaying = 2, outgoing = 1 (unchanged)
 * New class `CrossSectionAccumulator` to keep track of Cross Section of the
   different subproccess
 * New template struct `Parameters` similar to old `Weights`
   - `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
     not change
   - `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
     will be removed in HEJ Version 2.3.0
 * Function to multiplication and division of `EventParameters.weight` by double
   - This can be combined with `Parameters`, e.g.
     `Parameters<EventParameters>*Weights`, see also `Events.parameters()`
   - Moved `EventParameters` to `Parameters.hh` header
 * Restructured `Event` class
   - `Event` can now only be build from a (new) `Event::EventData` class
   - Removed default constructor for `Event`
   - `Event::EventData` replaces the old `UnclusteredEvent` struct.
   - `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
     2.3.0
+  - Removed `Event.unclustered()` function
   - Added new member function `Events.parameters()`, to directly access
     (underlying) `Parameters<EventParameters>`
 
 ## 2.0.4
 
 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
 * Correctly include rivet headers
 * New `EXCLUDE_package` variable in `cmake` to not interface to specific
   packages
 * Consistent search and include for side packages in `cmake`
 
 ## 2.0.3
 
 * Fixed parsing of (numerical factor) * (base scale) in configuration
 * Don't change scale names, but sanitise Rivet output file names instead
 
 ## 2.0.2
 
 * Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
   `"/"` and `"*"` before)
 
 ## 2.0.1
 
 * Fixed name of fixed-order generator in error message.
diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 3fe740a..a1084fe 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,370 +1,363 @@
 /** \file
  *  \brief Declares the Event class and helpers
  *
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 
 #pragma once
 
 #include <array>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
 #include "HEJ/event_types.hh"
 #include "HEJ/Parameters.hh"
 #include "HEJ/Particle.hh"
 
 #include "fastjet/ClusterSequence.hh"
 
 namespace LHEF{
   class HEPEUP;
   class HEPRUP;
 }
 
 namespace fastjet{
   class JetDefinition;
 }
 
 namespace HEJ{
 
   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.
     *
     * \note Use EventData to build this class.
     *       There is no other constructor available.
     */
   class Event{
   public:
     class EventData;
     //! No default Constructor
     Event() = delete;
     //! Event Constructor adding jet clustering to an unclustered event
     //! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.3.0
     [[deprecated("UnclusteredEvent will be replaced by EventData")]]
     Event(
       UnclusteredEvent const & ev,
       fastjet::JetDefinition const & jet_def, double min_jet_pt
     );
 
 
     //! The jets formed by the outgoing partons
     std::vector<fastjet::PseudoJet> jets() const;
 
-    //! The corresponding event before jet clustering
-    // [[deprecated]]
-    // UnclusteredEvent unclustered() const {
-      // return UnclusteredEvent(ev_);
-      // TODO what to do with this?
-    // }
-
     //! Incoming particles
     std::array<Particle, 2> const &  incoming() const{
       return incoming_;
     }
     //! Outgoing particles
     std::vector<Particle> const &  outgoing() const{
       return outgoing_;
     }
     //! Particle decays
     /**
      *  The key in the returned map corresponds to the index in the
      *  vector returned by outgoing()
      */
     std::unordered_map<size_t, std::vector<Particle>> const &  decays() const{
       return decays_;
     }
 
     //! 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(size_t i) const{
       return parameters_.variations[i];
     }
     //! Parameter (scale) variation
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters & variations(size_t i){
       return parameters_.variations[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);
     }
 
     //! 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_;
     }
 
   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<size_t, std::vector<Particle>> && decays,
       Parameters<EventParameters> && parameters,
       fastjet::JetDefinition const & jet_def,
       double const min_jet_pt
     ): incoming_{std::move(incoming)},
        outgoing_{std::move(outgoing)},
        decays_{std::move(decays)},
        parameters_{std::move(parameters)},
        cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
        min_jet_pt_{min_jet_pt}
     {};
 
     std::array<Particle, 2> incoming_;
     std::vector<Particle> outgoing_;
     std::unordered_map<size_t, std::vector<Particle>> decays_;
     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> const & incoming,
       std::vector<Particle> const & outgoing,
       std::unordered_map<size_t, std::vector<Particle>> const & decays,
       Parameters<EventParameters> && parameters
     ):
       incoming_(incoming), outgoing_(outgoing),
       decays_(decays), parameters_(std::move(parameters))
     {};
     //! Move Constructor with all values given
     EventData(
       std::array<Particle, 2> && incoming,
       std::vector<Particle> && outgoing,
       std::unordered_map<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 const 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();
 
     //! Get Incoming Particles
     std::array<Particle, 2> const & get_incoming() const{
       return incoming_;
     }
     //! Set Incoming Particles
     void set_incoming(std::array<Particle, 2> const & in){
       incoming_ = in;
     }
     //! Set Incoming Particles by index
     void set_incoming(Particle const & in, size_t i){
       if(i>1) throw std::out_of_range("Can only set incoming particle 0 or 1.");
       incoming_[i] = in;
     }
 
     //! Get Outgoing Particles
     std::vector<Particle> const & get_outgoing() const{
       return outgoing_;
     }
     //! Set Outgoing Particles
     void set_outgoing(std::vector<Particle> const & out){
       outgoing_ = out;
     }
     //! Replace Outgoing Particles at index i
     void replace_outgoing(Particle const & out, size_t i){
       if(i>=outgoing_.size()) throw std::out_of_range(
         "Can not replace outgoing particle; not set before.");
       outgoing_[i] = out;
     }
     //! Add an Outgoing Particle to the end
     void add_outgoing(Particle const & out){
       outgoing_.push_back(out);
     }
     //! Add an Outgoing Particle to the end (move version)
     void add_outgoing(Particle && out){
       outgoing_.emplace_back(std::move(out));
     }
 
     //! Get Particle decays in the format {outgoing index, decay products}
     std::unordered_map<size_t, std::vector<Particle>> const & get_decays() const{
       return decays_;
     }
     //! Set Particle decays in the format {outgoing index, decay products}
     void set_decays(
       std::unordered_map<size_t, std::vector<Particle>> const & decays
     ){
       decays_ = decays;
     }
     //! Add a Particle decay
     //! @note this will replace previously set decays if out_idx is already used
     void add_decay(
       size_t const out_idx, std::vector<Particle> const & products
     ){
       decays_[out_idx] = products;
     }
 
     //! Get all parameters, i.e. scale choices
     Parameters<EventParameters> const & get_parameters() const{
       return parameters_;
     }
     //! Set all parameters, i.e. scale choices
     void set_parameters(Parameters<EventParameters> const & parameters){
       parameters_ = parameters;
     }
 
     //! Get central parameter choice
     EventParameters const & get_central() const{
       return parameters_.central;
     }
     //! Set central parameter choice
     void set_central(EventParameters const & central){
       parameters_.central = central;
     }
 
     //! Get parameter (scale) variations
     std::vector<EventParameters> const & get_variations() const{
       return parameters_.variations;
     }
     //! Set parameter (scale) variations
     void set_variations(std::vector<EventParameters> const & variations){
       parameters_.variations = variations;
     }
     //! Replace parameter (scale) variation at index i
     void replace_variation(EventParameters const & variation, size_t i){
       if(i>=parameters_.variations.size()) throw std::out_of_range(
         "Can not replace variation; not set before.");
       parameters_.variations[i] = variation;
     }
     //! Add a parameter (scale) variation to the end
     void add_variation(EventParameters const & variation){
       parameters_.variations.push_back(variation);
     }
     //! Add a parameter (scale) variation to the end (move version)
     void add_variation(EventParameters && variation){
       parameters_.variations.push_back(variation);
     }
 
   private:
     std::array<Particle, 2> incoming_;
     std::vector<Particle> outgoing_;
     std::unordered_map<size_t, std::vector<Particle>> decays_;
     Parameters<EventParameters> parameters_;
   }; // end class EventData
 
   //! 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 *);
 
   // 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.3.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<size_t, std::vector<Particle>> decays;
     //! Central parameter (e.g. scale) choice
     EventParameters central;
     std::vector<EventParameters> variations;    /**< For parameter variation */
   };
 
 }