diff --git a/.clang-tidy b/.clang-tidy
index 0dc056c..f96b147 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -1,38 +1,39 @@
 Checks: '
   *,
   -cert-dcl21-cpp,
   -clang-analyzer-optin.cplusplus.VirtualCall,
   -cppcoreguidelines-avoid-magic-numbers,
   -cppcoreguidelines-pro-bounds-array-to-pointer-decay,
   -cppcoreguidelines-pro-bounds-constant-array-index,
   -fuchsia*,
   -google-build-using-namespace,
   -google-explicit-constructor,
   -google-readability*,
   -google-runtime-int,
   -google-runtime-references,
   -hicpp*,
   -llvm-header-guard,
   -modernize-use-trailing-return-type,
   -readability-braces-around-statements,
   -readability-uppercase-literal-suffix,
+  -bugprone-macro-parentheses,
   '
 CheckOptions:
   - key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor
     value: '1'
   - key: llvm-namespace-comment.ShortNamespaceLines
     value: '10'
   - key: misc-non-private-member-variables-in-classes.IgnoreClassesWithAllMemberVariablesBeingPublic
     value: '1'
   - key: performance-move-const-arg.CheckTriviallyCopyableMove
     value: '0'
   - key: readability-identifier-naming.GlobalConstantCase
     value: UPPER_CASE
   - key: readability-identifier-naming.GlobalVariableCase
     value: UPPER_CASE
   - key: readability-identifier-naming.PrivateMemberSuffix
     value: '_'
   - key: readability-implicit-bool-conversion.AllowIntegerConditions
     value: '1'
   - key: readability-magic-numbers.IgnoreAllFloatingPointValues
     value: '1'
diff --git a/include/HEJ/PDF.hh b/include/HEJ/PDF.hh
index 44a09b2..a1c2e39 100644
--- a/include/HEJ/PDF.hh
+++ b/include/HEJ/PDF.hh
@@ -1,81 +1,81 @@
 /** \file
  *
  * \brief Contains all the necessary classes and functions for interaction with PDFs.
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <array>
 #include <cstddef>
 #include <memory>
 
 #include "HEJ/PDG_codes.hh"
 
 namespace LHAPDF {
   class PDF;
 }
 
 namespace HEJ {
   //! Class for interaction with a PDF set
   class PDF {
   public:
     /**
      * \brief PDF Constructor
      * @param id         Particle ID according to PDG
      * @param beam1      Particle ID of particle in beam 1
      * @param beam2      Particle ID of particle in beam 2
      */
     PDF(int id, ParticleID beam1, ParticleID beam2);
 
     /**
      * \brief Calculate the pdf value x*f(x, q)
      * @param beam_idx  Beam number (0 or 1)
      * @param x         Momentum fraction
      * @param q         Energy scale
      * @param id        PDG particle id
      * @returns         x*f(x, q)
      *
      * Returns 0 if x or q are outside the range covered by the PDF set
      */
     double pdfpt(std::size_t beam_idx, double x, double q, ParticleID id) const;
 
     /**
      * \brief Value of the strong coupling \f$\alpha_s(q)\f$ at the given scale
      * @param q         Renormalisation scale
      * @returns         Value of the strong coupling constant
      */
     double Halphas(double q) const;
 
     //! Check if the energy scale is within the range covered by the PDF set
     /**
      * @param q        Energy Scale
      * @returns        true if q is within the covered range, false otherwise
      */
     bool inRangeQ(double q) const;
 
     //! Check if the momentum fraction is within the range covered by the PDF set
     /**
      * @param x        Momentum Fraction
      * @returns        true if x is within the covered range, false otherwise
      */
     bool inRangeX(double x) const;
 
     //! PDF id of the current set
     int id() const;
 
     PDF(PDF const & other) = delete; //!< Copy not possible
     PDF & operator=(PDF const & other) = delete; //!< Copy not possible
-    PDF(PDF && other); //!< Moving allowed
-    PDF & operator=(PDF && other); //!< Moving allowed
+    PDF(PDF && other) noexcept; //!< Moving allowed
+    PDF & operator=(PDF && other) noexcept; //!< Moving allowed
     ~PDF();
 
   private:
     //! @internal unique_ptr does not allow copy
     std::unique_ptr<LHAPDF::PDF> pdf_;
 
-    std::array<int, 2> beamtype_;
+    std::array<int, 2> beamtype_{};
   };
 } // namespace HEJ
diff --git a/src/PDF.cc b/src/PDF.cc
index 4cdd837..34f2e10 100644
--- a/src/PDF.cc
+++ b/src/PDF.cc
@@ -1,76 +1,76 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/PDF.hh"
 
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
 #include <stdexcept>
 #include <string>
 
 #include "LHAPDF/LHAPDF.h"
 
 namespace HEJ {
 
   namespace {
 
     int to_beam(ParticleID id){
       if(std::abs(id) == pid::proton){
         return (id > 0)?1:-1;
       }
       throw std::invalid_argument(
           "unknown beam type: " + std::to_string(id)
       );
     }
 
   } // namespace
 
   PDF::PDF(int id, ParticleID beam1, ParticleID beam2):
     pdf_{LHAPDF::mkPDF(id)},
     beamtype_{{to_beam(beam1), to_beam(beam2)}}
   {}
 
   double PDF::pdfpt(size_t beam_idx, double x, double q, ParticleID id) const{
     if(!(inRangeQ(q) && inRangeX(x))) return 0.;
 
     if(id == pid::gluon){
       return pdf_->xfxQ(pid::gluon,x,q);
     }
     if(std::abs(id) <= pid::top){
       return pdf_->xfxQ(id*beamtype_[beam_idx],x,q);
     }
     std::cerr << "particle type unknown: "<< id << std::endl;
     return 0.0;
 
   }
 
   double PDF::Halphas(double q) const{
     double as = pdf_->alphasQ(q);
     if (std::isnan(as) || as > 0.5) {
       as = 0.5;
     }
     return as;
   }
 
   int PDF::id() const{
     return pdf_->lhapdfID();
   }
 
   bool PDF::inRangeQ(double q) const{
     return pdf_->inRangeQ(q);
   }
 
   bool PDF::inRangeX(double x) const{
     return pdf_->inRangeX(x);
   }
 
   // can't be in header since we forward declare LHAPDF::PDF
   PDF::~PDF() = default;
-  PDF::PDF(PDF && other) = default;
-  PDF & PDF::operator=(PDF && other) = default;
+  PDF::PDF(PDF && other) noexcept = default;
+  PDF & PDF::operator=(PDF && other) noexcept = default;
 
 
 } // namespace HEJ
diff --git a/src/jets.cc b/src/jets.cc
index 265ea9b..88efa02 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,357 +1,355 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/jets.hh"
 
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 
 #include "HEJ/Constants.hh"
 
 // generated headers
 #include "HEJ/currents/j_j.hh"
 #include "HEJ/currents/j_jqqbar.hh"
 #include "HEJ/currents/j_qqbar_j.hh"
 #include "HEJ/currents/juno_j.hh"
 
 namespace {
   // short hand for math functions
-  using std::abs;
-  using std::conj;
-  using std::sqrt;
-} // Anonymous Namespace
+   using std::conj;
+  } // Anonymous Namespace
 
 namespace HEJ {
 namespace currents {
 
   // Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
   // @TODO: this is not a current and should be moved somewhere else
   double K_g(double p1minus, double paminus) {
     return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
   }
   double K_g(
     HLV const & pout,
     HLV const & pin
     ) {
     if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
     return K_g(pout.minus(), pin.minus());
   }
 
   namespace {
   //@{
   /**
     * @brief Pure Jet FKL Contributions, function to handle all incoming types.
     * @param p1out             Outgoing Particle 1.
     * @param p1in              Incoming Particle 1.
     * @param p2out             Outgoing Particle 2
     * @param p2in              Incoming Particle 2
     *
     * Calculates j_\mu    j^\mu.
     * Handles all possible incoming states. Helicity doesn't matter since we sum
     * over all of them. In addition, we only have to distinguish between the two
     * possibilities of contracting same-helicity or opposite-helicity currents.
     */
     double j_j(HLV const & p1out, HLV const & p1in,
         HLV const & p2out, HLV const & p2in
     ){
       using helicity::plus;
       using helicity::minus;
 
       COM const Mp = HEJ::j_j<minus, plus>(p1in, p1out, p2in, p2out);
       COM const Mm = HEJ::j_j<minus, minus>(p1in, p1out, p2in, p2out);
 
       double const sst = std::norm(Mm) + std::norm(Mp);
 
       HLV const q1=p1in-p1out;
       HLV const q2=-(p2in-p2out);
 
       // Multiply by Cf^2 (colour) * 2 (helicities)
       return 2.*C_F*C_F*(sst)/(q1.m2()*q2.m2());
     }
   } // Anonymous Namespace
   double ME_qQ(HLV const & p1out, HLV const & p1in,
                HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in);
   }
 
   double ME_qQbar(HLV const & p1out, HLV const & p1in,
                   HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in);
   }
 
   double ME_qbarQbar(HLV const & p1out, HLV const & p1in,
                      HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in);
   }
 
   double ME_qg(HLV const & p1out, HLV const & p1in,
                HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/C_F;
   }
 
   double ME_qbarg(HLV const & p1out, HLV const & p1in,
                   HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(C_F);
   }
 
   double ME_gg(HLV const & p1out, HLV const & p1in,
                HLV const & p2out, HLV const & p2in
   ){
     return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(C_F*C_F);
   }
   //@}
 
   namespace {
     template<Helicity h1, Helicity hg>
     double amp_juno_j(
       HLV const & pa, HLV const & pb,
       HLV const & pg, HLV const & p1, HLV const & p2
     ) {
       // TODO: code duplication with Wjets
       const COM u1 = U1_j<h1, hg>(pa,p1,pb,p2,pg);
       const COM u2 = U2_j<h1, hg>(pa,p1,pb,p2,pg);
       const COM l  = L_j<h1, hg>(pa,p1,pb,p2,pg);
 
       const COM x = u1 - l;
       const COM y = u2 + l;
 
       return C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
     }
 
     double juno_j(HLV const & pg, HLV const & p1out,
       HLV const & p1in, HLV const & p2out, HLV const & p2in
     ){
       using helicity::minus;
       using helicity::plus;
 
       const HLV q2 = p2out - p2in;       // Bottom End
       const HLV qg = p1in - p1out - pg;  // Extra bit post-gluon
 
 
       // only calculate half of the helicity amplitudes,
       // the remaining ones follow from CP symmetry
       const double amm = amp_juno_j<minus, minus>(p1in, p2in, pg, p1out, p2out);
       const double amp = amp_juno_j<minus,  plus>(p1in, p2in, pg, p1out, p2out);
       const double apm = amp_juno_j< plus, minus>(p1in, p2in, pg, p1out, p2out);
       const double app = amp_juno_j< plus,  plus>(p1in, p2in, pg, p1out, p2out);
 
       constexpr double hel_fac = 2.;
 
       return hel_fac/(4.*C_A*C_A)*(amm + amp + apm + app)/(q2.m2()*qg.m2());
     }
   } // Anonymous Namespace
 
   //Unordered bits for pure jet
   double ME_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
                     HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in);
   }
 
   double ME_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
                        HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in);
   }
 
   double ME_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
                        HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in);
   }
 
   double ME_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
                           HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in);
   }
 
   double ME_unob_qg( HLV const & pg, HLV const & p1out, HLV const & p1in,
                      HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
   }
 
   double ME_unob_qbarg(HLV const & pg, HLV const & p1out, HLV const & p1in,
                        HLV const & p2out, HLV const & p2in
   ){
     return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
   }
 
   namespace {
     // helicity amplitudes for j jqqbar contraction
     template<Helicity h1, Helicity hg>
     double amp_j_jqqbar(
       HLV const & pa, HLV const & pb,
       HLV const & p1, HLV const & p2, HLV const & p3
     ) {
       // TODO: code duplication with Wjets.cc
 
       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 = j_qggm1<h1,hg>(pb,p2,p3,pa,p1);
       const COM m2 = j_qggm2<h1,hg>(pb,p2,p3,pa,p1);
       const COM m3 = j_qggm3<h1,hg>(pb,p2,p3,pa,p1);
 
       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))
         ;
     }
 
     //Now the function to give helicity/colour sum/average
     double MqgtqQQ(HLV const & pa, HLV const & pb,
                    HLV const & p1, HLV const & p2, HLV const & p3
     ){
       using helicity::minus;
       using helicity::plus;
 
       const double Mmmm = amp_j_jqqbar<minus, minus>(pa, pb, p1, p2, p3);
       const double Mmmp = amp_j_jqqbar<minus,  plus>(pa, pb, p1, p2, p3);
       const double Mpmm = amp_j_jqqbar< plus, minus>(pa, pb, p1, p2, p3);
       const double Mpmp = amp_j_jqqbar< plus,  plus>(pa, pb, p1, p2, p3);
 
       // Factor of 2 for the helicity for conjugate configurations
       return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2();
     }
   } // Anonymous Namespace
 
   // Extremal qqx
   double ME_Exqqx_qbarqQ(HLV const & pgin, HLV const &  pqout,
                          HLV const & pqbarout, HLV const & p2out, HLV const & p2in
   ){
     return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
   }
 
   double ME_Exqqx_qqbarQ(HLV const & pgin, HLV const & pqout,
                          HLV const & pqbarout, HLV const & p2out, HLV const & p2in
   ){
     return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
   }
 
   double ME_Exqqx_qbarqg(HLV const & pgin, HLV const & pqout,
                          HLV const & pqbarout, HLV const & p2out, HLV const & p2in
   ){
     return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/C_F;
   }
 
   double ME_Exqqx_qqbarg(HLV const & pgin, HLV const & pqout,
                          HLV const & pqbarout, HLV const & p2out, HLV const & p2in
   ){
     return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/C_F;
   }
 
   namespace {
     // helicity amplitudes for matrix element with central qqbar
     template<Helicity h1a, Helicity h4b>
     double amp_Cenqqx_qq(
       HLV const & pa, HLV const &  p1,
       HLV const & pb, HLV const &  p4,
       HLV const &  pq, HLV const &  pqbar,
       HLV const &  q11, HLV const &  q12
     ){
       using std::norm;
 
       const COM sym = M_sym<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, q11, q12
       );
       const COM cross = M_cross<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, q11, q12
       );
       const COM uncross = M_qbar<h1a, h4b>(
         pa, p1, pb, p4, pq, pqbar, 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
 
   double ME_Cenqqx_qq(
     HLV const & pa, HLV const & pb,
     std::vector<HLV> const & partons,
     bool /* aqlinepa */, bool /* aqlinepb */,
     const bool qqxmarker, const std::size_t nabove
   ){
     using helicity::plus;
     using helicity::minus;
 
     CLHEP::HepLorentzVector q1 = pa;
     for(std::size_t 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 = partons[nabove+1];
     HLV pqbar = partons[nabove+2];
     if(qqxmarker) std::swap(pq, pqbar);
 
     const HLV p1 = partons.front();
     const HLV pn = partons.back();
 
     // 8  helicity choices, but only 4 independent ones
     //(complex conjugation related).
     // In principle, the proper helicity labeling depends on
     // whether we have antiquark lines (aqlinea and aqlineb).
     // However, this only corresponds to a relabeling,
     // which doesn't change the sum over all helicities below.
     const double amp_mm = amp_Cenqqx_qq<minus, minus>(
       pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
     );
     const double amp_mp = amp_Cenqqx_qq<minus, plus>(
       pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
     );
     const double amp_pm = amp_Cenqqx_qq<plus, minus>(
       pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
     );
     const double amp_pp = amp_Cenqqx_qq<plus, plus>(
       pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
     );
 
     //Result (averaged, without coupling or t-channel props). Factor of
     //2 for the 4 helicity configurations I didn't work out explicitly
     const HLV q3 = q1 - pq - pqbar;
 
     return (2.*(amp_mm+amp_mp+amp_pm+amp_pp)/9./4.) /
     ((pa-p1).m2()*(pb-pn).m2()*q1.m2()*q3.m2());
   }
 
 } // namespace currents
 } // namespace HEJ