diff --git a/Core/include/Acts/Definitions/Units.hpp b/Core/include/Acts/Definitions/Units.hpp
index 65d7cee775d..9f6c2f8d187 100644
--- a/Core/include/Acts/Definitions/Units.hpp
+++ b/Core/include/Acts/Definitions/Units.hpp
@@ -255,13 +255,15 @@ ACTS_DEFINE_UNIT_LITERAL(mol)
 ///
 /// Unit constants are intentionally not listed.
 namespace PhysicalConstants {
-// Speed of light
+/// Speed of light
 constexpr double c = 1.0;
 /// Reduced Planck constant h/2*pi.
 ///
 /// Computed from CODATA 2018 constants to double precision.
 constexpr double hbar =
     6.582119569509066e-25 * UnitConstants::GeV * UnitConstants::s;
+/// Avogadro constant
+constexpr double kAvogadro = 6.02214076e23 / Acts::UnitConstants::mol;
 }  // namespace PhysicalConstants
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Material/AccumulatedMaterialSlab.hpp b/Core/include/Acts/Material/AccumulatedMaterialSlab.hpp
index 5f0013bcc3b..35d41da7378 100644
--- a/Core/include/Acts/Material/AccumulatedMaterialSlab.hpp
+++ b/Core/include/Acts/Material/AccumulatedMaterialSlab.hpp
@@ -46,7 +46,7 @@ class AccumulatedMaterialSlab {
   ///
   ///  Vacuum steps with a non-zero thickness can be added to account for holes
   ///  in material structures.
-  void accumulate(MaterialSlab slabAlongTrack, float pathCorrection = 1);
+  void accumulate(MaterialSlab slabAlongTrack, double pathCorrection = 1);
 
   /// Use the accumulated material to update the material variance
   ///
@@ -89,7 +89,7 @@ class AccumulatedMaterialSlab {
   /// If there have been additional calls to `.accumulate(...)` afterwards, the
   /// information is not part of the total average. The number of tracks is only
   /// updated on the call of `.trackAverage(...)`
-  std::pair<float, unsigned int> totalVariance() const;
+  std::pair<double, unsigned int> totalVariance() const;
 
  private:
   /// Averaged properties for a single track.
@@ -97,9 +97,9 @@ class AccumulatedMaterialSlab {
   /// Averaged properties over multiple tracks.
   MaterialSlab m_totalAverage;
   /// Averaged variance over multiple tracks.
-  float m_totalVariance = 0.0;
+  double m_totalVariance = 0;
   // Number of tracks contributing to the total average.
-  unsigned int m_totalCount = 0u;
+  unsigned int m_totalCount = 0;
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Material/Interactions.hpp b/Core/include/Acts/Material/Interactions.hpp
index 07ccbd7541f..daeccf73be7 100644
--- a/Core/include/Acts/Material/Interactions.hpp
+++ b/Core/include/Acts/Material/Interactions.hpp
@@ -9,11 +9,8 @@
 #pragma once
 
 #include "Acts/Definitions/PdgParticle.hpp"
-#include "Acts/Definitions/Units.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
 
-#include <cmath>
-
 namespace Acts {
 
 /// Compute the mean energy loss due to ionisation and excitation.
@@ -30,13 +27,13 @@ namespace Acts {
 ///
 /// where -dE/dx is given by the Bethe formula. The computations are valid
 /// for intermediate particle energies.
-float computeEnergyLossBethe(const MaterialSlab& slab, float m, float qOverP,
-                             float absQ);
+double computeEnergyLossBethe(const MaterialSlab& slab, double m, double qOverP,
+                              double absQ);
 /// Derivative of the Bethe energy loss with respect to q/p.
 ///
 /// @copydoc computeEnergyLossBethe
-float deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
-                                  float qOverP, float absQ);
+double deriveEnergyLossBetheQOverP(const MaterialSlab& slab, double m,
+                                   double qOverP, double absQ);
 
 /// Compute the most propable energy loss due to ionisation and excitation.
 ///
@@ -46,13 +43,13 @@ float deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
 /// the given properties and thickness as described by the mode of the
 /// Landau-Vavilov-Bichsel distribution. The computations are valid
 /// for intermediate particle energies.
-float computeEnergyLossLandau(const MaterialSlab& slab, float m, float qOverP,
-                              float absQ);
+double computeEnergyLossLandau(const MaterialSlab& slab, double m,
+                               double qOverP, double absQ);
 /// Derivative of the most probable ionisation energy loss with respect to q/p.
 ///
 /// @copydoc computeEnergyLossBethe
-float deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
-                                   float qOverP, float absQ);
+double deriveEnergyLossLandauQOverP(const MaterialSlab& slab, double m,
+                                    double qOverP, double absQ);
 
 /// Compute the Gaussian-equivalent sigma for the ionisation loss fluctuations.
 ///
@@ -61,20 +58,20 @@ float deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
 /// This is the sigma parameter of a Gaussian distribution with the same
 /// full-width-half-maximum as the Landau-Vavilov-Bichsel distribution. The
 /// computations are valid for intermediate particle energies.
-float computeEnergyLossLandauSigma(const MaterialSlab& slab, float m,
-                                   float qOverP, float absQ);
+double computeEnergyLossLandauSigma(const MaterialSlab& slab, double m,
+                                    double qOverP, double absQ);
 
 /// Compute the full with half maximum of landau energy loss distribution
 ///
 /// @see computeEnergyLossBethe for parameters description
-float computeEnergyLossLandauFwhm(const MaterialSlab& slab, float m,
-                                  float qOverP, float absQ);
+double computeEnergyLossLandauFwhm(const MaterialSlab& slab, double m,
+                                   double qOverP, double absQ);
 
 /// Compute q/p Gaussian-equivalent sigma due to ionisation loss fluctuations.
 ///
 /// @copydoc computeEnergyLossBethe
-float computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab, float m,
-                                         float qOverP, float absQ);
+double computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab, double m,
+                                          double qOverP, double absQ);
 
 /// Compute the mean energy loss due to radiative effects at high energies.
 ///
@@ -87,14 +84,14 @@ float computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab, float m,
 /// This computes the mean energy loss -dE(x) using an approximative formula.
 /// Bremsstrahlung is always included; direct e+e- pair production and
 /// photo-nuclear interactions only for muons.
-float computeEnergyLossRadiative(const MaterialSlab& slab, PdgParticle absPdg,
-                                 float m, float qOverP, float absQ);
+double computeEnergyLossRadiative(const MaterialSlab& slab, PdgParticle absPdg,
+                                  double m, double qOverP, double absQ);
 /// Derivative of the mean radiative energy loss with respect to q/p.
 ///
 /// @copydoc computeEnergyLossRadiative
-float deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
-                                      PdgParticle absPdg, float m, float qOverP,
-                                      float absQ);
+double deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
+                                       PdgParticle absPdg, double m,
+                                       double qOverP, double absQ);
 
 /// Compute the combined mean energy loss.
 ///
@@ -107,24 +104,24 @@ float deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
 /// This computes the combined mean energy loss -dE(x) including ionisation and
 /// radiative effects. The computations are valid over a wide range of particle
 /// energies.
-float computeEnergyLossMean(const MaterialSlab& slab, PdgParticle absPdg,
-                            float m, float qOverP, float absQ);
+double computeEnergyLossMean(const MaterialSlab& slab, PdgParticle absPdg,
+                             double m, double qOverP, double absQ);
 /// Derivative of the combined mean energy loss with respect to q/p.
 ///
 /// @copydoc computeEnergyLossMean
-float deriveEnergyLossMeanQOverP(const MaterialSlab& slab, PdgParticle absPdg,
-                                 float m, float qOverP, float absQ);
+double deriveEnergyLossMeanQOverP(const MaterialSlab& slab, PdgParticle absPdg,
+                                  double m, double qOverP, double absQ);
 
 /// Compute the combined most probably energy loss.
 ///
 /// @copydoc computeEnergyLossMean
-float computeEnergyLossMode(const MaterialSlab& slab, PdgParticle absPdg,
-                            float m, float qOverP, float absQ);
+double computeEnergyLossMode(const MaterialSlab& slab, PdgParticle absPdg,
+                             double m, double qOverP, double absQ);
 /// Derivative of the combined most probable energy loss with respect to q/p.
 ///
 /// @copydoc computeEnergyLossMean
-float deriveEnergyLossModeQOverP(const MaterialSlab& slab, PdgParticle absPdg,
-                                 float m, float qOverP, float absQ);
+double deriveEnergyLossModeQOverP(const MaterialSlab& slab, PdgParticle absPdg,
+                                  double m, double qOverP, double absQ);
 
 /// Compute the core width of the projected planar scattering distribution.
 ///
@@ -133,8 +130,8 @@ float deriveEnergyLossModeQOverP(const MaterialSlab& slab, PdgParticle absPdg,
 /// @param m         Particle mass
 /// @param qOverP    Particle charge divided by absolute momentum
 /// @param absQ      Absolute particle charge
-float computeMultipleScatteringTheta0(const MaterialSlab& slab,
-                                      PdgParticle absPdg, float m, float qOverP,
-                                      float absQ);
+double computeMultipleScatteringTheta0(const MaterialSlab& slab,
+                                       PdgParticle absPdg, double m,
+                                       double qOverP, double absQ);
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Material/Material.hpp b/Core/include/Acts/Material/Material.hpp
index 1caf59e97ec..b798d6a51f7 100644
--- a/Core/include/Acts/Material/Material.hpp
+++ b/Core/include/Acts/Material/Material.hpp
@@ -40,13 +40,35 @@ class Material {
  public:
   using ParametersVector = Eigen::Matrix<float, 5, 1>;
 
-  // Both mass and molar density are stored as a float and can thus not be
+  /// Return vacuum material parameters.
+  static ParametersVector vacuumParameters() {
+    return (ParametersVector() << std::numeric_limits<float>::infinity(),
+            std::numeric_limits<float>::infinity(), 0., 0., 0.)
+        .finished();
+  }
+
+  /// Return almost vacuum material parameters.
+  static ParametersVector almostVacuumParameters() {
+    return (ParametersVector() << std::numeric_limits<float>::max(),
+            std::numeric_limits<float>::max(), 0., 0., 1.)
+        .finished();
+  }
+
+  enum ParametersIndices {
+    eRadiationLength = 0,
+    eInteractionLength = 1,
+    eRelativeAtomicMass = 2,
+    eNuclearCharge = 3,
+    eMolarDensity = 4,
+  };
+
+  // Both mass and molar density are stored as a double and can thus not be
   // distinguished by their types. Just changing the last element in the
-  // previously existing constructor that took five floats as input to represent
-  // molar density instead of mass density could have lead to significant
-  // confusion compared to the previous behaviour. To avoid any ambiguity,
-  // construction from separate material parameters must happen through the
-  // following named constructors.
+  // previously existing constructor that took five doubles as input to
+  // represent molar density instead of mass density could have lead to
+  // significant confusion compared to the previous behaviour. To avoid any
+  // ambiguity, construction from separate material parameters must happen
+  // through the following named constructors.
 
   /// Construct from material parameters using the molar density.
   ///
@@ -55,8 +77,8 @@ class Material {
   /// @param ar is the relative atomic mass
   /// @param z is the nuclear charge number
   /// @param molarRho is the molar density
-  static Material fromMolarDensity(float x0, float l0, float ar, float z,
-                                   float molarRho);
+  static Material fromMolarDensity(double x0, double l0, double ar, double z,
+                                   double molarRho);
   /// Construct from material parameters using the mass density.
   ///
   /// @param x0 is the radiation length
@@ -68,8 +90,8 @@ class Material {
   /// @warning Due to the choice of native mass units, using the mass density
   ///   can lead to numerical problems. Typical mass densities lead to
   ///   computations with values differing by 20+ orders of magnitude.
-  static Material fromMassDensity(float x0, float l0, float ar, float z,
-                                  float massRho);
+  static Material fromMassDensity(double x0, double l0, double ar, double z,
+                                  double massRho);
   /// Construct a vacuum representation.
   Material() = default;
   /// Construct from an encoded parameters vector.
@@ -82,34 +104,37 @@ class Material {
   Material& operator=(const Material& mat) = default;
 
   /// Check if the material is valid, i.e. it is not vacuum.
-  bool isValid() const { return 0.0f < m_ar; }
+  bool isValid() const { return m_ar > 0; }
+
+  /// Check if the material is valid, i.e. it is not vacuum.
+  bool isVacuum() const { return m_ar > 0; }
 
   /// Return the radiation length. Infinity in case of vacuum.
-  constexpr float X0() const { return m_x0; }
+  constexpr double X0() const { return m_x0; }
   /// Return the nuclear interaction length. Infinity in case of vacuum.
-  constexpr float L0() const { return m_l0; }
+  constexpr double L0() const { return m_l0; }
   /// Return the relative atomic mass.
-  constexpr float Ar() const { return m_ar; }
+  constexpr double Ar() const { return m_ar; }
   /// Return the nuclear charge number.
-  constexpr float Z() const { return m_z; }
+  constexpr double Z() const { return m_z; }
   /// Return the molar density.
-  constexpr float molarDensity() const { return m_molarRho; }
+  constexpr double molarDensity() const { return m_molarRho; }
   /// Return the molar electron density.
-  constexpr float molarElectronDensity() const { return m_z * m_molarRho; }
+  constexpr double molarElectronDensity() const { return m_z * m_molarRho; }
   /// Return the mass density.
-  float massDensity() const;
+  double massDensity() const;
   /// Return the mean electron excitation energy.
-  float meanExcitationEnergy() const;
+  double meanExcitationEnergy() const;
 
   /// Encode the properties into an opaque parameters vector.
   ParametersVector parameters() const;
 
  private:
-  float m_x0 = std::numeric_limits<float>::infinity();
-  float m_l0 = std::numeric_limits<float>::infinity();
-  float m_ar = 0.0f;
-  float m_z = 0.0f;
-  float m_molarRho = 0.0f;
+  double m_x0 = std::numeric_limits<double>::infinity();
+  double m_l0 = std::numeric_limits<double>::infinity();
+  double m_ar = 0;
+  double m_z = 0;
+  double m_molarRho = 0;
 
   /// @brief Check if two materials are exactly equal.
   ///
diff --git a/Core/include/Acts/Material/MaterialComposition.hpp b/Core/include/Acts/Material/MaterialComposition.hpp
index 843736ebd58..25dd885a1e2 100644
--- a/Core/include/Acts/Material/MaterialComposition.hpp
+++ b/Core/include/Acts/Material/MaterialComposition.hpp
@@ -37,7 +37,7 @@ class ElementFraction {
   ///
   /// @param e is the atomic number of the element
   /// @param f is the relative fraction and must be a value in [0,1]
-  constexpr ElementFraction(unsigned int e, float f)
+  constexpr ElementFraction(unsigned int e, double f)
       : m_element(static_cast<std::uint8_t>(e)),
         m_fraction(static_cast<std::uint8_t>(
             f * std::numeric_limits<std::uint8_t>::max())) {
@@ -55,19 +55,11 @@ class ElementFraction {
     assert((w < 256u) && "Integer weight must be in [0,256)");
   }
 
-  /// Must always be created with valid data.
-  ElementFraction() = delete;
-  ElementFraction(ElementFraction&&) = default;
-  ElementFraction(const ElementFraction&) = default;
-  ~ElementFraction() = default;
-  ElementFraction& operator=(ElementFraction&&) = default;
-  ElementFraction& operator=(const ElementFraction&) = default;
-
   /// The element atomic number.
   constexpr std::uint8_t element() const { return m_element; }
   /// The relative fraction of this element.
-  constexpr float fraction() const {
-    return static_cast<float>(m_fraction) /
+  constexpr double fraction() const {
+    return static_cast<double>(m_fraction) /
            std::numeric_limits<std::uint8_t>::max();
   }
 
@@ -107,19 +99,13 @@ class MaterialComposition {
       total += element.m_fraction;
     }
     // compute scale factor into the [0, 256) range
-    float scale = float{std::numeric_limits<std::uint8_t>::max()} / total;
+    double scale = double{std::numeric_limits<std::uint8_t>::max()} / total;
     for (auto& element : m_elements) {
       element.m_fraction =
           static_cast<std::uint8_t>(element.m_fraction * scale);
     }
   }
 
-  MaterialComposition(MaterialComposition&&) = default;
-  MaterialComposition(const MaterialComposition&) = default;
-  ~MaterialComposition() = default;
-  MaterialComposition& operator=(MaterialComposition&&) = default;
-  MaterialComposition& operator=(const MaterialComposition&) = default;
-
   // Support range-based iteration over contained elements.
   auto begin() const { return m_elements.begin(); }
   auto end() const { return m_elements.end(); }
diff --git a/Core/include/Acts/Material/MaterialSlab.hpp b/Core/include/Acts/Material/MaterialSlab.hpp
index 63fd7378e8b..660b6cceaa3 100644
--- a/Core/include/Acts/Material/MaterialSlab.hpp
+++ b/Core/include/Acts/Material/MaterialSlab.hpp
@@ -45,12 +45,12 @@ class MaterialSlab {
   /// Construct vacuum without thickness.
   MaterialSlab() = default;
   /// Construct vacuum with thickness.
-  explicit MaterialSlab(float thickness);
+  explicit MaterialSlab(double thickness);
   /// Construct from material description.
   ///
   /// @param material  is the material description
   /// @param thickness is the thickness of the material
-  MaterialSlab(const Material& material, float thickness);
+  MaterialSlab(const Material& material, double thickness);
   ~MaterialSlab() = default;
 
   MaterialSlab(MaterialSlab&&) = default;
@@ -59,25 +59,25 @@ class MaterialSlab {
   MaterialSlab& operator=(const MaterialSlab&) = default;
 
   /// Scale the material thickness by the given factor.
-  void scaleThickness(float scale);
+  void scaleThickness(double scale);
 
   /// Check if the material is valid, i.e. it is finite and not vacuum.
-  bool isValid() const { return m_material.isValid() && (0.0f < m_thickness); }
+  bool isValid() const { return m_material.isValid() && (m_thickness > 0); }
 
   /// Access the (average) material parameters.
   constexpr const Material& material() const { return m_material; }
   /// Return the thickness.
-  constexpr float thickness() const { return m_thickness; }
+  constexpr double thickness() const { return m_thickness; }
   /// Return the radiation length fraction.
-  constexpr float thicknessInX0() const { return m_thicknessInX0; }
+  constexpr double thicknessInX0() const { return m_thicknessInX0; }
   /// Return the nuclear interaction length fraction.
-  constexpr float thicknessInL0() const { return m_thicknessInL0; }
+  constexpr double thicknessInL0() const { return m_thicknessInL0; }
 
  private:
   Material m_material;
-  float m_thickness = 0.0f;
-  float m_thicknessInX0 = 0.0f;
-  float m_thicknessInL0 = 0.0f;
+  double m_thickness = 0;
+  double m_thicknessInX0 = 0;
+  double m_thicknessInL0 = 0;
 
   /// @brief Check if two materials are exactly equal.
   ///
diff --git a/Core/include/Acts/Material/VolumeMaterialMapper.hpp b/Core/include/Acts/Material/VolumeMaterialMapper.hpp
index baf29d7b3d1..9ec33b56a64 100644
--- a/Core/include/Acts/Material/VolumeMaterialMapper.hpp
+++ b/Core/include/Acts/Material/VolumeMaterialMapper.hpp
@@ -20,7 +20,6 @@
 #include "Acts/Material/MaterialGridHelper.hpp"
 #include "Acts/Material/MaterialInteraction.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
-#include "Acts/Propagator/MaterialInteractor.hpp"
 #include "Acts/Propagator/Navigator.hpp"
 #include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Propagator/StraightLineStepper.hpp"
@@ -32,7 +31,6 @@
 #include <map>
 #include <memory>
 #include <utility>
-#include <vector>
 
 namespace Acts {
 
@@ -73,7 +71,7 @@ class VolumeMaterialMapper {
   /// Nested Configuration struct for the material mapper
   struct Config {
     /// Size of the step for the step extrapolation
-    float mappingStep = 1.;
+    double mappingStep = 1;
   };
 
   /// @struct State
diff --git a/Core/include/Acts/Utilities/Interpolation.hpp b/Core/include/Acts/Utilities/Interpolation.hpp
index 27b71e10358..4e9ce20e0c5 100644
--- a/Core/include/Acts/Utilities/Interpolation.hpp
+++ b/Core/include/Acts/Utilities/Interpolation.hpp
@@ -8,11 +8,9 @@
 
 #pragma once
 
-#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Utilities/detail/interpolation_impl.hpp"
 
 #include <array>
-#include <type_traits>
 
 namespace Acts {
 
diff --git a/Core/src/Material/AccumulatedMaterialSlab.cpp b/Core/src/Material/AccumulatedMaterialSlab.cpp
index b8b7ebf88fe..a7db07dff84 100644
--- a/Core/src/Material/AccumulatedMaterialSlab.cpp
+++ b/Core/src/Material/AccumulatedMaterialSlab.cpp
@@ -11,22 +11,24 @@
 #include "Acts/Material/Material.hpp"
 #include "Acts/Material/detail/AverageMaterials.hpp"
 
-void Acts::AccumulatedMaterialSlab::accumulate(MaterialSlab slab,
-                                               float pathCorrection) {
+namespace Acts {
+
+void AccumulatedMaterialSlab::accumulate(MaterialSlab slab,
+                                         double pathCorrection) {
   // scale the recorded material to the equivalence contribution along the
   // surface normal
   slab.scaleThickness(1 / pathCorrection);
   m_trackAverage = detail::combineSlabs(m_trackAverage, slab);
 }
 
-void Acts::AccumulatedMaterialSlab::trackVariance(MaterialSlab slabReference,
-                                                  bool useEmptyTrack) {
+void AccumulatedMaterialSlab::trackVariance(MaterialSlab slabReference,
+                                            bool useEmptyTrack) {
   // Only use real tracks or if empty tracks are allowed.
   if (useEmptyTrack || (0 < m_trackAverage.thickness())) {
-    float variance = ((1 / m_trackAverage.material().X0()) -
-                      (1 / slabReference.material().X0())) *
-                     ((1 / m_trackAverage.material().X0()) -
-                      (1 / slabReference.material().X0()));
+    double variance = ((1 / m_trackAverage.material().X0()) -
+                       (1 / slabReference.material().X0())) *
+                      ((1 / m_trackAverage.material().X0()) -
+                       (1 / slabReference.material().X0()));
     if (m_totalCount == 0u) {
       m_totalVariance = variance;
     } else {
@@ -37,7 +39,7 @@ void Acts::AccumulatedMaterialSlab::trackVariance(MaterialSlab slabReference,
   }
 }
 
-void Acts::AccumulatedMaterialSlab::trackAverage(bool useEmptyTrack) {
+void AccumulatedMaterialSlab::trackAverage(bool useEmptyTrack) {
   // average only real tracks or if empty tracks are allowed.
   if (useEmptyTrack || (0 < m_trackAverage.thickness())) {
     if (m_totalCount == 0u) {
@@ -58,12 +60,13 @@ void Acts::AccumulatedMaterialSlab::trackAverage(bool useEmptyTrack) {
   m_trackAverage = MaterialSlab();
 }
 
-std::pair<Acts::MaterialSlab, unsigned int>
-Acts::AccumulatedMaterialSlab::totalAverage() const {
+std::pair<MaterialSlab, unsigned int> AccumulatedMaterialSlab::totalAverage()
+    const {
   return {m_totalAverage, m_totalCount};
 }
 
-std::pair<float, unsigned int> Acts::AccumulatedMaterialSlab::totalVariance()
-    const {
+std::pair<double, unsigned int> AccumulatedMaterialSlab::totalVariance() const {
   return {m_totalVariance, m_totalCount};
 }
+
+}  // namespace Acts
diff --git a/Core/src/Material/AverageMaterials.cpp b/Core/src/Material/AverageMaterials.cpp
index e9770e910bf..13e2831c3bb 100644
--- a/Core/src/Material/AverageMaterials.cpp
+++ b/Core/src/Material/AverageMaterials.cpp
@@ -12,8 +12,17 @@
 
 #include <cmath>
 
-Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
-                                              const MaterialSlab& slab2) {
+namespace Acts {
+
+MaterialSlab detail::combineSlabs(const MaterialSlab& slab1,
+                                  const MaterialSlab& slab2) {
+  if (slab1.thickness() <= 0) {
+    return slab2;
+  }
+  if (slab2.thickness() <= 0) {
+    return slab1;
+  }
+
   const auto& mat1 = slab1.material();
   const auto& mat2 = slab2.material();
 
@@ -28,43 +37,36 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
   // In the case of the atomic number, the averaging is based on the log
   // to properly account for the energy loss in multiple material.
 
-  // use double for (intermediate) computations to avoid precision loss
-  const double thickness1 = static_cast<double>(slab1.thickness());
-  const double thickness2 = static_cast<double>(slab2.thickness());
-
   // the thickness properties are purely additive
-  const double thickness = thickness1 + thickness2;
+  const double thickness = slab1.thickness() + slab2.thickness();
 
   // if the two materials are the same there is no need for additional
   // computation
   if (mat1 == mat2) {
-    return {mat1, static_cast<float>(thickness)};
+    return {mat1, thickness};
   }
 
   // molar amount-of-substance assuming a unit area, i.e. volume = thickness*1*1
-  const double molarDensity1 = static_cast<double>(mat1.molarDensity());
-  const double molarDensity2 = static_cast<double>(mat2.molarDensity());
-
-  const double molarAmount1 = molarDensity1 * thickness1;
-  const double molarAmount2 = molarDensity2 * thickness2;
+  const double molarAmount1 = mat1.molarDensity() * slab1.thickness();
+  const double molarAmount2 = mat2.molarDensity() * slab2.thickness();
   const double molarAmount = molarAmount1 + molarAmount2;
 
   // handle vacuum specially
-  if (!(0.0 < molarAmount)) {
-    return {Material(), static_cast<float>(thickness)};
+  if (molarAmount <= 0) {
+    return {Material(), thickness};
   }
 
   // radiation/interaction length follows from consistency argument
-  const double thicknessX01 = static_cast<double>(slab1.thicknessInX0());
-  const double thicknessX02 = static_cast<double>(slab2.thicknessInX0());
-  const double thicknessL01 = static_cast<double>(slab1.thicknessInL0());
-  const double thicknessL02 = static_cast<double>(slab2.thicknessInL0());
+  const double thicknessX01 = slab1.thicknessInX0();
+  const double thicknessX02 = slab2.thicknessInX0();
+  const double thicknessL01 = slab1.thicknessInL0();
+  const double thicknessL02 = slab2.thicknessInL0();
 
   const double thicknessInX0 = thicknessX01 + thicknessX02;
   const double thicknessInL0 = thicknessL01 + thicknessL02;
 
-  const float x0 = static_cast<float>(thickness / thicknessInX0);
-  const float l0 = static_cast<float>(thickness / thicknessInL0);
+  const double x0 = thickness / thicknessInX0;
+  const double l0 = thickness / thicknessInL0;
 
   // assume two slabs of materials with N1,N2 atoms/molecules each with atomic
   // masses A1,A2 and nuclear charges. We have a total of N = N1 + N2
@@ -83,12 +85,9 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
   //        = (Vi*rhoi) / (V1*rho1 + V2*rho2)
   //
   // which can be computed from the molar amount-of-substance above.
-  const double ar1 = static_cast<double>(mat1.Ar());
-  const double ar2 = static_cast<double>(mat2.Ar());
-
   const double molarWeight1 = molarAmount1 / molarAmount;
   const double molarWeight2 = molarAmount2 / molarAmount;
-  const float ar = static_cast<float>(molarWeight1 * ar1 + molarWeight2 * ar2);
+  const double ar = molarWeight1 * mat1.Ar() + molarWeight2 * mat2.Ar();
 
   // In the case of the atomic number, its main use is the computation
   // of the mean excitation energy approximated in ATL-SOFT-PUB-2008-003 as :
@@ -103,23 +102,21 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
   // To respect this the average atomic number thus need to be defined as :
   //     ln(Z)*t = ln(Z1)*t1 + ln(Z2)*t2
   //           Z = Exp( ln(Z1)*t1/t + ln(Z2)*t2/t )
-  const double z1 = static_cast<double>(mat1.Z());
-  const double z2 = static_cast<double>(mat2.Z());
-
-  const double thicknessWeight1 = thickness1 / thickness;
-  const double thicknessWeight2 = thickness2 / thickness;
-  float z = 0.f;
-  if (z1 > 0. && z2 > 0.) {
-    z = static_cast<float>(std::exp(thicknessWeight1 * std::log(z1) +
-                                    thicknessWeight2 * std::log(z2)));
+  const double thicknessWeight1 = slab1.thickness() / thickness;
+  const double thicknessWeight2 = slab2.thickness() / thickness;
+  double z = 0;
+  if (mat1.Z() > 0 && mat2.Z() > 0) {
+    z = std::exp(thicknessWeight1 * std::log(mat1.Z()) +
+                 thicknessWeight2 * std::log(mat2.Z()));
   } else {
-    z = static_cast<float>(thicknessWeight1 * z1 + thicknessWeight2 * z2);
+    z = thicknessWeight1 * mat1.Z() + thicknessWeight2 * mat2.Z();
   }
 
   // compute average molar density by dividing the total amount-of-substance by
   // the total volume for the same unit area, i.e. volume = totalThickness*1*1
-  const float molarDensity = static_cast<float>(molarAmount / thickness);
+  const double molarDensity = molarAmount / thickness;
 
-  return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity),
-          static_cast<float>(thickness)};
+  return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity), thickness};
 }
+
+}  // namespace Acts
diff --git a/Core/src/Material/Interactions.cpp b/Core/src/Material/Interactions.cpp
index a0612e4d0f8..0b217049986 100644
--- a/Core/src/Material/Interactions.cpp
+++ b/Core/src/Material/Interactions.cpp
@@ -9,6 +9,7 @@
 #include "Acts/Material/Interactions.hpp"
 
 #include "Acts/Definitions/PdgParticle.hpp"
+#include "Acts/Definitions/Units.hpp"
 #include "Acts/Material/Material.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 
@@ -21,20 +22,20 @@ namespace {
 
 // values from RPP2018 table 33.1
 // electron mass
-constexpr float Me = 0.5109989461_MeV;
+constexpr double Me = 0.5109989461_MeV;
 // Bethe formular prefactor. 1/mol unit is just a factor 1 here.
-constexpr float K = 0.307075_MeV * 1_cm * 1_cm;
+constexpr double K = 0.307075_MeV * 1_cm * 1_cm;
 // Energy scale for plasma energy.
-constexpr float PlasmaEnergyScale = 28.816_eV;
+constexpr double PlasmaEnergyScale = 28.816_eV;
 
 /// Additional derived relativistic quantities.
 struct RelativisticQuantities {
-  float q2OverBeta2 = 0.0f;
-  float beta2 = 0.0f;
-  float betaGamma = 0.0f;
-  float gamma = 0.0f;
+  double q2OverBeta2 = 0;
+  double beta2 = 0;
+  double betaGamma = 0;
+  double gamma = 0;
 
-  RelativisticQuantities(float mass, float qOverP, float absQ) {
+  RelativisticQuantities(double mass, double qOverP, double absQ) {
     assert((0 < mass) && "Mass must be positive");
     assert((qOverP != 0) && "q/p must be non-zero");
     assert((absQ > 0) && "Absolute charge must be non-zero and positive");
@@ -43,10 +44,10 @@ struct RelativisticQuantities {
     q2OverBeta2 = absQ * absQ + (mass * qOverP) * (mass * qOverP);
     assert((q2OverBeta2 >= 0) && "Negative q2OverBeta2");
     // 1/p = q/(qp) = (q/p)/q
-    const float mOverP = mass * std::abs(qOverP / absQ);
-    const float pOverM = 1.0f / mOverP;
+    const double mOverP = mass * std::abs(qOverP / absQ);
+    const double pOverM = 1 / mOverP;
     // beta² = p²/E² = p²/(m² + p²) = 1/(1 + (m/p)²)
-    beta2 = 1.0f / (1.0f + mOverP * mOverP);
+    beta2 = 1 / (1 + mOverP * mOverP);
     assert((beta2 >= 0) && "Negative beta2");
     // beta*gamma = (p/sqrt(m² + p²))*(sqrt(m² + p²)/m) = p/m
     betaGamma = pOverM;
@@ -57,17 +58,17 @@ struct RelativisticQuantities {
 };
 
 /// Compute q/p derivative of beta².
-inline float deriveBeta2(float qOverP, const RelativisticQuantities& rq) {
+inline double deriveBeta2(double qOverP, const RelativisticQuantities& rq) {
   return -2 / (qOverP * rq.gamma * rq.gamma);
 }
 
 /// Compute the 2 * mass * (beta * gamma)² mass term.
-inline float computeMassTerm(float mass, const RelativisticQuantities& rq) {
+inline double computeMassTerm(double mass, const RelativisticQuantities& rq) {
   return 2 * mass * rq.betaGamma * rq.betaGamma;
 }
 
 /// Compute mass term logarithmic derivative w/ respect to q/p.
-inline float logDeriveMassTerm(float qOverP) {
+inline double logDeriveMassTerm(double qOverP) {
   // only need to compute d((beta*gamma)²)/(beta*gamma)²; rest cancels.
   return -2 / qOverP;
 }
@@ -75,22 +76,22 @@ inline float logDeriveMassTerm(float qOverP) {
 /// Compute the maximum energy transfer in a single collision.
 ///
 /// Uses RPP2018 eq. 33.4.
-inline float computeWMax(float mass, const RelativisticQuantities& rq) {
-  const float mfrac = Me / mass;
-  const float nominator = 2 * Me * rq.betaGamma * rq.betaGamma;
-  const float denonimator = 1.0f + 2 * rq.gamma * mfrac + mfrac * mfrac;
+inline double computeWMax(double mass, const RelativisticQuantities& rq) {
+  const double mfrac = Me / mass;
+  const double nominator = 2 * Me * rq.betaGamma * rq.betaGamma;
+  const double denonimator = 1 + 2 * rq.gamma * mfrac + mfrac * mfrac;
   return nominator / denonimator;
 }
 
 /// Compute WMax logarithmic derivative w/ respect to q/p.
-inline float logDeriveWMax(float mass, float qOverP,
-                           const RelativisticQuantities& rq) {
+inline double logDeriveWMax(double mass, double qOverP,
+                            const RelativisticQuantities& rq) {
   // this is (q/p) * (beta/q).
   // both quantities have the same sign and the product must always be
   // positive. we can thus reuse the known (unsigned) quantity (q/beta)².
-  const float a = std::abs(qOverP / std::sqrt(rq.q2OverBeta2));
+  const double a = std::abs(qOverP / std::sqrt(rq.q2OverBeta2));
   // (m² + me²) / me = me (1 + (m/me)²)
-  const float b = Me * (1.0f + (mass / Me) * (mass / Me));
+  const double b = Me * (1 + (mass / Me) * (mass / Me));
   return -2 * (a * b - 2 + rq.beta2) / (qOverP * (a * b + 2));
 }
 
@@ -102,40 +103,41 @@ inline float logDeriveWMax(float mass, float qOverP,
 ///
 /// where (Z/A)*rho is the electron density in the material and x is the
 /// traversed length (thickness) of the material.
-inline float computeEpsilon(float molarElectronDensity, float thickness,
-                            const RelativisticQuantities& rq) {
-  return 0.5f * K * molarElectronDensity * thickness * rq.q2OverBeta2;
+inline double computeEpsilon(double molarElectronDensity, double thickness,
+                             const RelativisticQuantities& rq) {
+  return 0.5 * K * molarElectronDensity * thickness * rq.q2OverBeta2;
 }
 
 /// Compute epsilon logarithmic derivative w/ respect to q/p.
-inline float logDeriveEpsilon(float qOverP, const RelativisticQuantities& rq) {
+inline double logDeriveEpsilon(double qOverP,
+                               const RelativisticQuantities& rq) {
   // only need to compute d(q²/beta²)/(q²/beta²); everything else cancels.
   return 2 / (qOverP * rq.gamma * rq.gamma);
 }
 
 /// Compute the density correction factor delta/2.
-inline float computeDeltaHalf(float meanExitationPotential,
-                              float molarElectronDensity,
-                              const RelativisticQuantities& rq) {
+inline double computeDeltaHalf(double meanExitationPotential,
+                               double molarElectronDensity,
+                               const RelativisticQuantities& rq) {
   /// Uses RPP2018 eq. 33.6 which is only valid for high energies.
   // only relevant for very high ernergies; use arbitrary cutoff
-  if (rq.betaGamma < 10.0f) {
-    return 0.0f;
+  if (rq.betaGamma < 10.0) {
+    return 0;
   }
   // pre-factor according to RPP2019 table 33.1
-  const float plasmaEnergy =
-      PlasmaEnergyScale * std::sqrt(1000.f * molarElectronDensity);
+  const double plasmaEnergy =
+      PlasmaEnergyScale * std::sqrt(1000 * molarElectronDensity);
   return std::log(rq.betaGamma) +
-         std::log(plasmaEnergy / meanExitationPotential) - 0.5f;
+         std::log(plasmaEnergy / meanExitationPotential) - 0.5;
 }
 
 /// Compute derivative w/ respect to q/p for the density correction.
-inline float deriveDeltaHalf(float qOverP, const RelativisticQuantities& rq) {
+inline double deriveDeltaHalf(double qOverP, const RelativisticQuantities& rq) {
   // original equation is of the form
   //     log(beta*gamma) + log(eplasma/I) - 1/2
   // which the resulting derivative as
   //     d(beta*gamma)/(beta*gamma)
-  return (rq.betaGamma < 10.0f) ? 0.0f : (-1.0f / qOverP);
+  return (rq.betaGamma < 10.0) ? 0 : (-1 / qOverP);
 }
 
 /// Convert Landau full-width-half-maximum to an equivalent Gaussian sigma,
@@ -145,22 +147,22 @@ inline float deriveDeltaHalf(float qOverP, const RelativisticQuantities& rq) {
 ///     fwhm = 2 * sqrt(2 * log(2)) * sigma
 /// -> sigma = fwhm / (2 * sqrt(2 * log(2)))
 ///
-inline float convertLandauFwhmToGaussianSigma(float fwhm) {
-  // return fwhm / (2 * std::sqrt(2 * std::log(2.0f)));
-  return fwhm * 0.42466090014400953f;
+inline double convertLandauFwhmToGaussianSigma(double fwhm) {
+  // return fwhm / (2 * std::sqrt(2 * std::log(2.0)));
+  return fwhm * 0.42466090014400953;
 }
 
 namespace detail {
 
-inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
-                                         const RelativisticQuantities& rq) {
+inline double computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
+                                          const RelativisticQuantities& rq) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
   // the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
   return 4 * computeEpsilon(Ne, thickness, rq);
 }
@@ -169,46 +171,46 @@ inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
 
 }  // namespace
 
-float Acts::computeEnergyLossBethe(const MaterialSlab& slab, float m,
-                                   float qOverP, float absQ) {
+double Acts::computeEnergyLossBethe(const MaterialSlab& slab, double m,
+                                    double qOverP, double absQ) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float I = slab.material().meanExcitationEnergy();
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
-  const float eps = computeEpsilon(Ne, thickness, rq);
-  const float dhalf = computeDeltaHalf(I, Ne, rq);
-  const float u = computeMassTerm(Me, rq);
-  const float wmax = computeWMax(m, rq);
+  const double I = slab.material().meanExcitationEnergy();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
+  const double eps = computeEpsilon(Ne, thickness, rq);
+  const double dhalf = computeDeltaHalf(I, Ne, rq);
+  const double u = computeMassTerm(Me, rq);
+  const double wmax = computeWMax(m, rq);
   // uses RPP2018 eq. 33.5 scaled from mass stopping power to linear stopping
   // power and multiplied with the material thickness to get a total energy loss
   // instead of an energy loss per length.
   // the required modification only change the prefactor which becomes
   // identical to the prefactor epsilon for the most probable value.
-  const float running =
-      std::log(u / I) + std::log(wmax / I) - 2.0f * rq.beta2 - 2.0f * dhalf;
+  const double running =
+      std::log(u / I) + std::log(wmax / I) - 2 * rq.beta2 - 2 * dhalf;
   return eps * running;
 }
 
-float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
-                                        float qOverP, float absQ) {
+double Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, double m,
+                                         double qOverP, double absQ) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float I = slab.material().meanExcitationEnergy();
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
-  const float eps = computeEpsilon(Ne, thickness, rq);
-  const float dhalf = computeDeltaHalf(I, Ne, rq);
-  const float u = computeMassTerm(Me, rq);
-  const float wmax = computeWMax(m, rq);
+  const double I = slab.material().meanExcitationEnergy();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
+  const double eps = computeEpsilon(Ne, thickness, rq);
+  const double dhalf = computeDeltaHalf(I, Ne, rq);
+  const double u = computeMassTerm(Me, rq);
+  const double wmax = computeWMax(m, rq);
   // original equation is of the form
   //
   //     eps * (log(u/I) + log(wmax/I) - 2 * beta² - delta)
@@ -220,51 +222,51 @@ float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
   //     + eps * (d(u)/(u) + d(wmax)/(wmax) - 2 * d(beta²) - d(delta))
   //
   // where we can use d(eps) = eps * (d(eps)/eps) for further simplification.
-  const float logDerEps = logDeriveEpsilon(qOverP, rq);
-  const float derDHalf = deriveDeltaHalf(qOverP, rq);
-  const float logDerU = logDeriveMassTerm(qOverP);
-  const float logDerWmax = logDeriveWMax(m, qOverP, rq);
-  const float derBeta2 = deriveBeta2(qOverP, rq);
-  const float rel = logDerEps * (std::log(u / I) + std::log(wmax / I) -
-                                 2.0f * rq.beta2 - 2.0f * dhalf) +
-                    logDerU + logDerWmax - 2.0f * derBeta2 - 2.0f * derDHalf;
+  const double logDerEps = logDeriveEpsilon(qOverP, rq);
+  const double derDHalf = deriveDeltaHalf(qOverP, rq);
+  const double logDerU = logDeriveMassTerm(qOverP);
+  const double logDerWmax = logDeriveWMax(m, qOverP, rq);
+  const double derBeta2 = deriveBeta2(qOverP, rq);
+  const double rel = logDerEps * (std::log(u / I) + std::log(wmax / I) -
+                                  2 * rq.beta2 - 2 * dhalf) +
+                     logDerU + logDerWmax - 2 * derBeta2 - 2 * derDHalf;
   return eps * rel;
 }
 
-float Acts::computeEnergyLossLandau(const MaterialSlab& slab, float m,
-                                    float qOverP, float absQ) {
+double Acts::computeEnergyLossLandau(const MaterialSlab& slab, double m,
+                                     double qOverP, double absQ) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float I = slab.material().meanExcitationEnergy();
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
-  const float eps = computeEpsilon(Ne, thickness, rq);
-  const float dhalf = computeDeltaHalf(I, Ne, rq);
-  const float u = computeMassTerm(Me, rq);
+  const double I = slab.material().meanExcitationEnergy();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
+  const double eps = computeEpsilon(Ne, thickness, rq);
+  const double dhalf = computeDeltaHalf(I, Ne, rq);
+  const double u = computeMassTerm(Me, rq);
   // uses RPP2018 eq. 33.12
-  const float running =
-      std::log(u / I) + std::log(eps / I) + 0.2f - rq.beta2 - 2 * dhalf;
+  const double running =
+      std::log(u / I) + std::log(eps / I) + 0.2 - rq.beta2 - 2 * dhalf;
   return eps * running;
 }
 
-float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
-                                         float qOverP, float absQ) {
+double Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, double m,
+                                          double qOverP, double absQ) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float I = slab.material().meanExcitationEnergy();
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
-  const float eps = computeEpsilon(Ne, thickness, rq);
-  const float dhalf = computeDeltaHalf(I, Ne, rq);
-  const float t = computeMassTerm(Me, rq);
+  const double I = slab.material().meanExcitationEnergy();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
+  const double eps = computeEpsilon(Ne, thickness, rq);
+  const double dhalf = computeDeltaHalf(I, Ne, rq);
+  const double t = computeMassTerm(Me, rq);
   // original equation is of the form
   //
   //     eps * (log(t/I) - log(eps/I) - 0.2 - beta² - delta)
@@ -276,43 +278,43 @@ float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
   //     + eps * (d(t)/t + d(eps)/eps - d(beta²) - 2*d(delta/2))
   //
   // where we can use d(eps) = eps * (d(eps)/eps) for further simplification.
-  const float logDerEps = logDeriveEpsilon(qOverP, rq);
-  const float derDHalf = deriveDeltaHalf(qOverP, rq);
-  const float logDerT = logDeriveMassTerm(qOverP);
-  const float derBeta2 = deriveBeta2(qOverP, rq);
-  const float rel = logDerEps * (std::log(t / I) + std::log(eps / I) - 0.2f -
-                                 rq.beta2 - 2 * dhalf) +
-                    logDerT + logDerEps - derBeta2 - 2 * derDHalf;
+  const double logDerEps = logDeriveEpsilon(qOverP, rq);
+  const double derDHalf = deriveDeltaHalf(qOverP, rq);
+  const double logDerT = logDeriveMassTerm(qOverP);
+  const double derBeta2 = deriveBeta2(qOverP, rq);
+  const double rel = logDerEps * (std::log(t / I) + std::log(eps / I) - 0.2 -
+                                  rq.beta2 - 2 * dhalf) +
+                     logDerT + logDerEps - derBeta2 - 2 * derDHalf;
   return eps * rel;
 }
 
-float Acts::computeEnergyLossLandauSigma(const MaterialSlab& slab, float m,
-                                         float qOverP, float absQ) {
+double Acts::computeEnergyLossLandauSigma(const MaterialSlab& slab, double m,
+                                          double qOverP, double absQ) {
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float Ne = slab.material().molarElectronDensity();
-  const float thickness = slab.thickness();
+  const double Ne = slab.material().molarElectronDensity();
+  const double thickness = slab.thickness();
   // the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
-  const float fwhm = 4 * computeEpsilon(Ne, thickness, rq);
+  const double fwhm = 4 * computeEpsilon(Ne, thickness, rq);
   return convertLandauFwhmToGaussianSigma(fwhm);
 }
 
-float Acts::computeEnergyLossLandauFwhm(const MaterialSlab& slab, float m,
-                                        float qOverP, float absQ) {
+double Acts::computeEnergyLossLandauFwhm(const MaterialSlab& slab, double m,
+                                         double qOverP, double absQ) {
   const RelativisticQuantities rq{m, qOverP, absQ};
   return detail::computeEnergyLossLandauFwhm(slab, rq);
 }
 
-float Acts::computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab,
-                                               float m, float qOverP,
-                                               float absQ) {
+double Acts::computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab,
+                                                double m, double qOverP,
+                                                double absQ) {
   const RelativisticQuantities rq{m, qOverP, absQ};
-  const float fwhm = detail::computeEnergyLossLandauFwhm(slab, rq);
-  const float sigmaE = convertLandauFwhmToGaussianSigma(fwhm);
+  const double fwhm = detail::computeEnergyLossLandauFwhm(slab, rq);
+  const double sigmaE = convertLandauFwhmToGaussianSigma(fwhm);
   //  var(q/p) = (d(q/p)/dE)² * var(E)
   // d(q/p)/dE = d/dE (q/sqrt(E²-m²))
   //           = q * -(1/2) * 1/p³ * 2E
@@ -320,20 +322,20 @@ float Acts::computeEnergyLossLandauSigmaQOverP(const MaterialSlab& slab,
   //  var(q/p) = (q/p)^4 * (q/beta)² * (1/q)^4 * var(E)
   //           = (1/p)^4 * (q/beta)² * var(E)
   // do not need to care about the sign since it is only used squared
-  const float pInv = qOverP / absQ;
-  const float qOverBeta = std::sqrt(rq.q2OverBeta2);
+  const double pInv = qOverP / absQ;
+  const double qOverBeta = std::sqrt(rq.q2OverBeta2);
   return qOverBeta * pInv * pInv * sigmaE;
 }
 
 namespace {
 
 /// Compute mean energy loss from bremsstrahlung per radiation length.
-inline float computeBremsstrahlungLossMean(float mass, float energy) {
+inline double computeBremsstrahlungLossMean(double mass, double energy) {
   return energy * (Me / mass) * (Me / mass);
 }
 
 /// Derivative of the bremsstrahlung loss per rad length with respect to energy.
-inline float deriveBremsstrahlungLossMeanE(float mass) {
+inline double deriveBremsstrahlungLossMeanE(double mass) {
   return (Me / mass) * (Me / mass);
 }
 
@@ -345,7 +347,7 @@ inline float deriveBremsstrahlungLossMeanE(float mass) {
 /// factored out and the coefficients must be scaled to the native units such
 /// that the evaluated expansion with terms E^n has dimension energy in
 /// native units.
-constexpr float MuonHighLowThreshold = 1_TeV;
+constexpr double MuonHighLowThreshold = 1_TeV;
 // [low0 / X0] = MeV / mm -> [low0] = MeV
 constexpr double MuonLow0 = -0.5345_MeV;
 // [low1 * E / X0] = MeV / mm -> [low1] = 1
@@ -360,7 +362,7 @@ constexpr double MuonHigh0 = -2.986_MeV;
 constexpr double MuonHigh1 = 9.253e-5;
 
 /// Compute additional radiation energy loss for muons per radiation length.
-inline float computeMuonDirectPairPhotoNuclearLossMean(double energy) {
+inline double computeMuonDirectPairPhotoNuclearLossMean(double energy) {
   if (energy < MuonHighLowThreshold) {
     return MuonLow0 +
            (MuonLow1 + (MuonLow2 + MuonLow3 * energy) * energy) * energy;
@@ -370,7 +372,7 @@ inline float computeMuonDirectPairPhotoNuclearLossMean(double energy) {
 }
 
 /// Derivative of the additional rad loss per rad length with respect to energy.
-inline float deriveMuonDirectPairPhotoNuclearLossMeanE(double energy) {
+inline double deriveMuonDirectPairPhotoNuclearLossMeanE(double energy) {
   if (energy < MuonHighLowThreshold) {
     return MuonLow1 + (2 * MuonLow2 + 3 * MuonLow3 * energy) * energy;
   } else {
@@ -380,25 +382,25 @@ inline float deriveMuonDirectPairPhotoNuclearLossMeanE(double energy) {
 
 }  // namespace
 
-float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
-                                       PdgParticle absPdg, float m,
-                                       float qOverP, float absQ) {
+double Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
+                                        PdgParticle absPdg, double m,
+                                        double qOverP, double absQ) {
   assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   // relative radiation length
-  const float x = slab.thicknessInX0();
+  const double x = slab.thicknessInX0();
   // particle momentum and energy
   // do not need to care about the sign since it is only used squared
-  const float momentum = absQ / qOverP;
-  const float energy = fastHypot(m, momentum);
+  const double momentum = absQ / qOverP;
+  const double energy = fastHypot(m, momentum);
 
-  float dEdx = computeBremsstrahlungLossMean(m, energy);
+  double dEdx = computeBremsstrahlungLossMean(m, energy);
 
   // muon- or muon+
   // TODO magic number 8_GeV
@@ -409,26 +411,26 @@ float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
   return dEdx * x;
 }
 
-float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
-                                            PdgParticle absPdg, float m,
-                                            float qOverP, float absQ) {
+double Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
+                                             PdgParticle absPdg, double m,
+                                             double qOverP, double absQ) {
   assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   // relative radiation length
-  const float x = slab.thicknessInX0();
+  const double x = slab.thicknessInX0();
   // particle momentum and energy
   // do not need to care about the sign since it is only used squared
-  const float momentum = absQ / qOverP;
-  const float energy = fastHypot(m, momentum);
+  const double momentum = absQ / qOverP;
+  const double energy = fastHypot(m, momentum);
 
   // compute derivative w/ respect to energy.
-  float derE = deriveBremsstrahlungLossMeanE(m);
+  double derE = deriveBremsstrahlungLossMeanE(m);
 
   // muon- or muon+
   // TODO magic number 8_GeV
@@ -441,80 +443,80 @@ float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
   //     E = sqrt(m² + p²) = sqrt(m² + q²/(q/p)²)
   // and the resulting derivative
   //     dE/d(q/p) = -q² / ((q/p)³ * E)
-  const float derQOverP = -(absQ * absQ) / (qOverP * qOverP * qOverP * energy);
+  const double derQOverP = -(absQ * absQ) / (qOverP * qOverP * qOverP * energy);
   return derE * derQOverP * x;
 }
 
-float Acts::computeEnergyLossMean(const MaterialSlab& slab, PdgParticle absPdg,
-                                  float m, float qOverP, float absQ) {
+double Acts::computeEnergyLossMean(const MaterialSlab& slab, PdgParticle absPdg,
+                                   double m, double qOverP, double absQ) {
   return computeEnergyLossBethe(slab, m, qOverP, absQ) +
          computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
 }
 
-float Acts::deriveEnergyLossMeanQOverP(const MaterialSlab& slab,
-                                       PdgParticle absPdg, float m,
-                                       float qOverP, float absQ) {
+double Acts::deriveEnergyLossMeanQOverP(const MaterialSlab& slab,
+                                        PdgParticle absPdg, double m,
+                                        double qOverP, double absQ) {
   return deriveEnergyLossBetheQOverP(slab, m, qOverP, absQ) +
          deriveEnergyLossRadiativeQOverP(slab, absPdg, m, qOverP, absQ);
 }
 
-float Acts::computeEnergyLossMode(const MaterialSlab& slab, PdgParticle absPdg,
-                                  float m, float qOverP, float absQ) {
+double Acts::computeEnergyLossMode(const MaterialSlab& slab, PdgParticle absPdg,
+                                   double m, double qOverP, double absQ) {
   // see ATL-SOFT-PUB-2008-003 section 3 for the relative fractions
   // TODO this is inconsistent with the text of the note
-  return 0.9f * computeEnergyLossLandau(slab, m, qOverP, absQ) +
-         0.15f * computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
+  return 0.9 * computeEnergyLossLandau(slab, m, qOverP, absQ) +
+         0.15 * computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
 }
 
-float Acts::deriveEnergyLossModeQOverP(const MaterialSlab& slab,
-                                       PdgParticle absPdg, float m,
-                                       float qOverP, float absQ) {
+double Acts::deriveEnergyLossModeQOverP(const MaterialSlab& slab,
+                                        PdgParticle absPdg, double m,
+                                        double qOverP, double absQ) {
   // see ATL-SOFT-PUB-2008-003 section 3 for the relative fractions
   // TODO this is inconsistent with the text of the note
-  return 0.9f * deriveEnergyLossLandauQOverP(slab, m, qOverP, absQ) +
-         0.15f * deriveEnergyLossRadiativeQOverP(slab, absPdg, m, qOverP, absQ);
+  return 0.9 * deriveEnergyLossLandauQOverP(slab, m, qOverP, absQ) +
+         0.15 * deriveEnergyLossRadiativeQOverP(slab, absPdg, m, qOverP, absQ);
 }
 
 namespace {
 
 /// Multiple scattering theta0 for minimum ionizing particles.
-inline float theta0Highland(float xOverX0, float momentumInv,
-                            float q2OverBeta2) {
+inline double theta0Highland(double xOverX0, double momentumInv,
+                             double q2OverBeta2) {
   // RPP2018 eq. 33.15 (treats beta and q² consistently)
-  const float t = std::sqrt(xOverX0 * q2OverBeta2);
+  const double t = std::sqrt(xOverX0 * q2OverBeta2);
   // log((x/X0) * (q²/beta²)) = log((sqrt(x/X0) * (q/beta))²)
   //                          = 2 * log(sqrt(x/X0) * (q/beta))
-  return 13.6_MeV * momentumInv * t * (1.0f + 0.038f * 2 * std::log(t));
+  return 13.6_MeV * momentumInv * t * (1 + 0.038 * 2 * std::log(t));
 }
 
 /// Multiple scattering theta0 for electrons.
-inline float theta0RossiGreisen(float xOverX0, float momentumInv,
-                                float q2OverBeta2) {
+inline double theta0RossiGreisen(double xOverX0, double momentumInv,
+                                 double q2OverBeta2) {
   // TODO add source paper/ resource
-  const float t = std::sqrt(xOverX0 * q2OverBeta2);
-  return 17.5_MeV * momentumInv * t *
-         (1.0f + 0.125f * std::log10(10.0f * xOverX0));
+  const double t = std::sqrt(xOverX0 * q2OverBeta2);
+  return 17.5_MeV * momentumInv * t * (1 + 0.125 * std::log10(10.0 * xOverX0));
 }
 
 }  // namespace
 
-float Acts::computeMultipleScatteringTheta0(const MaterialSlab& slab,
-                                            PdgParticle absPdg, float m,
-                                            float qOverP, float absQ) {
+double Acts::computeMultipleScatteringTheta0(const MaterialSlab& slab,
+                                             PdgParticle absPdg, double m,
+                                             double qOverP, double absQ) {
   assert((absPdg == Acts::makeAbsolutePdgParticle(absPdg)) &&
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
   if (!slab.isValid()) {
-    return 0.0f;
+    return 0;
   }
 
   // relative radiation length
-  const float xOverX0 = slab.thicknessInX0();
+  const double xOverX0 = slab.thicknessInX0();
   // 1/p = q/(pq) = (q/p)/q
-  const float momentumInv = std::abs(qOverP / absQ);
+  const double momentumInv = std::abs(qOverP / absQ);
   // q²/beta²; a smart compiler should be able to remove the unused computations
-  const float q2OverBeta2 = RelativisticQuantities(m, qOverP, absQ).q2OverBeta2;
+  const double q2OverBeta2 =
+      RelativisticQuantities(m, qOverP, absQ).q2OverBeta2;
 
   // electron or positron
   if (absPdg == PdgParticle::eElectron) {
diff --git a/Core/src/Material/Material.cpp b/Core/src/Material/Material.cpp
index 29d9502fde9..4e7b2f81d6c 100644
--- a/Core/src/Material/Material.cpp
+++ b/Core/src/Material/Material.cpp
@@ -15,22 +15,10 @@
 
 namespace Acts {
 
-namespace {
-enum MaterialClassificationNumberIndices {
-  eRadiationLength = 0,
-  eInteractionLength = 1,
-  eRelativeAtomicMass = 2,
-  eNuclearCharge = 3,
-  eMolarDensity = 4,
-};
-
-// Avogadro constant
-constexpr double kAvogadro = 6.02214076e23 / UnitConstants::mol;
-}  // namespace
-
-Material Material::fromMassDensity(float x0, float l0, float ar, float z,
-                                   float massRho) {
+Material Material::fromMassDensity(double x0, double l0, double ar, double z,
+                                   double massRho) {
   using namespace UnitLiterals;
+  using namespace PhysicalConstants;
 
   Material mat;
   mat.m_x0 = x0;
@@ -53,8 +41,8 @@ Material Material::fromMassDensity(float x0, float l0, float ar, float z,
   return mat;
 }
 
-Material Material::fromMolarDensity(float x0, float l0, float ar, float z,
-                                    float molarRho) {
+Material Material::fromMolarDensity(double x0, double l0, double ar, double z,
+                                    double molarRho) {
   Material mat;
   mat.m_x0 = x0;
   mat.m_l0 = l0;
@@ -71,8 +59,9 @@ Material::Material(const ParametersVector& parameters)
       m_z(parameters[eNuclearCharge]),
       m_molarRho(parameters[eMolarDensity]) {}
 
-float Material::massDensity() const {
+double Material::massDensity() const {
   using namespace UnitLiterals;
+  using namespace PhysicalConstants;
 
   // perform computations in double precision to avoid loss of precision
   const double atomicMass = static_cast<double>(m_ar) * 1_u;
@@ -80,7 +69,7 @@ float Material::massDensity() const {
   return atomicMass * numberDensity;
 }
 
-float Material::meanExcitationEnergy() const {
+double Material::meanExcitationEnergy() const {
   using namespace UnitLiterals;
 
   // use approximative computation as defined in ATL-SOFT-PUB-2008-003
diff --git a/Core/src/Material/MaterialMapUtils.cpp b/Core/src/Material/MaterialMapUtils.cpp
index b51b1a3e223..f6adfd8d11a 100644
--- a/Core/src/Material/MaterialMapUtils.cpp
+++ b/Core/src/Material/MaterialMapUtils.cpp
@@ -16,10 +16,7 @@
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
-#include <initializer_list>
 #include <limits>
-#include <set>
-#include <tuple>
 #include <utility>
 
 using Acts::VectorHelpers::perp;
@@ -90,10 +87,9 @@ auto Acts::materialMapperRZ(
           materialVectorToGridMapper({{i - 1, j - 1}}, nIndices));
     }
   }
-  Material::ParametersVector vec;
-  vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
-      0., 0., 0.;
-  grid.setExteriorBins(vec);
+
+  // use almost vacuum so the interpolation can handle the edges correctly
+  grid.setExteriorBins(Material::almostVacuumParameters());
 
   // [4] Create the transformation for the position
   // map (x,y,z) -> (r,z)
@@ -186,10 +182,9 @@ auto Acts::materialMapperXYZ(
       }
     }
   }
-  Material::ParametersVector vec;
-  vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
-      0., 0., 0.;
-  grid.setExteriorBins(vec);
+
+  // use almost vacuum so the interpolation can handle the edges correctly
+  grid.setExteriorBins(Material::almostVacuumParameters());
 
   // [4] Create the transformation for the position
   // map (x,y,z) -> (r,z)
diff --git a/Core/src/Material/MaterialSlab.cpp b/Core/src/Material/MaterialSlab.cpp
index 4fad731c6c1..5ae36b06095 100644
--- a/Core/src/Material/MaterialSlab.cpp
+++ b/Core/src/Material/MaterialSlab.cpp
@@ -17,12 +17,12 @@
 namespace Acts {
 
 namespace {
-static constexpr auto eps = 2 * std::numeric_limits<float>::epsilon();
+static constexpr auto eps = 2 * std::numeric_limits<double>::epsilon();
 }
 
-MaterialSlab::MaterialSlab(float thickness) : m_thickness(thickness) {}
+MaterialSlab::MaterialSlab(double thickness) : m_thickness(thickness) {}
 
-MaterialSlab::MaterialSlab(const Material& material, float thickness)
+MaterialSlab::MaterialSlab(const Material& material, double thickness)
     : m_material(material),
       m_thickness(thickness),
       m_thicknessInX0((eps < material.X0()) ? (thickness / material.X0()) : 0),
@@ -54,7 +54,7 @@ MaterialSlab MaterialSlab::averageLayers(
   return result;
 }
 
-void MaterialSlab::scaleThickness(float scale) {
+void MaterialSlab::scaleThickness(double scale) {
   if (scale < 0) {
     throw std::runtime_error("scale < 0");
   }
diff --git a/Core/src/Material/SurfaceMaterialMapper.cpp b/Core/src/Material/SurfaceMaterialMapper.cpp
index f4ad7ecc6f9..c905b5c457d 100644
--- a/Core/src/Material/SurfaceMaterialMapper.cpp
+++ b/Core/src/Material/SurfaceMaterialMapper.cpp
@@ -277,7 +277,7 @@ void Acts::SurfaceMaterialMapper::mapInteraction(
   GeometryIdentifier lastID = GeometryIdentifier();
   GeometryIdentifier currentID = GeometryIdentifier();
   Vector3 currentPos(0., 0., 0);
-  float currentPathCorrection = 1.;
+  double currentPathCorrection = 1;
   auto currentAccMaterial = mState.accumulatedMaterial.end();
 
   // To remember the bins of this event
diff --git a/Core/src/Material/VolumeMaterialMapper.cpp b/Core/src/Material/VolumeMaterialMapper.cpp
index ca7c00b5730..21b60b100db 100644
--- a/Core/src/Material/VolumeMaterialMapper.cpp
+++ b/Core/src/Material/VolumeMaterialMapper.cpp
@@ -253,7 +253,7 @@ void Acts::VolumeMaterialMapper::createExtraHits(
   // Computing the extra hits properties based on the mappingStep length
   int volumeStep =
       static_cast<int>(std::floor(properties.thickness() / m_cfg.mappingStep));
-  float remainder = properties.thickness() - m_cfg.mappingStep * volumeStep;
+  double remainder = properties.thickness() - m_cfg.mappingStep * volumeStep;
   properties.scaleThickness(m_cfg.mappingStep / properties.thickness());
   direction = direction * (m_cfg.mappingStep / direction.norm());
 
@@ -452,7 +452,7 @@ void Acts::VolumeMaterialMapper::mapMaterialTrack(
       if (currentBinning != mState.materialBin.end() &&
           rmIter->materialSlab.thickness() > 0) {
         // check if there is vacuum between this material point and the last one
-        float vacuumThickness = (rmIter->position - lastPositionEnd).norm();
+        double vacuumThickness = (rmIter->position - lastPositionEnd).norm();
         if (vacuumThickness > s_epsilon) {
           auto properties = Acts::MaterialSlab(vacuumThickness);
           // creat vacuum hits
@@ -482,7 +482,7 @@ void Acts::VolumeMaterialMapper::mapMaterialTrack(
             double distVol = (volIter->position - mTrack.first.first).norm();
             double distSur = (sfIter->position - mTrack.first.first).norm();
             if (distSur - distVol > s_epsilon) {
-              float vacuumThickness =
+              double vacuumThickness =
                   (sfIter->position - lastPositionEnd).norm();
               // if the last material slab stop before the boundary surface
               // create vacuum hits
diff --git a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
index b18fd373038..35efefbcdb1 100644
--- a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
+++ b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
@@ -28,7 +28,7 @@ void PointwiseMaterialInteraction::covarianceContributions(
   // Compute contributions from interactions
   if (multipleScattering) {
     // TODO use momentum before or after energy loss in backward mode?
-    const float theta0 =
+    const double theta0 =
         computeMultipleScatteringTheta0(slab, absPdg, mass, qOverP, absQ);
     // sigmaPhi = theta0 / sin(theta)
     const auto sigmaPhi = theta0 * (dir.norm() / VectorHelpers::perp(dir));
@@ -38,7 +38,7 @@ void PointwiseMaterialInteraction::covarianceContributions(
   }
   // TODO just ionisation loss or full energy loss?
   if (energyLoss) {
-    const float sigmaQoverP =
+    const double sigmaQoverP =
         computeEnergyLossLandauSigmaQOverP(slab, mass, qOverP, absQ);
     varianceQoverP = sigmaQoverP * sigmaQoverP;
   }
diff --git a/Tests/UnitTests/Core/Material/AccumulatedMaterialSlabTests.cpp b/Tests/UnitTests/Core/Material/AccumulatedMaterialSlabTests.cpp
index 34f9fe76f16..cc8bd03facc 100644
--- a/Tests/UnitTests/Core/Material/AccumulatedMaterialSlabTests.cpp
+++ b/Tests/UnitTests/Core/Material/AccumulatedMaterialSlabTests.cpp
@@ -14,7 +14,6 @@
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
 
-#include <limits>
 #include <utility>
 
 namespace {
@@ -25,7 +24,7 @@ using Acts::MaterialSlab;
 using Acts::Test::makeSilicon;
 using Acts::Test::makeUnitSlab;
 
-constexpr auto eps = std::numeric_limits<float>::epsilon();
+constexpr double eps = 1e-10;
 
 }  // namespace
 
@@ -74,8 +73,8 @@ BOOST_AUTO_TEST_CASE(MultipleIdenticalThicknessTrackSteps) {
     BOOST_CHECK_EQUAL(trackCount, 1u);
     BOOST_CHECK_EQUAL(average.material(), unit.material());
     BOOST_CHECK_EQUAL(average.thickness(), 3 * unit.thickness());
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0f);
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0f);
+    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0);
+    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0);
   }
   // accumulate three identical steps for one additional track
   {
@@ -88,8 +87,8 @@ BOOST_AUTO_TEST_CASE(MultipleIdenticalThicknessTrackSteps) {
     // averages must stay the same since we added the same material again
     BOOST_CHECK_EQUAL(average.material(), unit.material());
     BOOST_CHECK_EQUAL(average.thickness(), 3 * unit.thickness());
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0f);
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0f);
+    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0);
+    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0);
   }
 }
 
@@ -108,8 +107,8 @@ BOOST_AUTO_TEST_CASE(MultipleDifferentThicknessTrackSteps) {
     BOOST_CHECK_EQUAL(trackCount, 1u);
     BOOST_CHECK_EQUAL(average.material(), unit.material());
     BOOST_CHECK_EQUAL(average.thickness(), 3 * unit.thickness());
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0f);
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0f);
+    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0);
+    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0);
   }
   // accumulate one step with thickness 1, one with thickness 2
   {
@@ -123,8 +122,8 @@ BOOST_AUTO_TEST_CASE(MultipleDifferentThicknessTrackSteps) {
     // averages must stay the same
     BOOST_CHECK_EQUAL(average.material(), unit.material());
     BOOST_CHECK_EQUAL(average.thickness(), 3 * unit.thickness());
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0f);
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0f);
+    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0);
+    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0);
   }
   // accumulate one step with thickness 3
   {
@@ -137,8 +136,8 @@ BOOST_AUTO_TEST_CASE(MultipleDifferentThicknessTrackSteps) {
     // averages must stay the same
     BOOST_CHECK_EQUAL(average.material(), unit.material());
     BOOST_CHECK_EQUAL(average.thickness(), 3 * unit.thickness());
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0f);
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0f);
+    BOOST_CHECK_EQUAL(average.thicknessInX0(), 3.0);
+    BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.0);
   }
 }
 
@@ -182,19 +181,19 @@ BOOST_AUTO_TEST_CASE(MultipleDifferentTracks) {
     auto [average, trackCount] = a.totalAverage();
     BOOST_CHECK_EQUAL(trackCount, 4u);
     // average material density halved
-    CHECK_CLOSE_REL(average.material().X0(), 2 * unit.material().X0(), eps);
-    CHECK_CLOSE_REL(average.material().L0(), 2 * unit.material().L0(), eps);
-    CHECK_CLOSE_REL(average.material().molarDensity(),
-                    0.5f * unit.material().molarDensity(), eps);
+    CHECK_CLOSE_ABS(average.material().X0(), 2 * unit.material().X0(), eps);
+    CHECK_CLOSE_ABS(average.material().L0(), 2 * unit.material().L0(), eps);
+    CHECK_CLOSE_ABS(average.material().molarDensity(),
+                    0.5 * unit.material().molarDensity(), eps);
     // average atom is still the same species
-    CHECK_CLOSE_REL(average.material().Ar(), unit.material().Ar(), eps);
+    CHECK_CLOSE_ABS(average.material().Ar(), unit.material().Ar(), eps);
     // average atomic number proportional to the thickness
-    CHECK_CLOSE_REL(average.material().Z(), 0.5 * unit.material().Z(), eps);
+    CHECK_CLOSE_ABS(average.material().Z(), 0.5 * unit.material().Z(), eps);
     // thickness in x0/l0 depends on density and thus halved as well
-    BOOST_CHECK_EQUAL(average.thicknessInX0(), 1 * unit.thicknessInX0());
-    BOOST_CHECK_EQUAL(average.thicknessInL0(), 1 * unit.thicknessInL0());
+    CHECK_CLOSE_ABS(average.thicknessInX0(), 1 * unit.thicknessInX0(), eps);
+    CHECK_CLOSE_ABS(average.thicknessInL0(), 1 * unit.thicknessInL0(), eps);
     // average real thickness stays the same
-    BOOST_CHECK_EQUAL(average.thickness(), 2 * unit.thickness());
+    CHECK_CLOSE_ABS(average.thickness(), 2 * unit.thickness(), eps);
   }
 }
 
diff --git a/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp b/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
index ee8a40d9ff3..a937a93f825 100644
--- a/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
+++ b/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
@@ -21,13 +21,13 @@ namespace {
 
 using Acts::detail::combineSlabs;
 
-constexpr auto eps = std::numeric_limits<float>::epsilon();
+constexpr auto eps = std::numeric_limits<double>::epsilon();
 
 // vacuum w/ different thickness
-const Acts::MaterialSlab zeroVacuum = Acts::MaterialSlab(0.0f);
-const Acts::MaterialSlab unitVacuum = Acts::MaterialSlab(1.0f);
+const Acts::MaterialSlab zeroVacuum = Acts::MaterialSlab(0);
+const Acts::MaterialSlab unitVacuum = Acts::MaterialSlab(1);
 // same material corresponding to 0%, 1% and 100% radiation/interaction length
-const Acts::MaterialSlab zero(Acts::Test::makeSilicon(), 0.0f);
+const Acts::MaterialSlab zero(Acts::Test::makeSilicon(), 0);
 const Acts::MaterialSlab percent = Acts::Test::makePercentSlab();
 const Acts::MaterialSlab unit = Acts::Test::makeUnitSlab();
 
@@ -43,18 +43,18 @@ BOOST_AUTO_TEST_CASE(CombineSlabsVacuum) {
     auto slab = combineSlabs(zeroVacuum, zeroVacuum);
     BOOST_CHECK(!slab.isValid());
     BOOST_CHECK(!slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 0.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 0.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 0.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 0.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 0.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 0.0);
   }
   // vacuum with unit thickness
   {
     auto slab = combineSlabs(unitVacuum, unitVacuum);
     BOOST_CHECK(!slab.isValid());
     BOOST_CHECK(!slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 2.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 0.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 0.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 2.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 0.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 0.0);
   }
 }
 
@@ -186,13 +186,13 @@ BOOST_AUTO_TEST_CASE(CombineSlabsUnitZero) {
 
 BOOST_AUTO_TEST_CASE(CombineSlabsEqualThicknessVacuum) {
   const auto mat = Acts::Test::makeSilicon();
-  const auto slabMat = Acts::MaterialSlab(mat, 1.0f);
-  const auto slabVac = Acts::MaterialSlab(Acts::Material(), 1.0f);
+  const auto slabMat = Acts::MaterialSlab(mat, 1);
+  const auto slabVac = Acts::MaterialSlab(Acts::Material(), 1);
   {
     auto slab = combineSlabs(slabMat, slabVac);
     BOOST_CHECK(slab.isValid());
     BOOST_CHECK(slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 2.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 2.0);
     BOOST_CHECK_EQUAL(slab.thicknessInX0(), slabMat.thicknessInX0());
     BOOST_CHECK_EQUAL(slab.thicknessInL0(), slabMat.thicknessInL0());
     // atomic mass and nuclear charge are per atom, adding any amount vacuum
@@ -200,18 +200,17 @@ BOOST_AUTO_TEST_CASE(CombineSlabsEqualThicknessVacuum) {
     BOOST_CHECK_EQUAL(slab.material().Ar(), mat.Ar());
     BOOST_CHECK_EQUAL(slab.material().Z(), 0.5 * mat.Z());
     // we have the same type of interactions just spread over twice the length
-    CHECK_CLOSE_REL(slab.material().X0(), 2.0f * mat.X0(), eps);
-    CHECK_CLOSE_REL(slab.material().L0(), 2.0f * mat.L0(), eps);
+    CHECK_CLOSE_REL(slab.material().X0(), 2 * mat.X0(), eps);
+    CHECK_CLOSE_REL(slab.material().L0(), 2 * mat.L0(), eps);
     // we have the same atoms just spread over more volume
-    BOOST_CHECK_EQUAL(slab.material().molarDensity(),
-                      0.5f * mat.molarDensity());
+    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 0.5 * mat.molarDensity());
   }
   // reverse input order
   {
     auto slab = combineSlabs(slabVac, slabMat);
     BOOST_CHECK(slab.isValid());
     BOOST_CHECK(slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 2.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 2.0);
     BOOST_CHECK_EQUAL(slab.thicknessInX0(), slabMat.thicknessInX0());
     BOOST_CHECK_EQUAL(slab.thicknessInL0(), slabMat.thicknessInL0());
     // atomic mass and nuclear charge are per atom, adding any amount vacuum
@@ -219,11 +218,10 @@ BOOST_AUTO_TEST_CASE(CombineSlabsEqualThicknessVacuum) {
     BOOST_CHECK_EQUAL(slab.material().Ar(), mat.Ar());
     BOOST_CHECK_EQUAL(slab.material().Z(), 0.5 * mat.Z());
     // we have the same type of interactions just spread over twice the length
-    CHECK_CLOSE_REL(slab.material().X0(), 2.0f * mat.X0(), eps);
-    CHECK_CLOSE_REL(slab.material().L0(), 2.0f * mat.L0(), eps);
+    CHECK_CLOSE_REL(slab.material().X0(), 2 * mat.X0(), eps);
+    CHECK_CLOSE_REL(slab.material().L0(), 2 * mat.L0(), eps);
     // we have the same atoms just spread over more volume
-    BOOST_CHECK_EQUAL(slab.material().molarDensity(),
-                      0.5f * mat.molarDensity());
+    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 0.5 * mat.molarDensity());
   }
 }
 
@@ -232,51 +230,53 @@ BOOST_AUTO_TEST_CASE(CombineSlabsEqualThicknessVacuum) {
 BOOST_AUTO_TEST_CASE(CombineSlabs) {
   const auto mat0 = Acts::Material::fromMolarDensity(1, 1, 8, 12, 2);
   const auto mat1 = Acts::Material::fromMolarDensity(2, 2, 2, 6, 5);
-  const auto slabMat0 = Acts::MaterialSlab(mat0, 0.5f);
-  const auto slabMat1 = Acts::MaterialSlab(mat1, 1.0f);
+  const auto slabMat0 = Acts::MaterialSlab(mat0, 0.5);
+  const auto slabMat1 = Acts::MaterialSlab(mat1, 1.0);
   // verify derived quantities for the input slabs. these tests are not really
   // needed, but to show the input values for the tests below.
   BOOST_CHECK(slabMat0.isValid());
   BOOST_CHECK(slabMat0.material().isValid());
-  BOOST_CHECK_EQUAL(slabMat0.thickness(), 0.5f);
-  BOOST_CHECK_EQUAL(slabMat0.thicknessInX0(), 0.5f);
-  BOOST_CHECK_EQUAL(slabMat0.thicknessInL0(), 0.5f);
+  BOOST_CHECK_EQUAL(slabMat0.thickness(), 0.5);
+  BOOST_CHECK_EQUAL(slabMat0.thicknessInX0(), 0.5);
+  BOOST_CHECK_EQUAL(slabMat0.thicknessInL0(), 0.5);
   BOOST_CHECK(slabMat1.isValid());
   BOOST_CHECK(slabMat1.material().isValid());
-  BOOST_CHECK_EQUAL(slabMat1.thickness(), 1.0f);
-  BOOST_CHECK_EQUAL(slabMat1.thicknessInX0(), 0.5f);
-  BOOST_CHECK_EQUAL(slabMat1.thicknessInL0(), 0.5f);
+  BOOST_CHECK_EQUAL(slabMat1.thickness(), 1.0);
+  BOOST_CHECK_EQUAL(slabMat1.thicknessInX0(), 0.5);
+  BOOST_CHECK_EQUAL(slabMat1.thicknessInL0(), 0.5);
   // check combined slabs
   {
     auto slab = combineSlabs(slabMat0, slabMat1);
     BOOST_CHECK(slab.isValid());
     BOOST_CHECK(slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 1.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 1.0f);
-    BOOST_CHECK_EQUAL(slab.material().X0(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.material().L0(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.material().Ar(), 3.0f);
-    BOOST_CHECK_EQUAL(slab.material().Z(),
-                      static_cast<float>(
-                          exp((0.5 / 1.5) * log(12.0) + (1.0 / 1.5) * log(6))));
-    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 4.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 1.5);
+    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 1.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 1.0);
+    BOOST_CHECK_EQUAL(slab.material().X0(), 1.5);
+    BOOST_CHECK_EQUAL(slab.material().L0(), 1.5);
+    BOOST_CHECK_EQUAL(slab.material().Ar(), 3.0);
+    CHECK_CLOSE_ABS(
+        slab.material().Z(),
+        std::exp((0.5 / 1.5) * std::log(12.0) + (1.0 / 1.5) * std::log(6.0)),
+        1e-6);
+    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 4.0);
   }
   // reverse input order
   {
     auto slab = combineSlabs(slabMat0, slabMat1);
     BOOST_CHECK(slab.isValid());
     BOOST_CHECK(slab.material().isValid());
-    BOOST_CHECK_EQUAL(slab.thickness(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 1.0f);
-    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 1.0f);
-    BOOST_CHECK_EQUAL(slab.material().X0(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.material().L0(), 1.5f);
-    BOOST_CHECK_EQUAL(slab.material().Ar(), 3.0f);
-    BOOST_CHECK_EQUAL(slab.material().Z(),
-                      static_cast<float>(
-                          exp((0.5 / 1.5) * log(12.0) + (1.0 / 1.5) * log(6))));
-    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 4.0f);
+    BOOST_CHECK_EQUAL(slab.thickness(), 1.5);
+    BOOST_CHECK_EQUAL(slab.thicknessInX0(), 1.0);
+    BOOST_CHECK_EQUAL(slab.thicknessInL0(), 1.0);
+    BOOST_CHECK_EQUAL(slab.material().X0(), 1.5);
+    BOOST_CHECK_EQUAL(slab.material().L0(), 1.5);
+    BOOST_CHECK_EQUAL(slab.material().Ar(), 3.0);
+    CHECK_CLOSE_ABS(
+        slab.material().Z(),
+        std::exp((0.5 / 1.5) * std::log(12.0) + (1.0 / 1.5) * std::log(6.0)),
+        1e-6);
+    BOOST_CHECK_EQUAL(slab.material().molarDensity(), 4.0);
   }
 }
 
diff --git a/Tests/UnitTests/Core/Material/MaterialTests.cpp b/Tests/UnitTests/Core/Material/MaterialTests.cpp
index 8e2e7736449..d36b6d3c626 100644
--- a/Tests/UnitTests/Core/Material/MaterialTests.cpp
+++ b/Tests/UnitTests/Core/Material/MaterialTests.cpp
@@ -78,7 +78,6 @@ BOOST_DATA_TEST_CASE(EncodingDecodingRoundtrip,
   // encode material again
   Acts::Material::ParametersVector numbers1 = fromNumbers.parameters();
 
-  BOOST_CHECK_EQUAL(material, fromNumbers);
   BOOST_CHECK_EQUAL(numbers0, numbers1);
 }
 
diff --git a/Tests/UnitTests/Core/Utilities/MaterialMapUtilsTests.cpp b/Tests/UnitTests/Core/Utilities/MaterialMapUtilsTests.cpp
index 8f02a12ba68..ef8436b52e8 100644
--- a/Tests/UnitTests/Core/Utilities/MaterialMapUtilsTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/MaterialMapUtilsTests.cpp
@@ -121,4 +121,5 @@ BOOST_AUTO_TEST_CASE(materialmap_creation) {
   CHECK_CLOSE_ABS(value1_xyz.parameters(), mat1_xyz.parameters(), 1e-9);
   CHECK_CLOSE_ABS(value2_xyz.parameters(), mat2_xyz.parameters(), 1e-9);
 }
+
 }  // namespace Acts::Test