diff --git a/Core/include/Acts/EventData/ChargeConcept.hpp b/Core/include/Acts/EventData/ChargeConcept.hpp
index 4b0540ac407..aea68f72ebd 100644
--- a/Core/include/Acts/EventData/ChargeConcept.hpp
+++ b/Core/include/Acts/EventData/ChargeConcept.hpp
@@ -20,14 +20,14 @@
 namespace Acts {
 
 template <typename C>
-concept ChargeConcept = requires(C c, C c2, float f, double d) {
+concept ChargeConcept = requires(C c, C c2, double f, double d) {
   { C{f} };
 
   { c == c2 } -> std::same_as<bool>;
   { c != c2 } -> std::same_as<bool>;
 
-  { c.absQ() } -> std::same_as<float>;
-  { c.extractCharge(d) } -> std::same_as<float>;
+  { c.absQ() } -> std::same_as<double>;
+  { c.extractCharge(d) } -> std::same_as<double>;
   { c.extractMomentum(d) } -> std::same_as<double>;
   { c.qOverP(d, d) } -> std::same_as<double>;
 };
diff --git a/Core/include/Acts/EventData/GenericParticleHypothesis.hpp b/Core/include/Acts/EventData/GenericParticleHypothesis.hpp
index 38c70409dc6..520aee1480a 100644
--- a/Core/include/Acts/EventData/GenericParticleHypothesis.hpp
+++ b/Core/include/Acts/EventData/GenericParticleHypothesis.hpp
@@ -33,7 +33,7 @@ class GenericParticleHypothesis {
   /// @param absPdg the absolute PDG
   /// @param mass the particle mass
   /// @param chargeType the type of charge
-  constexpr GenericParticleHypothesis(PdgParticle absPdg, float mass,
+  constexpr GenericParticleHypothesis(PdgParticle absPdg, double mass,
                                       ChargeType chargeType)
       : m_absPdg{absPdg}, m_mass{mass}, m_chargeType{std::move(chargeType)} {
     assert(absPdg == makeAbsolutePdgParticle(absPdg) &&
@@ -67,10 +67,10 @@ class GenericParticleHypothesis {
   constexpr PdgParticle absolutePdg() const noexcept { return m_absPdg; }
 
   /// Get the hypothesized mass.
-  constexpr float mass() const noexcept { return m_mass; }
+  constexpr double mass() const noexcept { return m_mass; }
 
   /// Get the hypothesized absolute charge.
-  constexpr float absoluteCharge() const noexcept {
+  constexpr double absoluteCharge() const noexcept {
     return m_chargeType.absQ();
   }
 
@@ -125,7 +125,7 @@ class GenericParticleHypothesis {
 
  private:
   PdgParticle m_absPdg;
-  float m_mass;
+  double m_mass;
   ChargeType m_chargeType;
 
   friend bool operator==(const GenericParticleHypothesis<ChargeType>& lhs,
diff --git a/Core/include/Acts/EventData/ParticleHypothesis.hpp b/Core/include/Acts/EventData/ParticleHypothesis.hpp
index 8fb7fe9ec03..c6593215f66 100644
--- a/Core/include/Acts/EventData/ParticleHypothesis.hpp
+++ b/Core/include/Acts/EventData/ParticleHypothesis.hpp
@@ -25,7 +25,7 @@ namespace Acts {
 class SinglyChargedParticleHypothesis
     : public GenericParticleHypothesis<SinglyCharged> {
  public:
-  constexpr SinglyChargedParticleHypothesis(PdgParticle absPdg, float mass)
+  constexpr SinglyChargedParticleHypothesis(PdgParticle absPdg, double mass)
       : GenericParticleHypothesis(absPdg, mass, {}) {}
   SinglyChargedParticleHypothesis(PdgParticle absPdg)
       : GenericParticleHypothesis(absPdg) {}
@@ -68,7 +68,7 @@ class SinglyChargedParticleHypothesis
 /// @note This serves as a factory for common neutral particles.
 class NeutralParticleHypothesis : public GenericParticleHypothesis<Neutral> {
  public:
-  constexpr NeutralParticleHypothesis(PdgParticle absPdg, float mass)
+  constexpr NeutralParticleHypothesis(PdgParticle absPdg, double mass)
       : GenericParticleHypothesis(absPdg, mass, {}) {}
   NeutralParticleHypothesis(PdgParticle absPdg)
       : GenericParticleHypothesis(absPdg) {}
@@ -99,7 +99,7 @@ class NeutralParticleHypothesis : public GenericParticleHypothesis<Neutral> {
 class NonNeutralChargedParticleHypothesis
     : public GenericParticleHypothesis<NonNeutralCharge> {
  public:
-  constexpr NonNeutralChargedParticleHypothesis(PdgParticle absPdg, float mass,
+  constexpr NonNeutralChargedParticleHypothesis(PdgParticle absPdg, double mass,
                                                 NonNeutralCharge chargeType)
       : GenericParticleHypothesis(absPdg, mass, chargeType) {}
   NonNeutralChargedParticleHypothesis(PdgParticle absPdg)
@@ -126,7 +126,7 @@ class NonNeutralChargedParticleHypothesis
     return SinglyChargedParticleHypothesis::proton();
   }
 
-  static NonNeutralChargedParticleHypothesis pionLike(float absQ) {
+  static NonNeutralChargedParticleHypothesis pionLike(double absQ) {
     return NonNeutralChargedParticleHypothesis(pion().absolutePdg(),
                                                pion().mass(), absQ);
   }
@@ -135,7 +135,7 @@ class NonNeutralChargedParticleHypothesis
     static const auto cache = chargedGeantino(Acts::UnitConstants::e);
     return cache;
   }
-  static NonNeutralChargedParticleHypothesis chargedGeantino(float absQ) {
+  static NonNeutralChargedParticleHypothesis chargedGeantino(double absQ) {
     return NonNeutralChargedParticleHypothesis(PdgParticle::eInvalid, 0, absQ);
   }
 };
@@ -145,7 +145,7 @@ class NonNeutralChargedParticleHypothesis
 /// @note This serves as a factory for common particles with any kind of charge.
 class ParticleHypothesis : public GenericParticleHypothesis<AnyCharge> {
  public:
-  constexpr ParticleHypothesis(PdgParticle absPdg, float mass,
+  constexpr ParticleHypothesis(PdgParticle absPdg, double mass,
                                AnyCharge chargeType)
       : GenericParticleHypothesis(absPdg, mass, chargeType) {}
   ParticleHypothesis(PdgParticle absPdg) : GenericParticleHypothesis(absPdg) {}
@@ -178,7 +178,7 @@ class ParticleHypothesis : public GenericParticleHypothesis<AnyCharge> {
     return NeutralParticleHypothesis::pion0();
   }
 
-  static ParticleHypothesis pionLike(float absQ) {
+  static ParticleHypothesis pionLike(double absQ) {
     return ParticleHypothesis(pion().absolutePdg(), pion().mass(), absQ);
   }
 
@@ -189,7 +189,7 @@ class ParticleHypothesis : public GenericParticleHypothesis<AnyCharge> {
     static const auto cache = chargedGeantino(Acts::UnitConstants::e);
     return cache;
   }
-  static ParticleHypothesis chargedGeantino(float absQ) {
+  static ParticleHypothesis chargedGeantino(double absQ) {
     return ParticleHypothesis(PdgParticle::eInvalid, 0, absQ);
   }
 };
diff --git a/Core/include/Acts/Material/Interactions.hpp b/Core/include/Acts/Material/Interactions.hpp
index 4ca87b1ca0c..57d95ac09ef 100644
--- a/Core/include/Acts/Material/Interactions.hpp
+++ b/Core/include/Acts/Material/Interactions.hpp
@@ -30,13 +30,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 +46,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 +61,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 +87,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 +107,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 +133,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/Propagator/detail/VolumeMaterialInteraction.hpp b/Core/include/Acts/Propagator/detail/VolumeMaterialInteraction.hpp
index e5a8efbd29f..63f0c9833e9 100644
--- a/Core/include/Acts/Propagator/detail/VolumeMaterialInteraction.hpp
+++ b/Core/include/Acts/Propagator/detail/VolumeMaterialInteraction.hpp
@@ -29,13 +29,13 @@ struct VolumeMaterialInteraction {
   /// The particle current direction
   const Vector3 dir = Vector3::Zero();
   /// The particle q/p at the interaction
-  const float qOverP = 0;
+  const double qOverP = 0;
   /// The absolute particle charge
-  const float absQ = 0;
+  const double absQ = 0;
   /// The particle momentum at the interaction
-  const float momentum = 0;
+  const double momentum = 0;
   /// The particle mass
-  const float mass = 0;
+  const double mass = 0;
   /// The particle pdg
   const PdgParticle absPdg = eInvalid;
   /// The covariance transport decision at the interaction
diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
index e87a52410d8..2e835855440 100644
--- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
+++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
@@ -665,12 +665,10 @@ class Gx2Fitter {
             const auto& particle =
                 parametersWithHypothesis->particleHypothesis();
 
-            const double sigma =
-                static_cast<double>(Acts::computeMultipleScatteringTheta0(
-                    interaction.slab, particle.absolutePdg(), particle.mass(),
-                    static_cast<float>(
-                        parametersWithHypothesis->parameters()[eBoundQOverP]),
-                    particle.absoluteCharge()));
+            const double sigma = Acts::computeMultipleScatteringTheta0(
+                interaction.slab, particle.absolutePdg(), particle.mass(),
+                parametersWithHypothesis->parameters()[eBoundQOverP],
+                particle.absoluteCharge());
             ACTS_VERBOSE(
                 "        The Highland formula gives sigma = " << sigma);
             invSigma2 = 1. / std::pow(sigma, 2);
diff --git a/Core/src/Definitions/ParticleData.cpp b/Core/src/Definitions/ParticleData.cpp
index 7254dd11dfe..4ed30e06042 100644
--- a/Core/src/Definitions/ParticleData.cpp
+++ b/Core/src/Definitions/ParticleData.cpp
@@ -51,18 +51,18 @@ static inline auto findByPdg(std::int32_t pdg, const ColumnContainer& column)
   return column[*index];
 }
 
-static constexpr inline float extractCharge(float value) {
+static constexpr inline double extractCharge(double value) {
   // convert three charge to regular charge in native units
   return (value / 3.0f) * Acts::UnitConstants::e;
 }
 
-static constexpr inline float extractMass(float value) {
+static constexpr inline double extractMass(double value) {
   return value * Acts::UnitConstants::MeV;
 }
 
 }  // namespace
 
-std::optional<float> Acts::findCharge(Acts::PdgParticle pdg) {
+std::optional<double> Acts::findCharge(Acts::PdgParticle pdg) {
   const auto charge =
       findByPdg(static_cast<std::int32_t>(pdg), kParticlesThreeCharge);
   if (!charge) {
@@ -71,7 +71,7 @@ std::optional<float> Acts::findCharge(Acts::PdgParticle pdg) {
   return extractCharge(*charge);
 }
 
-std::optional<float> Acts::findMass(Acts::PdgParticle pdg) {
+std::optional<double> Acts::findMass(Acts::PdgParticle pdg) {
   const auto mass =
       findByPdg(static_cast<std::int32_t>(pdg), kParticlesMassMeV);
   if (!mass) {
diff --git a/Core/src/Material/AccumulatedMaterialSlab.cpp b/Core/src/Material/AccumulatedMaterialSlab.cpp
index f1099578b46..f740277fea1 100644
--- a/Core/src/Material/AccumulatedMaterialSlab.cpp
+++ b/Core/src/Material/AccumulatedMaterialSlab.cpp
@@ -15,7 +15,7 @@ void Acts::AccumulatedMaterialSlab::accumulate(MaterialSlab slab,
                                                float pathCorrection) {
   // scale the recorded material to the equivalence contribution along the
   // surface normal
-  slab.scaleThickness(1 / pathCorrection);
+  slab.scaleThickness(1. / static_cast<double>(pathCorrection));
   m_trackAverage = detail::combineSlabs(m_trackAverage, slab);
 }
 
diff --git a/Core/src/Material/Interactions.cpp b/Core/src/Material/Interactions.cpp
index b6a725cb4ca..367c3ec8cb8 100644
--- a/Core/src/Material/Interactions.cpp
+++ b/Core/src/Material/Interactions.cpp
@@ -20,53 +20,58 @@ 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");
+
     // beta²/q² = (p/E)²/q² = p²/(q²m² + q²p²) = 1/(q² + (m²(q/p)²)
     // q²/beta² = q² + m²(q/p)²
     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;
     assert((betaGamma >= 0) && "Negative betaGamma");
+
     // gamma = sqrt(m² + p²)/m = sqrt(1 + (p/m)²)
-    gamma = std::sqrt(1.0f + pOverM * pOverM);
+    gamma = std::sqrt(1. + pOverM * pOverM);
   }
 };
 
 /// Compute q/p derivative of beta².
-inline float deriveBeta2(float qOverP, const RelativisticQuantities& rq) {
-  return -2 / (qOverP * rq.gamma * rq.gamma);
+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;
 }
@@ -74,22 +79,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));
 }
 
@@ -101,40 +106,42 @@ 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.) {
+    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. : (-1. / qOverP);
 }
 
 /// Convert Landau full-width-half-maximum to an equivalent Gaussian sigma,
@@ -144,22 +151,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.)));
+  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);
 }
@@ -168,46 +175,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)
@@ -219,51 +226,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)
@@ -275,43 +282,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
@@ -319,20 +326,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);
 }
 
@@ -344,7 +351,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
@@ -359,7 +366,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;
@@ -369,7 +376,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 {
@@ -379,25 +386,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 = std::hypot(m, momentum);
+  const double momentum = absQ / qOverP;
+  const double energy = std::hypot(m, momentum);
 
-  float dEdx = computeBremsstrahlungLossMean(m, energy);
+  double dEdx = computeBremsstrahlungLossMean(m, energy);
 
   // muon- or muon+
   // TODO magic number 8_GeV
@@ -408,26 +415,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 = std::hypot(m, momentum);
+  const double momentum = absQ / qOverP;
+  const double energy = std::hypot(m, momentum);
 
   // compute derivative w/ respect to energy.
-  float derE = deriveBremsstrahlungLossMeanE(m);
+  double derE = deriveBremsstrahlungLossMeanE(m);
 
   // muon- or muon+
   // TODO magic number 8_GeV
@@ -440,80 +447,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. * 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/VolumeMaterialMapper.cpp b/Core/src/Material/VolumeMaterialMapper.cpp
index 2bb9113f4a9..6b7452c56af 100644
--- a/Core/src/Material/VolumeMaterialMapper.cpp
+++ b/Core/src/Material/VolumeMaterialMapper.cpp
@@ -254,7 +254,8 @@ 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;
+  const double remainder =
+      properties.thickness() - m_cfg.mappingStep * volumeStep;
   properties.scaleThickness(m_cfg.mappingStep / properties.thickness());
   direction = direction * (m_cfg.mappingStep / direction.norm());
 
@@ -453,7 +454,8 @@ 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();
+        const double vacuumThickness =
+            (rmIter->position - lastPositionEnd).norm();
         if (vacuumThickness > s_epsilon) {
           auto properties = Acts::MaterialSlab(vacuumThickness);
           // creat vacuum hits
@@ -483,7 +485,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 =
+              const 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 8314003111e..2be7d333012 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/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
index 1e5deef90b5..2af9e7915fc 100644
--- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
+++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
@@ -267,7 +267,7 @@ struct SimulationActor {
   void interact(const Acts::MaterialSlab &slab, result_type &result) const {
     // run the continuous processes over a fraction of the material. returns
     // true on break condition (same as the underlying physics lists).
-    auto runContinuousPartial = [&, this](float fraction) {
+    auto runContinuousPartial = [&, this](double fraction) {
       Acts::MaterialSlab partialSlab = slab;
       partialSlab.scaleThickness(fraction);
       // material after passing this slab
@@ -313,22 +313,20 @@ struct SimulationActor {
     //   material. simulation is limited to the continuous processes.
     //
     // `clamp` ensures a valid range in all cases.
-    const float fracX0 =
-        std::clamp(static_cast<float>(x0Dist / slabX0), 0.0f, 1.0f);
-    const float fracL0 =
-        std::clamp(static_cast<float>(l0Dist / slabL0), 0.0f, 1.0f);
+    const double fracX0 = std::clamp(x0Dist / slabX0, 0., 1.);
+    const double fracL0 = std::clamp(l0Dist / slabL0, 0., 1.);
     // fraction of the material where the first point-like interaction occurs
-    const float frac = std::min(fracX0, fracL0);
+    const double frac = std::min(fracX0, fracL0);
 
     // do not run if there is zero material before the point-like interaction
-    if (0.0f < frac) {
+    if (0. < frac) {
       // simulate continuous processes before the point-like interaction
       if (runContinuousPartial(frac)) {
         return;
       }
     }
     // do not run if there is no point-like interaction
-    if (frac < 1.0f) {
+    if (frac < 1.) {
       // select which process to simulate
       const std::size_t process =
           (fracX0 < fracL0) ? result.x0Process : result.l0Process;
diff --git a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp
index 08c6a8a42bb..933785cfa64 100644
--- a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp
+++ b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp
@@ -41,14 +41,14 @@ struct BetheBloch {
                                      const Acts::MaterialSlab &slab,
                                      Particle &particle) const {
     // compute energy loss distribution parameters
-    const float m = particle.mass();
-    const float qOverP = particle.qOverP();
-    const float absQ = particle.absoluteCharge();
+    const double m = particle.mass();
+    const double qOverP = particle.qOverP();
+    const double absQ = particle.absoluteCharge();
     // most probable value
     const double energyLoss =
         Acts::computeEnergyLossLandau(slab, m, qOverP, absQ);
     // Gaussian-equivalent sigma
-    const float energyLossSigma =
+    const double energyLossSigma =
         Acts::computeEnergyLossLandauSigma(slab, m, qOverP, absQ);
 
     // Simulate the energy loss
diff --git a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp
index a408b1ff97f..5b1d54da038 100644
--- a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp
+++ b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp
@@ -183,7 +183,7 @@ struct NuclearInteraction {
   const detail::NuclearInteractionParameters& findParameters(
       double rnd,
       const detail::NuclearInteractionParametrisation& parametrisation,
-      float particleMomentum) const;
+      double particleMomentum) const;
 
   /// Estimates the interaction type
   ///
@@ -191,7 +191,7 @@ struct NuclearInteraction {
   /// @param [in] probability The probability for a soft interaction
   ///
   /// @return True if a soft interaction occurs
-  inline bool softInteraction(double rnd, float probability) const {
+  inline bool softInteraction(double rnd, double probability) const {
     return rnd <= probability;
   }
 
@@ -248,7 +248,7 @@ struct NuclearInteraction {
       generator_t& generator,
       const detail::NuclearInteractionParameters::
           ParametersWithFixedMultiplicity& parametrisation,
-      float initialMomentum) const;
+      double initialMomentum) const;
 
   /// Tests whether the final state momenta and invariant masses are
   /// matching to each other to allow the evaluation of particle directions.
@@ -261,7 +261,7 @@ struct NuclearInteraction {
   /// not.
   bool match(const Acts::ActsDynamicVector& momenta,
              const Acts::ActsDynamicVector& invariantMasses,
-             float parametrizedMomentum) const;
+             double parametrizedMomentum) const;
 
   /// This method samples the kinematics of the final state particles
   ///
@@ -276,7 +276,7 @@ struct NuclearInteraction {
   sampleKinematics(generator_t& generator,
                    const detail::NuclearInteractionParameters::
                        ParametersWithFixedMultiplicity& parameters,
-                   float momentum) const;
+                   double momentum) const;
 
   /// Converts relative angles to absolute angles wrt the global
   /// coordinate system.
@@ -293,8 +293,8 @@ struct NuclearInteraction {
   /// coordinate system
   std::pair<ActsFatras::Particle::Scalar, ActsFatras::Particle::Scalar>
   globalAngle(ActsFatras::Particle::Scalar phi1,
-              ActsFatras::Particle::Scalar theta1, float phi2,
-              float theta2) const;
+              ActsFatras::Particle::Scalar theta1, double phi2,
+              double theta2) const;
 
   /// Converter from sampled numbers to a vector of particles
   ///
@@ -313,7 +313,7 @@ struct NuclearInteraction {
       generator_t& generator, const std::vector<int>& pdgId,
       const Acts::ActsDynamicVector& momenta,
       const Acts::ActsDynamicVector& invariantMasses, Particle& initialParticle,
-      float parametrizedMomentum, bool soft) const;
+      double parametrizedMomentum, bool soft) const;
 
   /// This function performs an inverse sampling to provide a discrete
   /// value from a distribution.
@@ -357,7 +357,7 @@ std::vector<int> NuclearInteraction::samplePdgIds(
   std::vector<int> pdgIds;
   pdgIds.reserve(multiplicity);
 
-  std::uniform_real_distribution<float> uniformDistribution{0., 1.};
+  std::uniform_real_distribution<double> uniformDistribution{0., 1.};
 
   // Find the producers probability distribution
   auto citProducer = pdgMap.cbegin();
@@ -365,19 +365,19 @@ std::vector<int> NuclearInteraction::samplePdgIds(
     citProducer++;
   }
 
-  const std::vector<std::pair<int, float>>& mapInitial = citProducer->second;
+  const std::vector<std::pair<int, double>>& mapInitial = citProducer->second;
   // Set the first particle depending on the interaction type
   if (soft) {
     // Store the initial particle if the interaction is soft
     pdgIds.push_back(particlePdg);
   } else {
     // Otherwise dice the particle
-    const float rndInitial = uniformDistribution(generator);
+    const double rndInitial = uniformDistribution(generator);
 
     pdgIds.push_back(
         std::lower_bound(mapInitial.begin(), mapInitial.end(), rndInitial,
-                         [](const std::pair<int, float>& element,
-                            float random) { return element.second < random; })
+                         [](const std::pair<int, double>& element,
+                            double random) { return element.second < random; })
             ->first);
   }
 
@@ -391,12 +391,12 @@ std::vector<int> NuclearInteraction::samplePdgIds(
     }
 
     // Set the next particle
-    const std::vector<std::pair<int, float>>& map = citProducer->second;
-    const float rnd = uniformDistribution(generator);
+    const std::vector<std::pair<int, double>>& map = citProducer->second;
+    const double rnd = uniformDistribution(generator);
     pdgIds.push_back(
         std::lower_bound(map.begin(), map.end(), rnd,
-                         [](const std::pair<int, float>& element,
-                            float random) { return element.second < random; })
+                         [](const std::pair<int, double>& element,
+                            double random) { return element.second < random; })
             ->first);
   }
   return pdgIds;
@@ -414,7 +414,7 @@ Acts::ActsDynamicVector NuclearInteraction::sampleInvariantMasses(
 
   // Sample in the eigenspace
   for (unsigned int i = 0; i < size; i++) {
-    float variance = parametrisation.eigenvaluesInvariantMass[i];
+    double variance = parametrisation.eigenvaluesInvariantMass[i];
     std::normal_distribution<Acts::ActsScalar> dist{
         parametrisation.meanInvariantMass[i], sqrtf(variance)};
     parameters[i] = dist(generator);
@@ -436,7 +436,7 @@ Acts::ActsDynamicVector NuclearInteraction::sampleMomenta(
     generator_t& generator,
     const detail::NuclearInteractionParameters::ParametersWithFixedMultiplicity&
         parametrisation,
-    float initialMomentum) const {
+    double initialMomentum) const {
   // The resulting vector
   Acts::ActsDynamicVector parameters;
   const unsigned int size = parametrisation.eigenvaluesMomentum.size();
@@ -444,7 +444,7 @@ Acts::ActsDynamicVector NuclearInteraction::sampleMomenta(
 
   // Sample in the eigenspace
   for (unsigned int i = 0; i < size; i++) {
-    float variance = parametrisation.eigenvaluesMomentum[i];
+    double variance = parametrisation.eigenvaluesMomentum[i];
     std::normal_distribution<Acts::ActsScalar> dist{
         parametrisation.meanMomentum[i], sqrtf(variance)};
     parameters[i] = dist(generator);
@@ -455,15 +455,15 @@ Acts::ActsDynamicVector NuclearInteraction::sampleMomenta(
 
   // Perform the inverse sampling from the distributions
   for (unsigned int i = 0; i < size; i++) {
-    const float cdf = (std::erff(parameters[i]) + 1) * 0.5;
+    const double cdf = (std::erff(parameters[i]) + 1) * 0.5;
     parameters[i] =
         sampleContinuousValues(cdf, parametrisation.momentumDistributions[i]);
   }
 
   // Scale the momenta
   Acts::ActsDynamicVector momenta = parameters.head(size - 1);
-  const float sum = momenta.sum();
-  const float scale = parameters.template tail<1>()(0, 0) / sum;
+  const double sum = momenta.sum();
+  const double scale = parameters.template tail<1>()(0, 0) / sum;
   momenta *= scale * initialMomentum;
   return momenta;
 }
@@ -474,7 +474,7 @@ NuclearInteraction::sampleKinematics(
     generator_t& generator,
     const detail::NuclearInteractionParameters::ParametersWithFixedMultiplicity&
         parameters,
-    float momentum) const {
+    double momentum) const {
   unsigned int trials = 0;
   Acts::ActsDynamicVector invariantMasses =
       sampleInvariantMasses(generator, parameters);
@@ -500,7 +500,7 @@ std::vector<Particle> NuclearInteraction::convertParametersToParticles(
     generator_t& generator, const std::vector<int>& pdgId,
     const Acts::ActsDynamicVector& momenta,
     const Acts::ActsDynamicVector& invariantMasses, Particle& initialParticle,
-    float parametrizedMomentum, bool soft) const {
+    double parametrizedMomentum, bool soft) const {
   std::uniform_real_distribution<double> uniformDistribution{0., 1.};
   const auto& initialDirection = initialParticle.direction();
   const double phi = Acts::VectorHelpers::phi(initialDirection);
@@ -512,10 +512,10 @@ std::vector<Particle> NuclearInteraction::convertParametersToParticles(
 
   // Build the particles
   for (unsigned int i = 0; i < size; i++) {
-    const float momentum = momenta[i];
-    const float invariantMass = invariantMasses[i];
-    const float p1p2 = 2. * momentum * parametrizedMomentum;
-    const float costheta = 1. - invariantMass * invariantMass / p1p2;
+    const double momentum = momenta[i];
+    const double invariantMass = invariantMasses[i];
+    const double p1p2 = 2. * momentum * parametrizedMomentum;
+    const double costheta = 1. - invariantMass * invariantMass / p1p2;
 
     const auto phiTheta =
         globalAngle(phi, theta, uniformDistribution(generator) * 2. * M_PI,
diff --git a/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp b/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp
index de9446a046e..168f697a319 100644
--- a/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp
+++ b/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp
@@ -20,7 +20,7 @@ namespace ActsFatras {
 const detail::NuclearInteractionParameters& NuclearInteraction::findParameters(
     double rnd,
     const detail::NuclearInteractionParametrisation& parametrisation,
-    float particleMomentum) const {
+    double particleMomentum) const {
   // Return lowest/highest if momentum outside the boundary
   if (particleMomentum <= parametrisation.front().first) {
     return parametrisation.front().second;
@@ -32,16 +32,16 @@ const detail::NuclearInteractionParameters& NuclearInteraction::findParameters(
   // Find the two neighbouring parametrisations
   const auto lowerBound = std::lower_bound(
       parametrisation.begin(), parametrisation.end(), particleMomentum,
-      [](const std::pair<const float,
+      [](const std::pair<const double,
                          ActsFatras::detail::NuclearInteractionParameters>&
              params,
-         const float mom) { return params.first < mom; });
-  const float momentumUpperNeighbour = lowerBound->first;
-  const float momentumLowerNeighbour = std::prev(lowerBound, 1)->first;
+         const double mom) { return params.first < mom; });
+  const double momentumUpperNeighbour = lowerBound->first;
+  const double momentumLowerNeighbour = std::prev(lowerBound, 1)->first;
 
   // Pick one randomly
-  const float weight = (momentumUpperNeighbour - particleMomentum) /
-                       (momentumUpperNeighbour - momentumLowerNeighbour);
+  const double weight = (momentumUpperNeighbour - particleMomentum) /
+                        (momentumUpperNeighbour - momentumLowerNeighbour);
   return (rnd < weight) ? std::prev(lowerBound, 1)->second : lowerBound->second;
 }  // namespace ActsFatras
 
@@ -113,8 +113,8 @@ unsigned int NuclearInteraction::finalStateMultiplicity(
 
 std::pair<ActsFatras::Particle::Scalar, ActsFatras::Particle::Scalar>
 NuclearInteraction::globalAngle(ActsFatras::Particle::Scalar phi1,
-                                ActsFatras::Particle::Scalar theta1, float phi2,
-                                float theta2) const {
+                                ActsFatras::Particle::Scalar theta1,
+                                double phi2, double theta2) const {
   // Rotation around the global y-axis
   Acts::SquareMatrix3 rotY = Acts::SquareMatrix3::Zero();
   rotY(0, 0) = std::cos(theta1);
@@ -138,22 +138,22 @@ NuclearInteraction::globalAngle(ActsFatras::Particle::Scalar phi1,
   const Acts::Vector3 vectorSum = rotZ * rotY * vector2;
 
   // Calculate the global angles
-  const float theta = std::acos(vectorSum.z() / vectorSum.norm());
-  const float phi = std::atan2(vectorSum.y(), vectorSum.x());
+  const double theta = std::acos(vectorSum.z() / vectorSum.norm());
+  const double phi = std::atan2(vectorSum.y(), vectorSum.x());
 
   return std::make_pair(phi, theta);
 }
 
 bool NuclearInteraction::match(const Acts::ActsDynamicVector& momenta,
                                const Acts::ActsDynamicVector& invariantMasses,
-                               float parametrizedMomentum) const {
+                               double parametrizedMomentum) const {
   const unsigned int size = momenta.size();
   for (unsigned int i = 0; i < size; i++) {
     // Calculate the invariant masses
-    const float momentum = momenta[i];
-    const float invariantMass = invariantMasses[i];
-    const float p1p2 = 2. * momentum * parametrizedMomentum;
-    const float costheta = 1. - invariantMass * invariantMass / p1p2;
+    const double momentum = momenta[i];
+    const double invariantMass = invariantMasses[i];
+    const double p1p2 = 2. * momentum * parametrizedMomentum;
+    const double costheta = 1. - invariantMass * invariantMass / p1p2;
 
     // Abort if an angle cannot be calculated
     if (std::abs(costheta) > 1) {
diff --git a/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp b/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
index 7004f490cc5..45dd6194e29 100644
--- a/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
+++ b/Tests/UnitTests/Core/Material/AverageMaterialsTests.cpp
@@ -274,7 +274,7 @@ BOOST_AUTO_TEST_CASE(CombineSlabs) {
   }
   // reverse input order
   {
-    auto slab = combineSlabs(slabMat0, slabMat1);
+    auto slab = combineSlabs(slabMat1, slabMat0);
     BOOST_CHECK(slab.isValid());
     BOOST_CHECK(slab.material().isValid());
     BOOST_CHECK_EQUAL(slab.thickness(), 1.5);
@@ -285,7 +285,7 @@ BOOST_AUTO_TEST_CASE(CombineSlabs) {
     BOOST_CHECK_EQUAL(slab.material().Ar(), 3.);
     BOOST_CHECK_EQUAL(
         slab.material().Z(),
-        std::exp((0.5 / 1.5) * std::log(12.) + (1. / 1.5) * std::log(6.)));
+        std::exp((1. / 1.5) * std::log(6.) + (0.5 / 1.5) * std::log(12.)));
     BOOST_CHECK_EQUAL(slab.material().molarDensity(), 4.);
   }
 }
diff --git a/Tests/UnitTests/Core/Material/InteractionsTests.cpp b/Tests/UnitTests/Core/Material/InteractionsTests.cpp
index 80134c43f30..0f10a263a06 100644
--- a/Tests/UnitTests/Core/Material/InteractionsTests.cpp
+++ b/Tests/UnitTests/Core/Material/InteractionsTests.cpp
@@ -146,10 +146,10 @@ BOOST_AUTO_TEST_CASE(silicon_landau) {
   const Acts::Material silicon = Acts::Test::makeSilicon();
   const auto thickness = 0.17_cm;
   const auto slab = Acts::MaterialSlab(silicon, thickness);
-  const float m = 105.7_MeV;
-  const float q = -1_e;
-  const float qOverP = q / 10_GeV;
-  const float absQ = std::abs(q);
+  const double m = 105.7_MeV;
+  const double q = -1_e;
+  const double qOverP = q / 10_GeV;
+  const double absQ = std::abs(q);
 
   // Difference is within 5% from PDG value
   const auto dE = computeEnergyLossLandau(slab, m, qOverP, absQ) / 1_MeV;