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