diff --git a/t/test_colours.cc b/t/test_colours.cc
index d7804ab..5582a2e 100644
--- a/t/test_colours.cc
+++ b/t/test_colours.cc
@@ -1,362 +1,363 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2022
  *  \copyright GPLv2 or later
  */
 #include "hej_test.hh"
 
 #include <cstdlib>
 #include <iostream>
 #include <utility>
 #include <vector>
 
 #include "HEJ/Constants.hh"
 #include "HEJ/Event.hh"
 #include "HEJ/event_types.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PDG_codes.hh"
 #include "HEJ/RNG.hh"
 
 #include "fastjet/JetDefinition.hh"
 
 namespace {
   /// biased RNG to connect always to colour
   class dum_rnd: public HEJ::RNG {
   public:
     dum_rnd() = default;
     double flat() override {
       return 0.;
     }
   };
 
   HEJ::Event::EventData decay_boson( HEJ::Event::EventData ev ){
     for( std::size_t i=0; i<ev.outgoing.size(); ++i ){
       if( std::abs(ev.outgoing[i].type) == HEJ::ParticleID::Wp){
         ev.decays[i] = decay_W(ev.outgoing[i]);
       }
     }
     return ev;
   }
 
 
   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()){
       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_parton(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 >= HEJ::COLOUR_OFFSET
         && anti_colour >= HEJ::COLOUR_OFFSET;
     if(HEJ::is_quark(part))
       return anti_colour == 0 && colour >= HEJ::COLOUR_OFFSET;
     return colour == 0 && anti_colour >= HEJ::COLOUR_OFFSET;
   }
 
   bool correct_colour(HEJ::Event const & ev){
     if(!ev.is_leading_colour())
       return false;
     // some of these additional checks are also in ev.is_leading_colour()
     for(auto const & part: ev.incoming()){
       if(!correct_colour(part))
         return false;
     }
     for(auto const & part: ev.outgoing()){
       if(!correct_colour(part))
         return false;
     }
     return true;
   }
 
   bool match_expected(
     HEJ::Event const & ev,
     std::vector<HEJ::Colour> const & expected
   ){
     ASSERT(ev.outgoing().size()+2==expected.size());
     for(std::size_t i=0; i<ev.incoming().size(); ++i){
       ASSERT(ev.incoming()[i].colour);
       if( *ev.incoming()[i].colour != expected[i])
         return false;
     }
     for(std::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::Event::EventData unc_ev,
     std::vector<HEJ::Colour> const & expected_colours
 
   ){
+    repair_momenta(unc_ev);
     unc_ev = decay_boson(std::move(unc_ev));
     shuffle_particles(unc_ev); // make sure incoming order doesn't matter
     HEJ::Event ev{unc_ev.cluster(
       fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.)
     };
     ASSERT(HEJ::event_type::is_resummable(ev.type()));
     dum_rnd rng;
     ASSERT(!ev.is_leading_colour());
     ASSERT(ev.generate_colours(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");
     }
   }
 
   HEJ::Event::EventData reset_colour(
     HEJ::Event::EventData ev, std::vector<HEJ::Colour> const & goal
   ){
     for(std::size_t i=0; i<2; ++i){
       ev.incoming[i].colour = goal[i];
     }
     for(std::size_t i=0; i<ev.outgoing.size(); ++i){
       auto const & col_goal{ goal[i+2] };
       if(col_goal.first == 0 && col_goal.second == 0)
         ev.outgoing[i].colour = std::optional<HEJ::Colour>{};
       else
         ev.outgoing[i].colour = col_goal;
     }
     return ev;
   }
 } // namespace
 
 int main() {
   HEJ::Event::EventData ev;
   std::vector<HEJ::Colour> expected_colours(7);
 
   /// pure gluon (they all have a mass of 4 GeV to allow decays)
   ev.incoming[0] =      { HEJ::ParticleID::gluon, {  0,   0, -205, 205}, {}};
   ev.incoming[1] =      { HEJ::ParticleID::gluon, {  0,   0,  279, 279}, {}};
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-15, -82,  -82, 117}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 68,  93,   20, 117}, {}});
-  ev.outgoing.push_back({ HEJ::ParticleID::higgs, {-30, -65,   22,  75}, {}});
+  ev.outgoing.push_back({ HEJ::ParticleID::higgs, {-30, -65,   42,  75}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-12,  92,   76, 120}, {}});
   ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-11, -38,   38,  55}, {}});
 
   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};
   // default colours is always forbidden!
   // default: swap last two (anti-)colour -> crossing
   ev=reset_colour(ev, expected_colours);
   std::swap(ev.outgoing[4].colour, ev.outgoing[3].colour);
   check_event(ev, expected_colours);
 
   /// last g to Qbar (=> gQbar -> g ... Qbar)
   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;
   // default: swap last two anti-colours -> last gluon colour singlet
   ev=reset_colour(ev, expected_colours);
   std::swap(ev.outgoing[4].colour->second, ev.outgoing[3].colour->second);
   check_event(ev, expected_colours);
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno forward (=> gQbar -> g ... Qbar 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
     // default: swap last two anti-colours -> crossing
     new_ev=reset_colour(new_ev, new_expected);
     std::swap(new_ev.outgoing[4].colour->second, new_ev.outgoing[3].colour->second);
     check_event(new_ev, new_expected);
   }
 
   /// swap Qbar <-> Q (=> gQ -> g ... Q)
   ev.incoming[1].type = HEJ::ParticleID::d;
   ev.outgoing[4].type = HEJ::ParticleID::d;
   // => swap: colour<->anti && initial<->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);
   // default: swap incoming <-> outgoing
   ev=reset_colour(ev, expected_colours);
   std::swap(ev.incoming[0].colour, ev.outgoing[0].colour);
   check_event(ev, expected_colours);
 
   /// first g to qbar (=> qbarQ -> qbar ... 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};
   // default: closed qbar->qbar g
   ev=reset_colour(ev, expected_colours);
   ev.outgoing[1].colour->first = ev.outgoing[0].colour->second;
   ev.outgoing[1].colour->second = ev.incoming[0].colour->second;
   ev.outgoing[4].colour->first = ev.outgoing[3].colour->second;
   check_event(ev, expected_colours);
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno backward (=> qbarQ -> g qbar ... 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};
     // default: Colourful Higgs
     new_ev=reset_colour(new_ev, new_expected);
     new_ev.outgoing[2].colour = std::make_pair(1,1);
     check_event(new_ev, new_expected);
 
     /// swap qbar <-> 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;
     // default: swap 2 quarks -> disconnected
     new_ev=reset_colour(new_ev, new_expected);
     std::swap(new_ev.outgoing[1].colour, new_ev.outgoing[4].colour);
     check_event(new_ev, new_expected);
   }
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// uno forward (=> qbarQ -> qbar ... 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};
     // default: no colour on last gluon
     new_ev=reset_colour(new_ev, new_expected);
     new_ev.incoming[1].colour->first = new_ev.outgoing[4].colour->second;
     new_ev.outgoing[4].colour = {};
     check_event(new_ev, new_expected);
   }
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// qqbar backward (=> gQ -> qbar q ... Q) with Wp
     // => swap: incoming q <-> outgoing gluon
     std::swap(new_ev.incoming[0].type, new_ev.outgoing[1].type);
     new_ev.outgoing[1].type=static_cast<HEJ::ParticleID>(
       -(new_ev.outgoing[1].type+1) );
     new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
     // incoming q -> outgoing q (colour<->anti)
     std::swap(new_expected[0], new_expected[3]);
     std::swap(new_expected[3].first, new_expected[3].second);
     new_expected[3].first+=2;
     new_expected[0].first-=1; // skip one index
     // couple first in to first out
     new_expected[2].second=new_expected[0].second;
     // default: swap qqbar <-> first g
     new_ev=reset_colour(new_ev, new_expected);
     std::swap(new_ev.outgoing[0].colour->second, new_ev.outgoing[3].colour->second);
     std::swap(new_ev.outgoing[1].colour->first, new_ev.outgoing[3].colour->first);
     check_event(new_ev, new_expected);
   }
 
   {
     // don't overwrite
     auto new_expected = expected_colours;
     auto new_ev = ev;
     /// qqbar forward (=> qbar g -> qbar ... Qbar Q) with Wp
     // => swap: incoming Q <-> outgoing gluon
     std::swap(new_ev.incoming[1].type, new_ev.outgoing[3].type);
     new_ev.outgoing[3].type=static_cast<HEJ::ParticleID>(
       -(new_ev.outgoing[3].type+1));
     new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
     // incoming q -> outgoing q (colour<->anti)
     std::swap(new_expected[1], new_expected[5]);
     std::swap(new_expected[5].first, new_expected[5].second);
     new_expected[5].second-=2;
     new_expected[1].second-=1; // skip one index
     // couple last in to last out
     new_expected[6].first=new_expected[1].first;
     // default: uncoloured quark
     new_ev=reset_colour(new_ev, new_expected);
     new_ev.outgoing[0].colour = {};
     check_event(new_ev, new_expected);
 
     // move Higgs to position 1 (=> qbar g -> qbar h g Qbar Q)
     std::swap(new_ev.outgoing[1].type, new_ev.outgoing[2].type);
     std::swap(new_expected[3], new_expected[4]); // trivial
     // default: incoming qbar wrong colour
     new_ev=reset_colour(new_ev, new_expected);
     new_ev.incoming[0].colour->first  = 1;
     check_event(new_ev, new_expected);
 
     // central qqbar (=> qbar g -> qbar h Q Qbar g)
     // => swap: Q <-> g
     std::swap(new_ev.outgoing[2].type, new_ev.outgoing[4].type);
     std::swap(new_expected[4], new_expected[6]);
     // gluon was connected on left side, i.e. doesn't matter for QQbar
     // => couple Q to out qbar
     new_expected[4].first = new_expected[2].second;
     // Qbar next in line
     new_expected[5].second = new_expected[4].first+2;
     // incoming g shifted by one position in line
     new_expected[1].first-=2;
     new_expected[1].second+=2;
     // default: wrong colour in last incoming
     new_ev=reset_colour(new_ev, new_expected);
     std::swap(new_ev.incoming[1].colour->first,
               new_ev.incoming[1].colour->second);
     check_event(new_ev, new_expected);
   }
 
   return EXIT_SUCCESS;
 }