diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc
index aa12ae4..5b033b6 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,205 +1,227 @@
 /**
  *  \brief     check that the PSP generates the all W+jet subleading processes
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
 #include <iomanip>
 #include <iostream>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
 #include "HEJ/Event.hh"
 #include "HEJ/event_types.hh"
 #include "HEJ/EWConstants.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/Mixmax.hh"
 #include "HEJ/PDF.hh"
 #include "HEJ/PDG_codes.hh"
 
 #include "fastjet/JetDefinition.hh"
 
 #include "Decay.hh"
 #include "JetParameters.hh"
 #include "PhaseSpacePoint.hh"
 #include "Process.hh"
 #include "Status.hh"
 #include "Subleading.hh"
 
 namespace {
   using namespace HEJFOG;
   using namespace HEJ;
 
   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};
   }
+
+  #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6)
+  // gcc version < 6 explicitly needs hash function for enum
+  // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970
+  using type_map = std::unordered_map<event_type::EventType, std::size_t, std::hash<std::size_t>>;
+  #else
+  using type_map = std::unordered_map<event_type::EventType, std::size_t>;
+  #endif
 }
 
 int main(){
   constexpr std::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 ParticlesDecayMap boson_decays{
       {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }},
       {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }}
     };
   const EWConstants ew_constants{246.2196508,
     ParticleProperties{80.385, 2.085},
     ParticleProperties{91.187, 2.495},
     ParticleProperties{125, 0.004165}
   };
   HEJ::Mixmax ran{};
 
   Subleading subl_channels(~0l);
   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, {}};
   std::size_t n_psp = n_psp_base;
 
-  #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6)
-  // gcc version < 6 explicitly needs hash function for enum
-  // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970
-  std::unordered_map<event_type::EventType, std::size_t, std::hash<std::size_t>> type_counter;
-  #else
-  std::unordered_map<event_type::EventType, std::size_t> type_counter;
-  #endif
+  type_map type_counter;
 
   for( std::size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==good){
       const Event ev{ to_EventData(psp).cluster(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()) {
+          == allowed_types.cend()) {
         bail_out(psp, "Found not allowed event of type " + name(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)
               << (name(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 " << name(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.reset();
   subl_channels.set(subleading::uno);
   allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
   type_counter.clear();
 
   for( std::size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==good){
       const Event ev{ to_EventData(psp).cluster(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 " + name(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)
               << (name(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 " << name(t) << std::endl;
       return EXIT_FAILURE;
     }
   }
 
   // Wm4j
   proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}};
   n_psp = n_psp_base;
 
   subl_channels.set();
   allowed_types = {event_type::FKL,
     event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
     event_type::qqxmid};
   type_counter.clear();
+  std::array<type_map,3> wpos_counter; // position of the W boson (back, central, forward)
 
   for( std::size_t i = 0; i<n_psp; ++i){
     const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==good){
       const Event ev{ to_EventData(psp).cluster(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()) {
+          == allowed_types.cend()) {
         bail_out(psp, "Found not allowed event of type " + name(ev.type()) );
       }
+      if(ev.outgoing().begin()->type==pid::Wm){
+        ++wpos_counter[0][ev.type()];
+      } else if(ev.outgoing().rbegin()->type==pid::Wm){
+        ++wpos_counter[2][ev.type()];
+      } else {
+        ++wpos_counter[1][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)
               << (name(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 " << name(t) << std::endl;
       return EXIT_FAILURE;
     }
   }
+  std::cout << "Stats by Wpos:\n";
+  for(std::size_t i=0; i<wpos_counter.size(); ++i){
+    std::cout << "position " << i << std::endl;;
+    for(auto const & t: allowed_types){
+      std::cout <<  name(t) << " " << wpos_counter[i][t] << std::endl;
+      if(wpos_counter[i][t] < 0.05 * type_counter[t]){
+        std::cerr << "Less than 5% of the events are of type " << name(t)
+          << " with W at position " << i << std::endl;
+        return EXIT_FAILURE;
+      }
+    }
+  }
 
   std::cout << "All processes passed." << std::endl;
   return EXIT_SUCCESS;
 }