diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh
index d8f0aa2..21cf292 100644
--- a/include/HEJ/Particle.hh
+++ b/include/HEJ/Particle.hh
@@ -1,225 +1,227 @@
 /**
  * \file Particle.hh
  * \brief Contains the particle struct
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <utility>
 #include <vector>
 
 #include "fastjet/PseudoJet.hh"
 
 #include "HEJ/PDG_codes.hh"
 #include "HEJ/optional.hh"
 
 namespace HEJ {
 
   using Colour = std::pair<int,int>;
 
   //! Class representing a particle
   struct Particle {
     //! particle type
     ParticleID type = pid::unspecified;
     //! particle momentum
     fastjet::PseudoJet p;
     //! (optional) colour & anti-colour
     optional<Colour> colour;
 
     //! get rapidity
     double rapidity() const{
       return p.rapidity();
     }
     //! get transverse momentum
     double perp() const{
       return p.perp();
     }
     //! get transverse momentum
     double pt() const{
       return perp();
     }
     //! get momentum in x direction
     double px() const{
       return p.px();
     }
     //! get momentum in y direction
     double py() const{
       return p.py();
     }
     //! get momentum in z direction
     double pz() const{
       return p.pz();
     }
     //! get energy
     double E() const{
       return p.E();
     }
     //! get mass
     double m() const{
       return p.m();
     }
   };
 
   //! Functor to compare rapidities
   /**
    *  This can be used whenever a rapidity comparison function is needed,
    *  for example in many standard library functions.
    *
    *  @see pz_less
    */
   struct rapidity_less{
     template<class FourVector>
     bool operator()(FourVector const & p1, FourVector const & p2){
       return p1.rapidity() < p2.rapidity();
     }
   };
 
   //! Functor to compare momenta in z direction
   /**
    *  This can be used whenever a pz comparison function is needed,
    *  for example in many standard library functions.
    *
    *  @see rapidity_less
    */
   struct pz_less{
     template<class FourVector>
     bool operator()(FourVector const & p1, FourVector const & p2){
       return p1.pz() < p2.pz();
     }
   };
 
   //! Convert a vector of Particles to a vector of particle momenta
   inline
   std::vector<fastjet::PseudoJet> to_PseudoJet(
       std::vector<Particle> const & v
   ){
     std::vector<fastjet::PseudoJet> result;
     result.reserve(v.size());
     for(auto && sp: v) result.emplace_back(sp.p);
     return result;
   }
 
   //! Check if a particle is a parton, i.e. quark, antiquark, or gluon
   inline
-  bool is_parton(Particle const & p){
+  constexpr bool is_parton(Particle const & p){
     return is_parton(p.type);
   }
 
   //! Check if a particle is a quark
   inline
-  bool is_quark(Particle const & p){
+  constexpr bool is_quark(Particle const & p){
     return is_quark(p.type);
   }
 
   //! Check if a particle is an anti-quark
   inline
-  bool is_antiquark(Particle const & p){
+  constexpr bool is_antiquark(Particle const & p){
     return is_antiquark(p.type);
   }
 
   //! Check if a particle is a quark or anit-quark
   inline
-  bool is_anyquark(Particle const & p){
+  constexpr bool is_anyquark(Particle const & p){
     return is_anyquark(p.type);
   }
 
   /**
    * \brief Function to determine if particle is a lepton
    * @param p         the particle
    * @returns         true if the particle is a lepton, false otherwise
    */
   inline
   constexpr bool is_lepton(Particle const & p){
     return is_lepton(p.type);
   }
 
   /**
    * \brief Function to determine if particle is an antilepton
    * @param p         the particle
    * @returns         true if the particle is an antilepton, false otherwise
    */
   inline
   constexpr bool is_antilepton(Particle const & p){
     return is_antilepton(p.type);
   }
 
   /**
    * \brief Function to determine if particle is an (anti-)lepton
    * @param p         the particle
    * @returns         true if the particle is a lepton or antilepton, false otherwise
    */
   inline
   constexpr bool is_anylepton(Particle const & p){
     return is_anylepton(p.type);
   }
 
   /**
    * \brief Function to determine if particle is a neutrino
    * @param p         the particle
    * @returns         true if the particle is a neutrino, false otherwise
    */
   inline
   constexpr bool is_neutrino(Particle const & p){
     return is_neutrino(p.type);
   }
 
   /**
    * \brief Function to determine if particle is an antineutrino
    * @param p         the particle
    * @returns         true if the particle is an antineutrino, false otherwise
    */
   inline
   constexpr bool is_antineutrino(Particle const & p){
     return is_antineutrino(p.type);
   }
 
   /**
    * \brief Function to determine if particle is an (anti-)neutrino
    * @param p         the particle
    * @returns         true if the particle is a neutrino or antineutrino, false otherwise
    */
   inline
   constexpr bool is_anyneutrino(Particle const & p){
     return is_anyneutrino(p.type);
   }
 
   //! Check if a particle is a photon, W or Z boson
-  inline bool is_AWZ_boson(Particle const & particle){
+  inline
+  constexpr bool is_AWZ_boson(Particle const & particle){
     return is_AWZ_boson(particle.type);
   }
 
   //! Check if a particle is a photon, W, Z, or Higgs boson
-  inline bool is_AWZH_boson(Particle const & particle){
+  inline
+  constexpr bool is_AWZH_boson(Particle const & particle){
     return is_AWZH_boson(particle.type);
   }
 
   //! Extract all partons from a vector of particles
   inline
   std::vector<Particle> filter_partons(
       std::vector<Particle> const & v
   ){
     std::vector<Particle> result;
     result.reserve(v.size());
     std::copy_if(
         begin(v), end(v), std::back_inserter(result),
         [](Particle const & p){ return is_parton(p); }
     );
     return result;
   }
 
   //! Extract all AWZH bosons from a vector of particles
   inline
   std::vector<Particle> filter_AWZH_bosons(
       std::vector<Particle> const & v
   ){
     std::vector<Particle> result;
     std::copy_if(
         begin(v), end(v), std::back_inserter(result),
         [](Particle const & p){ return is_AWZH_boson(p); }
     );
     return result;
   }
 } // namespace HEJ
diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
index 9723543..b53ba36 100644
--- a/include/HEJ/event_types.hh
+++ b/include/HEJ/event_types.hh
@@ -1,108 +1,109 @@
 /** \file
  *  \brief Define different types of events.
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <string>
 
 #include "HEJ/exceptions.hh"
 
 namespace HEJ {
 
   //! Namespace for event types
   namespace event_type {
     //! Possible event types
     enum EventType: std::size_t {
       non_resummable = 0,     //!< event configuration not covered by All Order resummation
       bad_final_state = 1,    //!< event with an unsupported final state
       no_2_jets = 2,          //!< event with less than two jets
       FKL = 4,                //!< FKL-type event
       unordered_backward = 8, //!< event with unordered backward emission
       unordered_forward = 16, //!< event with unordered forward emission
       extremal_qqxb = 32,     //!< event with a backward extremal qqbar
       extremal_qqxf = 64,     //!< event with a forward extremal qqbar
       central_qqx = 128,      //!< event with a central qqbar
       unob = unordered_backward,   //!< alias for unordered_backward
       unof = unordered_forward,    //!< alias for unordered_forward
       qqxexb = extremal_qqxb,      //!< alias for extremal_qqxb
       qqxexf = extremal_qqxf,      //!< alias for extremal_qqxf
       qqxmid = central_qqx,        //!< alias for central_qqx
       first_type = non_resummable, //!< alias for non_resummable
       last_type = central_qqx      //!< alias for central_qqx
     };
 
     //! Event type names
     /**
      * For example, name(FKL) is the string "FKL"
      */
-    inline std::string name(EventType type) {
+    inline
+    std::string name(EventType type) {
       switch(type) {
       case FKL:
         return "FKL";
       case unordered_backward:
         return "unordered backward";
       case unordered_forward:
         return "unordered forward";
       case extremal_qqxb:
         return "extremal qqbar backward";
       case extremal_qqxf:
         return "extremal qqbar forward";
       case central_qqx:
         return "central qqbar";
       case non_resummable:
         return "non-resummable";
       case no_2_jets:
         return "no 2 jets";
       case bad_final_state:
         return "bad final state";
       default:
         throw std::logic_error{"Unreachable"};
       }
     }
 
     //! Returns True for a HEJ \ref event_type::EventType "EventType"
     inline
-    bool is_resummable(EventType type) {
+    constexpr bool is_resummable(EventType type) {
       switch(type) {
       case FKL:
       case unordered_backward:
       case unordered_forward:
       case extremal_qqxb:
       case extremal_qqxf:
       case central_qqx:
         return true;
       default:
         return false;
       }
     }
 
     //! Returns True for an unordered \ref event_type::EventType "EventType"
     inline
-    bool is_uno(EventType type) {
+    constexpr bool is_uno(EventType type) {
       return type == unordered_backward || type == unordered_forward;
     }
 
     //! Returns True for an extremal_qqx \ref event_type::EventType "EventType"
     inline
-    bool is_ex_qqx(EventType type) {
+    constexpr bool is_ex_qqx(EventType type) {
       return type == extremal_qqxb || type == extremal_qqxf;
     }
 
     //! Returns True for an central_qqx \ref event_type::EventType "EventType"
     inline
-    bool is_mid_qqx(EventType type) {
+    constexpr bool is_mid_qqx(EventType type) {
       return type == central_qqx;
     }
 
     //! Returns True for any qqx event \ref event_type::EventType "EventType"
     inline
-    bool is_qqx(EventType type) {
+    constexpr bool is_qqx(EventType type) {
       return is_ex_qqx(type) || is_mid_qqx(type);
     }
 
   } // namespace event_type
 } // namespace HEJ
diff --git a/include/HEJ/utility.hh b/include/HEJ/utility.hh
index 6f94205..997dce8 100644
--- a/include/HEJ/utility.hh
+++ b/include/HEJ/utility.hh
@@ -1,103 +1,104 @@
 /**
  * \file
  * \brief Contains various utilities
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <array>
 #include <memory>
 #include <string>
 
 #include "boost/core/demangle.hpp"
 
 #include "fastjet/PseudoJet.hh"
 
 namespace HEJ {
 
   inline
   std::string join(
       std::string const & /* delim */
   ){
     return "";
   }
 
   inline
   std::string join(
       std::string const & /* delim */, std::string const & str
   ){
     return str;
   }
 
   //! Join strings with a delimiter
   /**
    *   @param delim      Delimiter to be put between consecutive strings
    *   @param first      First string
    *   @param second     Second string
    *   @param rest       Remaining strings
    */
   template<typename... Strings>
   std::string join(
       std::string const & delim,
       std::string const & first, std::string const & second,
       Strings&&... rest
   ){
     return join(delim, first + delim + second, std::forward<Strings>(rest)...);
   }
 
   //! Return the name of the argument's type
   template<typename T>
   std::string type_string(T&& /*unused*/){
     return boost::core::demangle(typeid(T).name());
   }
 
   //! Eliminate compiler warnings for unused variables
   template<typename... T>
   constexpr void ignore(T&&... /*unused*/) {}
 
   //! Check whether two doubles are closer than ep > 0 to each other
   inline
-  bool nearby_ep(double a, double b, double ep){
+  constexpr bool nearby_ep(double a, double b, double ep){
     assert(ep > 0);
     return std::abs(a-b) < ep;
   }
 
   //! Check whether all components of two PseudoJets are closer than ep to each other
   inline
   bool nearby_ep(
       fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb,
       double ep
   ){
     assert(ep > 0);
     for(size_t i = 0; i < 4; ++i){
       if(!nearby_ep(pa[i], pb[i], ep)) return false;
     }
     return true;
   }
 
   inline
   bool nearby(
-      fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb, double const norm = 1.
+      fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb,
+      double const norm = 1.
   ){
     return nearby_ep(pa, pb, 1e-7*norm);
   }
 
   namespace detail {
     template<typename T, std::size_t N, std::size_t... Ns>
     struct ArrayTag{
       using type = typename ArrayTag<std::array<T, N>, Ns...>::type;
     };
     template<typename T, std::size_t N>
     struct ArrayTag<T, N> {
       using type = std::array<T, N>;
     };
   }
   // helper for multidimensional std::array, for example
   // MultiArray<T, N1, N2> = std::array<std::array<T, N1>, N2>
   template<typename T, std::size_t N, std::size_t... Ns>
   using MultiArray = typename detail::ArrayTag<T, N, Ns...>::type;
 
 } // namespace HEJ