diff --git a/t/test_classify.cc b/t/test_classify.cc
index b09c6a0..a46dedc 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,641 +1,641 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "hej_test.hh"
 
 #include <array>
 #include <cstdlib>
 #include <iostream>
 #include <random>
 #include <string>
 #include <vector>
 
 #include "fastjet/JetDefinition.hh"
 
 #include "HEJ/Event.hh"
 #include "HEJ/event_types.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/PDG_codes.hh"
 
 namespace {
   const fastjet::JetDefinition JET_DEF{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
   const double MIN_JET_PT{30.};
 
   const std::vector<std::string> ALL_QUARKS{"-4","-1","1","2","3","4"};
   const std::vector<std::string> ALL_PARTONS{"g","-2","-1","1","2","3","4"};
   const std::vector<std::string> ALL_BOSONS{"h", "Wp", "Wm", "Z_photon_mix"};
   const std::vector<std::string> ALL_G_Z{"photon", "Z"};
   const std::vector<std::string> ALL_W{"W+", "W-"};
   const std::size_t MAX_MULTI = 6;
 
   std::mt19937_64 RAN{0};
 
   //! Upon clustering does EventData match the expected classification
   bool match_expectation(HEJ::event_type::EventType expected,
                          HEJ::Event::EventData evd
   ){
     HEJ::Event ev { evd.cluster(JET_DEF, MIN_JET_PT) };
     if(ev.type() != expected){
       std::cerr << "Expected type " << name(expected)
                 << " but found " << name(ev.type()) << "\n" << ev;
       auto jet_idx{ ev.particle_jet_indices() };
       std::cout << "Particle Jet indices: ";
       for(int const i: jet_idx)
         std::cout << i << " ";
       std::cout << std::endl;
       return false;
     }
     return true;
   }
 
   //! Does match_expectation for EventData from parse_configuration
   bool match_expectation( HEJ::event_type::EventType expected,
     std::array<std::string,2> const & in, std::vector<std::string> const & out,
       int const overwrite_boson = 0
   ){
     return match_expectation(expected, parse_configuration( in,out,overwrite_boson ));
   }
 
   //! test FKL configurations
   //! if implemented==false : check processes that are not in HEJ yet
   bool check_fkl( bool const implemented=true ){
     using namespace HEJ;
     auto const type{ implemented?event_type::FKL:event_type::non_resummable };
     std::vector<std::string> bosons;
     if(implemented)
       bosons = ALL_BOSONS;
     else {
       bosons = ALL_G_Z;
     }
     for(std::string const & first: ALL_PARTONS)       // all quark flavours
       for(std::string const & last: ALL_PARTONS){
         for(std::size_t njet=1; njet<=MAX_MULTI; ++njet){  // all multiplicities
           if(njet==5) continue;
           std::array<std::string,2> base_in{first,last};
           std::vector<std::string> base_out(njet, "g");
           if(njet>1){ // pure jets only for >=2j
             base_out.front() = first;
             base_out.back() = last;
             if(implemented && !match_expectation(type, base_in, base_out))
               return false;
           }
           for(auto const & boson: bosons){         // any boson
             for(std::size_t pos=0; pos<=njet; ++pos){    // at any position
               if(njet == 1){
                 // special case: one jet only for g->h
                 if(boson != "h" || ((pos==0?first:last) != "g")){
                   continue;
                 } else {
                   base_out.front() = pos==0?last:first;
                 }
               }
               auto const & in{base_in};
               auto out{base_out};
               // change quark flavours for W
               const bool couple_idx
                 = std::uniform_int_distribution<int>{0,1}(RAN) != 0;
               if(!couple_quark(boson, couple_idx?out.back():out.front()))
                 continue;
               out.insert(out.begin()+pos, boson);
               if(!match_expectation(type, in, out))
                 return false;
             }
           }
         }
       }
     return true;
   }
 
   //! test unordered configurations
   //! if implemented==false : check processes that are not in HEJ yet
   bool check_uno( bool const implemented=true ){
     using namespace HEJ;
     auto const b{ implemented?event_type::unob:event_type::non_resummable };
     auto const f{ implemented?event_type::unof:event_type::non_resummable };
     std::vector<std::string> bosons;
     if(implemented) {
       bosons = ALL_BOSONS;
     } else {
       bosons = ALL_G_Z;
     }
     for(std::string const & uno: ALL_QUARKS)          // all quark flavours
       for(std::string const & fkl: ALL_PARTONS){
         for(std::size_t njet=3; njet<=MAX_MULTI; ++njet){  // all multiplicities >2
           if(njet==5) continue;
           for(std::size_t i=0; i<2; ++i){                 // forward & backwards
             std::array<std::string,2> base_in;
             std::vector<std::string> base_out(njet, "g");
             const std::size_t uno_pos = i?1:(njet-2);
             const std::size_t fkl_pos = i?(njet-1):0;
             base_in[i?0:1] = uno;
             base_in[i?1:0] = fkl;
             base_out[uno_pos] = uno;
             base_out[fkl_pos] = fkl;
             auto expectation{ i?b:f };
             if( implemented
                 && !match_expectation(expectation, base_in, base_out) )
               return false;
             for(auto const & boson: bosons){      // any boson
               // at any position (higgs only inside FKL chain)
               std::size_t start = 0;
               std::size_t end = njet;
               if(to_ParticleID(boson) == pid::higgs){
                 start = i?(uno_pos+1):fkl_pos;
                 end = i?(fkl_pos+1):uno_pos;
               }
               for(std::size_t pos=start; pos<=end; ++pos){
                 auto const & in{base_in};
                 auto out{base_out};
                 // change quark flavours for W
                 const bool couple_idx
                   = std::uniform_int_distribution<int>{0,1}(RAN) != 0;
                 if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos]))
                   continue;
                 out.insert(out.begin()+pos, boson);
                 if(!match_expectation(expectation, in, out))
                   return false;
               }
             }
           }
         }
       }
     return true;
   }
 
   //! test extremal qqbar configurations
   //! if implemented==false : check processes that are not in HEJ yet
   bool check_extremal_qqbar( bool const implemented=true ){
     using namespace HEJ;
     auto const b{ implemented?event_type::qqbar_exb:event_type::non_resummable };
     auto const f{ implemented?event_type::qqbar_exf:event_type::non_resummable };
     std::vector<std::string> bosons;
     if(implemented)
       bosons = ALL_W;
     else {
       bosons = ALL_G_Z;
       bosons.emplace_back("h");
       bosons.emplace_back("Z_photon_mix");
     }
     for(std::string const & qqbar: ALL_QUARKS)          // all quark flavours
       for(std::string const & fkl: ALL_PARTONS){
         std::string const qqbar2{ std::to_string(HEJ::to_ParticleID(qqbar)*-1) };
         for(std::size_t njet=3; njet<=MAX_MULTI; ++njet){  // all multiplicities >2
           if(njet==5) continue;
           for(std::size_t i=0; i<2; ++i){             // forward & backwards
             std::array<std::string,2> base_in;
             std::vector<std::string> base_out(njet, "g");
             const std::size_t qqbar_pos = i?0:(njet-2);
             const std::size_t fkl_pos = i?(njet-1):0;
             base_in[i?0:1] = "g";
             base_in[i?1:0] = fkl;
             base_out[fkl_pos]   = fkl;
             base_out[qqbar_pos]   = qqbar;
             base_out[qqbar_pos+1] = qqbar2;
             auto expectation{ i?b:f };
             if( implemented
                 && !match_expectation(expectation, base_in, base_out) )
               return false;
             for(auto const & boson: bosons){ // all bosons
               // at any position (higgs only inside FKL chain)
               std::size_t start = 0;
               std::size_t end = njet;
               if(to_ParticleID(boson) == pid::higgs){
                 start = i?(qqbar_pos+2):fkl_pos;
                 end = i?(fkl_pos+1):qqbar_pos;
               }
               for(std::size_t pos=start; pos<=end; ++pos){
                 auto const & in{base_in};
                 auto out{base_out};
                 // change quark flavours for W
                 const bool couple_idx
                   = std::uniform_int_distribution<int>{0,1}(RAN) != 0;
                 if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){
                   // (randomly) try couple to FKL, else fall-back to qqbar
                   if(!couple_quark(boson, out[qqbar_pos]))
                     couple_quark(boson, out[qqbar_pos+1]);
                 }
                 out.insert(out.begin()+pos, boson);
                 if(!match_expectation(expectation, in, out))
                   return false;
               }
             }
           }
         }
         // test allowed jet configurations
         if( implemented){
           if( !( match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -3)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -4)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -5)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -5)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -6)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -7)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -7)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -8)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -8)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -9)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -10)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -11)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -11)
               && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqbar,qqbar2}, -12)
               && match_expectation(b,{"g",fkl},{qqbar,qqbar2,"g","g","g","g",fkl}, -12)
             ))
             return false;
          if (fkl == "2") {
           if( !( match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqbar,qqbar2}, -3)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","Wp","g","g","1"}, -4)
               && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqbar,qqbar2}, -5)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","Wp","g","g","1"}, -5)
               && match_expectation(f,{"2","g"},{"1","g","Wp","g","g",qqbar,qqbar2}, -6)
               && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqbar,qqbar2}, -7)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","g","g","Wp","1"}, -7)
               && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqbar,qqbar2}, -8)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"Wp","g","g","g","1"}, -8)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","Wp","g","g","1"}, -9)
               && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqbar,qqbar2}, -10)
               && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqbar,qqbar2}, -11)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","g","g","Wp","1"}, -11)
               && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqbar,qqbar2}, -12)
               && match_expectation(b,{"g","2"},{qqbar,qqbar2,"g","Wp","g","g","1"}, -12)
             ))
             return false;
          }
         }
       }
     return true;
   }
 
   //! test central qqbar configurations
   //! if implemented==false : check processes that are not in HEJ yet
   bool check_central_qqbar(bool const implemented=true){
     using namespace HEJ;
     auto const t{ implemented?event_type::qqbar_mid:event_type::non_resummable };
     std::vector<std::string> bosons;
     if(implemented)
       bosons = ALL_W;
     else {
       bosons = ALL_G_Z;
       bosons.emplace_back("h");
       bosons.emplace_back("Z_photon_mix");
     }
     for(std::string const & qqbar: ALL_QUARKS)        // all quark flavours
       for(std::string const & fkl1: ALL_PARTONS)
         for(std::string const & fkl2: ALL_PARTONS){
           std::string const qqbar2{ std::to_string(HEJ::to_ParticleID(qqbar)*-1) };
           for(std::size_t njet=4; njet<=MAX_MULTI; ++njet){                 // all multiplicities >3
             if(njet==5) continue;
             for(std::size_t qqbar_pos=1; qqbar_pos<njet-2; ++qqbar_pos){ // any qqbar position
               std::array<std::string,2> base_in;
               std::vector<std::string> base_out(njet, "g");
               base_in[0] = fkl1;
               base_in[1] = fkl2;
               base_out.front()    = fkl1;
               base_out.back()     = fkl2;
               base_out[qqbar_pos]   = qqbar;
               base_out[qqbar_pos+1] = qqbar2;
               if( implemented && !match_expectation(t, base_in, base_out) )
                 return false;
               for(auto const & boson: bosons)         // any boson
                 for(std::size_t pos=0; pos<=njet; ++pos){         // at any position
                   if( to_ParticleID(boson) == pid::higgs
                     && (pos==qqbar_pos || pos==qqbar_pos+1) )
                     continue;
                   auto const & in{base_in};
                   auto out{base_out};
                   // change quark flavours for W
                   const int couple_idx{ std::uniform_int_distribution<int>{0,2}(RAN) };
                   // (randomly) try couple to FKL, else fall-back to qqbar
                   if( couple_idx == 0 && couple_quark(boson, out.front()) ){}
                   else if( couple_idx == 1 && couple_quark(boson, out.back()) ){}
                   else {
                     if(!couple_quark(boson, out[qqbar_pos]))
                       couple_quark(boson, out[qqbar_pos+1]);
                   }
                   out.insert(out.begin()+pos, boson);
                   if(!match_expectation(t, in, out))
                     return false;
                 }
             }
           }
         }
     return true;
   }
 
   // this checks a (non excessive) list of non-resummable states
   bool check_non_resummable(){
     auto type{ HEJ::event_type::non_resummable};
     return
       // 1j - crossing lines
          match_expectation(type, {"g","2"},  {"2","h"})
       && match_expectation(type, {"-1","g"}, {"h","-1"})
       // 2j - crossing lines
       && match_expectation(type, {"g","2"},  {"2","g"})
       && match_expectation(type, {"-1","g"}, {"g","-1"})
       && match_expectation(type, {"1","-1"}, {"-1","1"})
       && match_expectation(type, {"g","2"},  {"2","g","h"})
       && match_expectation(type, {"1","2"},  {"2","h","1"})
       && match_expectation(type, {"1","-1"}, {"h","-1","1"})
       && match_expectation(type, {"g","2"},  {"Wp","1","g"})
       && match_expectation(type, {"1","-1"}, {"-2","Wp","1"})
       && match_expectation(type, {"4","g"},  {"g","3","Wp"})
       && match_expectation(type, {"1","-2"}, {"-1","Wm","1"})
       && match_expectation(type, {"g","3"},  {"4","g","Wm"})
       && match_expectation(type, {"1","3"},  {"Wm","4","1"})
       && match_expectation(type, {"g","2"},  {"Z_photon_mix","2","g"})
       && match_expectation(type, {"1","-1"}, {"-1","Z_photon_mix","1"})
       && match_expectation(type, {"4","g"},  {"g","4","Z_photon_mix"})
       // 2j - qqbar
       && match_expectation(type, {"g","g"},  {"1","-1"})
       && match_expectation(type, {"g","g"},  {"-2","2","h"})
       && match_expectation(type, {"g","g"},  {"-4","Wp","3"})
       && match_expectation(type, {"g","g"},  {"Wm","-1","2"})
       && match_expectation(type, {"g","g"},  {"-3","Z_photon_mix","3"})
       // 3j - crossing lines
       && match_expectation(type, {"g","4"},  {"4","g","g"})
       && match_expectation(type, {"-1","g"}, {"g","g","-1"})
       && match_expectation(type, {"1","3"},  {"3","g","1"})
       && match_expectation(type, {"-2","2"}, {"2","g","-2","h"})
       && match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"})
       && match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"})
       && match_expectation(type, {"-1","g"}, {"1","-1","-1"})
       && match_expectation(type, {"1","-4"}, {"Z_photon_mix","-4","g","1"})
       // higgs inside uno
       && match_expectation(type, {"-1","g"}, {"g","h","-1","g"})
       && match_expectation(type, {"-1","1"}, {"g","h","-1","1"})
       && match_expectation(type, {"g","2"},  {"g","2","h","g"})
       && match_expectation(type, {"-1","1"}, {"-1","1","h","g"})
       // higgs outside uno
       && match_expectation(type, {"-1","g"}, {"h","g","-1","g"})
       && match_expectation(type, {"-1","1"}, {"-1","1","g","h"})
       // higgs inside qqbar
       && match_expectation(type, {"g","g"}, {"-1","h","1","g","g"})
       && match_expectation(type, {"g","g"}, {"g","-1","h","1","g"})
       && match_expectation(type, {"g","g"}, {"g","g","2","h","-2"})
       // higgs outside qqbar
       && match_expectation(type, {"g","g"}, {"h","-1","1","g","g"})
       && match_expectation(type, {"g","g"}, {"g","g","2","-2","h"})
       // 4j - two uno
       && match_expectation(type, {"-2","2"}, {"g","-2","2","g"})
       && match_expectation(type, {"1","3"},  {"g","1","h","3","g"})
       && match_expectation(type, {"1","2"},  {"g","1","3","Wp","g"})
       && match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"})
       && match_expectation(type, {"3","2"},  {"g","3","Z_photon_mix","2","g"})
       // 4j - two gluon outside
       && match_expectation(type, {"g","4"},  {"g","4","g","g"})
       && match_expectation(type, {"1","3"},  {"1","3","h","g","g"})
       && match_expectation(type, {"1","2"},  {"1","3","g","Wp","g"})
       && match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"})
       && match_expectation(type, {"-1","g"}, {"g","g","-1","g"})
       && match_expectation(type, {"1","3"},  {"g","g","1","3","h"})
       && match_expectation(type, {"1","2"},  {"g","g","1","Wp","3"})
       && match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"})
       && match_expectation(type, {"-1","2"}, {"g","g","-1","Z_photon_mix","2"})
       // 4j - ggx+uno
       && match_expectation(type, {"g","4"},  {"1","-1","4","g"})
       && match_expectation(type, {"2","g"},  {"g","2","-3","3"})
       && match_expectation(type, {"g","4"},  {"1","-1","h","4","g"})
       && match_expectation(type, {"2","g"},  {"g","2","-3","3","h"})
       && match_expectation(type, {"g","4"},  {"Wp","1","-1","3","g"})
       && match_expectation(type, {"2","g"},  {"g","2","-4","Wp","3"})
       && match_expectation(type, {"g","4"},  {"2","Wm","-1","4","g"})
       && match_expectation(type, {"2","g"},  {"g","2","Wp","-3","4"})
       && match_expectation(type, {"-4","g"}, {"g","-4","-3","3","Z_photon_mix"})
       // 3j - crossing+uno
       && match_expectation(type, {"1","4"},  {"g","4","1"})
       && match_expectation(type, {"1","4"},  {"4","1","g"})
       && match_expectation(type, {"1","4"},  {"g","h","4","1"})
       && match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"})
       && match_expectation(type, {"1","4"},  {"3","1","Wp","g"})
       && match_expectation(type, {"1","4"},  {"3","1","g","h"})
       && match_expectation(type, {"2","3"},  {"3","2","Z_photon_mix","g"})
       // 3j - crossing+qqbar
       && match_expectation(type, {"1","g"},  {"-1","1","g","1"})
       && match_expectation(type, {"1","g"},  {"-1","1","1","g"})
       && match_expectation(type, {"g","1"},  {"1","g","1","-1"})
       && match_expectation(type, {"g","1"},  {"g","1","1","-1"})
       && match_expectation(type, {"1","g"},  {"2","-2","g","1"})
       && match_expectation(type, {"1","g"},  {"2","-2","1","g"})
       && match_expectation(type, {"g","1"},  {"1","g","-2","2"})
       && match_expectation(type, {"g","1"},  {"g","1","-2","2"})
       && match_expectation(type, {"1","g"},  {"-1","1","h","g","1"})
       && match_expectation(type, {"1","g"},  {"-1","h","1","1","g"})
       && match_expectation(type, {"g","1"},  {"1","g","1","h","-1"})
       && match_expectation(type, {"g","1"},  {"h","g","1","1","-1"})
       && match_expectation(type, {"1","g"},  {"2","-2","1","g","h"})
       && match_expectation(type, {"g","1"},  {"g","h","1","-2","2"})
       && match_expectation(type, {"1","g"},  {"Wp","3","-4","g","1"})
       && match_expectation(type, {"3","g"},  {"-2","Wm","1","3","g"})
       && match_expectation(type, {"g","1"},  {"1","g","Wm","-3","4"})
       && match_expectation(type, {"g","-3"},  {"g","-3","-1","Wp","2"})
       && match_expectation(type, {"g","2"},  {"2","g","Z_photon_mix","4","-4"})
       // 4j- gluon in qqbar
       && match_expectation(type, {"g","1"},  {"1","g","-1","1"})
       && match_expectation(type, {"1","g"},  {"1","-1","g","1"})
       && match_expectation(type, {"g","1"},  {"1","g","Wm","-2","1"})
       && match_expectation(type, {"2","g"},  {"2","-2","g","Wp","1"})
       && match_expectation(type, {"g","g"},  {"Wp","3","g","-4","g"})
       && match_expectation(type, {"1","g"},  {"1","h","-1","g","1"})
       && match_expectation(type, {"3","g"},  {"3","1","g","Z_photon_mix","-1"})
       // 6j - two qqbar
       && match_expectation(type, {"g","g"},  {"1","-1","g","g","1","-1"})
       && match_expectation(type, {"g","g"},  {"1","-1","g","1","-1","g"})
       && match_expectation(type, {"g","g"},  {"g","1","-1","g","1","-1"})
       && match_expectation(type, {"g","g"},  {"g","1","-1","1","-1","g"})
       && match_expectation(type, {"g","g"},  {"g","1","1","-1","-1","g"})
       && match_expectation(type, {"g","g"},  {"h","1","-1","g","g","1","-1"})
       && match_expectation(type, {"g","g"},  {"1","Wp","-2","g","1","-1","g"})
       && match_expectation(type, {"g","g"},  {"g","1","Wp","-1","g","1","-2"})
       && match_expectation(type, {"g","g"},  {"g","1","-1","Wm","2","-1","g"})
       && match_expectation(type, {"g","g"},  {"g","1","2","-1","Wm","-1","g"})
       && match_expectation(type, {"g","g"},  {"2","-2","g","-1","1","Z_photon_mix","g"})
       // random stuff (can be non-physical)
       && match_expectation(type, {"g","g"},  {"1","-2","2","-1"}) // != 2 qqbar
       && match_expectation(type, {"g","g"},  {"1","-2","2","g"})  // could be qqbar
       && match_expectation(type, {"e+","e-"},{"1","-1"})          // bad initial state
       && match_expectation(type, {"1","e-"}, {"g","1","Wm"})      // bad initial state
       && match_expectation(type, {"h","g"},  {"g","g"})           // bad initial state
       && match_expectation(type, {"-1","g"}, {"-1","1","1"})      // bad qqbar
       && match_expectation(type, {"-1","g"}, {"1","1","-1"})      // crossing in bad qqbar
       && match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqbar
       && match_expectation(type, {"1","2"},  {"1","-1","g","g","g","2"})  // bad qqbar
       && match_expectation(type, {"1","2"},  {"1","-1","-2","g","g","2"}) // gluon in bad qqbar
       && match_expectation(type, {"g","g"},  {"-1","2","g","g"}) // wrong back qqbar
       && match_expectation(type, {"g","g"},  {"g","g","2","1"}) // wrong forward qqbar
       && match_expectation(type, {"g","g"},  {"g","-2","1","g"}) // wrong central qqbar
       && match_expectation(type, {"1","g"},  {"1","-2","g","g","Wp"}) // extra quark
       && match_expectation(type, {"g","1"},  {"g","g","-2","1","Wp"}) // extra quark
       && match_expectation(type, {"g","1"},  {"g","g","Wp","-2","1"}) // extra quark
       && match_expectation(type, {"g","1"},  {"g","-2","1","g","Wp"}) // extra quark
       && match_expectation(type, {"g","g"},  {"g","g","g","-2","1","-1","Wp"}) // extra quark
       && match_expectation(type, {"1","g"},  {"g","Wp","1","-2","g"}) // extra quark
       && match_expectation(type, {"g","g"},  {"1","-1","-2","g","g","g","Wp"}) // extra quark
       ;
   }
 
   // Two boson states, that are currently not implemented
   // Check for supported final states
   bool check_bad_FS(){
     auto type{ HEJ::event_type::bad_final_state};
     return
          match_expectation(type, {"g","g"},  {"g","p","p","g"}) // protons
       && match_expectation(type, {"-4","-1"},{"-4","g","11","-11","-2"}) // leptons should be in decay
       && match_expectation(type, {"-4","-1"},{"-4","g","-13","g","-2"}) // leptons should be in decay
     ;
   }
 
   // not enough jets
   bool check_not_enough_jets(){
     auto type{ HEJ::event_type::not_enough_jets};
     return
          match_expectation(type, {"g","g"},  {})
       && match_expectation(type, {"1","-1"}, {})
       && match_expectation(type, {"g","-1"}, {"-1"})
       && match_expectation(type, {"g","g"},  {"g"})
       && match_expectation(type, {"g","g"},  {"h"})
       && match_expectation(type, {"g","2"},  {"Wp","1"})
       && match_expectation(type, {"g","2"},  {"Z","2"})
     ;
   }
 
   // 2 boson final-states
   bool check_2_boson(bool const implemented=true){
     auto type = (implemented) ? HEJ::event_type::FKL : HEJ::event_type::non_resummable;
 
     // Implemented
     if ( implemented ) {
       //Checks WpWp and WmWm separatelty to use couple_quark method.
       for(std::string const & q1 : ALL_QUARKS){
         for(std::string const & q2 : ALL_QUARKS){
           for(size_t njet = 2; njet < 6; ++njet){
             // Event setup
             std::array<std::string,2> in {q1, q2};
             std::vector<std::string>  out(njet, "g");
             out.front() = q1;
             out.back()  = q2;
             // Adjust quark flavours for the boson, if possible
             // WpWp
             if( couple_quark("Wp", out.front()) &&
                 couple_quark("Wp", out.back() )
               ) {
 
                 std::unordered_map< size_t, std::vector<std::string> > decays;
                 for(size_t i = 0; i < 2; ++i){
                   out.emplace_back("Wp");
                   decays.emplace( out.size()-1,
                               std::initializer_list<std::string>{"e+", "nu_e"} );
                 }
                 // Check classification
                 if( !match_expectation(type, rapidity_order_ps(in, out, false, decays)) ){return false;}
             }
           }
         }
       }
 
       for(std::string const & q1 : ALL_QUARKS){
         for(std::string const & q2 : ALL_QUARKS){
           for(size_t njet = 2; njet < 6; ++njet){
             // Event setup
             std::array<std::string,2> in {q1, q2};
             std::vector<std::string>  out(njet, "g");
             out.front() = q1;
             out.back()  = q2;
             // Adjust quark flavours for the boson, if possible
             // WmWm
             if( couple_quark("Wm", out.front()) &&
                 couple_quark("Wm", out.back() )
               ) {
 
                 std::unordered_map< size_t, std::vector<std::string> > decays;
                 for(size_t i = 0; i < 2; ++i){
                   out.emplace_back("Wm");
                   decays.emplace( out.size()-1,
                               std::initializer_list<std::string>{"e-", "nu_e_bar"} );
                 }
                 // Check classification
                 if( !match_expectation(type, rapidity_order_ps(in, out, false, decays)) ){return false;}
             }
           }
         }
       }
       return true;
     }
 
     // Not Implemented (non_resummable)
     else {
       return
           // h h
           match_expectation(type, {"g","g"}, {"h","g","h","g"})
           // Wp Wm
           && match_expectation(type,
               rapidity_order_ps(
                 {"u","d"},
                 {"e+","nu_e","d","mu-","nu_mu_bar","u"},
                 true)
              )
           // Wm h
           && match_expectation(type,
               rapidity_order_ps(
                 {"u","d"},
                 {"e-","nu_e_bar","u","h","u"},
                 true)
              )
           // Wp h
           && match_expectation(type,
               rapidity_order_ps(
                 {"u","d"},
-                {"e+","nu_e","d","h","u"},
+                {"e+","nu_e","d","h","d"},
                 true)
              )
       ;
     }
   }
 
   // not implemented processes
   bool check_not_implemented(){
     return check_fkl(false)
         && check_uno(false)
 
         && check_extremal_qqbar(false)
         && check_central_qqbar(false)
         && check_2_boson(false);
   }
 
   // h + single jet
   bool check_hj(){
     using namespace HEJ::event_type;
     return
       match_expectation(FKL, {"g","g"},  {"h", "g"})
       && match_expectation(FKL, {"g","g"},  {"g", "h"})
       && match_expectation(FKL, {"u","g"},  {"u", "h"})
       && match_expectation(FKL, {"g","u"},  {"h", "u"})
       && match_expectation(non_resummable, {"g","u"},  {"u", "h"})
       && match_expectation(non_resummable, {"u","u"},  {"u", "h"})
     ;
 
   }
 
 } // namespace
 
 
 int main() {
   // tests for "no false negatives"
   // i.e. all HEJ-configurations get classified correctly
   if(!check_fkl()) return EXIT_FAILURE;
   if(!check_uno()) return EXIT_FAILURE;
 
 
   if(!check_extremal_qqbar()) return EXIT_FAILURE;
   if(!check_central_qqbar()) return EXIT_FAILURE;
   if(!check_2_boson()) return EXIT_FAILURE;
 
   // test for "no false positive"
   // i.e. non-resummable gives non-resummable
   if(!check_non_resummable()) return EXIT_FAILURE;
   if(!check_bad_FS()) return EXIT_FAILURE;
   if(!check_not_enough_jets()) return EXIT_FAILURE;
   if(!check_not_implemented()) return EXIT_FAILURE;
   if(!check_hj()) return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
 }