diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 2e89f8c..fe82359 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,353 +1,393 @@
 /** \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 <memory>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
 #include "HEJ/event_types.hh"
 #include "HEJ/Particle.hh"
 
 #include "fastjet/ClusterSequence.hh"
 
 namespace LHEF{
   class HEPEUP;
   class HEPRUP;
 }
 
 namespace fastjet{
   class JetDefinition;
 }
 
 namespace HEJ{
 
 
   struct ParameterDescription;
 
   //! Event parameters
   struct EventParameters{
     double mur;              /**< Value of the Renormalisation Scale */
     double muf;              /**< Value of the Factorisation Scale */
     double weight;           /**< Event Weight */
     //! Optional description
     std::shared_ptr<ParameterDescription> description = nullptr;
   };
 
   //! Description of event parameters
   struct ParameterDescription {
     //! Name of central scale choice (e.g. "H_T/2")
     std::string scale_name;
     //! Actual renormalisation scale divided by central scale
     double mur_factor;
     //! Actual factorisation scale divided by central scale
     double muf_factor;
 
     ParameterDescription() = default;
     ParameterDescription(
-        std::string scale_name, double mur_factor, double muf_factor
+      std::string scale_name, double mur_factor, double muf_factor
     ):
       scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor}
     {};
   };
 
   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_;
     }
 
     //! Central parameter choice (const version)
     EventParameters const & central() const{
       return central_;
     }
     //! Central parameter choice
     EventParameters & central(){
       return central_;
     }
 
     //! Parameter (scale) variations (const version)
     std::vector<EventParameters> const & variations() const{
       return variations_;
     }
     //! Parameter (scale) variations
     std::vector<EventParameters> & variations(){
       return variations_;
     }
 
     //! Parameter (scale) variation (const version)
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters const & variations(size_t i) const{
       return variations_[i];
     }
     //! Parameter (scale) variation
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters & variations(size_t i){
       return 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,
-        EventParameters && central,
-        std::vector<EventParameters> && variations,
-        fastjet::JetDefinition const & jet_def,
-        double const min_jet_pt
+      std::array<Particle, 2> && incoming,
+      std::vector<Particle> && outgoing,
+      std::unordered_map<size_t, std::vector<Particle>> && decays,
+      EventParameters && central,
+      std::vector<EventParameters> && variations,
+      fastjet::JetDefinition const & jet_def,
+      double const min_jet_pt
     ): incoming_{incoming},
       outgoing_{outgoing},
       decays_{decays},
       central_{central},
       variations_{variations},
       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_;
     //! @TODO replace this by "ParameterVariations"
     EventParameters central_;
     //! @TODO replace this by "ParameterVariations"
     std::vector<EventParameters> variations_;
     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,
       EventParameters const & central,
       std::vector<EventParameters> const & variations
     ):
       incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
       decays_(std::move(decays)),
       central_(std::move(central)), variations_(std::move(variations))
     {};
     //! 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,
       EventParameters && central,
       std::vector<EventParameters> && variations
     ):
       incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
       decays_(std::move(decays)),
       central_(std::move(central)), variations_(std::move(variations))
     {};
 
     //! 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_.push_back(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 Central parameter (e.g. scale) choice
     EventParameters const & get_central() const{
       return central_;
     }
     //! Set Central parameter (e.g. scale) choice
     void set_central(EventParameters const & central){
       central_ = central;
     }
 
     //! Get parameter variation
     std::vector<EventParameters> const & get_variations() const{
       return variations_;
     }
     //! Set parameter variation
     void set_variations(std::vector<EventParameters> const & variations){
       variations_ = variations;
     }
+    //! Replace parameter variation at index i
+    void replace_variation(EventParameters const & variation, size_t i){
+      if(i>=outgoing_.size()) throw std::out_of_range(
+        "Can not replace variation; not set before.");
+      variations_[i] = variation;
+    }
+    //! Add a parameter variation to the end
+    void add_variation(EventParameters const & variation){
+      variations_.push_back(variation);
+    }
+    //! Add a parameter variation to the end (move version)
+    void add_variation(EventParameters && variation){
+      variations_.push_back(variation);
+    }
 
   private:
     std::array<Particle, 2> incoming_;
     std::vector<Particle> outgoing_;
     std::unordered_map<size_t, std::vector<Particle>> decays_;
     //! @TODO replace this by "ParameterVariations"
     EventParameters central_;
     //! @TODO replace this by "ParameterVariations"
     std::vector<EventParameters> variations_;
   }; // 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 */
   };
 
 }
diff --git a/t/test_descriptions.cc b/t/test_descriptions.cc
index 186027b..f84960b 100644
--- a/t/test_descriptions.cc
+++ b/t/test_descriptions.cc
@@ -1,61 +1,60 @@
 #include <iostream>
 #include <cstddef>
 
 #include "HEJ/Event.hh"
 #include "HEJ/EventReweighter.hh"
 #include "HEJ/ScaleFunction.hh"
 
 #define ASSERT(x) if(!(x)) { \
     std::cerr << "Assertion '" #x "' failed.\n"; \
     return EXIT_FAILURE; \
   }
 
 int main() {
   constexpr double mu = 125.;
   HEJ::ScaleFunction fun{"125", HEJ::FixedScale{mu}};
   ASSERT(fun.name() == "125");
 
   HEJ::ScaleGenerator scale_gen{
     {std::move(fun)}, {0.5, 1, 2.}, 2.1
   };
   HEJ::Event::EventData tmp;
-  tmp.set_outgoing({
-      {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)},
-      {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)}
-    }
-  );
+  tmp.add_outgoing(
+      {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)});
+  tmp.add_outgoing(
+      {HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)});
   HEJ::Event ev{
     tmp.cluster(
       fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20.
     )
   };
 
   auto rescaled = scale_gen(std::move(ev));
   ASSERT(rescaled.central().description->scale_name == "125");
   for(auto const & var: rescaled.variations()) {
     ASSERT(var.description->scale_name == "125");
   }
   ASSERT(rescaled.central().description->mur_factor == 1.);
   ASSERT(rescaled.central().description->muf_factor == 1.);
 
   ASSERT(rescaled.variations(0).description->mur_factor == 1.);
   ASSERT(rescaled.variations(0).description->muf_factor == 1.);
 
   ASSERT(rescaled.variations(1).description->mur_factor == 0.5);
   ASSERT(rescaled.variations(1).description->muf_factor == 0.5);
 
   ASSERT(rescaled.variations(2).description->mur_factor == 0.5);
   ASSERT(rescaled.variations(2).description->muf_factor == 1.);
 
   ASSERT(rescaled.variations(3).description->mur_factor == 1.);
   ASSERT(rescaled.variations(3).description->muf_factor == 0.5);
 
   ASSERT(rescaled.variations(4).description->mur_factor == 1.);
   ASSERT(rescaled.variations(4).description->muf_factor == 2.);
 
   ASSERT(rescaled.variations(5).description->mur_factor == 2.);
   ASSERT(rescaled.variations(5).description->muf_factor == 1.);
 
   ASSERT(rescaled.variations(6).description->mur_factor == 2.);
   ASSERT(rescaled.variations(6).description->muf_factor == 2.);
 }