diff --git a/src/Wjets.cc b/src/Wjets.cc
index 609290c..c471e6a 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,648 +1,593 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/Wjets.hh"
 
 #include <algorithm>
 #include <array>
 #include <cassert>
 #include <cmath>
 #include <iostream>
 
 #include "HEJ/Constants.hh"
 #include "HEJ/EWConstants.hh"
 #include "HEJ/LorentzVector.hh"
 #include "HEJ/exceptions.hh"
 #include "HEJ/jets.hh"
 
 // generated headers
 #include "HEJ/currents/jV_j.hh"
 #include "HEJ/currents/jV_juno.hh"
 #include "HEJ/currents/jVuno_j.hh"
 #include "HEJ/currents/jW_jqqbar.hh"
 #include "HEJ/currents/jW_qqbar_j.hh"
 #include "HEJ/currents/jWqqbar_j.hh"
 #include "HEJ/currents/j_Wqqbar_j.hh"
 
 namespace HEJ {
 namespace currents {
 
   namespace {
     using COM = std::complex<double>;
 
     // --- Helper Functions ---
 
     // FKL W Helper Functions
     double WProp(const HLV & plbar, const HLV & pl,
                  ParticleProperties const & wprop
     ){
       COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
                               + COM(0.,1.)*wprop.mass*wprop.width);
       double PropFactor=(propW*conj(propW)).real();
       return PropFactor;
     }
 
     /**
      * @brief  Contraction of W + unordered jet current with FKL current
      *
      * See eq:wunocurrent in the developer manual for the definition
      * of the W + unordered jet current
      */
     template<Helicity h1, Helicity h2, Helicity pol>
     double jM2Wuno(
         HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
         HLV const & pa, HLV const & p2, HLV const & pb,
         ParticleProperties const & wprop
     ){
       using helicity::minus;
       const COM u1 = U1<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
       const COM u2 = U2<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
       const COM l = L<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
 
       const COM x = u1 - l;
       const COM y = u2 + l;
 
       double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
       amp *= WProp(plbar, pl, wprop);
 
       const HLV q1g = pa-pl-plbar-p1-pg;
       const HLV q2 = p2-pb;
 
       const double t1 = q1g.m2();
       const double t2 = q2.m2();
 
       amp /= (t1*t2);
       return amp;
     }
 
     /**
      * @brief Contraction of W + extremal qqbar current with FKL current
      *
      * See eq:crossed in the developer manual for the definition
      * of the W + extremal qqbar current.
      *
      */
     template<Helicity h1, Helicity h2, Helicity hg>
     double jM2_W_extremal_qqbar(
       HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
       HLV const & pqbar, HLV const & p3, HLV const & pb,
       ParticleProperties const & wprop
     ){
 
       const COM u1Xcontr = U1X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
       const COM u2Xcontr = U2X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
       const COM lXcontr = LX<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
 
       const COM x = u1Xcontr - lXcontr;
       const COM y = u2Xcontr + lXcontr;
 
       double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
       amp *= WProp(plbar, pl, wprop);
 
       const HLV q1 = pg-pl-plbar-pq-pqbar;
       const HLV q2 = p3-pb;
 
       const double t1 = q1.m2();
       const double t2 = q2.m2();
 
       amp /= (t1*t2);
       return amp;
     }
   } // Anonymous Namespace
 
   //! W+Jets FKL Contributions
   double ME_W_qQ(
     HLV const & p1out,
     HLV const & plbar, HLV const & pl,
     HLV const & p1in,
     HLV const & p2out, HLV const & p2in,
     ParticleProperties const & wprop
   ){
     using helicity::minus;
     using helicity::plus;
 
     const COM ampm = jV_j<minus, minus, minus>(p1in,p1out,p2in,p2out,pl,plbar);
     const COM ampp = jV_j<minus, minus, plus>(p1in,p1out,p2in,p2out,pl,plbar);
 
     const double Msqr = std::norm(ampm) + std::norm(ampp);
 
     return WProp(plbar, pl, wprop) * Msqr;
   }
 
   namespace {
     // helicity amplitude squares for contraction of W current with unordered
     // current
     template<Helicity h2, Helicity hg>
     double ampsq_jW_juno(
       HLV const & pa, HLV const & p1,
       HLV const & pb, HLV const & p2,
       HLV const & pg,
       HLV const & pl, HLV const & plbar
     ){
       using helicity::minus;
 
       const COM u1 = U1_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
       const COM u2 = U2_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
       const COM l  = L_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
 
       return 2.*C_F*std::real((l-u1)*std::conj(l+u2))
         + 2.*C_F*C_F/3.*std::norm(u1+u2)
         ;
     }
-
-    /**
-     * @brief W+Jets Unordered Contributions, function to handle all incoming types.
-     * @param p1out             Outgoing Particle 1. (W emission)
-     * @param plbar             Outgoing election momenta
-     * @param pl                Outgoing neutrino momenta
-     * @param p1in              Incoming Particle 1. (W emission)
-     * @param p2out             Outgoing Particle 2 (Quark, unordered emission this side.)
-     * @param p2in              Incoming Particle 2 (Quark, unordered emission this side.)
-     * @param pg                Unordered Gluon momenta
-     * @param aqlineb           Bool. Is Backwards quark line an anti-quark line?
-     * @param aqlinef           Bool. Is Forwards quark line an anti-quark line?
-     *
-     * Calculates j_W ^\mu    j_{uno}_\mu. Ie, unordered with W emission
-     * opposite side. Handles all possible incoming states.
-     *
-     * @TODO use const & for pl plbar
-     */
-    double jW_juno(
-      HLV const & p1out, HLV plbar, HLV pl, HLV const & p1in,
-      HLV const & p2out, HLV const & p2in, HLV const & pg,
-      bool /* aqlineb */, bool /* aqlinef */,
-      ParticleProperties const & wprop
-    ){
-      using helicity::minus;
-      using helicity::plus;
-
-      // helicity assignments assume quarks
-      // in the end, this is irrelevant since we sum over all helicities
-      const double ampsq =
-        + ampsq_jW_juno<minus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
-        + ampsq_jW_juno<minus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
-        + ampsq_jW_juno<plus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
-        + ampsq_jW_juno<plus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
-        ;
-
-      const HLV q1=p1in-p1out-plbar-pl;
-      const HLV q2=-(p2in-p2out-pg);
-
-      //! @TODO explain magic 16
-      return WProp(plbar, pl, wprop)*ampsq/(16.*q2.m2()*q1.m2());
-    }
   } // Anonymous Namespace
 
   double ME_W_unob_qQ(
       HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
       HLV const & pg, HLV const & plbar, HLV const & pl, ParticleProperties
       const & wprop
   ){
-    return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
+    using helicity::minus;
+    using helicity::plus;
+
+    // helicity assignments assume quarks
+    // in the end, this is irrelevant since we sum over all helicities
+    const double ampsq =
+      + ampsq_jW_juno<minus, minus>(p2in,p2out,p1in,p1out,pg,pl,plbar)
+      + ampsq_jW_juno<minus, plus> (p2in,p2out,p1in,p1out,pg,pl,plbar)
+      + ampsq_jW_juno<plus, minus> (p2in,p2out,p1in,p1out,pg,pl,plbar)
+      + ampsq_jW_juno<plus, plus>  (p2in,p2out,p1in,p1out,pg,pl,plbar)
+      ;
+
+    const double t1 = (p1in - p1out - pg).m2();
+    const double tn = (p2in - p2out - plbar - pl).m2();
+
+    //! @TODO explain magic 16
+    return WProp(plbar, pl, wprop)*ampsq/(16.*t1*tn);
   }
 
   namespace {
 
     // helicity sum helper for jWuno_j(...)
     template<Helicity h1>
     double jWuno_j_helsum(
         HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
         HLV const & pa, HLV const & p2, HLV const & pb,
         ParticleProperties const & wprop
     ){
       using helicity::minus;
       using helicity::plus;
 
       const double ME2h1pp = jM2Wuno<h1, plus, plus>(
         pg, p1, plbar, pl, pa, p2, pb, wprop
       );
 
       const double ME2h1pm = jM2Wuno<h1, plus, minus>(
         pg, p1, plbar, pl, pa, p2, pb, wprop
       );
 
       const double ME2h1mp = jM2Wuno<h1, minus, plus>(
         pg, p1, plbar, pl, pa, p2, pb, wprop
       );
 
       const double ME2h1mm = jM2Wuno<h1, minus, minus>(
         pg, p1, plbar, pl, pa, p2, pb, wprop
       );
 
       return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
     }
 
-  /**
-   * @brief W+Jets Unordered Contributions, function to handle all incoming
-   * @types.
-   * @param pg                Unordered Gluon momenta
-   * @param p1out             Outgoing Particle 1. (Quark - W and Uno emission)
-   * @param plbar             Outgoing election momenta
-   * @param pl                Outgoing neutrino momenta
-   * @param p1in              Incoming Particle 1. (Quark - W and Uno emission)
-   * @param p2out             Outgoing Particle 2
-   * @param p2in              Incoming Particle 2
-   * @param aqlineb           Bool. Is Backwards quark line an anti-quark line?
-   *
-   * Calculates j_W_{uno} ^\mu    j_\mu. Ie, unordered with W emission same
-   * side. Handles all possible incoming states. Note this handles both forward
-   * and back- -ward Wuno emission. For forward, ensure p1out is the uno and W
-   * emission parton.
-   * @TODO: Include separate wrapper functions for forward and backward to clean
-up *        ME_W_unof_current in `MatrixElement.cc`.
-   */
-    double jWuno_j(
-        HLV const & pg, HLV const & p1out, HLV const & plbar, HLV const & pl,
-        HLV const & p1in, HLV const & p2out, HLV const & p2in, bool aqlineb,
-        ParticleProperties const & wprop
-      ){
-      const auto helsum = aqlineb?
-        jWuno_j_helsum<helicity::plus>:
-        jWuno_j_helsum<helicity::minus>;
-
-      //Helicity sum and average over initial states
-      return helsum(pg, p1out,plbar,pl,p1in,p2out,p2in, wprop)/(4.*C_A*C_A);
-    }
   } // Anonymous Namespace
 
   double ME_Wuno_qQ(
       HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
       HLV const & pg, HLV const & plbar, HLV const & pl,
       ParticleProperties const & wprop
   ){
-    return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
+    return jWuno_j_helsum<helicity::minus>(
+      pg, p1out, plbar, pl, p1in, p2out, p2in, wprop
+    )/(4.*C_A*C_A);
   }
 
   // helicity sum helper for jWqqbar_j(...)
   template<Helicity h1>
   double jWqqbar_j_helsum(
     HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
     HLV const & pqbar, HLV const & p3, HLV const & pb,
     ParticleProperties const & wprop
   ){
     using helicity::minus;
     using helicity::plus;
 
     const double amp_h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
       pg, pq, plbar, pl, pqbar, p3, pb, wprop
     );
     const double amp_h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
       pg, pq, plbar, pl, pqbar, p3, pb, wprop
     );
     const double amp_h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
       pg, pq, plbar, pl, pqbar, p3, pb, wprop
     );
     const double amp_h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
       pg, pq, plbar, pl, pqbar, p3, pb, wprop
     );
 
     return amp_h1pp + amp_h1pm + amp_h1mp + amp_h1mm;
   }
 
   /**
    * @brief W+Jets Extremal qqbar Contributions, function to handle all incoming
    * @types.
    * @param pgin              Incoming gluon which will split into qqbar.
    * @param pqout             Quark of extremal qqbar outgoing (W-Emission).
    * @param plbar             Outgoing anti-lepton momenta
    * @param pl                Outgoing lepton momenta
    * @param pqbarout          Anti-quark of extremal qqbar pair. (W-Emission)
    * @param pout              Outgoing Particle 2 (end of FKL chain)
    * @param p2in              Incoming Particle 2
    * @param aqlinef           Bool. Is Forwards quark line an anti-quark line?
    *
    * Calculates j_W_{qqbar} ^\mu    j_\mu. Ie, Ex-QQbar with W emission same side.
    * Handles all possible incoming states. Calculated via crossing symmetry from
    * jWuno_j.
    */
   double jWqqbar_j(
       HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
       HLV const & pqbarout, HLV const & p2out, HLV const & p2in, bool aqlinef,
       ParticleProperties const & wprop
   ){
     const auto helsum = aqlinef?
       jWqqbar_j_helsum<helicity::plus>:
       jWqqbar_j_helsum<helicity::minus>;
 
     //Helicity sum and average over initial states.
     double ME2 = helsum(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, wprop)/
       (4.*C_A*C_A);
 
     //Correct colour averaging after crossing:
     ME2*=(3.0/8.0);
 
     return ME2;
   }
 
   double ME_WExqqbar_qbarqQ(
       HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
       HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
       ParticleProperties const & wprop
   ){
     return jWqqbar_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
   }
 
   double ME_WExqqbar_qqbarQ(
       HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
       HLV const & pqout, HLV const & p2out, HLV const & p2in,
       ParticleProperties const & wprop
   ){
     return jWqqbar_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
   }
 
   double ME_WExqqbar_qbarqg(
       HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
       HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
       ParticleProperties const & wprop
   ){
     return jWqqbar_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
             *K_g(p2out,p2in)/C_F;
   }
 
   double ME_WExqqbar_qqbarg(
       HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
       HLV const & pqout, HLV const & p2out, HLV const & p2in,
       ParticleProperties const & wprop
   ){
     return jWqqbar_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
             *K_g(p2out,p2in)/C_F;
   }
 
   namespace {
     template<Helicity h1, Helicity hg>
     double jW_jqqbar(
         HLV const & pb,
         HLV const & p2,
         HLV const & p3,
         HLV const & pa,
         HLV const & p1,
         HLV const & pl,
         HLV const & plbar
     ){
       using std::norm;
 
       static constexpr double cm1m1 = 8./3.;
       static constexpr double cm2m2 = 8./3.;
       static constexpr double cm3m3 = 6.;
       static constexpr double cm1m2 =-1./3.;
       static constexpr COM cm1m3 = COM{0.,-3.};
       static constexpr COM cm2m3 = COM{0.,3.};
 
       const COM m1 = jW_qggm1<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
       const COM m2 = jW_qggm2<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
       const COM m3 = jW_qggm3<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
 
       return
         + cm1m1*norm(m1)
         + cm2m2*norm(m2)
         + cm3m3*norm(m3)
         + 2.*real(cm1m2*m1*conj(m2))
         + 2.*real(cm1m3*m1*conj(m3))
         + 2.*real(cm2m3*m2*conj(m3)) ;
     }
   } // Anonymous Namespace
 
   // contraction of W current with extremal qqbar on the other side
   double ME_W_Exqqbar_QQq(
       HLV const & pa, HLV const & pb,
       HLV const & p1,  HLV const & p2,
       HLV const & p3, HLV const & plbar, HLV const & pl, bool aqlinepa,
       ParticleProperties const & wprop
   ){
     using helicity::minus;
     using helicity::plus;
 
     const double wPropfact = WProp(plbar, pl, wprop);
     const double prefactor = 2.*wPropfact
                 /24./4./((pa-p1-pl-plbar).m2()*(p2+p3-pb).m2());
     if(aqlinepa) {
       return prefactor*(
         jW_jqqbar<plus, minus>(pb,p2,p3,pa,p1,pl,plbar)
         + jW_jqqbar<plus, plus>(pb,p2,p3,pa,p1,pl,plbar)
       );
     }
     return prefactor*(
       jW_jqqbar<minus, minus>(pb,p2,p3,pa,p1,pl,plbar)
       + jW_jqqbar<minus, plus>(pb,p2,p3,pa,p1,pl,plbar)
     );
   }
 
   namespace {
     // helper function for matrix element for W + Jets with central qqbar
     template<Helicity h1a, Helicity h4b>
     double amp_WCenqqbar_qq(
         HLV const & pa, HLV const &  p1,
         HLV const & pb, HLV const &  p4,
         HLV const &  pq, HLV const &  pqbar,
         HLV const &  pl, HLV const &  plbar,
         HLV const &  q11, HLV const &  q12
     ){
       using std::norm;
 
       const COM sym = M_sym_W<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
       );
       const COM cross = M_cross_W<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
       );
       const COM uncross = M_uncross_W<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
       );
 
       // Colour factors
       static constexpr double cmsms = 3.;
       static constexpr double cmumu = 4./3.;
       static constexpr double cmcmc = 4./3.;
       static constexpr COM cmsmu = COM{0., 3./2.};
       static constexpr COM cmsmc = COM{0.,-3./2.};
       static constexpr double cmumc = -1./6.;
 
       return
         +cmsms*norm(sym)
         +cmumu*norm(uncross)
         +cmcmc*norm(cross)
         +2.*real(cmsmu*sym*conj(uncross))
         +2.*real(cmsmc*sym*conj(cross))
         +2.*real(cmumc*uncross*conj(cross))
         ;
     }
   } // Anonymous Namespace
 
   // matrix element for W + Jets with W emission off central qqbar
   double ME_WCenqqbar_qq(
       HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar,
       std::vector<HLV> const & partons, bool /* aqlinepa */, bool /* aqlinepb */,
       bool qqbar_marker, int nabove, ParticleProperties const & wprop
   ){
     using helicity::plus;
     using helicity::minus;
 
     CLHEP::HepLorentzVector q1 = pa;
     for(int i = 0; i <= nabove; ++i){
       q1 -= partons[i];
     }
     const auto qq = split_into_lightlike(q1);
     // since q1.m2() < 0 the following assertion is always true
     // see documentation for split_into_lightlike
     assert(qq.second.e() < 0);
 
     HLV pq;
     HLV pqbar;
     if (qqbar_marker){
       pqbar = partons[nabove+1];
       pq = partons[nabove+2];}
     else{
       pq = partons[nabove+1];
       pqbar = partons[nabove+2];
     }
     const HLV p1 = partons.front();
     const HLV p4 = partons.back();
 
     // 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
     // also the choice in qqbar helicity.
     // the first helicity label is for aqlinepa == true,
     // the second one for aqlinepb == true
     // In principle the corresponding helicity should be flipped
     // if either of them is false, but we sum up all helicities,
     // so this has no effect in the end.
     const double amp_mm = amp_WCenqqbar_qq<minus, minus>(
       pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_mp = amp_WCenqqbar_qq<minus, plus>(
       pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_pm = amp_WCenqqbar_qq<plus, minus>(
       pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_pp = amp_WCenqqbar_qq<plus, plus>(
       pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
 
     // Divide by t channels, extremal + adjacent central vertex
     const double ta = (pa-p1).m2();
     const double t1 = q1.m2();
     const double t3 = (q1-pl-plbar-pq-pqbar).m2();
     const double tb = (p4-pb).m2();
 
     const double amp = WProp(plbar, pl, wprop)*(
       amp_mm+amp_mp+amp_pm+amp_pp
     )/(9.*4.*ta*t1*t3*tb);
 
     return amp;
   }
 
   // helper function for matrix element for W + Jets with central qqbar
   // W emitted off extremal parton
   // @TODO: code duplication with amp_WCenqqbar_qq
   template<Helicity h1a, Helicity hqqbar>
   double amp_W_Cenqqbar_qq(
       HLV const & pa, HLV const &  p1,
       HLV const & pb, HLV const &  pn,
       HLV const &  pq, HLV const &  pqbar,
       HLV const &  pl, HLV const &  plbar,
       HLV const &  q11, HLV const &  q12
   ){
     using std::norm;
 
     const COM crossed = M_cross<h1a, hqqbar>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
     );
     const COM uncrossed = M_qbar<h1a, hqqbar>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
     );
     const COM sym = M_sym<h1a, hqqbar>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
     );
 
     //Colour factors:
     static constexpr double cmsms = 3.;
     static constexpr double cmumu = 4./3.;
     static constexpr double cmcmc = 4./3.;
     static constexpr COM cmsmu = COM{0.,3./2.};
     static constexpr COM cmsmc = COM{0.,-3./2.};
     static constexpr double cmumc = -1./6.;
 
     return
       +cmsms*norm(sym)
       +cmumu*norm(uncrossed)
       +cmcmc*norm(crossed)
       +2.*real(cmsmu*sym*conj(uncrossed))
       +2.*real(cmsmc*sym*conj(crossed))
       +2.*real(cmumc*uncrossed*conj(crossed))
       ;
   }
 
   // matrix element for W + Jets with W emission *not* off central qqbar
   double ME_W_Cenqqbar_qq(
       HLV pa, HLV pb, HLV pl,HLV plbar,
       std::vector<HLV> partons, bool aqlinepa, bool aqlinepb,
       bool qqbar_marker, int nabove, int nbelow, bool forwards,
       ParticleProperties const & wprop
   ){
     using helicity::minus;
     using helicity::plus;
 
     if (!forwards){ //If Emission from Leg a instead, flip process.
       std::swap(pa, pb);
       std::reverse(partons.begin(),partons.end());
       std::swap(aqlinepa, aqlinepb);
       qqbar_marker = !qqbar_marker;
       std::swap(nabove,nbelow);
     }
 
     HLV q1=pa;
     for(int i=0;i<nabove+1;++i){
       q1-=partons.at(i);
     }
     const auto qq = split_into_lightlike(q1);
 
     HLV pq;
     HLV pqbar;
     if (qqbar_marker){
       pqbar = partons[nabove+1];
       pq = partons[nabove+2];}
     else{
       pq = partons[nabove+1];
       pqbar = partons[nabove+2];}
 
     // we assume that the W is emitted off a quark line
     // if this is not the case, we have to apply CP conjugation,
     // which is equivalent to swapping lepton and antilepton momenta
     if(aqlinepb) std::swap(pl, plbar);
 
     const HLV p1 = partons.front();
     const HLV pn = partons.back();
 
     // helicity labels are for aqlinepa == aqlineb == false,
     // In principle the helicities should be flipped
     // if either of them is true, but we sum up all helicities,
     // so this has no effect in the end.
     const double amp_mm = amp_W_Cenqqbar_qq<minus, minus>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_mp = amp_W_Cenqqbar_qq<minus, plus>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_pm = amp_W_Cenqqbar_qq<plus, minus>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
     const double amp_pp = amp_W_Cenqqbar_qq<plus, plus>(
       pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
     );
 
     // Divide by t channels, extremal + adjacent central vertex
     const double ta = (pa-p1).m2();
     const double t1 = q1.m2();
     const double t3 = (q1 - pq - pqbar).m2();
     const double tb = (pn+pl+plbar-pb).m2();
 
     const double amp= WProp(plbar, pl, wprop)*(
       amp_mm+amp_mp+amp_pm+amp_pp
     )/(9.*4.*ta*t1*t3*tb);
 
     return amp;
   }
 } // namespace currents
 } // namespace HEJ