diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh
index 7f529eb..c00890d 100644
--- a/include/HEJ/Particle.hh
+++ b/include/HEJ/Particle.hh
@@ -1,118 +1,125 @@
 /**
  * \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
+    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 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;
   }
 }