refactor!: Use double for material

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

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,17 +89,17 @@ 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.
   MaterialSlab m_trackAverage;
   /// 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

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

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.
   ///

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(); }