fix!: make material validity checks and construction explicit (#3494)

By having the validity check explicit and forbid construction with type conversation, we avoid common pitfalls for the user[1].

The `AverageMaterials` had to be reworked a bit, since we were doing a lot of implicit type conversions in there. On different compilers we got different results. Now it should be stable.

Thanks also to @<!---->paulgessinger for sending me in the correct direction.

[1] The user is me, trying to debug for hours why I get always vacuum material, which could have been caught earlier, by not allowing implicit conversions.

Core/include/Acts/Material/ISurfaceMaterial.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -139,7 +139,7 @@ inline MaterialSlab ISurfaceMaterial::materialSlab(
   // The plain material properties associated to this bin
   MaterialSlab plainMatProp = materialSlab(lp);
   // Scale if you have material to scale
-  if (plainMatProp) {
+  if (plainMatProp.isValid()) {
     double scaleFactor = factor(pDir, mStage);
     if (scaleFactor == 0.) {
       return MaterialSlab();
@@ -154,7 +154,7 @@ inline MaterialSlab ISurfaceMaterial::materialSlab(
   // The plain material properties associated to this bin
   MaterialSlab plainMatProp = materialSlab(gp);
   // Scale if you have material to scale
-  if (plainMatProp) {
+  if (plainMatProp.isValid()) {
     double scaleFactor = factor(pDir, mStage);
     if (scaleFactor == 0.) {
       return MaterialSlab();

Core/include/Acts/Material/Material.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -74,7 +74,7 @@ class Material {
   /// Construct a vacuum representation.
   Material() = default;
   /// Construct from an encoded parameters vector.
-  Material(const ParametersVector& parameters);
+  explicit Material(const ParametersVector& parameters);
 
   Material(Material&& mat) = default;
   Material(const Material& mat) = default;
@@ -83,9 +83,9 @@ class Material {
   Material& operator=(const Material& mat) = default;
 
   /// Check if the material is valid, i.e. it is not vacuum.
-  constexpr operator bool() const { return 0.0f < m_ar; }
+  bool isValid() const { return 0.0f < m_ar; }
 
-  /// Return the radition length. Infinity in case of vacuum.
+  /// Return the radiation length. Infinity in case of vacuum.
   constexpr float X0() const { return m_x0; }
   /// Return the nuclear interaction length. Infinity in case of vacuum.
   constexpr float L0() const { return m_l0; }

Core/include/Acts/Material/MaterialSlab.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -45,7 +45,7 @@ class MaterialSlab {
   /// Construct vacuum without thickness.
   MaterialSlab() = default;
   /// Construct vacuum with thickness.
-  MaterialSlab(float thickness);
+  explicit MaterialSlab(float thickness);
   /// Construct from material description.
   ///
   /// @param material  is the material description
@@ -62,7 +62,7 @@ class MaterialSlab {
   void scaleThickness(float scale);
 
   /// Check if the material is valid, i.e. it is finite and not vacuum.
-  constexpr operator bool() const { return m_material && (0.0f < m_thickness); }
+  bool isValid() const { return m_material.isValid() && (0.0f < m_thickness); }
 
   /// Access the (average) material parameters.
   constexpr const Material& material() const { return m_material; }

Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -50,7 +50,7 @@ struct PointwiseMaterialInteraction {
   const Direction navDir;
 
   /// The effective, passed material properties including the path correction.
-  MaterialSlab slab;
+  MaterialSlab slab = MaterialSlab(0.);
   /// The path correction factor due to non-zero incidence on the surface.
   double pathCorrection = 0.;
   /// Expected phi variance due to the interactions.
@@ -89,6 +89,7 @@ struct PointwiseMaterialInteraction {
         navDir(state.options.direction) {}
 
   /// @brief This function evaluates the material properties to interact with
+  /// This updates the slab and then returns, if the resulting slab is valid
   ///
   /// @tparam propagator_state_t Type of the propagator state
   /// @tparam navigator_t Type of the navigator
@@ -119,8 +120,8 @@ struct PointwiseMaterialInteraction {
     pathCorrection = surface->pathCorrection(state.geoContext, pos, dir);
     slab.scaleThickness(pathCorrection);
 
-    // Get the surface material & properties from them
-    return slab;
+    // Check if the evaluated material is valid
+    return slab.isValid();
   }
 
   /// @brief This function evaluate the material effects

Core/include/Acts/Propagator/detail/VolumeMaterialInteraction.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -119,7 +119,7 @@ struct VolumeMaterialInteraction {
     } else {
       slab = MaterialSlab();
     }
-    return slab;
+    return slab.isValid();
   }
 };
 

Core/src/Material/AverageMaterials.cpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2020 CERN for the benefit of the Acts project
+// Copyright (C) 2020-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -29,41 +29,42 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
   // 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
-  double thickness = static_cast<double>(slab1.thickness()) +
-                     static_cast<double>(slab2.thickness());
+  const double thickness = thickness1 + thickness2;
 
   // if the two materials are the same there is no need for additional
   // computation
   if (mat1 == mat2) {
     return {mat1, static_cast<float>(thickness)};
   }
 
-  double thicknessInX0 = static_cast<double>(slab1.thicknessInX0()) +
-                         static_cast<double>(slab2.thicknessInX0());
-  double thicknessInL0 = static_cast<double>(slab1.thicknessInL0()) +
-                         static_cast<double>(slab2.thicknessInL0());
-
-  // radiation/interaction length follows from consistency argument
-  float x0 = thickness / thicknessInX0;
-  float l0 = thickness / thicknessInL0;
-
   // molar amount-of-substance assuming a unit area, i.e. volume = thickness*1*1
-  double molarAmount1 = static_cast<double>(mat1.molarDensity()) *
-                        static_cast<double>(slab1.thickness());
-  double molarAmount2 = static_cast<double>(mat2.molarDensity()) *
-                        static_cast<double>(slab2.thickness());
-  double molarAmount = molarAmount1 + molarAmount2;
+  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 molarAmount = molarAmount1 + molarAmount2;
 
   // handle vacuum specially
   if (!(0.0 < molarAmount)) {
     return {Material(), static_cast<float>(thickness)};
   }
 
-  // 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
-  float molarDensity = molarAmount / 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 thicknessInX0 = thicknessX01 + thicknessX02;
+  const double thicknessInL0 = thicknessL01 + thicknessL02;
+
+  const float x0 = static_cast<float>(thickness / thicknessInX0);
+  const float l0 = static_cast<float>(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
@@ -82,10 +83,12 @@ 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());
 
-  double molarWeight1 = molarAmount1 / molarAmount;
-  double molarWeight2 = molarAmount2 / molarAmount;
-  float ar = molarWeight1 * mat1.Ar() + molarWeight2 * mat2.Ar();
+  const double molarWeight1 = molarAmount1 / molarAmount;
+  const double molarWeight2 = molarAmount2 / molarAmount;
+  const float ar = static_cast<float>(molarWeight1 * ar1 + molarWeight2 * ar2);
 
   // 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 :
@@ -100,17 +103,23 @@ 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 )
-
-  double thicknessWeight1 = slab1.thickness() / thickness;
-  double thicknessWeight2 = slab2.thickness() / thickness;
-  float z = 0;
-  if (mat1.Z() != 0 && mat2.Z() != 0) {
-    z = exp(thicknessWeight1 * log(mat1.Z()) +
-            thicknessWeight2 * log(mat2.Z()));
+  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)));
   } else {
-    z = thicknessWeight1 * mat1.Z() + thicknessWeight2 * mat2.Z();
+    z = static_cast<float>(thicknessWeight1 * z1 + thicknessWeight2 * z2);
   }
 
+  // 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);
+
   return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity),
           static_cast<float>(thickness)};
 }

Core/src/Material/Interactions.cpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -154,7 +154,7 @@ namespace detail {
 inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
                                          const RelativisticQuantities& rq) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -171,7 +171,7 @@ inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
 float Acts::computeEnergyLossBethe(const MaterialSlab& slab, float m,
                                    float qOverP, float absQ) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -196,7 +196,7 @@ float Acts::computeEnergyLossBethe(const MaterialSlab& slab, float m,
 float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
                                         float qOverP, float absQ) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -233,7 +233,7 @@ float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
 float Acts::computeEnergyLossLandau(const MaterialSlab& slab, float m,
                                     float qOverP, float absQ) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -253,7 +253,7 @@ float Acts::computeEnergyLossLandau(const MaterialSlab& slab, float m,
 float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
                                          float qOverP, float absQ) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -288,7 +288,7 @@ float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
 float Acts::computeEnergyLossLandauSigma(const MaterialSlab& slab, float m,
                                          float qOverP, float absQ) {
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -386,7 +386,7 @@ float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -415,7 +415,7 @@ float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 
@@ -504,7 +504,7 @@ float Acts::computeMultipleScatteringTheta0(const MaterialSlab& slab,
          "pdg is not absolute");
 
   // return early in case of vacuum or zero thickness
-  if (!slab) {
+  if (!slab.isValid()) {
     return 0.0f;
   }
 

Core/src/Material/Material.cpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -96,7 +96,7 @@ Acts::Material::ParametersVector Acts::Material::parameters() const {
 }
 
 std::ostream& Acts::operator<<(std::ostream& os, const Material& material) {
-  if (!material) {
+  if (!material.isValid()) {
     os << "vacuum";
   } else {
     os << "x0=" << material.X0();

Examples/Io/Root/src/RootMaterialWriter.cpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2017-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -297,7 +297,7 @@ void ActsExamples::RootMaterialWriter::writeMaterial(
         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
         for (std::size_t point = 0; point < points; point++) {
           auto mat = Acts::Material(grid.at(point));
-          if (mat) {
+          if (mat.isValid()) {
             x0->SetBinContent(point + 1, mat.X0());
             l0->SetBinContent(point + 1, mat.L0());
             A->SetBinContent(point + 1, mat.Ar());
@@ -311,7 +311,7 @@ void ActsExamples::RootMaterialWriter::writeMaterial(
         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
         for (std::size_t point = 0; point < points; point++) {
           auto mat = Acts::Material(grid.at(point));
-          if (mat) {
+          if (mat.isValid()) {
             x0->SetBinContent(point + 1, mat.X0());
             l0->SetBinContent(point + 1, mat.L0());
             A->SetBinContent(point + 1, mat.Ar());