Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor(material)!: update slabs from float to ActsScalar #3609

Closed
wants to merge 9 commits into from
Closed

refactor(material)!: update slabs from float to ActsScalar #3609

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
10be219
function+tests
AJPfleger Sep 12, 2024
4dba75e
fix missing condition
AJPfleger Sep 12, 2024
c05aec5
outside usage
AJPfleger Sep 12, 2024
9f3907f
missed things
AJPfleger Sep 12, 2024
ec65174
switch to ActsScalar
AJPfleger Sep 12, 2024
162592b
Merge branch 'main' into material-slow
AJPfleger Sep 12, 2024
8514f74
Merge branch 'main' into material-slow
AJPfleger Oct 2, 2024
29c86d3
Merge branch 'main' into material-slow
AJPfleger Oct 8, 2024
e62ddeb
Merge branch 'main' into material-slow
AJPfleger Nov 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions Core/include/Acts/Material/MaterialSlab.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,12 +45,12 @@ class MaterialSlab {
/// Construct vacuum without thickness.
MaterialSlab() = default;
/// Construct vacuum with thickness.
explicit MaterialSlab(float thickness);
explicit MaterialSlab(ActsScalar 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, ActsScalar thickness);
~MaterialSlab() = default;

MaterialSlab(MaterialSlab&&) = default;
Expand All @@ -59,25 +59,25 @@ class MaterialSlab {
MaterialSlab& operator=(const MaterialSlab&) = default;

/// Scale the material thickness by the given factor.
void scaleThickness(float scale);
void scaleThickness(ActsScalar 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() && (0. < m_thickness); }

/// Access the (average) material parameters.
constexpr const Material& material() const { return m_material; }
/// Return the thickness.
constexpr float thickness() const { return m_thickness; }
constexpr ActsScalar thickness() const { return m_thickness; }
/// Return the radiation length fraction.
constexpr float thicknessInX0() const { return m_thicknessInX0; }
constexpr ActsScalar thicknessInX0() const { return m_thicknessInX0; }
/// Return the nuclear interaction length fraction.
constexpr float thicknessInL0() const { return m_thicknessInL0; }
constexpr ActsScalar thicknessInL0() const { return m_thicknessInL0; }

private:
Material m_material;
float m_thickness = 0.0f;
float m_thicknessInX0 = 0.0f;
float m_thicknessInL0 = 0.0f;
ActsScalar m_thickness = 0.;
ActsScalar m_thicknessInX0 = 0.;
ActsScalar m_thicknessInL0 = 0.;

friend constexpr bool operator==(const MaterialSlab& lhs,
const MaterialSlab& rhs) {
Expand Down
54 changes: 25 additions & 29 deletions Core/src/Material/AverageMaterials.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,40 +28,37 @@ 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<double>(slab1.thickness());
const double thickness2 = static_cast<double>(slab2.thickness());
const ActsScalar thickness1 = slab1.thickness();
const ActsScalar thickness2 = slab2.thickness();

// the thickness properties are purely additive
const double thickness = thickness1 + thickness2;
const ActsScalar thickness = thickness1 + thickness2;

// if the two materials are the same there is no need for additional
// computation
if (mat1 == mat2) {
return {mat1, static_cast<float>(thickness)};
return {mat1, thickness};
}

// molar amount-of-substance assuming a unit area, i.e. volume = thickness*1*1
const double molarDensity1 = static_cast<double>(mat1.molarDensity());
const double molarDensity2 = static_cast<double>(mat2.molarDensity());
// use ActsScalar for (intermediate) computations to avoid precision loss
const ActsScalar molarDensity1 = static_cast<ActsScalar>(mat1.molarDensity());
const ActsScalar molarDensity2 = static_cast<ActsScalar>(mat2.molarDensity());

const double molarAmount1 = molarDensity1 * thickness1;
const double molarAmount2 = molarDensity2 * thickness2;
const double molarAmount = molarAmount1 + molarAmount2;
const ActsScalar molarAmount1 = molarDensity1 * thickness1;
const ActsScalar molarAmount2 = molarDensity2 * thickness2;
const ActsScalar molarAmount = molarAmount1 + molarAmount2;

// handle vacuum specially
if (!(0.0 < molarAmount)) {
return {Material(), static_cast<float>(thickness)};
return {Material(), thickness};
}

// radiation/interaction length follows from consistency argument
const double thicknessX01 = static_cast<double>(slab1.thicknessInX0());
const double thicknessX02 = static_cast<double>(slab2.thicknessInX0());
const double thicknessL01 = static_cast<double>(slab1.thicknessInL0());
const double thicknessL02 = static_cast<double>(slab2.thicknessInL0());

const double thicknessInX0 = thicknessX01 + thicknessX02;
const double thicknessInL0 = thicknessL01 + thicknessL02;
const ActsScalar thicknessInX0 =
slab1.thicknessInX0() + slab2.thicknessInX0();
const ActsScalar thicknessInL0 =
slab1.thicknessInL0() + slab2.thicknessInL0();

const float x0 = static_cast<float>(thickness / thicknessInX0);
const float l0 = static_cast<float>(thickness / thicknessInL0);
Expand All @@ -83,11 +80,11 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
// = (Vi*rhoi) / (V1*rho1 + V2*rho2)
//
// which can be computed from the molar amount-of-substance above.
const double ar1 = static_cast<double>(mat1.Ar());
const double ar2 = static_cast<double>(mat2.Ar());
const ActsScalar ar1 = static_cast<ActsScalar>(mat1.Ar());
const ActsScalar ar2 = static_cast<ActsScalar>(mat2.Ar());

const double molarWeight1 = molarAmount1 / molarAmount;
const double molarWeight2 = molarAmount2 / molarAmount;
const ActsScalar molarWeight1 = molarAmount1 / molarAmount;
const ActsScalar molarWeight2 = molarAmount2 / molarAmount;
const float ar = static_cast<float>(molarWeight1 * ar1 + molarWeight2 * ar2);

// In the case of the atomic number, its main use is the computation
Expand All @@ -103,13 +100,13 @@ 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<double>(mat1.Z());
const double z2 = static_cast<double>(mat2.Z());
const ActsScalar z1 = static_cast<ActsScalar>(mat1.Z());
const ActsScalar z2 = static_cast<ActsScalar>(mat2.Z());

const double thicknessWeight1 = thickness1 / thickness;
const double thicknessWeight2 = thickness2 / thickness;
const ActsScalar thicknessWeight1 = thickness1 / thickness;
const ActsScalar thicknessWeight2 = thickness2 / thickness;
float z = 0.f;
if (z1 > 0. && z2 > 0.) {
if (z1 > 0. && z2 > 0. && thicknessWeight1 > 0. && thicknessWeight2 > 0.) {
z = static_cast<float>(std::exp(thicknessWeight1 * std::log(z1) +
thicknessWeight2 * std::log(z2)));
} else {
Expand All @@ -120,6 +117,5 @@ Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
// the total volume for the same unit area, i.e. volume = totalThickness*1*1
const float molarDensity = static_cast<float>(molarAmount / thickness);

return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity),
static_cast<float>(thickness)};
return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity), thickness};
}
22 changes: 12 additions & 10 deletions Core/src/Material/Interactions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,8 +102,10 @@ 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,
inline float computeEpsilon(float molarElectronDensity,
Acts::ActsScalar thickness_,
const RelativisticQuantities& rq) {
const float thickness = static_cast<float>(thickness_);
return 0.5f * K * molarElectronDensity * thickness * rq.q2OverBeta2;
}

Expand Down Expand Up @@ -160,7 +162,7 @@ inline float computeEnergyLossLandauFwhm(const Acts::MaterialSlab& slab,
}

const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const Acts::ActsScalar thickness = slab.thickness();
// the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
return 4 * computeEpsilon(Ne, thickness, rq);
}
Expand All @@ -179,7 +181,7 @@ float Acts::computeEnergyLossBethe(const MaterialSlab& slab, float m,
const RelativisticQuantities rq{m, qOverP, absQ};
const float I = slab.material().meanExcitationEnergy();
const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const ActsScalar thickness = slab.thickness();
const float eps = computeEpsilon(Ne, thickness, rq);
const float dhalf = computeDeltaHalf(I, Ne, rq);
const float u = computeMassTerm(Me, rq);
Expand All @@ -204,7 +206,7 @@ float Acts::deriveEnergyLossBetheQOverP(const MaterialSlab& slab, float m,
const RelativisticQuantities rq{m, qOverP, absQ};
const float I = slab.material().meanExcitationEnergy();
const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const ActsScalar thickness = slab.thickness();
const float eps = computeEpsilon(Ne, thickness, rq);
const float dhalf = computeDeltaHalf(I, Ne, rq);
const float u = computeMassTerm(Me, rq);
Expand Down Expand Up @@ -241,7 +243,7 @@ float Acts::computeEnergyLossLandau(const MaterialSlab& slab, float m,
const RelativisticQuantities rq{m, qOverP, absQ};
const float I = slab.material().meanExcitationEnergy();
const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const ActsScalar thickness = slab.thickness();
const float eps = computeEpsilon(Ne, thickness, rq);
const float dhalf = computeDeltaHalf(I, Ne, rq);
const float u = computeMassTerm(Me, rq);
Expand All @@ -261,7 +263,7 @@ float Acts::deriveEnergyLossLandauQOverP(const MaterialSlab& slab, float m,
const RelativisticQuantities rq{m, qOverP, absQ};
const float I = slab.material().meanExcitationEnergy();
const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const ActsScalar thickness = slab.thickness();
const float eps = computeEpsilon(Ne, thickness, rq);
const float dhalf = computeDeltaHalf(I, Ne, rq);
const float t = computeMassTerm(Me, rq);
Expand Down Expand Up @@ -295,7 +297,7 @@ float Acts::computeEnergyLossLandauSigma(const MaterialSlab& slab, float m,

const RelativisticQuantities rq{m, qOverP, absQ};
const float Ne = slab.material().molarElectronDensity();
const float thickness = slab.thickness();
const ActsScalar thickness = slab.thickness();
// the Landau-Vavilov fwhm is 4*eps (see RPP2018 fig. 33.7)
const float fwhm = 4 * computeEpsilon(Ne, thickness, rq);
return convertLandauFwhmToGaussianSigma(fwhm);
Expand Down Expand Up @@ -392,7 +394,7 @@ float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
}

// relative radiation length
const float x = slab.thicknessInX0();
const float x = static_cast<float>(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;
Expand Down Expand Up @@ -421,7 +423,7 @@ float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
}

// relative radiation length
const float x = slab.thicknessInX0();
const float x = static_cast<float>(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;
Expand Down Expand Up @@ -510,7 +512,7 @@ float Acts::computeMultipleScatteringTheta0(const MaterialSlab& slab,
}

// relative radiation length
const float xOverX0 = slab.thicknessInX0();
const float xOverX0 = static_cast<float>(slab.thicknessInX0());
// 1/p = q/(pq) = (q/p)/q
const float momentumInv = std::abs(qOverP / absQ);
// q²/beta²; a smart compiler should be able to remove the unused computations
Expand Down
17 changes: 9 additions & 8 deletions Core/src/Material/MaterialSlab.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,18 @@
namespace Acts {

namespace {
static constexpr auto eps = 2 * std::numeric_limits<float>::epsilon();
static constexpr auto eps = 2 * std::numeric_limits<ActsScalar>::epsilon();
}

MaterialSlab::MaterialSlab(float thickness) : m_thickness(thickness) {}
MaterialSlab::MaterialSlab(ActsScalar thickness) : m_thickness(thickness) {}

MaterialSlab::MaterialSlab(const Material& material, float thickness)
MaterialSlab::MaterialSlab(const Material& material, ActsScalar thickness)
: m_material(material),
m_thickness(thickness),
m_thicknessInX0((eps < material.X0()) ? (thickness / material.X0()) : 0),
m_thicknessInL0((eps < material.L0()) ? (thickness / material.L0()) : 0) {
if (thickness < 0) {
m_thicknessInX0((eps < material.X0()) ? (thickness / material.X0()) : 0.),
m_thicknessInL0((eps < material.L0()) ? (thickness / material.L0())
: 0.) {
if (thickness < 0.) {
throw std::runtime_error("thickness < 0");
}
}
Expand All @@ -55,8 +56,8 @@ MaterialSlab MaterialSlab::averageLayers(
return result;
}

void MaterialSlab::scaleThickness(float scale) {
if (scale < 0) {
void MaterialSlab::scaleThickness(ActsScalar scale) {
if (scale < 0.) {
throw std::runtime_error("scale < 0");
}

Expand Down
39 changes: 20 additions & 19 deletions Tests/UnitTests/Core/Material/AccumulatedMaterialSlabTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@ using Acts::MaterialSlab;
using Acts::Test::makeSilicon;
using Acts::Test::makeUnitSlab;

constexpr auto eps = std::numeric_limits<float>::epsilon();
constexpr auto epsF = std::numeric_limits<float>::epsilon();
constexpr auto eps = std::numeric_limits<Acts::ActsScalar>::epsilon();

} // namespace

Expand Down Expand Up @@ -74,8 +75,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.);
BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.);
}
// accumulate three identical steps for one additional track
{
Expand All @@ -88,8 +89,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.);
BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.);
}
}

Expand All @@ -108,8 +109,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.);
BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.);
}
// accumulate one step with thickness 1, one with thickness 2
{
Expand All @@ -123,8 +124,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.);
BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.);
}
// accumulate one step with thickness 3
{
Expand All @@ -137,8 +138,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.);
BOOST_CHECK_EQUAL(average.thicknessInL0(), 3.);
}
}

Expand Down Expand Up @@ -182,19 +183,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().X0(), 2 * unit.material().X0(), epsF);
CHECK_CLOSE_REL(average.material().L0(), 2 * unit.material().L0(), epsF);
CHECK_CLOSE_REL(average.material().molarDensity(),
0.5f * unit.material().molarDensity(), eps);
0.5f * unit.material().molarDensity(), epsF);
// average atom is still the same species
CHECK_CLOSE_REL(average.material().Ar(), unit.material().Ar(), eps);
CHECK_CLOSE_REL(average.material().Ar(), unit.material().Ar(), epsF);
// average atomic number proportional to the thickness
CHECK_CLOSE_REL(average.material().Z(), 0.5 * unit.material().Z(), eps);
CHECK_CLOSE_REL(average.material().Z(), 0.5f * unit.material().Z(), epsF);
// 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);
}
}

Expand Down
Loading
Loading