diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index a110609..a037a82 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,203 +1,210 @@
 /** \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 "HEJ/RNG.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
     ):
       scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor}
     {};
   };
 
   //! An event before jet clustering & classification
   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 */
-    //! @brief Generate the particle colour leading in the MRK limit, see \cite Andersen:2011zd
-    //! @note This will overwrite all existing colours
-    void generate_colour_flow(HEJ::RNG &);
-    void sort(); /**< Sort Particles in rapidity **/
   };
 
   /** 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:
     //! Default Event Constructor
     Event() = default;
     //! Event Constructor adding jet clustering to an unclustered event
     Event(
         UnclusteredEvent 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
     UnclusteredEvent const & unclustered() const {
       return ev_;
     }
 
     //! Central parameter choice
     EventParameters const & central() const{
       return ev_.central;
     }
 
     //! Central parameter choice
     EventParameters & central(){
       return ev_.central;
     }
 
     //! Incoming particles
     std::array<Particle, 2> const &  incoming() const{
       return ev_.incoming;
     }
 
     //! Outgoing particles
     std::vector<Particle> const &  outgoing() const{
       return ev_.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 ev_.decays;
     }
 
     //! Parameter (scale) variations
     std::vector<EventParameters> const & variations() const{
       return ev_.variations;
     }
 
     //! Parameter (scale) variations
     std::vector<EventParameters> & variations(){
       return ev_.variations;
     }
 
     //! Parameter (scale) variation
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters const & variations(size_t i) const{
       return ev_.variations[i];
     }
 
     //! Parameter (scale) variation
     /**
      *  @param i   Index of the requested variation
      */
     EventParameters & variations(size_t i){
       return ev_.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_;
     }
 
+    //! Assign colours to each particle
+    /**
+     * @details Colour ordering is done according to leading colour in the MRK
+     *          limit, see \cite Andersen:2011zd. This only affects \ref
+     *          is_HEJ() "HEJ" configurations, all other \ref event_type
+     *          "EventTypes" will be ignored.
+     * @note    This overwrites all previous set colours.
+     */
+    void generate_colour_flow(HEJ::RNG &);
+
   private:
+    void sort(); /**< Sort Particles in rapidity **/
     UnclusteredEvent ev_;
     fastjet::ClusterSequence cs_;
     double min_jet_pt_;
     event_type::EventType type_;
   };
 
   //! 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 *);
 
 }
diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
index 1b145a8..df60728 100644
--- a/include/HEJ/event_types.hh
+++ b/include/HEJ/event_types.hh
@@ -1,63 +1,65 @@
 /** \file
  *  \brief Define different types of events.
  *
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 
 #pragma once
 
 #include "HEJ/utility.hh"
 
 namespace HEJ{
 
   //! Namespace for event types
   namespace event_type{
     //! Possible event types
     enum EventType: size_t{
       FKL,                        /**< FKL-type event */
       unordered_backward,         /**< event with unordered backward emission */
       unordered_forward,          /**< event with unordered forward emission */
       nonHEJ,                     /**< event configuration not covered by HEJ */
       no_2_jets,                  /**< event with less than two jets */
       bad_final_state,            /**< event with an unsupported final state */
       unob = unordered_backward,
       unof = unordered_forward,
       first_type = FKL,
       last_type = bad_final_state
     };
 
     //! Event type names
     /**
      * For example, names[FKL] is the string "FKL"
      */
     static constexpr auto names = make_array(
       "FKL",
       "unordered backward",
       "unordered forward",
       "nonHEJ",
       "no 2 jets",
       "bad final state"
     );
 
+    //! Returns True for a HEJ \ref event_type::EventType "EventType"
     inline
     bool is_HEJ(EventType type) {
       switch(type) {
       case FKL:
       case unordered_backward:
       case unordered_forward:
         return true;
       default:
         return false;
       }
     }
 
+    //! Returns True for an unordered \ref event_type::EventType "EventType"
     inline
     bool is_uno(EventType type) {
       return type == unordered_backward || type == unordered_forward;
     }
 
   }
 
 }
diff --git a/src/Event.cc b/src/Event.cc
index 7df221a..31245b4 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,486 +1,490 @@
 /**
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 #include "HEJ/Event.hh"
 
 #include <algorithm>
 #include <assert.h>
 #include <numeric>
 #include <utility>
 
 #include "LHEF/LHEF.h"
 
 #include "fastjet/JetDefinition.hh"
 
 #include "HEJ/Constants.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/PDG_codes.hh"
 
 namespace HEJ{
 
   namespace{
     constexpr int status_in = -1;
     constexpr int status_decayed = 2;
     constexpr int status_out = 1;
 
     /// @name helper functions to determine event type
     //@{
 
     /**
      * \brief check if final state valid for HEJ
      *
      * check if there is at most one photon, W, H, Z in the final state
      * and all the rest are quarks or gluons
      */
     bool final_state_ok(std::vector<Particle> const & outgoing){
       bool has_AWZH_boson = false;
       for(auto const & out: outgoing){
         if(is_AWZH_boson(out.type)){
           if(has_AWZH_boson) return false;
           has_AWZH_boson = true;
         }
         else if(! is_parton(out.type)) return false;
       }
       return true;
     }
 
     template<class Iterator>
     Iterator remove_AWZH(Iterator begin, Iterator end){
       return std::remove_if(
           begin, end, [](Particle const & p){return is_AWZH_boson(p);}
       );
     }
 
     template<class Iterator>
     bool valid_outgoing(Iterator begin, Iterator end){
       return std::distance(begin, end) >= 2
         && std::is_sorted(begin, end, rapidity_less{})
         && std::count_if(
             begin, end, [](Particle const & s){return is_AWZH_boson(s);}
         ) < 2;
     }
 
     /// @note that this changes the outgoing range!
     template<class ConstIterator, class Iterator>
     bool is_FKL(
         ConstIterator begin_incoming, ConstIterator end_incoming,
         Iterator begin_outgoing, Iterator end_outgoing
     ){
       assert(std::distance(begin_incoming, end_incoming) == 2);
       assert(std::distance(begin_outgoing, end_outgoing) >= 2);
 
       // One photon, W, H, Z in the final state is allowed.
       // Remove it for remaining tests,
       end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
 
       // Test if this is a standard FKL configuration.
       return
         (begin_incoming->type == begin_outgoing->type)
         && ((end_incoming-1)->type == (end_outgoing-1)->type)
         && std::all_of(
             begin_outgoing + 1, end_outgoing - 1,
             [](Particle const & p){ return p.type == pid::gluon; }
         );
     }
 
     bool is_FKL(
         std::array<Particle, 2> const & incoming,
         std::vector<Particle> outgoing
     ){
       assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
       assert(valid_outgoing(begin(outgoing), end(outgoing)));
 
       return is_FKL(
           begin(incoming), end(incoming),
           begin(outgoing), end(outgoing)
       );
     }
 
     bool has_2_jets(Event const & event){
       return event.jets().size() >= 2;
     }
 
     /**
      * \brief Checks whether event is unordered backwards
      * @param ev        Event
      * @returns         Is Event Unordered Backwards
      *
      * - Checks there is more than 3 constuents in the final state
      * - Checks there is more than 3 jets
      * - Checks the most backwards parton is a gluon
      * - Checks the most forwards jet is not a gluon
      * - Checks the rest of the event is FKL
      * - Checks the second most backwards is not a different boson
      * - Checks the unordered gluon actually forms a jet
      */
     bool is_unordered_backward(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
       if(in.front().type == pid::gluon) return false;
       if(out.front().type != pid::gluon) return false;
       // When skipping the unordered emission
       // the remainder should be a regular FKL event,
       // except that the (new) first outgoing particle must not be a A,W,Z,H.
       const auto FKL_begin = next(begin(out));
       if(is_AWZH_boson(*FKL_begin)) return false;
       if(!is_FKL(in, {FKL_begin, end(out)})) return false;
       // check that the unordered gluon forms an extra jet
       const auto jets = sorted_by_rapidity(ev.jets());
       const auto indices = ev.particle_jet_indices({jets.front()});
       return indices[0] >= 0 && indices[1] == -1;
     }
 
     /**
      * \brief Checks for a forward unordered gluon emission
      * @param ev          Event
      * @returns           Is the event a forward unordered emission
      *
      * \see is_unordered_backward
      */
     bool is_unordered_forward(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
       if(in.back().type == pid::gluon) return false;
       if(out.back().type != pid::gluon) return false;
       // When skipping the unordered emission
       // the remainder should be a regular FKL event,
       // except that the (new) last outgoing particle must not be a A,W,Z,H.
       const auto FKL_end = prev(end(out));
       if(is_AWZH_boson(*prev(FKL_end))) return false;
       if(!is_FKL(in, {begin(out), FKL_end})) return false;
       // check that the unordered gluon forms an extra jet
       const auto jets = sorted_by_rapidity(ev.jets());
       const auto indices = ev.particle_jet_indices({jets.back()});
       return indices.back() >= 0 && indices[indices.size()-2] == -1;
     }
 
     using event_type::EventType;
 
     EventType classify(Event const & ev){
       if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state;
       if(! has_2_jets(ev)) return EventType::no_2_jets;
       if(is_FKL(ev.incoming(), ev.outgoing())) return EventType::FKL;
       if(is_unordered_backward(ev)){
         return EventType::unordered_backward;
       }
       if(is_unordered_forward(ev)){
         return EventType::unordered_forward;
       }
       return EventType::nonHEJ;
     }
     //@}
 
     Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
       const ParticleID id = static_cast<ParticleID>(hepeup.IDUP[i]);
       const fastjet::PseudoJet momentum{
         hepeup.PUP[i][0], hepeup.PUP[i][1],
         hepeup.PUP[i][2], hepeup.PUP[i][3]
       };
       if(is_parton(id))
         return Particle{ id, std::move(momentum), hepeup.ICOLUP[i] };
       return Particle{ id, std::move(momentum), {} };
     }
 
     bool is_decay_product(std::pair<int, int> const & mothers){
       if(mothers.first == 0) return false;
       return mothers.second == 0 || mothers.first == mothers.second;
     }
 
   } // namespace anonymous
 
   UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup):
     central(EventParameters{
       hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
     })
   {
     size_t in_idx = 0;
     for (int i = 0; i < hepeup.NUP; ++i) {
       // skip decay products
       // we will add them later on, but we have to ensure that
       // the decayed particle is added before
       if(is_decay_product(hepeup.MOTHUP[i])) continue;
 
       auto particle = extract_particle(hepeup, i);
       // needed to identify mother particles for decay products
       particle.p.set_user_index(i+1);
 
       if(hepeup.ISTUP[i] == status_in){
          if(in_idx > incoming.size()) {
            throw std::invalid_argument{
              "Event has too many incoming particles"
            };
          }
          incoming[in_idx++] = std::move(particle);
       }
       else outgoing.emplace_back(std::move(particle));
     }
 
     // add decay products
     for (int i = 0; i < hepeup.NUP; ++i) {
       if(!is_decay_product(hepeup.MOTHUP[i])) continue;
       const int mother_id = hepeup.MOTHUP[i].first;
       const auto mother = std::find_if(
           begin(outgoing), end(outgoing),
           [mother_id](Particle const & particle){
             return particle.p.user_index() == mother_id;
           }
       );
       if(mother == end(outgoing)){
         throw std::invalid_argument{"invalid decay product parent"};
       }
       const int mother_idx = std::distance(begin(outgoing), mother);
       assert(mother_idx >= 0);
       decays[mother_idx].emplace_back(extract_particle(hepeup, i));
     }
   }
 
   namespace {
     int connect_incoming(Particle & part, int & colour, int & anti_colour){
-      // TODO something is wrong here
       part.colour = std::make_pair(anti_colour, colour);
       // gluon
       if(part.type == pid::gluon) {
         return 0;
       }
       // q
       if(part.type > 0){
         part.colour->second = 0;
         colour*=-1;
         return -1;
       }
       // qx
       part.colour->first = 0;
       anti_colour*=-1;
       return +1;
     }
   }
 
-  void UnclusteredEvent::generate_colour_flow(RNG & ran){
-    sort();
+  void Event::generate_colour_flow(RNG & ran){
+    // only sort HEJ events
+    if(!event_type::is_HEJ(type()))
+      return;
+
+    assert(std::is_sorted(
+      begin(outgoing()), end(outgoing()), rapidity_less{}));
+    assert(incoming()[0].pz() < incoming()[1].pz());
 
     // initialise first
     int colour = COLOUR_OFFSET;
     int anti_colour = colour+1;
     // -1 anti, 0 no, +1 colour
-    int is_quark_line = connect_incoming(incoming[0], colour, anti_colour);
+    int is_quark_line = connect_incoming(ev_.incoming[0], colour, anti_colour);
 
-    for(auto & part: outgoing){
+    for(auto & part: ev_.outgoing){
       if(part.type == ParticleID::gluon){
         // gluon
         if(is_quark_line == 0){
           // on g line => connect to colour OR anti-colour (random)
           if(ran.flat() < 0.5){
             part.colour = std::make_pair(colour+2,colour);
             colour+=2;
           } else {
             part.colour = std::make_pair(anti_colour, anti_colour+2);
             anti_colour+=2;
           }
         } else if(is_quark_line > 0){
           // on q line => connect to available colour
             part.colour = std::make_pair(colour+2, colour);
             colour+=2;
         } else {
           // on qx line => connect to available anti-colour
           part.colour = std::make_pair(anti_colour, anti_colour+2);
           anti_colour+=2;
         }
       } else if(is_quark(part)) {
         // quark
         if(is_quark_line == 0){
           // on g line => connect and remove colour
           is_quark_line = 1;
           part.colour = std::make_pair(anti_colour, 0);
           anti_colour+=2;
           anti_colour*=-1;
         } else if(is_quark_line < 0){
           // on qx line => new colour
           is_quark_line = 0;
           colour*=-1;
           part.colour = std::make_pair(colour,0);
 
         } else {
           throw std::logic_error("Colour connection impossible: Expected anti-quark but found quark.");
         }
       } else if(is_antiquark(part)) {
         // anti-quark
         if(is_quark_line == 0){
           // on g line => connect and remove anti-colour
           is_quark_line = -1;
           part.colour = std::make_pair(0, colour);
           colour+=2;
           colour*=-1;
         } else if(is_quark_line > 0){
           // on q line => new anti-colour
           is_quark_line = 0;
           anti_colour*=-1;
           part.colour = std::make_pair(0,anti_colour);
         } else {
           throw std::logic_error("Colour connection impossible: Expected quark but found anti-quark.");
         }
       }
     }
     // Connect last
-    connect_incoming(incoming[1], anti_colour, colour);
+    connect_incoming(ev_.incoming[1], anti_colour, colour);
 
   }
 
-  void UnclusteredEvent::sort(){
+  void Event::sort(){
     // sort incoming
     std::sort(
-        begin(incoming), end(incoming),
+        begin(ev_.incoming), end(ev_.incoming),
         [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
     );
 
     //sort outgoing
-    if(std::is_sorted(begin(outgoing), end(outgoing), rapidity_less{}))
+    if(std::is_sorted(begin(ev_.outgoing), end(ev_.outgoing), rapidity_less{}))
       return;
-    auto old_outgoing = std::move(outgoing);
+    auto old_outgoing = std::move(ev_.outgoing);
     std::vector<size_t> idx(old_outgoing.size());
     std::iota(idx.begin(), idx.end(), 0);
     std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
       return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
     });
-    outgoing.clear();
-    outgoing.reserve(old_outgoing.size());
+    ev_.outgoing.clear();
     for(size_t i: idx) {
-      outgoing.emplace_back(std::move(old_outgoing[i]));
+      ev_.outgoing.emplace_back(std::move(old_outgoing[i]));
     }
 
     // find decays again
-    if(!decays.empty()){
-      auto old_decays = std::move(decays);
-      decays.clear();
+    if(!ev_.decays.empty()){
+      auto old_decays = std::move(ev_.decays);
+      ev_.decays.clear();
       for(size_t i=0; i<idx.size(); ++i) {
         auto decay = old_decays.find(idx[i]);
         if(decay != old_decays.end())
-          decays.emplace(i, std::move(decay->second));
+          ev_.decays.emplace(i, std::move(decay->second));
       }
-      assert(old_decays.size() == decays.size());
+      assert(old_decays.size() == ev_.decays.size());
     }
 
   }
 
   Event::Event(
       UnclusteredEvent ev,
       fastjet::JetDefinition const & jet_def, double min_jet_pt
   ):
     ev_{std::move(ev)},
     cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
     min_jet_pt_{min_jet_pt}
   {
     // sort particles
-    ev_.sort();
+    sort();
 
     // classify event
     type_ = classify(*this);
 
     assert(std::is_sorted(begin(outgoing()), end(outgoing()), rapidity_less{}));
   }
 
   std::vector<fastjet::PseudoJet> Event::jets() const{
     return cs_.inclusive_jets(min_jet_pt_);
   }
 
   double shat(Event const & ev){
     return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
   }
 
   namespace{
     // colour flow according to Les Houches standard
     // TODO: stub
     std::vector<std::pair<int, int>> colour_flow(
         std::array<Particle, 2> const & incoming,
         std::vector<Particle> const & outgoing
     ){
       std::vector<std::pair<int, int>> result(
           incoming.size() + outgoing.size()
       );
       for(auto & col: result){
         col = std::make_pair(-1, -1);
       }
       return result;
     }
   }
 
   LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
     LHEF::HEPEUP result;
     result.heprup = heprup;
     result.weights = {{event.central().weight, nullptr}};
     for(auto const & var: event.variations()){
       result.weights.emplace_back(var.weight, nullptr);
     }
     size_t num_particles = event.incoming().size() + event.outgoing().size();
     for(auto const & decay: event.decays()) num_particles += decay.second.size();
     result.NUP = num_particles;
     // the following entries are pretty much meaningless
     result.IDPRUP = event.type()+1;  // event ID
     result.AQEDUP = 1./128.;  // alpha_EW
     //result.AQCDUP = 0.118 // alpha_QCD
     // end meaningless part
     result.XWGTUP = event.central().weight;
     result.SCALUP = event.central().muf;
     result.scales.muf = event.central().muf;
     result.scales.mur = event.central().mur;
     result.scales.SCALUP = event.central().muf;
     result.pdfinfo.p1 = event.incoming().front().type;
     result.pdfinfo.p2 = event.incoming().back().type;
     result.pdfinfo.scale = event.central().muf;
     for(Particle const & in: event.incoming()){
       result.IDUP.emplace_back(in.type);
       result.ISTUP.emplace_back(status_in);
       result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
       result.MOTHUP.emplace_back(0, 0);
     }
     for(size_t i = 0; i < event.outgoing().size(); ++i){
       Particle const & out = event.outgoing()[i];
       result.IDUP.emplace_back(out.type);
       const int status = event.decays().count(i)?status_decayed:status_out;
       result.ISTUP.emplace_back(status);
       result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
       result.MOTHUP.emplace_back(1, 2);
     }
     result.ICOLUP = colour_flow(
         event.incoming(), filter_partons(event.outgoing())
     );
     if(result.ICOLUP.size() < num_particles){
       const size_t AWZH_boson_idx = std::find_if(
           begin(event.outgoing()), end(event.outgoing()),
           [](Particle const & s){ return is_AWZH_boson(s); }
       ) - begin(event.outgoing()) + event.incoming().size();
       assert(AWZH_boson_idx <= result.ICOLUP.size());
       result.ICOLUP.insert(
           begin(result.ICOLUP) + AWZH_boson_idx,
           std::make_pair(0,0)
       );
     }
     for(auto const & decay: event.decays()){
       for(auto const out: decay.second){
         result.IDUP.emplace_back(out.type);
         result.ISTUP.emplace_back(status_out);
         result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
         const size_t mother_idx = 1 + event.incoming().size() + decay.first;
         result.MOTHUP.emplace_back(mother_idx, mother_idx);
         result.ICOLUP.emplace_back(0,0);
       }
     }
     assert(result.ICOLUP.size() == num_particles);
     static constexpr double unknown_spin = 9.;     //per Les Houches accord
     result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
     result.SPINUP = result.VTIMUP;
     return result;
   }
 
 }
diff --git a/t/test_colour_flow.cc b/t/test_colour_flow.cc
index e76cd1f..44156c6 100644
--- a/t/test_colour_flow.cc
+++ b/t/test_colour_flow.cc
@@ -1,205 +1,209 @@
 #include <stdexcept>
 #include <utility>
 
 #include "HEJ/Event.hh"
 #include "HEJ/RNG.hh"
 
 /// biased RNG to connect always to colour
 class dum_rnd: public HEJ::DefaultRNG {
 public:
   dum_rnd() = default;
   double flat() override {
     return 0.;
   };
 };
 
-void dump_event(HEJ::UnclusteredEvent const & ev){
-  for(auto const & in: ev.incoming){
+void dump_event(HEJ::Event const & ev){
+  for(auto const & in: ev.incoming()){
     std::cerr << "in type=" << in.type
       << ", colour={" << (*in.colour).first
       << ", " << (*in.colour).second << "}\n";
   }
-  for(auto const & out: ev.outgoing){
+  for(auto const & out: ev.outgoing()){
     std::cerr << "out type=" << out.type << ", colour={";
     if(out.colour)
       std::cerr << (*out.colour).first << ", " << (*out.colour).second;
     else
       std::cerr << "non, non";
     std::cerr << "}\n";
   }
 }
 
 /// true if colour is allowed for particle
 bool correct_colour(HEJ::Particle const & part){
   if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
   if(!part.colour) return false;
   int const colour = part.colour->first;
   int const anti_colour = part.colour->second;
   if(part.type == HEJ::ParticleID::gluon)
     return colour != anti_colour && colour > 0 && anti_colour > 0;
   if(HEJ::is_quark(part))
     return anti_colour == 0 && colour > 0;
   return colour == 0 && anti_colour > 0;
 }
 
-bool correct_colour(HEJ::UnclusteredEvent const & ev){
-  for(auto const & part: ev.incoming){
+bool correct_colour(HEJ::Event const & ev){
+  for(auto const & part: ev.incoming()){
     if(!correct_colour(part))
       return false;
   }
-  for(auto const & part: ev.outgoing){
+  for(auto const & part: ev.outgoing()){
     if(!correct_colour(part))
       return false;
   }
   return true;
 }
 
 bool match_expected(
-  HEJ::UnclusteredEvent const & ev,
+  HEJ::Event const & ev,
   std::vector<HEJ::Colour> const & expected
 ){
-  assert(ev.outgoing.size()+2==expected.size());
-  for(size_t i=0; i<ev.incoming.size(); ++i){
-    assert(ev.incoming[i].colour);
-    if( *ev.incoming[i].colour != expected[i])
+  assert(ev.outgoing().size()+2==expected.size());
+  for(size_t i=0; i<ev.incoming().size(); ++i){
+    assert(ev.incoming()[i].colour);
+    if( *ev.incoming()[i].colour != expected[i])
       return false;
   }
-  for(size_t i=2; i<ev.outgoing.size()+2; ++i){
-    if( ev.outgoing[i-2].colour ){
-      if( *ev.outgoing[i-2].colour != expected[i] )
+  for(size_t i=2; i<ev.outgoing().size()+2; ++i){
+    if( ev.outgoing()[i-2].colour ){
+      if( *ev.outgoing()[i-2].colour != expected[i] )
         return false;
     } else if( expected[i].first != 0 || expected[i].second != 0)
       return false;
   }
   return true;
 }
 
 void check_event(
-  HEJ::UnclusteredEvent ev, std::vector<HEJ::Colour> const & expected_colours
+  HEJ::UnclusteredEvent unc_ev, std::vector<HEJ::Colour> const & expected_colours
 ){
+  HEJ::Event ev(
+    unc_ev,
+    fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.);
+  assert(HEJ::event_type::is_HEJ(ev.type()));
   dum_rnd rng;
   ev.generate_colour_flow(rng);
   if(!correct_colour(ev)){
     std::cerr << "Found illegal colours for event\n";
     dump_event(ev);
     throw std::invalid_argument("Illegal colour set");
   }
   if(!match_expected(ev, expected_colours)){
     std::cerr << "Colours didn't match expectation. Found\n";
     dump_event(ev);
     std::cerr << "but expected\n";
     for(auto const & col: expected_colours){
       std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
     }
     throw std::logic_error("Colours did not match expectation");
   }
 }
 
 int main() {
   HEJ::UnclusteredEvent ev;
   std::vector<HEJ::Colour> expected_colours(7);
 
   /// pure gluon
   ev.incoming[0] =      { HEJ::ParticleID::gluon, {   0,   0,-427, 427}, {}};
   ev.incoming[1] =      { HEJ::ParticleID::gluon, {   0,   0, 851, 851}, {}};
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 196, 124, -82, 246}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-167,-184,  16, 249}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::higgs, { 197, 180, 168, 339}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-190, -57, 126, 235}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, { -36, -63, 196, 209}, {}});
 
   expected_colours[0] = {502, 501};
   expected_colours[1] = {509, 502};
   expected_colours[2] = {503, 501};
   expected_colours[3] = {505, 503};
   expected_colours[4] = {  0,   0};
   expected_colours[5] = {507, 505};
   expected_colours[6] = {509, 507};
   check_event(ev, expected_colours);
 
   /// last g to Qx (=> gQx -> g ... Qx)
   ev.incoming[1].type = HEJ::ParticleID::d_bar;
   ev.outgoing[4].type = HEJ::ParticleID::d_bar;
   // => only end changes
   expected_colours[1].first = 0;
   expected_colours[6].first = 0;
   check_event(ev, expected_colours);
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno forward (=> gQx -> g ... Qx g)
     std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
     // => uno quarks eats colour and gluon connects to anti-colour
     new_expected[5] = {0, expected_colours[3].first};
     new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
     new_expected[1].second += 2; // one more anti-colour in line
     check_event(new_ev, new_expected);
   }
 
   /// swap Qx <-> Q (=> gQx -> g ... Qx)
   ev.incoming[1].type = HEJ::ParticleID::d;
   ev.outgoing[4].type = HEJ::ParticleID::d;
   // => swap: colour<->anti && inital<->final
   std::swap(expected_colours[1], expected_colours[6]);
   std::swap(expected_colours[1].first, expected_colours[1].second);
   std::swap(expected_colours[6].first, expected_colours[6].second);
   check_event(ev, expected_colours);
 
   /// first g to q (=> qxQ -> qx ... Q)
   ev.incoming[0].type = HEJ::ParticleID::u_bar;
   ev.outgoing[0].type = HEJ::ParticleID::u_bar;
   expected_colours[0] = {  0, 501};
   // => shift anti-colour index one up
   expected_colours[1].first -= 2;
   expected_colours[5] = expected_colours[3];
   expected_colours[3] = expected_colours[2];
   expected_colours[2] = {  0, 502};
   check_event(ev, expected_colours);
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno backward (=> qxQ -> g qx ... Q)
     std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
     // => uno gluon connects to quark colour
     new_expected[3] = expected_colours[2];
     new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
     check_event(new_ev, new_expected);
 
     /// swap qx <-> q (=> qQ -> g q ... Q)
     new_ev.incoming[0].type = HEJ::ParticleID::u;
     new_ev.outgoing[1].type = HEJ::ParticleID::u;
     // => swap: colour<->anti && inital<->final
     std::swap(new_expected[0], new_expected[3]);
     std::swap(new_expected[0].first, new_expected[0].second);
     std::swap(new_expected[3].first, new_expected[3].second);
     // => & connect first gluon with remaining anti-colour
     new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
     // shift colour line one down
     new_expected[1].first-=2;
     new_expected[5].first-=2;
     new_expected[5].second-=2;
     // shift anti-colour line one up
     new_expected[6].first+=2;
     check_event(new_ev, new_expected);
   }
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno forward (=> qxQ -> qx ... Q g)
     std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
     // => uno gluon connects to remaining colour
     new_expected[5] = expected_colours[6];
     new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
     check_event(new_ev, new_expected);
   }
 
   /// @TODO add qqx test when implemented (it should work)
 
   return EXIT_SUCCESS;
 }