diff --git a/src/Event.cc b/src/Event.cc
index 0c6e97f..3b50745 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,611 +1,611 @@
 
 #include "RHEJ/Event.hh"
 
 #include "RHEJ/utility.hh"
 
 #include "RHEJ/qqx.hh"
 
 namespace RHEJ{
 
   namespace{
     constexpr int status_in = -1;
     constexpr int status_decayed = 2;
     constexpr int status_out = 1;
 
     // helper functions to determine event type
 
     // check if there is at most one photon, W, H, Z in the final state
     // and all the rest are quarks or gluons
     bool final_state_ok(std::vector<Particle> const & outgoing){
       bool has_AWZH_boson = false;
       for(auto const & out: outgoing){
         if(is_AWZH_boson(out.type)){
           if(has_AWZH_boson) return false;
           has_AWZH_boson = true;
         }
         else if(! is_parton(out.type)) return false;
       }
       return true;
     }
 
     template<class Iterator>
     Iterator remove_AWZH(Iterator begin, Iterator end){
       return std::remove_if(
           begin, end, [](Particle const & p){return is_AWZH_boson(p);}
       );
     }
 
     template<class Iterator>
     bool valid_outgoing(Iterator begin, Iterator end){
       return std::distance(begin, end) >= 2
         && std::is_sorted(begin, end, rapidity_less{})
         && std::count_if(
             begin, end, [](Particle const & s){return is_AWZH_boson(s);}
         ) < 2;
     }
 
     /**
      * \brief function which determines if type change is consistent with W emission.
      * @param in                      incoming Particle
      * @param out                     outgoing Particle
      *
      * Ensures that change type of quark line is possible by a flavour changing
      * W emission.
      */
     bool is_W_Current(ParticleID in, ParticleID out){
       if((in==1 && out==2)||(in==2 && out==1)){
         return true;
       }
       else if((in==-1 && out==-2)||(in==-2 && out==-1)){
         return true;
       }
       else if((in==3 && out==4)||(in==4 && out==3)){
         return true;
       }
       else if((in==-3 && out==-4)||(in==-4 && out==-3)){
         return true;
       }
       else{
         return false;
       }
     }
 
     /**
      * \brief checks if particle type remains same from incoming to outgoing
      * @param in                      incoming Particle
      * @param out                     outgoing Particle
      */
     bool is_Pure_Current(ParticleID in, ParticleID out){
       if(abs(in)<=6 || in==21) return (in==out);
       else return false;
     }
 
 
     // Note that this changes the outgoing range!
     template<class ConstIterator, class Iterator>
     bool is_FKL(
         ConstIterator begin_incoming, ConstIterator end_incoming,
         Iterator begin_outgoing, Iterator end_outgoing
     ){
       assert(std::distance(begin_incoming, end_incoming) == 2);
       assert(std::distance(begin_outgoing, end_outgoing) >= 2);
 
       // One photon, W, H, Z in the final state is allowed.
       // Remove it for remaining tests,
       end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
 
       if(std::all_of(
             begin_outgoing + 1, end_outgoing - 1,
             [](Particle const & p){ return p.type == pid::gluon; })
       ){
         // Test if this is a standard FKL configuration.
         if (is_Pure_Current(begin_incoming->type, begin_outgoing->type)
             && is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
           return true;
         }
         else if(is_W_Current(begin_incoming->type, begin_outgoing->type)
             && is_Pure_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
           return true;
         }
         else if(is_Pure_Current(begin_incoming->type, begin_outgoing->type)
             && is_W_Current((end_incoming-1)->type, (end_outgoing-1)->type)){
           return true;
         }
       }
       return false;
     }
 
     bool is_FKL(
         std::array<Particle, 2> const & incoming,
         std::vector<Particle> outgoing
     ){
       assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
       assert(valid_outgoing(begin(outgoing), end(outgoing)));
 
       return is_FKL(
           begin(incoming), end(incoming),
           begin(outgoing), end(outgoing)
       );
     }
 
     bool has_2_jets(Event const & event){
       return event.jets().size() >= 2;
     }
 
     /**
      * \brief Checks whether event is unordered backwards
      * @param ev        Event
      * @returns         Is Event Unordered Backwards
      *
      * Checks there is more than 3 constuents in the final state
      * Checks there is more than 3 jets
      * Checks the most backwards parton is a gluon
      * Checks the most forwards jet is not a gluon
      * Checks the rest of the event is FKL
      * Checks the second most backwards is not a different boson
      * Checks the unordered gluon actually forms a jet
      */
     bool is_unordered_backward(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
       if(in.front().type == pid::gluon) return false;
       if(out.front().type != pid::gluon) return false;
       // When skipping the unordered emission
       // the remainder should be a regular FKL event,
       // except that the (new) first outgoing particle must not be a A,W,Z,H.
       const auto FKL_begin = next(begin(out));
       if(is_AWZH_boson(*FKL_begin)) return false;
       if(!is_FKL(in, {FKL_begin, end(out)})) return false;
       // check that the unordered gluon forms an extra jet
       const auto jets = sorted_by_rapidity(ev.jets());
       const auto indices = ev.particle_jet_indices({jets.front()});
       return indices[0] >= 0 && indices[1] == -1;
     }
 
     /**
      * \brief Checks for a forward unordered gluon emission
      * @param ev          Event
      * @returns           Is the event a forward unordered emission
      *
      * \see is_unordered_backward
      */
     bool is_unordered_forward(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
       if(in.back().type == pid::gluon) return false;
       if(out.back().type != pid::gluon) return false;
       // When skipping the unordered emission
       // the remainder should be a regular FKL event,
       // except that the (new) last outgoing particle must not be a A,W,Z,H.
       const auto FKL_end = prev(end(out));
       if(is_AWZH_boson(*prev(FKL_end))) return false;
       if(!is_FKL(in, {begin(out), FKL_end})) return false;
       // check that the unordered gluon forms an extra jet
       const auto jets = sorted_by_rapidity(ev.jets());
       const auto indices = ev.particle_jet_indices({jets.back()});
       return indices.back() >= 0 && indices[indices.size()-2] == -1;
     }
 
 
     /**
      * \brief Checks for a forward extremal qqx
      * @param ev          Event
      * @returns           Is the event a forward extremal qqx event
      *
      * Checks there is 3 or more than 3 constituents in the final state
      * Checks there is 3 or more than 3 jets
      * Checks most forwards incoming is gluon
      * Checks most extremal particle is not a Higgs (either direction)
      * Checks the second most forwards particle is not Higgs boson
      * Checks the most forwards parton is a either quark or anti-quark.
      * Checks the second most forwards parton is anti-quark or quark.
      */
     bool is_Ex_qqxf(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       int fkl_end=2;
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
-      if(in.back().type == pid::gluon) return false;
+      if(in.back().type != pid::gluon) return false;
       if(out.back().type == pid::Higgs || out.front().type == pid::Higgs
                              || out.rbegin()[1].type == pid::Higgs) return false;
 
       // if extremal AWZ
       if(is_AWZ_boson(out.back().type)){ // if extremal AWZ
         fkl_end++;
         if (is_quark(out.rbegin()[1].type)){ //if second quark
           if (!(is_antiquark(out.rbegin()[2].type))) return false;// third must be anti-quark
         } // end second quark
 
         else if (is_antiquark(out.rbegin()[1].type)){ //if second anti-quark
           if (!(is_quark(out.rbegin()[2].type))) return false;// third must be quark
           }
         else return false;
       } // end extremal AWZ
 
       // if extremal quark
         else if (is_quark(out.rbegin()[0].type)){ //if extremal quark
           if(is_AWZ_boson(out.rbegin()[1].type)){ // if second AWZ
             fkl_end++;
             if (!(is_antiquark(out.rbegin()[2].type))) return false;// third must be anti-quark
         } //end second AWZ
         else if (!(is_antiquark(out.rbegin()[1].type))) return false;// second must be anti-quark
         } // end extremal quark
 
         // if extremal anti-quark
         else if (is_antiquark(out.rbegin()[0].type)){ //if extremal anti-quark
           if(is_AWZ_boson(out.rbegin()[1].type)){ // if second AWZ
             fkl_end++;
             if (!(is_quark(out.rbegin()[2].type))) return false;// third must be quark
         } //end second AWZ
         else if (!(is_quark(out.rbegin()[1].type))) return false;// second must be quark
         } // end extremal antiquark
 
       // When skipping the qqbar
       // New last outgoing particle must not be a Higgs
       if (out.rbegin()[fkl_end].type == pid::Higgs) return false;
 
       // Opposite current should be logical to process
       if (is_AWZ_boson(out.front().type)){
         return (is_Pure_Current(in.front().type, out[1].type)
                 || is_W_Current(in.front().type,out[1].type));
       }
       else
         return (is_Pure_Current(in.front().type, out[0].type)
                 || is_W_Current(in.front().type,out[0].type));
     }
 
     /**
      * \brief Checks for a backward extremal qqx
      * @param ev          Event
      * @returns           Is the event a backward extremal qqx event
      *
      * Checks there is 3 or more than 3 constituents in the final state
      * Checks there is 3 or more than 3 jets
      * Checks most backwards incoming is gluon
      * Checks most extremal particle is not a Higgs (either direction) y
      * Checks the second most backwards particle is not Higgs boson y
      * Checks the most backwards parton is a either quark or anti-quark. y
      * Checks the second most backwards parton is anti-quark or quark. y
      */
     bool is_Ex_qqxb(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       int fkl_start=2;
 
       if(out.size() < 3) return false;
       if(ev.jets().size() < 3) return false;
-      if(in.front().type == pid::gluon) return false;
+      if(in.front().type != pid::gluon) return false;
       if(out.back().type == pid::Higgs || out.front().type == pid::Higgs
                             || out[1].type == pid::Higgs) return false;
 
       // if extremal AWZ
       if(is_AWZ_boson(out.front().type)){ // if extremal AWZ
         fkl_start++;
         if (is_quark(out[1].type)){ //if second quark
           if (!(is_antiquark(out[2].type))) return false;// third must be anti-quark
         } // end second quark
 
         else if (is_antiquark(out[1].type)){ //if second anti-quark
           if (!(is_quark(out[2].type))) return false;// third must be quark
           }
         else return false;
       } // end extremal AWZ
 
       // if extremal quark
         else if (is_quark(out[0].type)){ // if extremal quark
           if(is_AWZ_boson(out[1].type)){ // if second AWZ
             fkl_start++;
             if (!(is_antiquark(out[2].type))) return false;// third must be anti-quark
         } //end second AWZ
         else if (!(is_antiquark(out[1].type))) return false;// second must be anti-quark
         } // end extremal quark
 
         // if extremal anti-quark
         else if (is_antiquark(out[0].type)){ //if extremal anti-quark
           if(is_AWZ_boson(out[1].type)){ // if second AWZ
             fkl_start++;
             if (!(is_quark(out[2].type))) return false;// third must be quark
         } //end second AWZ
         else if (!(is_quark(out[1].type))) return false;// second must be quark
         } // end extremal antiquark
 
       // When skipping the qqbar
       // New last outgoing particle must not be a Higgs.
       if (out[fkl_start].type == pid::Higgs) return false;
 
       // Other current should be logical to process
       if (is_AWZ_boson(out.back().type)){
         return (is_Pure_Current(in.back().type, out.rbegin()[1].type)
                 || is_W_Current(in.back().type,out.rbegin()[1].type));
       }
       else
         return (is_Pure_Current(in.back().type, out.rbegin()[0].type)
                 || is_W_Current(in.back().type, out.rbegin()[0].type));
     }
 
 
     /**
      * \brief Checks for a central qqx
      * @param ev          Event
      * @returns           Is the event a central extremal qqx event
      *
      * Checks there is 4 or more than 4 constuents in the final state
      * Checks there is 4 or more than 4 jets
      * Checks most extremal particle is not a Higgs (either direction) y
      * Checks for a central quark in the outgoing states
      * Checks for adjacent anti-quark parton. (allowing for AWZ boson emission between)
      * Checks external currents are logically sound.
      */
     bool is_Mid_qqx(Event const & ev){
       auto const & in = ev.incoming();
       auto const & out = ev.outgoing();
       assert(std::is_sorted(begin(in), end(in), pz_less{}));
       assert(valid_outgoing(begin(out), end(out)));
 
       if(out.size() < 4) return false;
       if(ev.jets().size() < 4) return false;
 
       if(out.back().type == pid::Higgs || out.front().type == pid::Higgs) return false;
 
       int start_FKL=0;
       int end_FKL=0;
 
       if (is_AWZ_boson(out.back())){
         end_FKL++;
       }
       if (is_AWZ_boson(out.front())){
         start_FKL++;
       }
 
       if ((is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
            && is_Pure_Current(in.front().type,out[start_FKL].type))){
         //nothing to do
       }
       else if (is_W_Current(in.back().type,out.rbegin()[end_FKL].type)
                && is_Pure_Current(in.front().type,out[start_FKL].type)){
         //nothing to do
       }
       else if (!(is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
                  && is_W_Current(in.front().type,out[start_FKL].type))){
         return false;
       }
 
       for(auto i=1+start_FKL; i<out.size()-1-end_FKL;i++){
           if (is_quark(out[i].type)){
             if ((is_antiquark(out[i-1].type) && i!=1)
                 || (is_antiquark(out[i+1]) && i!=out.size()-1-end_FKL)){
               return true;}
             else if (is_AWZ_boson(out[i-1]) && (is_antiquark(out[i-2].type) && i!=2) ){
               return true;}
             else if (is_AWZ_boson(out[i+1]) && (is_antiquark(out[i+2].type) && i!=out.size()-2) ){
               return true;}
           }
       }
       return false;
   }
 
     using event_type::EventType;
 
     EventType classify(Event const & ev){
       if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state;
       if(! has_2_jets(ev)) return EventType::no_2_jets;
       if(is_FKL(ev.incoming(), ev.outgoing())) {
         return EventType::FKL;
       }
       if(is_unordered_backward(ev)){
         return EventType::unordered_backward;
       }
       if(is_unordered_forward(ev)){
         return EventType::unordered_forward;
       }
       if(is_Ex_qqxb(ev)){
         return EventType::extremal_qqxb;
       }
       if(is_Ex_qqxf(ev)){
         return EventType::extremal_qqxf;
       }
       if(is_Mid_qqx(ev)){
         return EventType::central_qqx;
       }
       return EventType::nonHEJ;
     }
 
     Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
       return Particle{
         static_cast<ParticleID>(hepeup.IDUP[i]),
         fastjet::PseudoJet{
           hepeup.PUP[i][0], hepeup.PUP[i][1],
           hepeup.PUP[i][2], hepeup.PUP[i][3]
         }
       };
     }
 
     bool is_decay_product(std::pair<int, int> const & mothers){
       if(mothers.first == 0) return false;
       return mothers.second == 0 || mothers.first == mothers.second;
     }
 
   }
 
   UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup):
     central(EventParameters{
       hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
     })
   {
     size_t in_idx = 0;
     for (int i = 0; i < hepeup.NUP; ++i) {
       // skip decay products
       // we will add them later on, but we have to ensure that
       // the decayed particle is added before
       if(is_decay_product(hepeup.MOTHUP[i])) continue;
 
       auto particle = extract_particle(hepeup, i);
       // needed to identify mother particles for decay products
       particle.p.set_user_index(i+1);
 
       if(hepeup.ISTUP[i] == status_in){
          if(in_idx > incoming.size()) {
            throw std::invalid_argument{
              "Event has too many incoming particles"
            };
          }
          incoming[in_idx++] = std::move(particle);
       }
       else outgoing.emplace_back(std::move(particle));
     }
     std::sort(
         begin(incoming), end(incoming),
         [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
     );
     std::sort(begin(outgoing), end(outgoing), rapidity_less{});
 
     // add decay products
     for (int i = 0; i < hepeup.NUP; ++i) {
       if(!is_decay_product(hepeup.MOTHUP[i])) continue;
       const int mother_id = hepeup.MOTHUP[i].first;
       const auto mother = std::find_if(
           begin(outgoing), end(outgoing),
           [mother_id](Particle const & particle){
             return particle.p.user_index() == mother_id;
           }
       );
       if(mother == end(outgoing)){
         throw std::invalid_argument{"invalid decay product parent"};
       }
       const int mother_idx = std::distance(begin(outgoing), mother);
       assert(mother_idx >= 0);
       decays[mother_idx].emplace_back(extract_particle(hepeup, i));
     }
   }
 
   Event::Event(
       UnclusteredEvent ev,
       fastjet::JetDefinition const & jet_def, double min_jet_pt
   ):
     ev_{std::move(ev)},
     cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
     min_jet_pt_{min_jet_pt}
   {
     type_ = classify(*this);
   }
 
   std::vector<fastjet::PseudoJet> Event::jets() const{
     return cs_.inclusive_jets(min_jet_pt_);
   }
 
   /**
    * \brief Returns the invarient mass of the event
    * @param ev         Event
    * @returns          s hat
    *
    * Makes use of the FastJet PseudoJet function m2().
    * Applies this function to the sum of the incoming partons.
    */
   double shat(Event const & ev){
     return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
   }
 
   namespace{
     // colour flow according to Les Houches standard
     // TODO: stub
     std::vector<std::pair<int, int>> colour_flow(
         std::array<Particle, 2> const & incoming,
         std::vector<Particle> const & outgoing
     ){
       std::vector<std::pair<int, int>> result(
           incoming.size() + outgoing.size()
       );
       for(auto & col: result){
         col = std::make_pair(-1, -1);
       }
       return result;
     }
   }
 
   LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
     LHEF::HEPEUP result;
     result.heprup = heprup;
     result.weights = {{event.central().weight, nullptr}};
     for(auto const & var: event.variations()){
       result.weights.emplace_back(var.weight, nullptr);
     }
     size_t num_particles = event.incoming().size() + event.outgoing().size();
     for(auto const & decay: event.decays()) num_particles += decay.second.size();
     result.NUP = num_particles;
     // the following entries are pretty much meaningless
     result.IDPRUP = event.type()+1;  // event ID
     result.AQEDUP = 1./128.;  // alpha_EW
     //result.AQCDUP = 0.118 // alpha_QCD
     // end meaningless part
     result.XWGTUP = event.central().weight;
     result.SCALUP = event.central().muf;
     result.scales.muf = event.central().muf;
     result.scales.mur = event.central().mur;
     result.scales.SCALUP = event.central().muf;
     result.pdfinfo.p1 = event.incoming().front().type;
     result.pdfinfo.p2 = event.incoming().back().type;
     result.pdfinfo.scale = event.central().muf;
     for(Particle const & in: event.incoming()){
       result.IDUP.emplace_back(in.type);
       result.ISTUP.emplace_back(status_in);
       result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
       result.MOTHUP.emplace_back(0, 0);
     }
     for(size_t i = 0; i < event.outgoing().size(); ++i){
       Particle const & out = event.outgoing()[i];
       result.IDUP.emplace_back(out.type);
       const int status = event.decays().count(i)?status_decayed:status_out;
       result.ISTUP.emplace_back(status);
       result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
       result.MOTHUP.emplace_back(1, 2);
     }
     result.ICOLUP = colour_flow(
         event.incoming(), filter_partons(event.outgoing())
     );
     if(result.ICOLUP.size() < num_particles){
       const size_t AWZH_boson_idx = std::find_if(
           begin(event.outgoing()), end(event.outgoing()),
           [](Particle const & s){ return is_AWZH_boson(s); }
       ) - begin(event.outgoing()) + event.incoming().size();
       assert(AWZH_boson_idx <= result.ICOLUP.size());
       result.ICOLUP.insert(
           begin(result.ICOLUP) + AWZH_boson_idx,
           std::make_pair(0,0)
       );
     }
     for(auto const & decay: event.decays()){
       for(auto const out: decay.second){
         result.IDUP.emplace_back(out.type);
         result.ISTUP.emplace_back(status_out);
         result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
         const int mother_idx = 1 + event.incoming().size() + decay.first;
         result.MOTHUP.emplace_back(mother_idx, mother_idx);
         result.ICOLUP.emplace_back(0,0);
       }
     }
     assert(result.ICOLUP.size() == num_particles);
     static constexpr double unknown_spin = 9.;     //per Les Houches accord
     result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
     result.SPINUP = result.VTIMUP;
     return result;
   }
 
 }