diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
index eb99958..cc35e45 100644
--- a/FixedOrderGen/src/EventGenerator.cc
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -1,80 +1,82 @@
 #include "EventGenerator.hh"
 
 #include "Process.hh"
 #include "Beam.hh"
 #include "JetParameters.hh"
 #include "PhaseSpacePoint.hh"
 
 #include "HEJ/Event.hh"
 #include "HEJ/config.hh"
 
 namespace HEJFOG{
   EventGenerator::EventGenerator(
       Process process,
       Beam beam,
       HEJ::ScaleGenerator scale_gen,
       JetParameters jets,
       int pdf_id,
       double subl_change,
       unsigned int subl_channels,
       ParticlesPropMap particles_properties,
       HEJ::HiggsCouplingSettings Higgs_coupling,
       HEJ::RNG & ran
   ):
     pdf_{pdf_id, beam.particles[0], beam.particles[1]},
     ME_{
       [this](double mu){ return pdf_.Halphas(mu); },
       HEJ::MatrixElementConfig{
         false,
         std::move(Higgs_coupling)
       }
     },
     scale_gen_{std::move(scale_gen)},
     process_{std::move(process)},
     jets_{std::move(jets)},
     beam_{std::move(beam)},
     subl_change_{subl_change},
     subl_channels_{subl_channels},
     particles_properties_{std::move(particles_properties)},
     ran_{ran}
   {
   }
 
   HEJ::Event EventGenerator::gen_event(){
     HEJFOG::PhaseSpacePoint psp{
       process_,
       jets_,
       pdf_, beam_.energy,
       subl_change_, subl_channels_,
       particles_properties_,
       ran_
     };
     status_ = psp.status();
     if(status_ != good) return {};
 
     HEJ::Event ev = scale_gen_(
         HEJ::Event{
           to_UnclusteredEvent(std::move(psp)),
           jets_.def, jets_.min_pt
         }
     );
 
+    ev.generate_colour_flow(ran_);
+
     const double shat = HEJ::shat(ev);
     const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
     const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);
 
     // evaluate matrix element on this point
     const auto ME_weights = ME_.tree(ev);
     ev.central().weight *= ME_weights.central/(shat*shat);
     ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
     ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
     for(size_t i = 0; i < ev.variations().size(); ++i){
       auto & var = ev.variations(i);
       var.weight *= ME_weights.variations[i]/(shat*shat);
       var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
       var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
     }
     return ev;
   }
 
 }
diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh
index 5b287a1..67af9db 100644
--- a/include/HEJ/Particle.hh
+++ b/include/HEJ/Particle.hh
@@ -1,144 +1,143 @@
 /**
  * \file Particle.hh
  * \brief Contains the particle struct
  *
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 
 #pragma once
 
 #include <utility>
 
 #include "fastjet/PseudoJet.hh"
 
 #include "HEJ/optional.hh"
 #include "HEJ/PDG_codes.hh"
 
 namespace HEJ {
 
   using Colour = std::pair<int,int>;
 
   //! Class representing a particle
   struct Particle {
     //! particle type
     ParticleID type;
     //! particle momentum
     fastjet::PseudoJet p;
     //! (optional) colour & anti-colour
-    //! @TODO resolve compiler warnings of "missing initializer for member"
     optional<Colour> colour;
 
     //! get rapidity
     double rapidity() const{
       return p.rapidity();
     }
     //! get transverse momentum
     double perp() const{
       return p.perp();
     }
     //! get momentum in x direction
     double px() const{
       return p.px();
     }
     //! get momentum in y direction
     double py() const{
       return p.py();
     }
     //! get momentum in z direction
     double pz() const{
       return p.pz();
     }
     //! get energy
     double E() const{
       return p.E();
     }
     //! get mass
     double m() const{
       return p.m();
     }
   };
 
   //! Functor to compare rapidities
   /**
    *  This can be used whenever a rapidity comparison function is needed,
    *  for example in many standard library functions.
    *
    *  @see pz_less
    */
   struct rapidity_less{
     template<class FourVector>
     bool operator()(FourVector const & p1, FourVector const & p2){
       return p1.rapidity() < p2.rapidity();
     }
   };
 
   //! Functor to compare momenta in z direction
   /**
    *  This can be used whenever a pz comparison function is needed,
    *  for example in many standard library functions.
    *
    *  @see rapidity_less
    */
   struct pz_less{
     template<class FourVector>
     bool operator()(FourVector const & p1, FourVector const & p2){
       return p1.pz() < p2.pz();
     }
   };
 
 
   //! Convert a vector of Particles to a vector of particle momenta
   inline
   std::vector<fastjet::PseudoJet> to_PseudoJet(
       std::vector<Particle> const & v
   ){
     std::vector<fastjet::PseudoJet> result;
     for(auto && sp: v) result.emplace_back(sp.p);
     return result;
   }
 
   //! Check if a particle is a parton, i.e. quark, antiquark, or gluon
   inline
   bool is_parton(Particle const & p){
     return is_parton(p.type);
   }
 
   //! Check if a particle is a quark
   inline
   bool is_quark(Particle const & p){
     return is_quark(p.type);
   }
 
   //! Check if a particle is an anti-quark
   inline
   bool is_antiquark(Particle const & p){
     return is_antiquark(p.type);
   }
 
   //! Check if a particle is a quark or anit-quark
   inline
   bool is_anyquark(Particle const & p){
     return is_anyquark(p.type);
   }
 
   //! Check if a particle is a photon, W, Z, or Higgs boson
   inline bool is_AWZH_boson(Particle const & particle){
     return is_AWZH_boson(particle.type);
   }
 
   //! Extract all partons from a vector of particles
   inline
   std::vector<Particle> filter_partons(
       std::vector<Particle> const & v
   ){
     std::vector<Particle> result;
     result.reserve(v.size());
     std::copy_if(
         begin(v), end(v), std::back_inserter(result),
         [](Particle const & p){ return is_parton(p); }
     );
     return result;
   }
 }
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 88ac45d..a8e6b00 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,275 +1,276 @@
 /**
  *  \authors   Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 #include "HEJ/EventReweighter.hh"
 
 #include <algorithm>
 #include <assert.h>
 #include <limits>
 #include <math.h>
 #include <stddef.h>
 #include <string>
 #include <unordered_map>
 
 #include "fastjet/ClusterSequence.hh"
 
 #include "LHEF/LHEF.h"
 
 #include "HEJ/Event.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PDG_codes.hh"
 #include "HEJ/PhaseSpacePoint.hh"
 #include "HEJ/Weights.hh"
 
 namespace HEJ{
   using EventType = event_type::EventType;
 
   namespace {
 
     static_assert(
         std::numeric_limits<double>::has_quiet_NaN,
         "no quiet NaN for double"
     );
     constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
 
     UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
       UnclusteredEvent result;
       result.incoming = psp.incoming();
       std::sort(
           begin(result.incoming), end(result.incoming),
           [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
       );
       assert(result.incoming.size() == 2);
       result.outgoing = psp.outgoing();
       assert(
           std::is_sorted(
               begin(result.outgoing), end(result.outgoing),
               rapidity_less{}
           )
       );
       assert(result.outgoing.size() >= 2);
       result.decays = psp.decays();
       result.central.mur = NaN;
       result.central.muf = NaN;
       result.central.weight = psp.weight();
       return result;
     }
 
   } // namespace anonymous
 
   EventReweighter::EventReweighter(
       LHEF::HEPRUP const & heprup,
       ScaleGenerator scale_gen,
       EventReweighterConfig conf,
       HEJ::RNG & ran
   ):
     EventReweighter{
       HEJ::Beam{
         heprup.EBMUP.first,
         {{
           static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
           static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
         }}
       },
       heprup.PDFSUP.first,
       std::move(scale_gen),
       std::move(conf),
       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 beam,
       int pdf_id,
       ScaleGenerator scale_gen,
       EventReweighterConfig conf,
       HEJ::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_{ran}
   {}
 
   PDF const & EventReweighter::pdf() const{
     return pdf_;
   }
 
   std::vector<Event> EventReweighter::reweight(
       Event const & input_ev, int 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_(event);
     return rescale(input_ev, std::move(res_events));
   }
 
   std::vector<Event> EventReweighter::gen_res_events(
       Event const & ev,
       int phase_space_points
   ){
     assert(ev.variations().empty());
 
     switch(param_.treat.at(ev.type())){
     case EventTreatment::discard: return {};
     case EventTreatment::keep:
       if(! jets_pass_resummation_cuts(ev)) return {};
       else return {ev};
     default:;
     }
     const double Born_shat = shat(ev);
 
     std::vector<Event> resummation_events;
     for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
       PhaseSpacePoint psp{ev, param_.psp_config, ran_};
       if(psp.weight() == 0.) continue;
       if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) continue;
 
       resummation_events.emplace_back(
           to_UnclusteredEvent(std::move(psp)),
           param_.jet_param.def, param_.jet_param.min_pt
       );
       auto & new_event = resummation_events.back();
+      new_event.generate_colour_flow(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.central().weight *= pdf.central*ME.central/(Born_pdf*Born_ME);
       for(size_t i = 0; i < cur_event.variations().size(); ++i){
         cur_event.variations(i).weight *=
           pdf.variations[i]*ME.variations[i]/(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 {
         switch(part.type){
         case pid::Higgs: {
           alpha_s_power += 2;
           break;
         }
         // TODO
         case pid::Wp:
         case pid::Wm:
         case pid::photon:
         case pid::Z:
         default:
           throw not_implemented("Emission of boson of unsupported type");
         }
       }
     }
 
     Weights result;
     result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
     for(auto const & var: ev.variations()){
       result.variations.emplace_back(
           pow(pdf_.Halphas(var.mur), alpha_s_power)
       );
     }
     return result;
   }
 
 } // namespace HEJ