diff --git a/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh
index f6b5934..8848eeb 100644
--- a/FixedOrderGen/include/EventGenerator.hh
+++ b/FixedOrderGen/include/EventGenerator.hh
@@ -1,68 +1,68 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <memory>
 
 #include "HEJ/EWConstants.hh"
 #include "HEJ/MatrixElement.hh"
 #include "HEJ/optional.hh"
 #include "HEJ/PDF.hh"
 #include "HEJ/ScaleFunction.hh"
 
 #include "Beam.hh"
 #include "Decay.hh"
 #include "JetParameters.hh"
 #include "Process.hh"
 #include "Status.hh"
 #include "Subleading.hh"
 
 namespace HEJ {
   class Event;
   class HiggsCouplingSettings;
   class RNG;
 }
 
 //! Namespace for HEJ Fixed Order Generator
 namespace HEJFOG {
   class EventGenerator{
   public:
     EventGenerator(
         Process process,
         Beam beam,
         HEJ::ScaleGenerator scale_gen,
         JetParameters jets,
         int pdf_id,
-        double subl_change,
+        double subl_chance,
         Subleading subl_channels,
         ParticlesDecayMap particle_decays,
         HEJ::HiggsCouplingSettings Higgs_coupling,
         HEJ::EWConstants ew_parameters,
         std::shared_ptr<HEJ::RNG> ran
     );
 
     HEJ::optional<HEJ::Event> gen_event();
 
     Status status() const {
       return status_;
     }
 
   private:
     HEJ::PDF pdf_;
     HEJ::MatrixElement ME_;
     HEJ::ScaleGenerator scale_gen_;
     Process process_;
     JetParameters jets_;
     Beam beam_;
     Status status_;
-    double subl_change_;
+    double subl_chance_;
     Subleading subl_channels_;
     ParticlesDecayMap particle_decays_;
     HEJ::EWConstants ew_parameters_;
     std::shared_ptr<HEJ::RNG> ran_;
   };
 
 }
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index bb1c45d..f35c164 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,282 +1,282 @@
 /** \file      PhaseSpacePoint.hh
  *  \brief     Contains the PhaseSpacePoint Class
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <array>
 #include <bitset>
 #include <cstddef>
 #include <unordered_map>
 #include <vector>
 
 #include "boost/iterator/filter_iterator.hpp"
 
 #include "HEJ/Event.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PDG_codes.hh"
 
 #include "fastjet/PseudoJet.hh"
 
 #include "Decay.hh"
 #include "Status.hh"
 #include "Subleading.hh"
 
 namespace HEJ {
   class EWConstants;
   class PDF;
   class RNG;
 }
 
 namespace HEJFOG {
   class JetParameters;
   class Process;
 
   //! A point in resummation phase space
   class PhaseSpacePoint{
   public:
 
     //! Iterator over partons
     using ConstPartonIterator = boost::filter_iterator<
       bool (*)(HEJ::Particle const &),
       std::vector<HEJ::Particle>::const_iterator
       >;
     //! Reverse Iterator over partons
     using ConstReversePartonIterator = std::reverse_iterator<
                                                 ConstPartonIterator>;
 
     //! Default PhaseSpacePoint Constructor
     PhaseSpacePoint() = delete;
 
     //! PhaseSpacePoint Constructor
     /**
      * @param proc                  The process to generate
      * @param jet_properties        Jet defintion & cuts
      * @param pdf                   The pdf set (used for sampling)
      * @param E_beam                Energie of the beam
      * @param subl_chance           Chance to turn a potentially unordered
      *                              emission into an actual one
      * @param subl_channels         Possible subleading channels.
      *                              see HEJFOG::Subleading
      * @param particle_properties   Properties of producted boson
      *
      * Initially, only FKL phase space points are generated. subl_chance gives
-     * the change of turning one emissions into a subleading configuration,
+     * the chance of turning one emissions into a subleading configuration,
      * i.e. either unordered or central quark/anti-quark pair. Unordered
      * emissions require that the most extremal emission in any direction is
      * a quark or anti-quark and the next emission is a gluon. Quark/anti-quark
      * pairs are only generated for W processes. At most one subleading
      * emission will be generated in this way.
      */
     PhaseSpacePoint(
       Process const & proc,
       JetParameters const & jet_properties,
       HEJ::PDF & pdf, double E_beam,
       double subl_chance,
       Subleading subl_channels,
       ParticlesDecayMap const & particle_decays,
       HEJ::EWConstants const & ew_parameters,
       HEJ::RNG & ran
     );
 
     //! Get Weight Function
     /**
      * @returns        Weight of Event
      */
     double weight() const{
       return weight_;
     }
 
     Status status() const{
       return status_;
     }
 
     //! Get Incoming Function
     /**
      * @returns        Incoming Particles
      */
     std::array<HEJ::Particle, 2> const & incoming() const{
       return incoming_;
     }
 
     //! Get Outgoing Function
     /**
      * @returns        Outgoing Particles
      */
     std::vector<HEJ::Particle> const & outgoing() const{
       return outgoing_;
     }
 
     std::unordered_map<std::size_t, std::vector<HEJ::Particle>> const & decays() const{
       return decays_;
     }
 
     //! Iterator to the first outgoing parton
     ConstPartonIterator begin_partons() const;
     //! Iterator to the first outgoing parton
     ConstPartonIterator cbegin_partons() const;
 
     //! Iterator to the end of the outgoing partons
     ConstPartonIterator end_partons() const;
     //! Iterator to the end of the outgoing partons
     ConstPartonIterator cend_partons() const;
 
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator rbegin_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator crbegin_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator rend_partons() const;
     //! Reverse Iterator to the first outgoing parton
     ConstReversePartonIterator crend_partons() const;
 
   private:
     friend HEJ::Event::EventData to_EventData(PhaseSpacePoint psp);
     /**
      * @internal
      * @brief Generate LO parton momentum
      *
      * @param count             Number of partons to generate
      * @param is_pure_jets      If true ensures momentum conservation in x and y
      * @param jet_param         Jet properties to fulfil
      * @param max_pt            max allowed pt for a parton (typically E_CMS)
      * @param ran               Random Number Generator
      *
      * @returns                 Momentum of partons
      *
      * Ensures that each parton is in its own jet.
      * Generation is independent of parton flavour. Output is sorted in rapidity.
      */
     std::vector<fastjet::PseudoJet> gen_LO_partons(
         int count, bool is_pure_jets,
         JetParameters const & jet_param,
         double max_pt,
         HEJ::RNG & ran
     );
     HEJ::Particle gen_boson(
         HEJ::ParticleID bosonid, double mass, double width,
         HEJ::RNG & ran
     );
     template<class ParticleMomenta>
     fastjet::PseudoJet gen_last_momentum(
         ParticleMomenta const & other_momenta,
         double mass_square, double y
     ) const;
 
     bool jets_ok(
         std::vector<fastjet::PseudoJet> const & Born_jets,
         std::vector<fastjet::PseudoJet> const & partons
     ) const;
     /**
      * @internal
      * @brief Generate incoming partons according to the PDF
      *
      * @param uf                Scale used in the PDF
      */
     void reconstruct_incoming(
         Process const & proc,
         double subl_chance, Subleading subl_channels,
         HEJ::PDF & pdf, double E_beam,
         double uf,
         HEJ::RNG & ran
     );
     HEJ::ParticleID generate_incoming_id(
         std::size_t beam_idx, double x, double uf, HEJ::PDF & pdf,
         std::bitset<11> allowed_partons, HEJ::RNG & ran
     );
     //! @brief Returns list of all allowed initial states partons
     std::array<std::bitset<11>,2> filter_partons(
         Process const & proc,
         double const subl_chance, Subleading const subl_channels,
         HEJ::RNG & ran
     );
 
     bool momentum_conserved(double ep) const;
 
     bool extremal_FKL_ok(
         std::vector<fastjet::PseudoJet> const & partons
     ) const;
     double random_normal(double stddev, HEJ::RNG & ran);
     /**
      * @internal
      * @brief Turns a FKL configuration into a subleading one
      *
-     * @param chance            Change to switch to subleading configuration
+     * @param chance            Chance to switch to subleading configuration
      * @param channels          Allowed channels for subleading process
      * @param proc              Process to decide which subleading
      *                          configurations are allowed
      *
      * With a chance of "chance" the FKL configuration is either turned into
      * a unordered configuration or, for A/W/Z bosons, a configuration with
      * a central quark/anti-quark pair.
      */
     void maybe_turn_to_subl(double chance, Subleading channels,
         Process const & proc, HEJ::RNG & ran);
     void turn_to_uno(bool can_be_uno_backward, bool can_be_uno_forward,
                      HEJ::RNG & ran);
     int select_qqx_flavour(bool allow_strange, HEJ::RNG & ran);
     void turn_to_qqx(bool allow_strange, HEJ::RNG & ran);
     void turn_to_cqqx(bool allow_strange, HEJ::RNG & ran);
     void turn_to_eqqx(bool allow_strange, HEJ::RNG & ran);
     //! decay where we select the decay channel
     std::vector<HEJ::Particle> decay_boson(
         HEJ::Particle const & parent,
         std::vector<Decay> const & decays,
         HEJ::RNG & ran
     );
     //! generate decay products of a boson
     std::vector<HEJ::Particle> decay_boson(
         HEJ::Particle const & parent,
         std::vector<HEJ::ParticleID> const & decays,
         HEJ::RNG & ran
     );
     /// @brief setup outgoing partons to ensure correct coupling to boson
     void couple_boson(HEJ::ParticleID boson, HEJ::RNG & ran);
     Decay select_decay_channel(
         std::vector<Decay> const & decays,
         HEJ::RNG & ran
     );
     double gen_hard_pt(
         int np, double ptmin, double ptmax, double y,
         HEJ::RNG & ran
     );
     double gen_soft_pt(int np, double ptmax, HEJ::RNG & ran);
     double gen_parton_pt(
         int count, JetParameters const & jet_param, double ptmax, double y,
         HEJ::RNG & ran
     );
 
     //! Iterator over partons (non-const)
     using PartonIterator = boost::filter_iterator<
       bool (*)(HEJ::Particle const &),
       std::vector<HEJ::Particle>::iterator
       >;
     //! Reverse Iterator over partons (non-const)
     using ReversePartonIterator = std::reverse_iterator<PartonIterator>;
 
     //! Iterator to the first outgoing parton (non-const)
     PartonIterator begin_partons();
     //! Iterator to the end of the outgoing partons (non-const)
     PartonIterator end_partons();
 
     //! Reverse Iterator to the first outgoing parton (non-const)
     ReversePartonIterator rbegin_partons();
     //! Reverse Iterator to the first outgoing parton (non-const)
     ReversePartonIterator rend_partons();
 
     double weight_;
 
     Status status_;
 
     std::array<HEJ::Particle, 2> incoming_;
     std::vector<HEJ::Particle> outgoing_;
     //! Particle decays in the format {outgoing index, decay products}
     std::unordered_map<std::size_t, std::vector<HEJ::Particle>> decays_;
   };
 
   //! Extract HEJ::Event::EventData from PhaseSpacePoint
   HEJ::Event::EventData to_EventData(PhaseSpacePoint psp);
 }
diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
index a7b5772..52bb69b 100644
--- a/FixedOrderGen/src/EventGenerator.cc
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -1,101 +1,101 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "EventGenerator.hh"
 
 #include <cstddef>
 #include <utility>
 
 #include "HEJ/Config.hh"
 #include "HEJ/Event.hh"
 #include "HEJ/event_types.hh"
 #include "HEJ/EWConstants.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/HiggsCouplingSettings.hh"
 
 #include "Process.hh"
 #include "Beam.hh"
 #include "JetParameters.hh"
 #include "PhaseSpacePoint.hh"
 
 namespace HEJFOG {
   EventGenerator::EventGenerator(
       Process process,
       Beam beam,
       HEJ::ScaleGenerator scale_gen,
       JetParameters jets,
       int pdf_id,
-      double subl_change,
+      double subl_chance,
       Subleading subl_channels,
       ParticlesDecayMap particle_decays,
       HEJ::HiggsCouplingSettings Higgs_coupling,
       HEJ::EWConstants ew_parameters,
       std::shared_ptr<HEJ::RNG> ran
   ):
     pdf_{pdf_id, beam.particles[0], beam.particles[1]},
     ME_{
       [this](double mu){ return pdf_.Halphas(mu); },
       HEJ::MatrixElementConfig{
         false,
         std::move(Higgs_coupling),
         ew_parameters
       }
     },
     scale_gen_{std::move(scale_gen)},
     process_{std::move(process)},
     jets_{std::move(jets)},
     beam_{std::move(beam)},
-    subl_change_{subl_change},
+    subl_chance_{subl_chance},
     subl_channels_{subl_channels},
     particle_decays_{std::move(particle_decays)},
     ew_parameters_{ew_parameters},
     ran_{std::move(ran)}
   {
   }
 
   HEJ::optional<HEJ::Event> EventGenerator::gen_event(){
     HEJFOG::PhaseSpacePoint psp{
       process_,
       jets_,
       pdf_, beam_.energy,
-      subl_change_, subl_channels_,
+      subl_chance_, subl_channels_,
       particle_decays_,
       ew_parameters_,
       *ran_
     };
     status_ = psp.status();
     if(status_ != Status::good) return {};
 
     HEJ::Event ev = scale_gen_(
         HEJ::Event{
           to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt)
         }
     );
     if(!is_resummable(ev.type())){
       std::cerr << ev << std::endl;
       throw HEJ::not_implemented("Tried to generate a event type, "
         "which is not yet implemented in HEJ.");
     }
 
     ev.generate_colours(*ran_);
 
     const double shat = HEJ::shat(ev);
     const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
     const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);
 
     // evaluate matrix element
     ev.parameters() *= ME_.tree(ev)/(shat*shat);
     // and PDFs
     ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
     ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
     for(std::size_t i = 0; i < ev.variations().size(); ++i){
       auto & var = ev.variations(i);
       var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
       var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
     }
     return ev;
   }
 
 }
diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc
index 5ea67a8..c9adfd2 100644
--- a/FixedOrderGen/t/W_2j_classify.cc
+++ b/FixedOrderGen/t/W_2j_classify.cc
@@ -1,165 +1,165 @@
 /**
  * \brief      check that the PSP generates only "valid" W + 2 jets events
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 
 #include <cstdlib>
 #include <iostream>
 #include <string>
 #include <vector>
 
 #include "fastjet/JetDefinition.hh"
 
 #include "HEJ/EWConstants.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/Mixmax.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PDF.hh"
 #include "HEJ/PDG_codes.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};
   }
 
   bool is_up_type(Particle const & part){
     return HEJ::is_anyquark(part) && !(std::abs(part.type)%2);
   }
   bool is_down_type(Particle const & part){
     return HEJ::is_anyquark(part) && std::abs(part.type)%2;
   }
 
   bool check_W2j(PhaseSpacePoint const & psp, ParticleID const W_type){
     bool found_quark = false;
     bool found_anti = false;
     std::vector<Particle> out_partons;
     std::vector<Particle> Wp;
     for(auto const & p: psp.outgoing()){
       if(p.type == W_type) Wp.push_back(p);
       else if(is_parton(p)) out_partons.push_back(p);
       else bail_out(psp, "Found particle with is not "
             +std::to_string(int(W_type))+" or parton");
     }
     if(Wp.size() != 1 || out_partons.size() != 2){
       bail_out(psp, "Found wrong number of outgoing partons");
     }
     for(std::size_t j=0; j<2; ++j){
       auto const & in = psp.incoming()[j];
       auto const & out = out_partons[j];
       if(is_quark(in) || is_antiquark(in)) {
         found_quark = true;
         if(in.type != out.type) { // switch in quark type -> Wp couples to it
           if(found_anti){ // already found qq for coupling to W
             bail_out(psp, "Found second up/down pair");
           } else if(std::abs(in.type)>4 || std::abs(out.type)>4){
             bail_out(psp, "Found bottom/top pair");
           }
           found_anti = true;
           if( is_up_type(in)) { // "up" in
             if(W_type > 0){
               // -> only allowed u -> Wp + d
               if(in.type < 0 || is_up_type(out) || out.type < 0)
                 bail_out(psp, "u -/> Wp + d");
 
             } else {
               // -> only allowed ux -> Wm + dx
               if(in.type > 0 || is_up_type(out) || out.type > 0)
                 bail_out(psp, "ux -/> Wm + dx");
             }
           } else { // "down" in
             if(W_type > 0){
               // -> only allowed dx -> Wp + ux
               if(in.type > 0 || is_down_type(out) || out.type > 0)
                 bail_out(psp, "dx -/> Wp + ux");
             } else {
               // -> only allowed d -> Wm + u
               if(in.type < 0 || is_down_type(out) || out.type < 0)
                 bail_out(psp, "d -/> Wm + u");
             }
           }
         }
       }
     }
     if(!found_quark) {
       bail_out(psp, "Found no initial quarks");
     } else if(!found_anti){
       bail_out(psp, "Found no up/down pair");
     }
     return true;
   }
 }
 
 int main(){
   constexpr std::size_t n_psp_base = 1337;
   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.5;
+  constexpr double subl_chance = 0.5;
   const Subleading subl_channels(~0l);
   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{};
 
   // Wp2j
   Process proc {{pid::proton,pid::proton}, 2, pid::Wp, {}};
   std::size_t n_psp = n_psp_base;
   for( std::size_t i = 0; i<n_psp; ++i){
-    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
+    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_chance,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==Status::good){
       check_W2j(psp, *proc.boson);
     } else { // bad process -> try again
       ++n_psp;
     }
   }
   std::cout << "Wp+2j: Took " << n_psp << " to generate "
     << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
   // Wm2j
   proc = Process{{pid::proton,pid::proton}, 2, pid::Wm, {}};
   n_psp = n_psp_base;
   for( std::size_t i = 0; i<n_psp; ++i){
-    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
+    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_chance,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==Status::good){
       check_W2j(psp, *proc.boson);
     } else { // bad process -> try again
       ++n_psp;
     }
   }
   std::cout << "Wm+2j: Took " << n_psp << " to generate "
     << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
 
   std::cout << "All processes passed." << std::endl;
   return EXIT_SUCCESS;
 }
diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc
index d0557f9..a27225c 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,228 +1,228 @@
 /**
  *  \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;
+  constexpr double subl_chance = 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;
 
   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,
+    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_chance,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==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 << "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,
+    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_chance,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==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();
   subl_channels.reset(subleading::qqx); //!< @TODO remove, temp fix to not double count
   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,
+    const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_chance,subl_channels,
       boson_decays, ew_constants, ran};
     if(psp.status()==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()) );
       }
       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;
 }