diff --git a/t/hej_test.cc b/t/hej_test.cc
index bc89eb2..c392c47 100644
--- a/t/hej_test.cc
+++ b/t/hej_test.cc
@@ -1,639 +1,639 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2022
  *  \copyright GPLv2 or later
  */
 #include "hej_test.hh"
 
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
 #include <cstdlib>
 #include <iterator>
 #include <memory>
 #include <numeric>
 #include <random>
 #include <utility>
 
 #include "HEJ/EWConstants.hh"
 #include "HEJ/Particle.hh"
 #include "HEJ/PDG_codes.hh"
 #include "HEJ/utility.hh"
 
 namespace {
   const HEJ::ParticleProperties Wprop{80.385, 2.085};
   const HEJ::ParticleProperties Zprop{91.187, 2.495};
   const HEJ::ParticleProperties Hprop{125,    0.004165};
   constexpr double vev = 246.2196508;
   const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop};
 }
 
 HEJ::Event::EventData get_process(int const njet, int const pos_boson){
   using namespace HEJ::pid;
   HEJ::Event::EventData ev;
   if(njet == 0){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, {   0,    0,   44,  132}, {}});
       return ev;
     default: // jet idx: -1 -1
       ev.outgoing.push_back({gluon, { -24,   12,  -57,   63}, {}});
       ev.outgoing.push_back({gluon, {  24,  -12,   41,   49}, {}});
       return ev;
     }
   }
   if(njet == 1){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, {  16,  -32,  -99,  163}, {}});
       ev.outgoing.push_back({gluon, { -16,   32,   76,   84}, {}});
       return ev;
     case 1:
       ev.outgoing.push_back({gluon, { -92,   84,  -57,  137}, {}});
       ev.outgoing.push_back({higgs, {  92,  -84,  -79,  193}, {}});
       return ev;
     default: // jet idx: 0 -1 -1
       ev.outgoing.push_back({gluon, {  23,   28,  -44,   57}, {}});
       ev.outgoing.push_back({gluon, { -11,  -24,  -12,   29}, {}});
       ev.outgoing.push_back({gluon, { -12,   -4,   39,   41}, {}});
       return ev;
     }
   }
   if(njet == 2){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, { 198,   33, -170,  291}, {}});
       ev.outgoing.push_back({gluon, {-154,   68,   44,  174}, {}});
       ev.outgoing.push_back({gluon, { -44, -101,   88,  141}, {}});
       return ev;
     case 1:
       ev.outgoing.push_back({gluon, {  -6,   82, -159,  179}, {}});
       ev.outgoing.push_back({higgs, { 195, -106,   74,  265}, {}});
       ev.outgoing.push_back({gluon, {-189,   24,  108,  219}, {}});
       return ev;
     case 2:
       ev.outgoing.push_back({gluon, { -80,  -80, -140,  180}, {}});
       ev.outgoing.push_back({gluon, { -60,  -32,    0,   68}, {}});
       ev.outgoing.push_back({higgs, { 140,  112,  177,  281}, {}});
       return ev;
     default:
       ev.outgoing.push_back({gluon, { -72,   24,   18,   78}, {}});
       ev.outgoing.push_back({gluon, {  72,  -24,   74,  106}, {}});
       return ev;
     }
   }
   if(njet == 3){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, { 152, -117,  -88,  245}, {}});
       ev.outgoing.push_back({gluon, {-146,   62,  -11,  159}, {}});
       ev.outgoing.push_back({gluon, { 126,  -72,   96,  174}, {}});
       ev.outgoing.push_back({gluon, {-132,  127,  144,  233}, {}});
       return ev;
     case 1:
       ev.outgoing.push_back({gluon, {-191,  188, -128,  297}, {}});
       ev.outgoing.push_back({higgs, { 199,   72,  -76,  257}, {}});
       ev.outgoing.push_back({gluon, { 184, -172,   -8,  252}, {}});
       ev.outgoing.push_back({gluon, {-192,  -88,   54,  218}, {}});
       return ev;
     case 2:
       ev.outgoing.push_back({gluon, { -42,   18,  -49,   67}, {}});
       ev.outgoing.push_back({gluon, { -12,  -54,  -28,   62}, {}});
       ev.outgoing.push_back({higgs, {  99,   32,  -16,  163}, {}});
       ev.outgoing.push_back({gluon, { -45,    4,   72,   85}, {}});
       return ev;
     case 3:
       ev.outgoing.push_back({gluon, { -65,  -32,  -76,  105}, {}});
       ev.outgoing.push_back({gluon, { -22,   31,  -34,   51}, {}});
       ev.outgoing.push_back({gluon, { -12,  -67,  -36,   77}, {}});
       ev.outgoing.push_back({higgs, {  99,   68,   -4,  173}, {}});
       return ev;
     default:
       ev.outgoing.push_back({gluon, { -90, -135,   30,  165}, {}});
       ev.outgoing.push_back({gluon, {-108,  198,   76,  238}, {}});
       ev.outgoing.push_back({gluon, { 198,  -63,  126,  243}, {}});
       return ev;
     }
   }
   if(njet == 4){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, { 199,   72,  -76,  257}, {}});
       ev.outgoing.push_back({gluon, {-200, -155,  -64,  261}, {}});
       ev.outgoing.push_back({gluon, { 198,  194,   57,  283}, {}});
       ev.outgoing.push_back({gluon, {   1,   32,    8,   33}, {}});
       ev.outgoing.push_back({gluon, {-198, -143,  186,  307}, {}});
       return ev;
     case 1:
       ev.outgoing.push_back({gluon, { 198,   61, -162,  263}, {}});
       ev.outgoing.push_back({higgs, { 199,   72,  -76,  257}, {}});
       ev.outgoing.push_back({gluon, {-200,  135,  144,  281}, {}});
       ev.outgoing.push_back({gluon, {-198, -186,  171,  321}, {}});
       ev.outgoing.push_back({gluon, {   1,  -82,  122,  147}, {}});
       return ev;
     case 2:
       ev.outgoing.push_back({gluon, {-180,  -27, -164,  245}, {}});
       ev.outgoing.push_back({gluon, {-108,   78,  -36,  138}, {}});
       ev.outgoing.push_back({higgs, { 196, -189,   68,  307}, {}});
       ev.outgoing.push_back({gluon, {-107,  136,   76,  189}, {}});
       ev.outgoing.push_back({gluon, { 199,    2,  178,  267}, {}});
       return ev;
     case 3:
       ev.outgoing.push_back({gluon, { -12,  -30,  -84,   90}, {}});
       ev.outgoing.push_back({gluon, { -72,   22,  -96,  122}, {}});
       ev.outgoing.push_back({gluon, {  68,    0,  -51,   85}, {}});
       ev.outgoing.push_back({higgs, {  64,   72,  -81,  177}, {}});
       ev.outgoing.push_back({gluon, { -48,  -64,   84,  116}, {}});
       return ev;
     case 4:
       ev.outgoing.push_back({gluon, { -72,  -49,  -72,  113}, {}});
       ev.outgoing.push_back({gluon, { -48,    0,  -36,   60}, {}});
       ev.outgoing.push_back({gluon, { -12,   54,  -36,   66}, {}});
       ev.outgoing.push_back({gluon, {  68,  -77,  -56,  117}, {}});
       ev.outgoing.push_back({higgs, {  64,   72,  -81,  177}, {}});
       return ev;
     default:
       ev.outgoing.push_back({gluon, { 248,  -56, -122,  282}, {}});
       ev.outgoing.push_back({gluon, { 249,   30,  -10,  251}, {}});
       ev.outgoing.push_back({gluon, {-249,  -18,   26,  251}, {}});
       ev.outgoing.push_back({gluon, {-248,   44,  199,  321}, {}});
       return ev;
     }
   }
   if(njet == 6){
     switch(pos_boson){
     case 0:
       ev.outgoing.push_back({higgs, { 349,  330,  -94,  505}, {}});
       ev.outgoing.push_back({gluon, {-315, -300,    0,  435}, {}});
       ev.outgoing.push_back({gluon, { 347,  306,   18,  463}, {}});
       ev.outgoing.push_back({gluon, {-249, -342,  162,  453}, {}});
       ev.outgoing.push_back({gluon, { 345,  312,  284,  545}, {}});
       ev.outgoing.push_back({gluon, {-324, -126,  292,  454}, {}});
       ev.outgoing.push_back({gluon, {-153, -180,  304,  385}, {}});
       return ev;
     case 1:
       ev.outgoing.push_back({gluon, { 242,  241, -182,  387}, {}});
       ev.outgoing.push_back({higgs, { 243,  238, -190,  409}, {}});
       ev.outgoing.push_back({gluon, {-218, -215,  -74,  315}, {}});
       ev.outgoing.push_back({gluon, {-224, -224,  112,  336}, {}});
       ev.outgoing.push_back({gluon, { 241,  182,  154,  339}, {}});
       ev.outgoing.push_back({gluon, { -53, -234,  126,  271}, {}});
       ev.outgoing.push_back({gluon, {-231,   12,  156,  279}, {}});
       return ev;
     case 2:
       ev.outgoing.push_back({gluon, { 151,  102,  -42,  187}, {}});
       ev.outgoing.push_back({gluon, { -86,  -46,  -17,   99}, {}});
       ev.outgoing.push_back({higgs, { 152,  153,    0,  249}, {}});
       ev.outgoing.push_back({gluon, { -60, -135,   64,  161}, {}});
       ev.outgoing.push_back({gluon, { 150,  123,  110,  223}, {}});
       ev.outgoing.push_back({gluon, {-154,  -49,   98,  189}, {}});
       ev.outgoing.push_back({gluon, {-153, -148,  144,  257}, {}});
       return ev;
     case 3:
       ev.outgoing.push_back({gluon, { 198,  197,  -66,  287}, {}});
       ev.outgoing.push_back({gluon, {-198, -189,  -54,  279}, {}});
       ev.outgoing.push_back({gluon, {-200,  -64,    2,  210}, {}});
       ev.outgoing.push_back({higgs, { 199,  158,    6,  283}, {}});
       ev.outgoing.push_back({gluon, {-199, -184,  172,  321}, {}});
       ev.outgoing.push_back({gluon, { 196,  168,  177,  313}, {}});
       ev.outgoing.push_back({gluon, {   4,  -86,   92,  126}, {}});
       return ev;
     case 4:
       ev.outgoing.push_back({gluon, { 151,  102,  -42,  187}, {}});
       ev.outgoing.push_back({gluon, { -86, -133,  -14,  159}, {}});
       ev.outgoing.push_back({gluon, {-154, -104,   -8,  186}, {}});
       ev.outgoing.push_back({gluon, { -60,   11,    0,   61}, {}});
       ev.outgoing.push_back({higgs, { 152,  153,    0,  249}, {}});
       ev.outgoing.push_back({gluon, { 150,  125,   90,  215}, {}});
       ev.outgoing.push_back({gluon, {-153, -154,  126,  251}, {}});
       return ev;
     case 5:
       ev.outgoing.push_back({gluon, { -15,  -90,  -94,  131}, {}});
       ev.outgoing.push_back({gluon, { -11,   82,  -74,  111}, {}});
       ev.outgoing.push_back({gluon, {  23,  -80,  -64,  105}, {}});
       ev.outgoing.push_back({gluon, { -48,  -25,  -36,   65}, {}});
       ev.outgoing.push_back({gluon, { -12,   99,  -16,  101}, {}});
       ev.outgoing.push_back({higgs, {  68,   92,  -18,  170}, {}});
       ev.outgoing.push_back({gluon, {  -5,  -78,   54,   95}, {}});
       return ev;
     case 6:
       ev.outgoing.push_back({gluon, { 198,  197,  -66,  287}, {}});
       ev.outgoing.push_back({gluon, {   4,  -84,  -18,   86}, {}});
       ev.outgoing.push_back({gluon, {-198,  -60,  -36,  210}, {}});
       ev.outgoing.push_back({gluon, { 196,  -78,  -36,  214}, {}});
       ev.outgoing.push_back({gluon, {-200,   45,    0,  205}, {}});
       ev.outgoing.push_back({gluon, {-199, -178,    2,  267}, {}});
       ev.outgoing.push_back({higgs, { 199,  158,    6,  283}, {}});
       return ev;
     default:
       ev.outgoing.push_back({gluon, {-350, -112, -280,  462}, {}});
       ev.outgoing.push_back({gluon, { 347,  266, -322,  543}, {}});
       ev.outgoing.push_back({gluon, {-349, -314,  -38,  471}, {}});
       ev.outgoing.push_back({gluon, { 349,  348,   12,  493}, {}});
       ev.outgoing.push_back({gluon, {-342,  -54,   23,  347}, {}});
       ev.outgoing.push_back({gluon, { 345, -134,  138,  395}, {}});
       return ev;
     }
   }
   if(njet == 7){
     switch(pos_boson){
       case -1: // jet idx: -1 0 1 2 3 4 5
         ev.outgoing.push_back({gluon, { -15,  -18,  -54,   59}, {}});
         ev.outgoing.push_back({gluon, { -11,   98,  -70,  121}, {}});
         ev.outgoing.push_back({gluon, {  23, -100,  -64,  121}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,  -20,  117}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -92,  -12,   93}, {}});
         ev.outgoing.push_back({gluon, { -48,  -76,   -2,   90}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         return ev;
       case -2: // jet idx: 0 1 2 3 4 -1 -1
         ev.outgoing.push_back({gluon, {  -5,  -86,  -82,  119}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -48,  -14,   20,   54}, {}});
         ev.outgoing.push_back({gluon, {  23,  -50,   26,   61}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         ev.outgoing.push_back({gluon, { -15,  -18,   54,   59}, {}});
         ev.outgoing.push_back({gluon, { -11,  -20,   88,   91}, {}});
         return ev;
       case -3: // jet idx: 0 0 1 2 2 3 4
         // jet pt fraction: 0.6 0.38 1 0.49 0.51 1 1
         ev.outgoing.push_back({gluon, {  23,  -94,  -62, 1.2e+02}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -62,  -34,   71}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20, 1.2e+02}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56, 1.1e+02}, {}});
         ev.outgoing.push_back({gluon, { -11,   98,   70, 1.2e+02}, {}});
         ev.outgoing.push_back({gluon, { -48, -1e+02,   82, 1.4e+02}, {}});
         ev.outgoing.push_back({gluon, { -15,  -30,   78,   85}, {}});
         return ev;
       case -4: // jet idx: 0 1 1 2 3 4 4
         // jet pt fraction: 1 0.51 0.49 1 1 0.25 0.75
         ev.outgoing.push_back({gluon, {  -5,  -88,  -64,  109}, {}});
         ev.outgoing.push_back({gluon, { -11,   98,  -70,  121}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, {  23,  -70,   22,   77}, {}});
         ev.outgoing.push_back({gluon, { -15,  -32,   16,   39}, {}});
         ev.outgoing.push_back({gluon, { -48,  -96,   75,  131}, {}});
         return ev;
       case -5: // jet idx: 0 1 -1 -1 2 3 4
         ev.outgoing.push_back({gluon, { -15,  -26,  -62,   69}, {}});
         ev.outgoing.push_back({gluon, { -48,  -60,  -54,   94}, {}});
         ev.outgoing.push_back({gluon, {  23,   10,  -14,   29}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -20,    0,   21}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -11,  -92,   40,  101}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         return ev;
       case -6: // jet idx: 0 1 1 2 -1 2 3
         // jet pt fraction: 1 0.63 0.36 0.49 1 0.51 1
         ev.outgoing.push_back({gluon, {  68,   93,  -20,  117}, {}});
         ev.outgoing.push_back({gluon, { -48, -100,   26,  114}, {}});
         ev.outgoing.push_back({gluon, { -15,  -62,   26,   69}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -28,   20,   35}, {}});
         ev.outgoing.push_back({gluon, { -11,   98,   70,  121}, {}});
         ev.outgoing.push_back({gluon, {  23,  -96,   92,  135}, {}});
         return ev;
       case -7: // jet idx: 0 1 2 2 3 3 4
         // jet pt fraction: 1 1 0.51 0.49 0.18 0.82 1
         ev.outgoing.push_back({gluon, { -15,  -80, -100,  129}, {}});
         ev.outgoing.push_back({gluon, {  23,  -96,  -92,  135}, {}});
         ev.outgoing.push_back({gluon, { -11,   98,  -70,  121}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -22,  -10,   25}, {}});
         ev.outgoing.push_back({gluon, { -48,  -88,  -31,  105}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         return ev;
       case -8: // jet idx: 0 1 2 2 2 3 4
         // jet pt fraction: 1 1 0.21 0.37 0.41 1 1
         ev.outgoing.push_back({gluon, { -48,  -44,  -62,   90}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -50,  -22,   55}, {}});
         ev.outgoing.push_back({gluon, {  23,  -90,  -34,   99}, {}});
         ev.outgoing.push_back({gluon, { -15, -100,  -28,  105}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,  -20,  117}, {}});
         ev.outgoing.push_back({gluon, { -11,   96,   76,  123}, {}});
         return ev;
       case -9: // jet idx: 0 1 2 1 3 0 4
         // jet pt fraction: 0.72 0.51 1 0.49 1 0.28 1
         ev.outgoing.push_back({gluon, { -15,  -98,  -62,  117}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  23,  -76,  -40,   89}, {}});
         ev.outgoing.push_back({gluon, { -11,   92,  -40,  101}, {}});
         ev.outgoing.push_back({gluon, { -48,  -68,  -34,   90}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -38,  -14,   41}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         return ev;
       case -10: // jet idx: 0 1 3 2 4 3 1
         // jet pt fraction: 1 0.33 0.51 1 1 0.49 0.67
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -48,   16,   51}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         ev.outgoing.push_back({gluon, {  23,  -76,   52,   95}, {}});
         ev.outgoing.push_back({gluon, { -48,  -60,   54,   94}, {}});
         ev.outgoing.push_back({gluon, { -11,   92,   68,  115}, {}});
         ev.outgoing.push_back({gluon, { -15,  -96,   72,  121}, {}});
         return ev;
       case -11: // jet idx: 0 1 2 3 4 -1 5
         // jet pt fraction: 1 1 1 1 1 1 1
         ev.outgoing.push_back({gluon, { -11,   98,  -70,  121}, {}});
         ev.outgoing.push_back({gluon, { -15,  -98,  -62,  117}, {}});
         ev.outgoing.push_back({gluon, {  23,  -90,   -2,   93}, {}});
         ev.outgoing.push_back({gluon, { -48,  -76,    2,   90}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -22,   10,   25}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         return ev;
       case -12: // jet idx: 0 1 -1 2 3 4 3
         // jet pt fraction: 1 1 1 1 0.35 1 0.65
         ev.outgoing.push_back({gluon, {  23,  -94,  -62,  115}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -28,    4,   29}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -15,  -58,   34,   69}, {}});
         ev.outgoing.push_back({gluon, { -11,   92,   68,  115}, {}});
         ev.outgoing.push_back({gluon, { -48, -100,   82,  138}, {}});
         return ev;
       case -13: // jet idx: 0 1 2 3 3 4 2
         // jet pt fraction: 1 1 0.5 0.35 0.65 1 0.5
         ev.outgoing.push_back({gluon, { -15,  -98,  -62,  117}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -28,   20,   35}, {}});
-        ev.outgoing.push_back({gluon, { -48,  -96,   75,  131}, {}});
+        ev.outgoing.push_back({gluon, { -48,  -96,   77,  131}, {}});
         ev.outgoing.push_back({gluon, {  23,  -62,   50,   83}, {}});
         ev.outgoing.push_back({gluon, { -11,   96,   76,  123}, {}});
         return ev;
       case -14: // jet idx: 0 1 2 3 3 4 2
         // jet pt fraction: 1 1 0.52 0.35 0.65 1 0.48
         ev.outgoing.push_back({gluon, {  23,  -94,  -62,  115}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,   56,  111}, {}});
         ev.outgoing.push_back({gluon, { -15,  -96,   72,  121}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -42,   38,   57}, {}});
         ev.outgoing.push_back({gluon, { -48,  -44,   62,   90}, {}});
         ev.outgoing.push_back({gluon, { -11,   88,   88,  125}, {}});
         return ev;
       case -15: // jet idx: 0 -1 1 0 2 3 4
         // jet pt fraction: 0.51 1 1 0.49 1 1 1
-        ev.outgoing.push_back({gluon, { -11,   98,  -70,  121}, {}});
+        ev.outgoing.push_back({gluon, { -11,   98,  -71,  121}, {}});
         ev.outgoing.push_back({gluon, {  -5,  -16,  -12,   21}, {}});
         ev.outgoing.push_back({gluon, {  23,  -94,  -62,  115}, {}});
         ev.outgoing.push_back({gluon, { -12,   95,  -56,  111}, {}});
         ev.outgoing.push_back({gluon, {  68,   93,   20,  117}, {}});
         ev.outgoing.push_back({gluon, { -48,  -76,   70,  114}, {}});
         ev.outgoing.push_back({gluon, { -15, -100,   80,  129}, {}});
         return ev;
     }
   }
   throw HEJ::unknown_option{"unknown process"};
 }
 
 double mass(HEJ::Particle const & p) {
   static constexpr double M_PROTON = 0.938;
   using namespace HEJ::pid;
   switch(p.type){
   case Higgs:
     return Hprop.mass;
   case Z:
   case Z_photon_mix:
     return Zprop.mass;
   case Wm:
   case Wp:
     return Wprop.mass;
   case proton:
     return M_PROTON;
   default:
     if(is_massless(p)) return 0.;
     throw std::invalid_argument{
       name(p.type) + " has unknown mass"
     };
   }
 }
 
 void repair_momenta(HEJ::Event::EventData & ev) {
   double E_sum = 0., pz_sum = 0.;
   for(auto & out: ev.outgoing) {
     const double m = mass(out);
     const double psq = out.px()*out.px() + out.py()*out.py() + out.pz()*out.pz();
     const double E = sqrt(psq + m * m);
     out.p.reset(out.px(), out.py(), out.pz(), E);
     assert(HEJ::nearby_ep(out.p.m(), m, 1e-7*E));
     E_sum += E;
     pz_sum += out.pz();
   }
   const double p0 = (E_sum - pz_sum)/2.;
   const double p1 = (E_sum + pz_sum)/2.;
   ev.incoming[0].p.reset(0, 0, -p0, p0);
   ev.incoming[1].p.reset(0, 0, +p1, p1);
 }
 
 HEJ::Event::EventData parse_configuration(
     std::array<std::string,2> const & in, std::vector<std::string> const & out,
     int const overwrite_boson
 ){
   auto boson = std::find_if(out.cbegin(), out.cend(),
     [](std::string const & id){
       return !HEJ::is_parton(HEJ::to_ParticleID(id)); });
   int const pos_boson = (overwrite_boson!=0)?overwrite_boson:
     ((boson == out.cend())?-1:std::distance(out.cbegin(), boson));
 
   std::size_t njets = out.size();
   if( (overwrite_boson == 0) && boson != out.cend()) --njets;
 
   HEJ::Event::EventData ev{get_process(njets, pos_boson)};
   ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs));
   for(std::size_t i=0; i<out.size(); ++i){
     ev.outgoing[i].type = HEJ::to_ParticleID(out[i]);
   }
   for(std::size_t i=0; i<in.size(); ++i){
     ev.incoming[i].type = HEJ::to_ParticleID(in[i]);
   }
   repair_momenta(ev);
   for(std::size_t i = 0; i < out.size(); ++i) {
     // decay W
     if( std::abs(ev.outgoing[i].type) == HEJ::ParticleID::Wp )
       ev.decays[i]=decay_W(ev.outgoing[i]);
     // decay Z
     if( ev.outgoing[i].type == HEJ::ParticleID::Z_photon_mix )
       ev.decays[i]=decay_Z(ev.outgoing[i]);
   }
   shuffle_particles(ev);
 
   return ev;
 }
 
 HEJ::Event::EventData rapidity_order_ps(
     std::array<std::string,2> const & in,
     std::vector<std::string> const & out,
     bool reconstruct /* = false */,
     std::unordered_map< size_t, std::vector<std::string> > decays /* = {} */,
     bool shuffle
   ) {
   using namespace HEJ::pid;
   HEJ::Event::EventData evd;
 
   const size_t n_out {out.size()};
 
   // Arbitrary non-zero mT2 shift in lieu of mass property of Z_photon_mix
   constexpr double arbitrary_mT2 = 80.*80.;
   // Arbitrary Parameters
   constexpr double pT {40.};
   constexpr double start_y {-4.5};
   constexpr double min_dy {0.5};  // chosen to be above typical jet-radius
   const double dy {std::max(-2.*start_y/(n_out-1), min_dy)};
 
   // Kinematic parameters
   std::array<double,2> cs_phi {1., 0.};
   std::array<double,2> cs_dphi {cos(2.*M_PI/n_out), sin(2.*M_PI/n_out)};
   double y {start_y};
 
   double sum_pE {0.};
   double sum_pz {0.};
   // Outgoing
   evd.outgoing.reserve(n_out);
   for(size_t i = 0; i < n_out; ++i) {
     auto pid = HEJ::to_ParticleID(out[i]);
     // Kinematics
     double mT2 = pT*pT;
     if(HEJ::is_AWZH_boson(pid)) {
       // Z_photon_mix has no properties
       if(pid == Z_photon_mix) { mT2 += arbitrary_mT2; }
       // regular bosons
       else {
         auto const ref = ew_parameters.prop(pid);
         mT2 += ref.mass*ref.mass;
       }
     }
     double pz {sqrt(mT2)*sinh(y)};
     double pE {sqrt(mT2)*cosh(y)};
     evd.outgoing.push_back({pid, {pT*cs_phi[0], pT*cs_phi[1], pz, pE}, {}});
     sum_pE+=pE;
     sum_pz+=pz;
     // Update cos,sin, y
     cs_phi = {cs_phi[0] * cs_dphi[0] - cs_phi[1] * cs_dphi[1],
               cs_phi[1] * cs_dphi[0] + cs_phi[0] * cs_dphi[1]};
     y += dy;
   }
 
   // Specified decays
   evd.decays.reserve(decays.size());
   for(auto const & decay : decays) {
     auto const parent = evd.outgoing.at(decay.first);
     std::vector<HEJ::Particle> progeny = decay_kinematics(parent);
     for(std::size_t i = 0; i < decay.second.size(); ++i) {
       progeny[i].type = HEJ::to_ParticleID(decay.second[i]);
     }
     evd.decays.emplace(decay.first, std::move(progeny));
   }
 
   // Reconstruct decays
   if(reconstruct) { evd.reconstruct_intermediate(ew_parameters); }
 
   // Incoming
   double p0 {0.5*(sum_pE - sum_pz)};
   evd.incoming[0] = {HEJ::to_ParticleID(in[0]),{0,0,-p0,p0},{}};
   double p1 {0.5*(sum_pE + sum_pz)};
   evd.incoming[1] = {HEJ::to_ParticleID(in[1]),{0,0,+p1,p1},{}};
 
   // Shuffle
   if(shuffle) { shuffle_particles(evd); }
   return evd;
 }
 
 namespace {
   std::mt19937_64 RAN{0};
 }
 
 void shuffle_particles(HEJ::Event::EventData & ev) {
   // incoming
   std::shuffle(ev.incoming.begin(), ev.incoming.end(), RAN);
   // outgoing (through index)
   auto old_outgoing = std::move(ev).outgoing;
   std::vector<std::size_t> idx(old_outgoing.size());
   std::iota(idx.begin(), idx.end(), 0);
   std::shuffle(idx.begin(), idx.end(), RAN);
   ev.outgoing.clear();
   ev.outgoing.reserve(old_outgoing.size());
   for(std::size_t i: idx) {
     ev.outgoing.emplace_back(std::move(old_outgoing[i]));
   }
   // find decays again
   if(!ev.decays.empty()){
     auto old_decays = std::move(ev).decays;
     ev.decays.clear();
     for(std::size_t i=0; i<idx.size(); ++i) {
       auto decay = old_decays.find(idx[i]);
       if(decay != old_decays.end())
         ev.decays.emplace(i, std::move(decay->second));
     }
     for(auto & decay: ev.decays){
       std::shuffle(decay.second.begin(), decay.second.end(), RAN);
     }
   }
 }
 
 bool couple_quark(std::string const & boson, std::string & quark){
   if(std::abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){
     auto qflav{ HEJ::to_ParticleID(quark) };
     if(!HEJ::is_anyquark(qflav)) return false;
     const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1;
     if(W_charge*qflav < 0 && !(std::abs(qflav)%2)) return false; // not anti-down
     if(W_charge*qflav > 0 &&  (std::abs(qflav)%2)) return false; // not up
     quark=std::to_string(qflav-W_charge);
   }
   if(HEJ::to_ParticleID(boson) == HEJ::ParticleID::Z_photon_mix){
     auto qflav{ HEJ::to_ParticleID(quark) };
     if(!HEJ::is_anyquark(qflav)) return false;
   }
   return true;
 }
 
 std::vector<HEJ::Particle> decay_kinematics( HEJ::Particle const & parent ) {
   std::vector<HEJ::Particle> decay_products(2);
   const double E = parent.m()/2;
   const double theta = 2.*M_PI*RAN()/static_cast<double>(RAN.max());
   const double cos_phi = 2.*RAN()/static_cast<double>(RAN.max())-1.;
   const double sin_phi = std::sqrt(1. - cos_phi*cos_phi);  // Know 0 < phi < pi
   const double px = E*std::cos(theta)*sin_phi;
   const double py = E*std::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;
 }
 
 std::vector<HEJ::Particle> decay_W( HEJ::Particle const & parent ){
   if(parent.m() == 0.) // we can't decay massless partons
     return {};
   std::array<HEJ::ParticleID, 2> decays{};
   if(parent.type==HEJ::ParticleID::Wp){
     // order matters: first particle, second anti
     decays[0] = HEJ::ParticleID::nu_e;
     decays[1] = HEJ::ParticleID::e_bar;
   } else {
     // this function is for testing: we don't check that parent==W boson
     decays[0] = HEJ::ParticleID::e;
     decays[1] = HEJ::ParticleID::nu_e_bar;
   }
   std::vector<HEJ::Particle> decay_products = decay_kinematics(parent);
   for(std::size_t i = 0; i < decay_products.size(); ++i){
     decay_products[i].type = decays[i];
   }
   return decay_products;
 }
 
 std::vector<HEJ::Particle> decay_Z( HEJ::Particle const & parent ){
   if(parent.m() == 0.) // we can't decay massless partons
     return {};
   std::array<HEJ::ParticleID, 2> decays{};
   // order matters: first particle, second anti
   decays[0] = HEJ::ParticleID::electron;
   decays[1] = HEJ::ParticleID::positron;
   std::vector<HEJ::Particle> decay_products = decay_kinematics(parent);
   for(std::size_t i = 0; i < decay_products.size(); ++i){
     decay_products[i].type = decays[i];
   }
   return decay_products;
 }