diff --git a/.github/workflows/lcg_linux_with_podio.yml b/.github/workflows/lcg_linux_with_podio.yml
index 22883f417..b4cb7c217 100644
--- a/.github/workflows/lcg_linux_with_podio.yml
+++ b/.github/workflows/lcg_linux_with_podio.yml
@@ -10,13 +10,9 @@ jobs:
       matrix:
         LCG: ["dev3/x86_64-el9-clang16-opt",
               "dev4/x86_64-el9-clang16-opt",
-              "dev4/x86_64-el9-gcc13-opt"]
+              "dev4/x86_64-el9-gcc13-opt",
+              "LCG_106/x86_64-el9-gcc13-opt"]
         CXX_STANDARD: [20]
-        include:
-          - LCG: "LCG_106/x86_64-el9-gcc13-opt"
-            CXX_STANDARD: 17
-          - LCG: "dev4/x86_64-ubuntu2004-gcc9-opt"
-            CXX_STANDARD: 17
     steps:
     - uses: actions/checkout@v3
     - uses: cvmfs-contrib/github-action-cvmfs@v3
diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml
index 11ee0a715..aaa6cb85a 100644
--- a/.github/workflows/pre-commit.yml
+++ b/.github/workflows/pre-commit.yml
@@ -19,7 +19,7 @@ jobs:
           cd podio
           mkdir build install
           cd build
-          cmake .. -DCMAKE_CXX_STANDARD=17 \
+          cmake .. -DCMAKE_CXX_STANDARD=20 \
             -DENABLE_SIO=ON \
             -DBUILD_TESTING=OFF \
             -DCMAKE_INSTALL_PREFIX=$(pwd)/../install \
@@ -40,7 +40,7 @@ jobs:
           mkdir build
           cd build
           cmake .. -DENABLE_SIO=ON \
-            -DCMAKE_CXX_STANDARD=17 \
+            -DCMAKE_CXX_STANDARD=20 \
             -DCMAKE_CXX_FLAGS=" -fdiagnostics-color=always -Werror "\
             -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \
             -DUSE_EXTERNAL_CATCH2=OFF
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1a1f8f8b1..3fd2b5a36 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -38,9 +38,9 @@ if(NOT CMAKE_CONFIGURATION_TYPES)
 endif()
 
 # - Define the C++ Standard to use (Simplest Possible)
-set(CMAKE_CXX_STANDARD 17 CACHE STRING "")
+set(CMAKE_CXX_STANDARD 20 CACHE STRING "")
 
-if(NOT CMAKE_CXX_STANDARD MATCHES "17|20")
+if(NOT CMAKE_CXX_STANDARD MATCHES "20|23")
   message(FATAL_ERROR "Unsupported C++ standard: ${CMAKE_CXX_STANDARD}")
 endif()
 
diff --git a/test/hepmc/edm4hep_testhepmc.cc b/test/hepmc/edm4hep_testhepmc.cc
index 353a554fe..4828cf49b 100644
--- a/test/hepmc/edm4hep_testhepmc.cc
+++ b/test/hepmc/edm4hep_testhepmc.cc
@@ -114,7 +114,7 @@ int main() {
     std::cout << "Converting particle with PDG ID: " << particle_i->pdg_id() << std::endl;
     std::cout << "\t and id: " << particle_i->id() << std::endl;
 
-    if (hepmcToEdmMap.find(particle_i->id()) == hepmcToEdmMap.end()) {
+    if (!hepmcToEdmMap.contains(particle_i->id())) {
       auto edm_particle = convert(particle_i);
       hepmcToEdmMap.insert({particle_i->id(), edm_particle});
     }
@@ -124,7 +124,7 @@ int main() {
     if (nullptr != prodvertex) {
 
       for (auto particle_mother : prodvertex->particles_in()) {
-        if (hepmcToEdmMap.find(particle_mother->id()) == hepmcToEdmMap.end()) {
+        if (!hepmcToEdmMap.contains(particle_mother->id())) {
           auto edm_particle = convert(particle_mother);
           hepmcToEdmMap.insert({particle_mother->id(), edm_particle});
         }
@@ -136,7 +136,7 @@ int main() {
     if (nullptr != prodvertex) {
 
       for (auto particle_daughter : prodvertex->particles_in()) {
-        if (hepmcToEdmMap.find(particle_daughter->id()) == hepmcToEdmMap.end()) {
+        if (!hepmcToEdmMap.contains(particle_daughter->id())) {
           auto edm_particle = convert(particle_daughter);
           hepmcToEdmMap.insert({particle_daughter->id(), edm_particle});
         }
diff --git a/test/utils/test_covmatrix_utils.cpp b/test/utils/test_covmatrix_utils.cpp
index 6fdb2f76b..4a387e766 100644
--- a/test/utils/test_covmatrix_utils.cpp
+++ b/test/utils/test_covmatrix_utils.cpp
@@ -14,10 +14,6 @@ TEST_CASE("CovarianceMatrix indexing", "[cov_matrix_utils]") {
 
   STATIC_REQUIRE(get_cov_dim(21) == 6);
   STATIC_REQUIRE(get_cov_dim(1) == 1);
-#if !__cpp_consteval
-  // This will fail to compile if we have consteval
-  REQUIRE_THROWS_AS(get_cov_dim(14), std::invalid_argument);
-#endif
 
   // clang-format off
   // For better interpretability of the tests below, these are the indices of a
diff --git a/utils/include/edm4hep/utils/cov_matrix_utils.h b/utils/include/edm4hep/utils/cov_matrix_utils.h
index 609cfd583..6b07c5f2e 100644
--- a/utils/include/edm4hep/utils/cov_matrix_utils.h
+++ b/utils/include/edm4hep/utils/cov_matrix_utils.h
@@ -49,11 +49,7 @@ namespace utils {
      *
      * @returns the dimension of the covariance matrix
      */
-#if cpp_consteval
     consteval std::size_t get_cov_dim(std::size_t N) {
-#else
-    constexpr std::size_t get_cov_dim(std::size_t N) {
-#endif
       switch (N) {
       case 21:
         return 6;
diff --git a/utils/include/edm4hep/utils/vector_utils.h b/utils/include/edm4hep/utils/vector_utils.h
index 48b5b6cf8..f2fccbc43 100644
--- a/utils/include/edm4hep/utils/vector_utils.h
+++ b/utils/include/edm4hep/utils/vector_utils.h
@@ -1,12 +1,6 @@
 #ifndef EDM4HEP_UTILS_VECTOR_HH
 #define EDM4HEP_UTILS_VECTOR_HH
 
-// These ultilies require concepts. If not available, use the fallback
-// vector_utils_legacy.h instead to capture most functionality.
-#if !__cpp_concepts
-#include <edm4hep/utils/vector_utils_legacy.h>
-#else
-
 #include <edm4hep/Vector3f.h>
 
 #include <cmath>
@@ -327,4 +321,3 @@ inline constexpr V operator/(const V& v, const double d) {
 } // namespace edm4hep
 
 #endif
-#endif
diff --git a/utils/include/edm4hep/utils/vector_utils_legacy.h b/utils/include/edm4hep/utils/vector_utils_legacy.h
deleted file mode 100644
index 3bb68fa1a..000000000
--- a/utils/include/edm4hep/utils/vector_utils_legacy.h
+++ /dev/null
@@ -1,305 +0,0 @@
-#ifndef EDM4HEP_UTILS_VECTOR_LEGACY_HH
-#define EDM4HEP_UTILS_VECTOR_LEGACY_HH
-
-// This is the legacy implementation of vector_utils. If possible, use
-// the better vector_utils.h instead (if concepts are available).
-#if __cpp_concepts
-#include <edm4hep/utils/vector_utils.h>
-#else
-
-#include <edm4hep/Vector2f.h>
-#include <edm4hep/Vector2i.h>
-#include <edm4hep/Vector3d.h>
-#include <edm4hep/Vector3f.h>
-#include <edm4hep/Vector4f.h>
-
-#include <cmath>
-#include <tuple>
-#include <type_traits>
-
-namespace edm4hep {
-
-namespace utils {
-
-  inline double etaToAngle(const double eta) {
-    return std::atan(std::exp(-eta)) * 2.;
-  }
-
-  inline double angleToEta(const double theta) {
-    return -std::log(std::tan(0.5 * theta));
-  }
-
-  // Utility getters to accomodate different vector types
-  template <class V>
-  constexpr auto vector_x(const V& v) {
-    return v.x;
-  }
-
-  template <class V>
-  constexpr auto vector_y(const V& v) {
-    return v.y;
-  }
-
-  template <class V>
-  constexpr auto vector_z(const V& v) {
-    return v.z;
-  }
-
-  template <class V>
-  constexpr auto vector_t(const V& v) {
-    return v.t;
-  }
-
-  // 2D vector uses a,b instead of x,y
-  template <>
-  inline constexpr auto vector_x<edm4hep::Vector2i>(const edm4hep::Vector2i& v) {
-    return v.a;
-  }
-
-  template <>
-  inline constexpr auto vector_y<edm4hep::Vector2i>(const edm4hep::Vector2i& v) {
-    return v.b;
-  }
-
-  template <>
-  inline constexpr auto vector_x<edm4hep::Vector2f>(const edm4hep::Vector2f& v) {
-    return v.a;
-  }
-
-  template <>
-  inline constexpr auto vector_y<edm4hep::Vector2f>(const edm4hep::Vector2f& v) {
-    return v.b;
-  }
-  // no z-component for 2D vectors
-  template <>
-  inline constexpr auto vector_z<edm4hep::Vector2f>(const edm4hep::Vector2f&) {
-    return 0;
-  }
-
-  template <>
-  inline constexpr auto vector_z<edm4hep::Vector2i>(const edm4hep::Vector2i&) {
-    return 0;
-  }
-
-  namespace detail {
-    /// Helper struct to determine the underlying type of edm4hep vector types
-    template <typename V>
-    struct ValueTypeHelper {
-      using type = decltype(vector_x(std::declval<V>()));
-    };
-
-    /// Helper struct to determine whether a type is in a tuple of other types
-    template <typename T, typename>
-    struct TypeInTupleHelper : std::false_type {};
-
-    template <typename T, typename... Ts>
-    struct TypeInTupleHelper<T, std::tuple<Ts...>> : std::bool_constant<(std::is_same_v<T, Ts> || ...)> {};
-
-    template <typename T, typename Tuple>
-    constexpr static bool isInTuple = TypeInTupleHelper<T, Tuple>::value;
-  } // namespace detail
-
-  /// Type alias that returns the underlying type of edm4hep
-  template <typename V>
-  using ValueType = typename detail::ValueTypeHelper<V>::type;
-
-  using EDM4hepVectorTypes = std::tuple<edm4hep::Vector3f, edm4hep::Vector3d, edm4hep::Vector2i, edm4hep::Vector2f>;
-  using EDM4hepVector4DTypes = std::tuple<edm4hep::Vector4f>;
-  using EDM4hepVector3DTypes = std::tuple<edm4hep::Vector3f, edm4hep::Vector3d>;
-  using EDM4hepVector2DTypes = std::tuple<edm4hep::Vector2f, edm4hep::Vector2i>;
-  using EDM4hepFloatVectorTypes = std::tuple<edm4hep::Vector2f, edm4hep::Vector3f, edm4hep::Vector3d>;
-
-  template <typename V>
-  using EnableIfEDM4hepVectorType = std::enable_if_t<detail::isInTuple<V, EDM4hepVectorTypes>, bool>;
-
-  template <typename V>
-  using EnableIfEDM4hepVector2DType = std::enable_if_t<detail::isInTuple<V, EDM4hepVector2DTypes>, bool>;
-
-  template <typename V>
-  using EnableIfEDM4hepVector3DType = std::enable_if_t<detail::isInTuple<V, EDM4hepVector3DTypes>, bool>;
-
-  template <typename V>
-  using EnableIfEdm4hepVector4DType = std::enable_if_t<detail::isInTuple<V, EDM4hepVector4DTypes>, bool>;
-  template <typename V>
-  using EnableIfEDM4hepFloatVectorType = std::enable_if_t<detail::isInTuple<V, EDM4hepFloatVectorTypes>, bool>;
-
-  // inline edm4hep::Vector2f VectorFromPolar(const double r, const double theta)
-  // {
-  //  return {r * sin(theta), r * cos(theta)};
-  //}
-
-  template <class V = edm4hep::Vector3f, typename = EnableIfEDM4hepVector3DType<V>>
-  V sphericalToVector(const double r, const double theta, const double phi) {
-    using FloatType = ValueType<V>;
-    const double sth = sin(theta);
-    const double cth = cos(theta);
-    const double sph = sin(phi);
-    const double cph = cos(phi);
-    const FloatType x = r * sth * cph;
-    const FloatType y = r * sth * sph;
-    const FloatType z = r * cth;
-    return {x, y, z};
-  }
-
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double anglePolar(const V& v) {
-    return std::atan2(std::hypot(vector_x(v), vector_y(v)), vector_z(v));
-  }
-
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double angleAzimuthal(const V& v) {
-    return std::atan2(vector_y(v), vector_x(v));
-  }
-
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double eta(const V& v) {
-    return angleToEta(anglePolar(v));
-  }
-
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double magnitude(const V& v) {
-    return std::hypot(vector_x(v), vector_y(v), vector_z(v));
-  }
-
-  template <class V, typename = EnableIfEDM4hepVector3DType<V>>
-  double magnitudeTransverse(const V& v) {
-    return std::hypot(vector_x(v), vector_y(v));
-  }
-
-  template <class V, typename = EnableIfEDM4hepVector3DType<V>>
-  double magnitudeLongitudinal(const V& v) {
-    return vector_z(v);
-  }
-
-  template <class V, typename = EnableIfEDM4hepFloatVectorType<V>>
-  V normalizeVector(const V& v, double norm = 1.) {
-    const double old = magnitude(v);
-    if (old == 0) {
-      return v;
-    }
-    return (norm > 0) ? v * norm / old : v * 0;
-  }
-
-  template <class V, typename = EnableIfEDM4hepVector3DType<V>>
-  constexpr V vectorTransverse(const V& v) {
-    return {vector_x(v), vector_y(v), 0};
-  }
-
-  template <class V, typename = EnableIfEDM4hepVector3DType<V>>
-  constexpr V vectorLongitudinal(const V& v) {
-    return {0, 0, vector_z(v)};
-  }
-
-  // Two vector functions
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double angleBetween(const V& v1, const V& v2) {
-    const double dot = v1 * v2;
-    if (dot == 0) {
-      return 0.;
-    }
-    return acos(dot / (magnitude(v1) * magnitude(v2)));
-  }
-
-  // Project v onto v1
-  template <class V, typename = EnableIfEDM4hepVectorType<V>>
-  double projection(const V& v, const V& v1) {
-    const double norm = magnitude(v1);
-    if (norm == 0) {
-      return magnitude(v);
-    }
-    return v * v1 / norm;
-  }
-
-} // namespace utils
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector2DType<V> = false>
-inline constexpr V operator+(const V& v1, const V& v2) {
-  const auto x = edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2);
-  const auto y = edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2);
-  return {x, y};
-}
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector3DType<V> = false>
-inline constexpr V operator+(const V& v1, const V& v2) {
-  const auto x = edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2);
-  const auto y = edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2);
-  const auto z = edm4hep::utils::vector_z(v1) + edm4hep::utils::vector_z(v2);
-  return {x, y, z};
-}
-
-template <typename V, edm4hep::utils::EnableIfEdm4hepVector4DType<V> = false>
-inline constexpr V operator+(const V& v1, const V& v2) {
-  const auto x = edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2);
-  const auto y = edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2);
-  const auto z = edm4hep::utils::vector_z(v1) + edm4hep::utils::vector_z(v2);
-  const auto t = edm4hep::utils::vector_t(v1) + edm4hep::utils::vector_t(v2);
-  return {x, y, z, t};
-}
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector2DType<V> = false>
-inline constexpr double operator*(const V& v1, const V& v2) {
-  return edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
-      edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2);
-}
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector3DType<V> = false>
-inline constexpr double operator*(const V& v1, const V& v2) {
-  return edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
-      edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2) +
-      edm4hep::utils::vector_z(v1) * edm4hep::utils::vector_z(v2);
-}
-
-template <typename V, edm4hep::utils::EnableIfEdm4hepVector4DType<V> = false>
-inline constexpr double operator*(const V& v1, const V& v2) {
-  return (edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
-          edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2) +
-          edm4hep::utils::vector_z(v1) * edm4hep::utils::vector_z(v2)) -
-      edm4hep::utils::vector_t(v1) * edm4hep::utils::vector_t(v2);
-}
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector2DType<V> = false>
-inline constexpr V operator*(const double d, const V& v) {
-  using VT = edm4hep::utils::ValueType<V>;
-  const VT x = d * edm4hep::utils::vector_x(v);
-  const VT y = d * edm4hep::utils::vector_y(v);
-  return {x, y};
-}
-
-template <typename V, edm4hep::utils::EnableIfEDM4hepVector3DType<V> = false>
-inline constexpr V operator*(const double d, const V& v) {
-  using VT = edm4hep::utils::ValueType<V>;
-  const VT x = d * edm4hep::utils::vector_x(v);
-  const VT y = d * edm4hep::utils::vector_y(v);
-  const VT z = d * edm4hep::utils::vector_z(v);
-  return {x, y, z};
-}
-
-template <typename V, edm4hep::utils::EnableIfEdm4hepVector4DType<V> = false>
-inline constexpr V operator*(const double d, const V& v) {
-  using VT = edm4hep::utils::ValueType<V>;
-  const VT x = d * edm4hep::utils::vector_x(v);
-  const VT y = d * edm4hep::utils::vector_y(v);
-  const VT z = d * edm4hep::utils::vector_z(v);
-  const VT t = d * edm4hep::utils::vector_t(v);
-  return {x, y, z, t};
-}
-
-template <typename V>
-inline constexpr V operator*(const V& v, const double d) {
-  return d * v;
-}
-
-template <typename V>
-inline constexpr V operator-(const V& v1, const V& v2) {
-  return v1 + (-1. * v2);
-}
-
-template <typename V>
-inline constexpr V operator/(const V& v, const double d) {
-  return (1. / d) * v;
-}
-
-} // namespace edm4hep
-
-#endif
-#endif // EDM4HEP_UTILS_VECTOR_LEGACY_HH
diff --git a/utils/src/ParticleIDUtils.cc b/utils/src/ParticleIDUtils.cc
index d4da00a36..2ec1b6952 100644
--- a/utils/src/ParticleIDUtils.cc
+++ b/utils/src/ParticleIDUtils.cc
@@ -29,7 +29,7 @@ ParticleIDMeta::ParticleIDMeta(const std::string& algName, const std::vector<std
 }
 
 std::optional<int> getParamIndex(const ParticleIDMeta& pidMetaInfo, const std::string& param) {
-  const auto nameIt = std::find(pidMetaInfo.paramNames.begin(), pidMetaInfo.paramNames.end(), param);
+  const auto nameIt = std::ranges::find(pidMetaInfo.paramNames, param);
   if (nameIt != pidMetaInfo.paramNames.end()) {
     return std::distance(pidMetaInfo.paramNames.begin(), nameIt);
   }