diff --git a/include/HEJ/utility.hh b/include/HEJ/utility.hh
index 997dce8..b57d7df 100644
--- a/include/HEJ/utility.hh
+++ b/include/HEJ/utility.hh
@@ -1,104 +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
+  //! Check whether two doubles are closer than ep >= 0 to each other
   inline
   constexpr bool nearby_ep(double a, double b, double ep){
-    assert(ep > 0);
+    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);
+    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.
   ){
     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