feat: Parameterize mean excitation energy in Material

Core/include/Acts/Definitions/Units.hpp

@@ -195,6 +195,8 @@ constexpr double Gauss = 1e-4 * T;
 constexpr double kGauss = 1e-1 * T;
 /// Amount of substance, native unit mol
 constexpr double mol = 1.0;
+/// Avogadro constant
+constexpr double Avogadro = 6.02214076e23 / mol;
 }  // namespace UnitConstants
 
 namespace UnitLiterals {
@@ -245,6 +247,7 @@ ACTS_DEFINE_UNIT_LITERAL(T)
 ACTS_DEFINE_UNIT_LITERAL(Gauss)
 ACTS_DEFINE_UNIT_LITERAL(kGauss)
 ACTS_DEFINE_UNIT_LITERAL(mol)
+ACTS_DEFINE_UNIT_LITERAL(Avogadro)
 // not needed anymore. undef to prevent littering the namespace
 #undef ACTS_DEFINE_UNIT_LITERAL
 }  // namespace UnitLiterals

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
@@ -12,7 +12,7 @@
 
 #include <iosfwd>
 #include <limits>
-#include <utility>
+#include <optional>
 
 namespace Acts {
 
@@ -26,6 +26,7 @@ namespace Acts {
 /// - relative atomic mass Ar (unitless number)
 /// - nuclear charge number Z (elementary charge e)
 /// - molar density (native amount-of-substance unit / (native length unit)³)
+/// - mean excitation energy I (native energy unit)
 ///
 /// The parameters can be effective or average parameters e.g. when a mixture
 /// of materials is described.
@@ -55,21 +56,27 @@ class Material {
   /// @param ar is the relative atomic mass
   /// @param z is the nuclear charge number
   /// @param molarRho is the molar density
+  /// @param I is the mean excitation energy
   static Material fromMolarDensity(float x0, float l0, float ar, float z,
-                                   float molarRho);
+                                   float molarRho,
+                                   std::optional<float> I = std::nullopt);
+
   /// Construct from material parameters using the mass density.
   ///
   /// @param x0 is the radiation length
   /// @param l0 is the nuclear interaction length
   /// @param ar is the relative atomic mass
   /// @param z is the nuclear charge number
   /// @param massRho is the mass density
+  /// @param I is the mean excitation energy
   ///
   /// @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);
+                                  float massRho,
+                                  std::optional<float> I = std::nullopt);
+
   /// Construct a vacuum representation.
   Material() = default;
   /// Construct from an encoded parameters vector.
@@ -99,7 +106,7 @@ class Material {
   /// Return the mass density.
   float massDensity() const;
   /// Return the mean electron excitation energy.
-  float meanExcitationEnergy() const;
+  constexpr float meanExcitationEnergy() const { return m_I; }
 
   /// Encode the properties into an opaque parameters vector.
   ParametersVector parameters() const;
@@ -110,6 +117,7 @@ class Material {
   float m_ar = 0.0f;
   float m_z = 0.0f;
   float m_molarRho = 0.0f;
+  float m_I = 0.0f;
 
   friend constexpr bool operator==(const Material& lhs, const Material& rhs) {
     return (lhs.m_x0 == rhs.m_x0) && (lhs.m_l0 == rhs.m_l0) &&

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
@@ -22,12 +22,17 @@ enum MaterialClassificationNumberIndices {
   eMolarDensity = 4,
 };
 
-// Avogadro constant
-constexpr double kAvogadro = 6.02214076e23 / Acts::UnitConstants::mol;
+float approximateMeanExcitationEnergy(float z) {
+  using namespace Acts::UnitLiterals;
+
+  // use approximative computation as defined in ATL-SOFT-PUB-2008-003
+  return 16_eV * std::pow(z, 0.9f);
+}
 }  // namespace
 
 Acts::Material Acts::Material::fromMassDensity(float x0, float l0, float ar,
-                                               float z, float massRho) {
+                                               float z, float massRho,
+                                               std::optional<float> I) {
   using namespace Acts::UnitLiterals;
 
   Material mat;
@@ -47,18 +52,34 @@ Acts::Material Acts::Material::fromMassDensity(float x0, float l0, float ar,
   //
   // perform computations in double precision to avoid loss of precision
   const double atomicMass = static_cast<double>(ar) * 1_u;
-  mat.m_molarRho = static_cast<double>(massRho) / (atomicMass * kAvogadro);
+  mat.m_molarRho =
+      static_cast<double>(massRho) / (atomicMass * UnitConstants::Avogadro);
+
+  if (I) {
+    mat.m_I = *I;
+  } else {
+    mat.m_I = approximateMeanExcitationEnergy(z);
+  }
+
   return mat;
 }
 
 Acts::Material Acts::Material::fromMolarDensity(float x0, float l0, float ar,
-                                                float z, float molarRho) {
+                                                float z, float molarRho,
+                                                std::optional<float> I) {
   Material mat;
   mat.m_x0 = x0;
   mat.m_l0 = l0;
   mat.m_ar = ar;
   mat.m_z = z;
   mat.m_molarRho = molarRho;
+
+  if (I) {
+    mat.m_I = *I;
+  } else {
+    mat.m_I = approximateMeanExcitationEnergy(z);
+  }
+
   return mat;
 }
 
@@ -67,24 +88,19 @@ Acts::Material::Material(const ParametersVector& parameters)
       m_l0(parameters[eInteractionLength]),
       m_ar(parameters[eRelativeAtomicMass]),
       m_z(parameters[eNuclearCharge]),
-      m_molarRho(parameters[eMolarDensity]) {}
+      m_molarRho(parameters[eMolarDensity]),
+      m_I(approximateMeanExcitationEnergy(m_z)) {}
 
 float Acts::Material::massDensity() const {
   using namespace Acts::UnitLiterals;
 
   // perform computations in double precision to avoid loss of precision
   const double atomicMass = static_cast<double>(m_ar) * 1_u;
-  const double numberDensity = static_cast<double>(m_molarRho) * kAvogadro;
+  const double numberDensity =
+      static_cast<double>(m_molarRho) * UnitConstants::Avogadro;
   return atomicMass * numberDensity;
 }
 
-float Acts::Material::meanExcitationEnergy() const {
-  using namespace Acts::UnitLiterals;
-
-  // use approximative computation as defined in ATL-SOFT-PUB-2008-003
-  return 16_eV * std::pow(m_z, 0.9f);
-}
-
 Acts::Material::ParametersVector Acts::Material::parameters() const {
   ParametersVector parameters;
   parameters[eRadiationLength] = m_x0;
@@ -104,6 +120,7 @@ std::ostream& Acts::operator<<(std::ostream& os, const Material& material) {
     os << "|ar=" << material.Ar();
     os << "|z=" << material.Z();
     os << "|molar_rho=" << material.molarDensity();
+    os << "|I=" << material.meanExcitationEnergy();
   }
   return os;
 }

Examples/Python/src/Base.cpp

@@ -74,6 +74,7 @@ void addUnits(Context& ctx) {
   UNIT(Gauss)
   UNIT(kGauss)
   UNIT(mol)
+  UNIT(Avogadro)
 
 #undef UNIT
 }