diff --git a/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml
index 26ccc5c..0633304 100644
--- a/FixedOrderGen/configFO.yml
+++ b/FixedOrderGen/configFO.yml
@@ -1,69 +1,68 @@
 # number of generated events
 events: 200
 
 jets:
   min pt: 20
   peak pt: 30
   algorithm: antikt
   R: 0.4
   max rapidity: 5
 
 beam:
   energy: 6500
   particles: [p, p]
 
 pdf: 230000
 
 process: p p => h 4j
 
 # fraction of events with two extremal emissions in one direction
 # that contain an subleading emission e.g. unordered emission
 subleading fraction: 0.05
 # Allow different subleading configurations
 # By default all implemented processes are allowed
 #
-# subleading channels: unordered
+# subleading channels: [unordered, qqx]
 
 scales: max jet pperp
 
 event output:
   - HEJFO.lhe
-#  - HEJ.lhe
-#  - HEJ_events.hepmc
+#  - HEJFO_events.hepmc
 
 particle properties:
   Higgs:
     mass: 125
     width: 0.004165
     decays: {into: [photon, photon], branching ratio: 0.0023568762400521404}
 
 random generator:
   name: mixmax
   # seed: 1
 
 # unweighting parameters
 # remove to obtain weighted events
 unweight:
   sample size: 200  # should be similar to "events:", but not more than ~10000
   max deviation: 0
 
 # to use a rivet analysis
 #
 # analysis:
 #   rivet: MC_XS # rivet analysis name
 #   output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
 #
 # to use a custom analysis
 #
 # analysis:
 #   plugin: /path/to/libmyanalysis.so
 #   my analysis parameter: some value
 
 # parameters for Higgs-gluon couplings
 # this requires compilation with qcdloop
 #
 # Higgs coupling:
 #   use impact factors: false
 #   mt: 174
 #   include bottom: true
 #   mb: 4.7
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index b471662..93ac98c 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,665 +1,663 @@
 #include "PhaseSpacePoint.hh"
 
 #include <random>
 #include <algorithm>
 
 #include "HEJ/Constants.hh"
 #include "HEJ/kinematics.hh"
 #include "HEJ/utility.hh"
 #include "HEJ/exceptions.hh"
 
 #include "Process.hh"
 #include "Subleading.hh"
 
 #include <CLHEP/Random/Randomize.h>
 #include <CLHEP/Random/RanluxEngine.h>
 
 using namespace HEJ;
 
 namespace HEJFOG{
 
   static_assert(
       std::numeric_limits<double>::has_quiet_NaN,
       "no quiet NaN for double"
   );
   constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
   HEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
     HEJ::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),
             HEJ::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{
     bool can_swap_to_uno(
         HEJ::Particle const & p1, HEJ::Particle const & p2
     ){
       return is_parton(p1)
         && p1.type != pid::gluon
         && p2.type == pid::gluon;
     }
 
     size_t count_gluons(std::vector<Particle>::const_iterator first,
       std::vector<Particle>::const_iterator last){
       return std::count_if(first, last, [](Particle const & p)
         {return p.type == pid::gluon;});
     }
 
     /** assumes FKL configurations between first and last,
      * else there can be a quark in a non-extreme position
      * e.g. uno configuration gqg would pass
      */
     bool can_change_to_qqx(
       std::vector<Particle>::const_iterator first,
       std::vector<Particle>::const_iterator last){
       return 1 < count_gluons(first,last);
     }
     bool is_AWZ_proccess(Process const & proc){
       return proc.boson && is_AWZ_boson(*proc.boson);
     }
 
     bool is_up_type(Particle const & part){
       return HEJ::is_anyquark(part) && !(abs(part.type)%2);
     }
     bool is_down_type(Particle const & part){
       return HEJ::is_anyquark(part) && abs(part.type)%2;
     }
     /// true iff parton can couple to a W
     bool can_couple_to_W(Particle const & part, int const sign_W){
       return abs(part.type)<5
         && ( (sign_W*part.type > 0 && is_up_type(part))
           || (sign_W*part.type < 0 && is_down_type(part)) );
     }
   }
 
   void PhaseSpacePoint::maybe_turn_to_subl(
       double chance,
       unsigned int const channels,
       Process const & proc,
       HEJ::RNG & ran
   ){
     if(proc.njets <= 2) return;
     assert(outgoing_.size() >= 2);
 
     // decide what kind of subleading process is allowed
     bool allow_uno = false;
     bool allow_strange = true;
     const size_t nout = outgoing_.size();
-    const bool can_be_uno_backward = channels&Subleading::uno?
-      can_swap_to_uno(outgoing_[0], outgoing_[1]):false;
-    const bool can_be_uno_forward = channels&Subleading::uno?
-      can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]):false;
+    const bool can_be_uno_backward = (channels&Subleading::uno)
+      && can_swap_to_uno(outgoing_[0], outgoing_[1]);
+    const bool can_be_uno_forward = (channels&Subleading::uno)
+      && can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]);
     allow_uno = can_be_uno_backward || can_be_uno_forward;
 
     bool allow_qqx = false;
     if(is_AWZ_proccess(proc)) {
-      allow_qqx = channels&Subleading::qqx?
-        can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend()):false;
+      allow_qqx = (channels&Subleading::qqx)
+        && can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend());
       const int sign_W = *proc.boson>0?1:-1;
       if(std::none_of(outgoing_.cbegin(), outgoing_.cend(),
           [sign_W](Particle const & p){ return can_couple_to_W(p, sign_W);})) {
-        if(!(channels&Subleading::qqx))
-          throw std::logic_error(
-            "Invalid configuration, require qqx pair but channel forbidden");
         // enforce qqx if A/W/Z can't couple somewhere else
+        assert(allow_qqx);
         allow_uno = false;
         chance = 1.;
         // strange not allowed for W
         if(abs(*proc.boson)== pid::Wp) allow_strange = false;
       }
     }
 
     if(!allow_uno && !allow_qqx) return;
     if(ran.flat() < chance){
       weight_ /= chance;
       if(allow_uno && !allow_qqx){
           turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
       } else if (!allow_uno && allow_qqx) {
         turn_to_qqx(allow_strange, ran);
       } else {
         assert( allow_uno && allow_qqx);
         if(ran.flat() < 0.5) turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
         else turn_to_qqx(allow_strange, ran);
         weight_ *= 2.;
       }
     } else weight_ /= 1 - chance;
 
   }
 
   void PhaseSpacePoint::turn_to_uno(
       const bool can_be_uno_backward, const bool can_be_uno_forward,
       HEJ::RNG & ran
   ){
     if(!can_be_uno_backward && !can_be_uno_forward) return;
     const size_t nout = outgoing_.size();
     if(can_be_uno_backward && can_be_uno_forward){
       if(ran.flat() < 0.5){
         std::swap(outgoing_[0].type, outgoing_[1].type);
       } else {
         std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
       }
       weight_ *= 2.;
     } else if(can_be_uno_backward){
       std::swap(outgoing_[0].type, outgoing_[1].type);
     } else {
       assert(can_be_uno_forward);
       std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
     }
   }
 
   void PhaseSpacePoint::turn_to_qqx(const bool allow_stange, HEJ::RNG & ran){
     /// find first and last gluon in FKL chain
     auto first = std::find_if(outgoing_.begin(), outgoing_.end(),
       [](Particle const & p){return p.type == pid::gluon;});
     std::vector<Particle*> FKL_gluons;
     for(auto p = first; p<outgoing_.end(); ++p){
       if((*p).type == pid::gluon) FKL_gluons.push_back(&*p);
       else if(is_quark(*p) || is_antiquark(*p)) break;
     }
     const size_t ng = FKL_gluons.size();
     if(ng < 2)
       throw std::logic_error("not enough gluons to create qqx");
     // select flavour of quark
     const double r1 = 2.*ran.flat()-1.;
     const double max_flavour = allow_stange?n_f:n_f-1;
     weight_ *= max_flavour*2;
     int flavour = pid::down;
     for (double sum = 1./max_flavour; sum < std::abs(r1); sum += 1./max_flavour)
       ++flavour;
     flavour*=r1<0.?-1:1;
     // select gluon for switch
     const size_t idx = floor((ng-1) * ran.flat());
     weight_ *= (ng-1);
     FKL_gluons[idx]->type = ParticleID(flavour);
     FKL_gluons[idx+1]->type = ParticleID(-flavour);
   }
 
   template<class ParticleMomenta>
   fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
       ParticleMomenta const & other_momenta,
       const double mass_square, const double y
   ) const {
     std::array<double,2> pt{0.,0.};
     for (auto const & p: other_momenta) {
       pt[0]-= p.px();
       pt[1]-= p.py();
     }
 
     const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
     const double pz=mperp*sinh(y);
     const double E=mperp*cosh(y);
 
     return {pt[0], pt[1], pz, E};
   }
 
   PhaseSpacePoint::PhaseSpacePoint(
     Process const & proc,
     JetParameters const & jet_param,
     HEJ::PDF & pdf, double E_beam,
     double const subl_chance,
     unsigned int const subl_channels,
     ParticlesPropMap const & particles_properties,
     HEJ::RNG & ran
   )
   {
     assert(proc.njets >= 2);
     if(proc.boson
       && particles_properties.find(*(proc.boson))
         == particles_properties.end())
       throw HEJ::missing_option("Boson "
         +std::to_string(*(proc.boson))+" can't be generated: missing properties");
     status_ = good;
     weight_ = 1;
 
     const int nout = proc.njets + (proc.boson?1:0);
     outgoing_.reserve(nout);
 
     // generate parton momenta
     const bool is_pure_jets = !proc.boson;
     auto partons = gen_LO_partons(
         proc.njets, is_pure_jets, jet_param, E_beam, ran
     );
     // pre fill flavour with gluons
     for(auto&& p_out: partons) {
       outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
     }
     if(status_ != good) return;
 
     // create boson
     if(proc.boson){
       const auto & boson_prop = particles_properties.at(*proc.boson);
       auto boson(gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran));
       const auto pos = std::upper_bound(
           begin(outgoing_),end(outgoing_),boson,rapidity_less{}
       );
       outgoing_.insert(pos, std::move(boson));
       if(! boson_prop.decays.empty()){
         const size_t boson_idx = std::distance(begin(outgoing_), pos);
         decays_.emplace(
             boson_idx,
             decay_boson(outgoing_[boson_idx], boson_prop.decays, ran)
         );
       }
     }
     // normalisation of momentum-conserving delta function
     weight_ *= pow(2*M_PI, 4);
 
     reconstruct_incoming(proc, subl_channels, pdf,E_beam,jet_param.min_pt, ran);
     if(status_ != good) return;
     // set outgoing states
     most_backward_FKL(outgoing_).type = incoming_[0].type;
     most_forward_FKL(outgoing_).type = incoming_[1].type;
 
     maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);
 
     if(proc.boson) couple_boson(*proc.boson, ran);
   }
 
   double PhaseSpacePoint::gen_hard_pt(
       int np , double ptmin, double ptmax, double y,
       HEJ::RNG & ran
   ) {
     // heuristic parameters for pt sampling
     const double ptpar = ptmin + np/5.;
     const double arg_small_y = atan((ptmax - ptmin)/ptpar);
     const double y_cut = 3.;
 
     const double r1 = ran.flat();
     if(y < y_cut){
       const double pt = ptmin + ptpar*tan(r1*arg_small_y);
       const double temp = cos(r1*arg_small_y);
       weight_ *= pt*ptpar*arg_small_y/(temp*temp);
       return pt;
     }
 
     const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
     const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
     const double pt = ptmin - ptpar2*std::log(1-r1*temp);
     weight_ *= pt*ptpar2*temp/(1-r1*temp);
     return pt;
   }
 
   double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) {
     constexpr double ptpar = 4.;
 
     const double r = ran.flat();
     const double pt = max_pt + ptpar/np*std::log(r);
     weight_ *= pt*ptpar/(np*r);
     return pt;
   }
 
   double PhaseSpacePoint::gen_parton_pt(
       int count, JetParameters const & jet_param, double max_pt, double y,
       HEJ::RNG & ran
   ) {
     constexpr double p_small_pt = 0.02;
 
     if(! jet_param.peak_pt) {
        return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran);
     }
     const double r = ran.flat();
     if(r > p_small_pt) {
       weight_ /= 1. - p_small_pt;
       return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran);
     }
     weight_ /= p_small_pt;
     const double pt =  gen_soft_pt(count, *jet_param.peak_pt, ran);
     if(pt < jet_param.min_pt) {
       weight_=0.0;
       status_ = not_enough_jets;
       return jet_param.min_pt;
     }
     return pt;
   }
 
   std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
       int np, bool is_pure_jets,
       JetParameters const & jet_param,
       double max_pt,
       HEJ::RNG & ran
   ){
     if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
 
     weight_ /= pow(16.*pow(M_PI,3),np);
     weight_ /= std::tgamma(np+1); //remove rapidity ordering
     std::vector<fastjet::PseudoJet> partons;
     partons.reserve(np);
     for(int i = 0; i < np; ++i){
       const double y = -jet_param.max_y + 2*jet_param.max_y*ran.flat();
       weight_ *= 2*jet_param.max_y;
 
       const bool is_last_parton = i+1 == np;
       if(is_pure_jets && is_last_parton) {
         constexpr double parton_mass_sq = 0.;
         partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y));
         break;
       }
 
       const double phi = 2*M_PI*ran.flat();
       weight_ *= 2.0*M_PI;
 
       const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran);
       if(weight_ == 0.0) return {};
 
       partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
       assert(jet_param.min_pt <= partons[i].pt());
       assert(partons[i].pt() <= max_pt+1e-5);
     }
     // Need to check that at LO, the number of jets = number of partons;
     fastjet::ClusterSequence cs(partons, jet_param.def);
     auto cluster_jets=cs.inclusive_jets(jet_param.min_pt);
     if (cluster_jets.size()!=unsigned(np)){
       weight_=0.0;
       status_ = not_enough_jets;
       return {};
     }
 
     std::sort(begin(partons), end(partons), rapidity_less{});
 
     return partons;
   }
 
   Particle PhaseSpacePoint::gen_boson(
       HEJ::ParticleID bosonid, double mass, double width,
       HEJ::RNG & ran
   ){
     // Usual phase space measure
     weight_ /= 16.*pow(M_PI, 3);
 
     // Generate a y Gaussian distributed around 0
     /// @TODO: magic number only for Higgs
     /// @TODO better sampling for W
     const double y = random_normal(1.6, ran);
     const double r1 = ran.flat();
     const double sH = mass*(
         mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width))
     );
 
     auto p = gen_last_momentum(outgoing_, sH, y);
 
     return Particle{bosonid, std::move(p)};
   }
 
   Particle const & PhaseSpacePoint::most_backward_FKL(
       std::vector<Particle> const & partons
   ) const{
     if(!HEJ::is_parton(partons[0])) return partons[1];
     return partons[0];
   }
 
   Particle const & PhaseSpacePoint::most_forward_FKL(
       std::vector<Particle> const & partons
   ) const{
     const size_t last_idx = partons.size() - 1;
     if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
     return partons[last_idx];
   }
 
   Particle & PhaseSpacePoint::most_backward_FKL(
       std::vector<Particle> & partons
   ) const{
     if(!HEJ::is_parton(partons[0])) return partons[1];
     return partons[0];
   }
 
   Particle & PhaseSpacePoint::most_forward_FKL(
       std::vector<Particle> & partons
   ) const{
     const size_t last_idx = partons.size() - 1;
     if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
     return partons[last_idx];
   }
 
   namespace {
     /// partons are ordered: even = anti, 0 = gluon
     ParticleID index_to_pid(size_t i){
       if(!i) return pid::gluon;
       return static_cast<ParticleID>(i%2?(i+1)/2:-i/2);
     }
 
     /// partons are ordered: even = anti, 0 = gluon
     size_t pid_to_index(ParticleID id){
       if(id==pid::gluon) return 0;
       return id>0?id*2-1:abs(id)*2;
     }
 
     std::bitset<11> init_allowed(ParticleID const id){
       if(abs(id) == pid::proton)
         return ~0;
       std::bitset<11> out = 0;
       if(is_parton(id))
         out[pid_to_index(id)] = 1;
       return out;
     }
 
     /// decides which "index" (see index_to_pid) are allowed for process
     std::bitset<11> allowed_quarks(ParticleID const boson){
       std::bitset<11> allowed = ~0;
       if(abs(boson) == pid::Wp){
         // special case W:
         // Wp: anti-down or up-type quark, no b/t -> 0001100110(1) = 205
         // Wm: down or anti-up-type quark, no b/t -> 0010011001(1) = 307
         allowed = boson>0?205:307;
       }
       return allowed;
     }
   }
 
   /**
    * checks which partons are allowed as initial state:
    * 1. only allow what is given in the Runcard (p -> all)
    * 2. A/W/Z require something to couple to
    *  a) no qqx => no incoming gluon
    *  b) 2j     => no incoming gluon
    *  c) 3j     => can couple OR is gluon => 2 gluons become qqx later
    */
   std::array<std::bitset<11>,2> PhaseSpacePoint::filter_partons(
       Process const & proc, unsigned int const subl_channels, HEJ::RNG & ran
   ){
     std::array<std::bitset<11>,2> allowed_partons{
         init_allowed(proc.incoming[0]),
         init_allowed(proc.incoming[1])
       };
     bool const allow_qqx = subl_channels&Subleading::qqx;
     // special case A/W/Z
     if(is_AWZ_proccess(proc) && ((proc.njets < 4) || !allow_qqx)){
       // all possible incoming states
       auto allowed(allowed_quarks(*proc.boson));
       if(proc.njets == 2 || !allow_qqx) allowed[0]=0;
 
       // possible states per leg
       std::array<std::bitset<11>,2> const maybe_partons{
         allowed_partons[0]&allowed, allowed_partons[1]&allowed};
 
       if(maybe_partons[0].any() && maybe_partons[1].any()){
         // two options to get allowed initial state => choose one at random
         const size_t idx = ran.flat() < 0.5;
         allowed_partons[idx] = maybe_partons[idx];
+        // else choose the possible
       } else if(maybe_partons[0].any()) {
-        // only first possible => choose
         allowed_partons[0] = maybe_partons[0];
       } else if(maybe_partons[1].any()) {
-        // only second possible => choose
         allowed_partons[1] = maybe_partons[1];
       } else{
         throw std::invalid_argument{"Incoming state not allowed."};
       }
     }
     return allowed_partons;
   }
 
   void PhaseSpacePoint::reconstruct_incoming(
       Process const & proc, unsigned int const subl_channels,
       HEJ::PDF & pdf, double E_beam,
       double uf,
       HEJ::RNG & ran
   ){
     std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
     // calculate xa, xb
     const double sqrts=2*E_beam;
     const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
     const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
     // abort if phase space point is outside of collider energy reach
     if (xa>1. || xb>1.){
       weight_=0;
       status_ = too_much_energy;
       return;
     }
     // pick pdfs
     /** @TODO
      *  ufa, ufb don't correspond to our final scale choice.
      *  The HEJ scale generators currently expect a full event as input,
      *  so fixing this is not completely trivial
      */
     auto const & ids = proc.incoming;
     std::array<std::bitset<11>,2> allowed_partons(
       filter_partons(proc, subl_channels, ran));
     for(size_t i = 0; i < 2; ++i){
       if(ids[i] == pid::proton || ids[i] == pid::p_bar){
         incoming_[i].type =
           generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran);
       } else {
         assert(allowed_partons[i][pid_to_index(ids[i])]);
         incoming_[i].type = ids[i];
       }
     }
     assert(momentum_conserved(1e-7));
   }
 
   HEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
       size_t const beam_idx, double const x, double const uf,
       HEJ::PDF & pdf, std::bitset<11> allowed_partons, HEJ::RNG & ran
   ){
     std::array<double,11> pdf_wt;
     pdf_wt[0] = allowed_partons[0]?fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.;
     double pdftot = pdf_wt[0];
     for(size_t i = 1; i < pdf_wt.size(); ++i){
       pdf_wt[i] = allowed_partons[i]?4./9.*fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0;
       pdftot += pdf_wt[i];
     }
     const double r1 = pdftot * ran.flat();
     double sum = 0;
     for(size_t i=0; i < pdf_wt.size(); ++i){
       if (r1 < (sum+=pdf_wt[i])){
         weight_*= pdftot/pdf_wt[i];
         return index_to_pid(i);
       }
     }
     std::cerr << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "
       <<sum<<" "<<pdftot<<" "<<r1<<std::endl;
     throw std::logic_error{"Failed to choose parton flavour"};
   }
 
   void PhaseSpacePoint::couple_boson(
     HEJ::ParticleID const boson, HEJ::RNG & ran
   ){
     if(abs(boson) != pid::Wp) return; // only matters for W
+    /// @TODO this could be use to sanity check gamma and Z
 
     // find all possible quarks
     const int sign_W = boson>0?1:-1;
     std::vector<Particle*> allowed_parts;
     for(auto & part: outgoing_){
       // Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
       if ( can_couple_to_W(part, sign_W) )
         allowed_parts.push_back(&part);
     }
     if(allowed_parts.size() == 0){
       throw std::logic_error{"Found no parton for coupling with boson"};
     }
 
     // select one and flip it
     size_t idx = 0;
     if(allowed_parts.size() > 1){
       /// @TODO more efficient sampling
       ///       old code: probability[i] = exp(parton[i].y - W.y)
       idx = floor(ran.flat()*allowed_parts.size());
       weight_ *= allowed_parts.size();
     }
     allowed_parts[idx]->type =
       static_cast<ParticleID>( allowed_parts[idx]->type - sign_W );
   }
 
   double PhaseSpacePoint::random_normal(
       double stddev,
       HEJ::RNG & ran
   ){
     const double r1 = ran.flat();
     const double r2 = ran.flat();
     const double lninvr1 = -log(r1);
     const double result = stddev*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
     weight_ *= exp(result*result/(2*stddev*stddev))*sqrt(2.*M_PI)*stddev;
     return result;
   }
 
   bool PhaseSpacePoint::momentum_conserved(double ep) const{
     fastjet::PseudoJet diff;
     for(auto const & in: incoming()) diff += in.p;
     for(auto const & out: outgoing()) diff -= out.p;
     return nearby_ep(diff, fastjet::PseudoJet{}, ep);
   }
 
   Decay PhaseSpacePoint::select_decay_channel(
       std::vector<Decay> const & decays,
       HEJ::RNG & ran
   ){
     double br_total = 0.;
     for(auto const & decay: decays) br_total += decay.branching_ratio;
     // adjust weight
     // this is given by (channel branching ratio)/(chance to pick channel)
     // where (chance to pick channel) =
     //       (channel branching ratio)/(total branching ratio)
     weight_ *= br_total;
     const double r1 = br_total*ran.flat();
     double br_sum = 0.;
     for(auto const & decay: decays){
       br_sum += decay.branching_ratio;
       if(r1 < br_sum) return decay;
     }
     throw std::logic_error{"unreachable"};
   }
 
   std::vector<Particle> PhaseSpacePoint::decay_boson(
       HEJ::Particle const & parent,
       std::vector<Decay> const & decays,
       HEJ::RNG & ran
   ){
     const auto channel = select_decay_channel(decays, ran);
     if(channel.products.size() != 2){
       throw HEJ::not_implemented{
         "only decays into two particles are implemented"
       };
     }
     std::vector<Particle> decay_products(channel.products.size());
     for(size_t i = 0; i < channel.products.size(); ++i){
       decay_products[i].type = channel.products[i];
     }
     // choose polar and azimuth angle in parent rest frame
     const double E = parent.m()/2;
     const double theta = 2.*M_PI*ran.flat();
     const double cos_phi = 2.*ran.flat()-1.;
     const double sin_phi = sqrt(1. - cos_phi*cos_phi);  // Know 0 < phi < pi
     const double px = E*cos(theta)*sin_phi;
     const double py = E*sin(theta)*sin_phi;
     const double pz = E*cos_phi;
     decay_products[0].p.reset(px, py, pz, E);
     decay_products[1].p.reset(-px, -py, -pz, E);
     for(auto & particle: decay_products) particle.p.boost(parent.p);
     return decay_products;
   }
 }
diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc
index b3cb317..06ff18f 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,177 +1,177 @@
 #ifdef NDEBUG
 #undef NDEBUG
 #endif
 
 #include <algorithm>
 
 #include "JetParameters.hh"
 #include "ParticleProperties.hh"
 #include "PhaseSpacePoint.hh"
 #include "Process.hh"
 #include "Subleading.hh"
 
 #include "HEJ/Event.hh"
 #include "HEJ/Mixmax.hh"
 #include "HEJ/PDF.hh"
 #include "HEJ/utility.hh"
 
 using namespace HEJFOG;
 using namespace HEJ;
 
 namespace {
   void print_psp(PhaseSpacePoint const & psp){
     std::cerr << "Process:\n"
       << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> ";
     for(auto const & out: psp.outgoing()){
       std::cerr << out.type << " ";
     }
     std::cerr << "\n";
   }
   void bail_out(PhaseSpacePoint const & psp, std::string msg){
     print_psp(psp);
     throw std::logic_error{msg};
   }
 }
 
 int main(){
-  constexpr size_t n_psp_base = 19375;
+  constexpr size_t n_psp_base = 10375;
   const JetParameters jet_para{
     fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30};
   PDF pdf(11000, pid::proton, pid::proton);
   constexpr double E_cms = 13000.;
   constexpr double subl_change = 0.8;
   const ParticlesPropMap boson_prop{
       {pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
       {pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
     };
   HEJ::Mixmax ran{};
 
   auto subl_channels = Subleading::all;
   std::vector<event_type::EventType> allowed_types{event_type::FKL,
     event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf};
 
   std::cout << "Wp3j" << std::endl;
   // Wp3j
   Process proc {{pid::proton,pid::proton}, 3, pid::Wp};
   size_t n_psp = n_psp_base;
   std::unordered_map<event_type::EventType, size_t> type_counter;
   for( size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_prop, ran};
     if(psp.status()==good){
       const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
       ++type_counter[ev.type()];
       if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
         == allowed_types.cend()) {
         bail_out(psp, "Found not allowed event of type "
           +std::string(event_type::names[ev.type()]));
       }
     } else { // bad process -> try again
       ++n_psp;
     }
   }
   std::cout << "Wp+3j: Took " << n_psp << " to generate "
     << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
   std::cout << "States by classification:\n";
   for(auto const & entry: type_counter){
     const double fraction = static_cast<double>(entry.second)/n_psp_base;
     const int percent = std::round(100*fraction);
     std::cout << std::left << std::setw(25)
               << (event_type::names[entry.first] + std::string(":"))
               << entry.second << " (" << percent << "%)\n";
 
   }
   for(auto const & t: allowed_types){
     if(type_counter[t] < 0.05 * n_psp_base){
       std::cerr << "Less than 5% of the events are of type " << event_type::names[t] << std::endl;
       return EXIT_FAILURE;
     }
   }
 
   // Wm3j - only uno
   proc = Process{{pid::proton,pid::proton}, 3, pid::Wm};
   n_psp = n_psp_base;
 
   subl_channels = Subleading::uno;
   allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
   type_counter.clear();
 
   for( size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_prop, ran};
     if(psp.status()==good){
       const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
       ++type_counter[ev.type()];
       if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
         == allowed_types.cend()) {
         bail_out(psp, "Found not allowed event of type "
           +std::string(event_type::names[ev.type()]));
       }
     } else { // bad process -> try again
       ++n_psp;
     }
   }
   std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate "
     << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
   std::cout << "States by classification:\n";
   for(auto const & entry: type_counter){
     const double fraction = static_cast<double>(entry.second)/n_psp_base;
     const int percent = std::round(100*fraction);
     std::cout << std::left << std::setw(25)
               << (event_type::names[entry.first] + std::string(":"))
               << entry.second << " (" << percent << "%)\n";
 
   }
   for(auto const & t: allowed_types){
     if(type_counter[t] < 0.05 * n_psp_base){
       std::cerr << "Less than 5% of the events are of type " << event_type::names[t] << std::endl;
       return EXIT_FAILURE;
     }
   }
 
   // Wm4j
   proc = Process{{pid::proton,pid::proton}, 4, pid::Wm};
   n_psp = n_psp_base;
 
   subl_channels = Subleading::all;
   allowed_types = {event_type::FKL,
     event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
     event_type::qqxmid};
   type_counter.clear();
 
   for( size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_prop, ran};
     if(psp.status()==good){
       const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
       ++type_counter[ev.type()];
       if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
         == allowed_types.cend()) {
         bail_out(psp, "Found not allowed event of type "
           +std::string(event_type::names[ev.type()]));
       }
     } else { // bad process -> try again
       ++n_psp;
     }
   }
   std::cout << "Wm+4j: Took " << n_psp << " to generate "
     << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
   std::cout << "States by classification:\n";
   for(auto const & entry: type_counter){
     const double fraction = static_cast<double>(entry.second)/n_psp_base;
     const int percent = std::round(100*fraction);
     std::cout << std::left << std::setw(25)
               << (event_type::names[entry.first] + std::string(":"))
               << entry.second << " (" << percent << "%)\n";
 
   }
   for(auto const & t: allowed_types){
     if(type_counter[t] < 0.03 * n_psp_base){
       std::cerr << "Less than 3% of the events are of type " << event_type::names[t] << std::endl;
       return EXIT_FAILURE;
     }
   }
 
   std::cout << "All processes passed." << std::endl;
   return EXIT_SUCCESS;
 }
diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
index 42c8af6..02c1e89 100644
--- a/doc/sphinx/HEJFOG.rst
+++ b/doc/sphinx/HEJFOG.rst
@@ -1,281 +1,291 @@
 The HEJ Fixed Order Generator
 =============================
 
 For high jet multiplicities event generation with standard fixed-order
 generators becomes increasingly cumbersome. For example, the
 leading-order production of a Higgs Boson with five or more jets is
 computationally prohibitively expensive.
 
 To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
 that allows to generate events with high jet multiplicities. To
 facilitate the computation the limit of Multi-Regge Kinematics with
 large invariant masses between all outgoing particles is assumed in the
 matrix elements. The typical use of the ``HEJFOG`` is to supplement
 low-multiplicity events from standard generators with high-multiplicity
 events before using the HEJ 2 program to add high-energy
 resummation.
 
 Installation
 ------------
 
 The ``HEJFOG`` comes bundled together with HEJ 2 and the
 installation is very similar. After downloading HEJ 2 and
 installing the prerequisites as described in :ref:`Installation` the
 ``HEJFOG`` can be installed with::
 
   cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory -DCMAKE_BUILD_TYPE=Release
   make install
 
 where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
 subdirectory in the HEJ 2 directory. If HEJ 2 was
 installed to a non-standard location, it may be necessary to specify the
 directory containing :file:`HEJ-config.cmake`. If the base installation
 directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
 found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
 installing the ``HEJFOG`` would read::
 
   cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory -DCMAKE_BUILD_TYPE=Release
   make install
 
 The installation can be tested with::
 
   make test
 
 provided that the NNPDF 3.0 PDF set is installed.
 
 Running the fixed-order generator
 ---------------------------------
 
 After installing the ``HEJFOG`` you can modify the provided
 configuration file :file:`configFO.yml` and run the generator with::
 
   HEJFOG configFO.yml
 
 The resulting event file, by default named :file:`HEJFO.lhe`, can then be
 fed into HEJ 2 like any event file generated from a standard
 fixed-order generator, see :ref:`Running HEJ 2`.
 
 Settings
 --------
 
 Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
 <http://yaml.org/>`_ configuration file. The settings are
 
 .. _`process`:
 
 **process**
    The scattering process for which events are being generated. The
    format is
 
    :code:`in1 in2 => out1 out2 ...`
 
    The incoming particles, :code:`in1`, :code:`in2` can be
 
    - quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
    - gluons: :code:`g`
    - protons :code:`p` or antiprotons :code:`p_bar`
 
    At most one of the outgoing particles can be a boson, the rest has to be
    partonic. At the moment only the Higgs boson :code:`h` is supported. All
    other outgoing particles are jets. Multiple jets can be grouped together, so
    :code:`p p => h j j` is the same as :code:`p p => h 2j`. There have to be at
    least two jets. Further decays of the boson can be added through the
    :ref:`particle properties<particle properties: particle: decays>`.
 
 .. _`events`:
 
 **events**
    Specifies the number of events to generate.
 
 .. _`jets`:
 
 **jets**
    Defines the properties of the generated jets.
 
    .. _`jets: min pt`:
 
    **min pt**
       Minimum jet transverse momentum in GeV.
 
    .. _`jets: peak pt`:
 
    **peak pt**
       Optional setting to specify the dominant jet transverse momentum
       in GeV. If the generated events are used as input for HEJ
       resummation, this should be set to the minimum transverse momentum
       of the resummation jets. The effect is that only a small fraction
       of jets will be generated with a transverse momentum below the
       value of this setting.
 
    .. _`jets: algorithm`:
 
    **algorithm**
       The algorithm used to define jets. Allowed settings are
       :code:`kt`, :code:`cambridge`, :code:`antikt`,
       :code:`cambridge for passive`. See the `FastJet
       <http://fastjet.fr/>`_ documentation for a description of these
       algorithms.
 
    .. _`jets: R`:
 
    **R**
       The R parameter used in the jet algorithm.
 
    .. _`jets: max rapidity`:
 
    **max rapidity**
       Maximum absolute value of the jet rapidity.
 
 .. _`beam`:
 
 **beam**
    Defines various properties of the collider beam.
 
    .. _`beam: energy`:
 
    **energy**
       The beam energy in GeV. For example, the 13
       TeV LHC corresponds to a value of 6500.
 
    .. _`beam: particles`:
 
    **particles**
       A list :code:`[p1, p2]` of two beam particles. The only supported
       entries are protons :code:`p` and antiprotons :code:`p_bar`.
 
 .. _`pdf`:
 
 **pdf**
    The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
    For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.
 
 .. _`subleading fraction`:
 
 **subleading fraction**
    This setting is related to the fraction of events, that are not a FKL
-   configuration. Currently only unordered emissions are implemented, and only
-   for Higgs plus multijet processes.
-   Typically, this value should be between 0.01 and 0.1.
+   configuration. All possible subleading process are listed in
+   :ref:`subleading channels<subleading channels>`. Typically, this value
+   should be between 0.01 and 0.1.
 
 .. _`subleading channels`:
 
 **subleading channels**
-   Optional parameters to select a specific unordered process, multiple
-   channels can be selected at once. Only matters if :code:`subleading fraction`
-   is non-zero. Currently three values are allowed:
+   Optional parameters to select a specific subleading process, multiple
+   channels can be selected at once. If multiple subleading configurations
+   are possible one will be selected at random for each event, thus each
+   final state will include at most one subleading process, e.g. either
+   :code:`unordered` or :code:`qqx`. The following values are allowed:
 
    - :code:`all`: All channels allowed, default if nothing else set
    - :code:`none`:  No subleading contribution, only FKL allowed
       Equivalent to :code:`subleading fraction: 0`
    - :code:`unordered`: Unordered emission allowed
       Unordered emission are any rapidity ordering, where exactly one gluon is
       emitted outside the rapidity ordering required in FKL events. More
       precisely, if at least one of the incoming particles is a quark or
       antiquark and there are more than two jets in the final state,
       :code:`subleading fraction` states the probability that the flavours
       of the outgoing particles are assigned in such a way that an unordered
-      configuration arises.
+      configuration arises. Unordered emissions are currently not implemented
+      for pure jet production.
+   - :code:`qqx`: Production of quark-antiquark pair
+      If allowed processes with at least three jets might include a
+      quark-antiquark pair, i.e. one quark (in rapidity) next to one antiquark.
+      The pair can be at any position in the FKL chain. Similar to above
+      :code:`subleading fraction` gives the probability of turning two gluons
+      into the quark-antiquark pair. Quark-antiquark pairs are currently only
+      implemented for W boson production.
 
 .. _`unweight`:
 
 **unweight**
    This setting defines the parameters for the partial unweighting of
    events. You can disable unweighting by removing this entry from the
    configuration file.
 
 .. _`unweight: sample size`:
 
    **sample size**
       The number of weighted events used to calibrate the unweighting.
       A good default is to set this to the number of target
       `events`_. If the number of `events`_ is large this can
       lead to significant memory consumption and a lower value should be
       chosen. Contrarily, for large multiplicities the unweighting
       efficiency becomes worse and the sample size should be increased.
 
 .. _`unweight: max deviation`:
 
    **max deviation**
      Controls the range of events to which unweighting is applied. A
      larger value means that a larger fraction of events are unweighted.
      Typical values are between -1 and 1.
 
 .. _`particle properties`:
 
 **particle properties**
    Specifies various properties of the different particles (Higgs, W or Z).
    This is only relevant if the chosen `process`_ is the production of the
    corresponding particles with jets. E.g. for the `process`_
    :code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
    :code:`decays` of the :code:`Higgs` boson are required, while all other
    particle properties will be ignored.
 
 .. _`particle properties: particle`:
 
    **Higgs, Wp, Wm or Z**
       Name of the boson. Can be any of the once above.
 
    .. _`particle properties: particle: mass`:
 
       **mass**
          The mass of the particle in GeV.
 
    .. _`particle properties: particle: width`:
 
       **width**
          The total decay width of the particle in GeV.
 
    .. _`particle properties: particle: decays`:
 
       **decays**
          Optional setting specifying the decays of the particle. Only the decay
          into two particles is implemented. Each decay has the form
 
          :code:`{into: [p1,p2], branching ratio: r}`
 
          where :code:`p1` and :code:`p2` are the particle names of the
          decay product (e.g. :code:`photon`) and :code:`r` is the branching
          ratio.
 
          Decays of a Higgs boson are treated as the production and subsequent
          decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
          supported.
 
 .. _`scales`:
 
 **scales**
    Specifies the renormalisation and factorisation scales for the output
    events. For details, see the corresponding entry in the HEJ 2
    :ref:`HEJ 2 settings`. Note that this should usually be a
    single value, as the weights resulting from additional scale choices
    will simply be ignored by HEJ 2.
 
 .. _`event output`:
 
 **event output**
    Specifies the name of a single event output file or a list of such
    files. See the corresponding entry in the HEJ 2
    :ref:`HEJ 2 settings` for details.
 
 .. _`RanLux init`:
 
 .. _`random generator`:
 
 **random generator**
    Sets parameters for random number generation. See the HEJ 2
    :ref:`HEJ 2 settings` for details.
 
 .. _`analysis`:
 
 **analysis**
    Specifies the name and settings for a custom analysis library. This
    can be useful to specify cuts at the fixed-order level. See the
    corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
    for details.
 
 .. _`Higgs coupling`:
 
 **Higgs coupling**
    This collects a number of settings concerning the effective coupling
    of the Higgs boson to gluons. See the corresponding entry in the
    HEJ 2 :ref:`HEJ 2 settings` for details.