cpu/material/materials.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Project include(s)
0010 #include "detray/definitions/units.hpp"
0011 #include "detray/geometry/mask.hpp"
0012 #include "detray/geometry/shapes/line.hpp"
0013 #include "detray/geometry/surface_descriptor.hpp"
0014 #include "detray/material/concepts.hpp"
0015 #include "detray/material/material.hpp"
0016 #include "detray/material/material_rod.hpp"
0017 #include "detray/material/material_slab.hpp"
0018 #include "detray/material/mixture.hpp"
0019 #include "detray/material/predefined_materials.hpp"
0020 #include "detray/navigation/intersection/ray_intersector.hpp"
0021
0022 // Detray test include(s)
0023 #include "detray/test/framework/types.hpp"
0024
0025 // GTest include(s)
0026 #include <gtest/gtest.h>
0027
0028 // System include(s)
0029 #include <limits>
0030
0031 using namespace detray;
0032
0033 using test_algebra = test::algebra;
0034 using scalar = test::scalar;
0035 using point2 = test::point2;
0036 using point3 = test::point3;
0037 using transform3 = test::transform3;
0038 using vector3 = test::vector3;
0039
0040 namespace {
0041
0042 constexpr scalar tol{1e-7f};
0043 constexpr auto max_val{std::numeric_limits<scalar>::max()};
0044
0045 }  // anonymous namespace
0046
0047 // This tests the material functionalities
0048 GTEST_TEST(detray_material, materials) {
0049   // vacuum
0050   constexpr vacuum<scalar> vac;
0051
0052   static_assert(concepts::material_params<decltype(vac)>);
0053   static_assert(concepts::homogeneous_material<decltype(vac)>);
0054   static_assert(!concepts::surface_material<decltype(vac)>);
0055   static_assert(concepts::volume_material<decltype(vac)>);
0056   static_assert(!concepts::material_map<decltype(vac)>);
0057
0058   EXPECT_EQ(vac.X0(), max_val);
0059   EXPECT_EQ(vac.L0(), max_val);
0060   EXPECT_EQ(vac.Ar(), scalar{0});
0061   EXPECT_EQ(vac.Z(), scalar{0});
0062   EXPECT_EQ(vac.molar_density(), scalar{0});
0063   EXPECT_EQ(vac.molar_electron_density(), scalar{0});
0064
0065   // beryllium
0066   constexpr beryllium_tml<scalar> beryll;
0067
0068   static_assert(concepts::material_params<decltype(beryll)>);
0069   static_assert(concepts::homogeneous_material<decltype(beryll)>);
0070   static_assert(!concepts::surface_material<decltype(beryll)>);
0071   static_assert(concepts::volume_material<decltype(beryll)>);
0072   static_assert(!concepts::material_map<decltype(beryll)>);
0073
0074   EXPECT_NEAR(beryll.X0(), static_cast<scalar>(352.8 * unit<double>::mm), 1e-4);
0075   EXPECT_NEAR(beryll.L0(), static_cast<scalar>(407.0 * unit<double>::mm), tol);
0076   EXPECT_NEAR(beryll.Ar(), 9.012f, tol);
0077   EXPECT_NEAR(beryll.Z(), 4.0f, tol);
0078
0079   // @note molar density is obtained by the following equation:
0080   // molar_density = mass_density [GeV/c²] / molar_mass [GeV/c²] * mol / cm³
0081   EXPECT_NEAR(beryll.molar_density(),
0082               static_cast<scalar>(1.848 / beryllium_tml<double>().Ar() *
0083                                   unit<double>::mol / unit<double>::cm3),
0084               tol);
0085
0086   // silicon
0087   constexpr silicon_tml<scalar> silicon;
0088
0089   static_assert(concepts::material_params<decltype(silicon)>);
0090   static_assert(concepts::homogeneous_material<decltype(silicon)>);
0091   static_assert(!concepts::surface_material<decltype(silicon)>);
0092   static_assert(concepts::volume_material<decltype(silicon)>);
0093   static_assert(!concepts::material_map<decltype(silicon)>);
0094
0095   EXPECT_NEAR(silicon.X0(), static_cast<scalar>(95.7 * unit<double>::mm), 1e-5);
0096   EXPECT_NEAR(silicon.L0(), static_cast<scalar>(465.2 * unit<double>::mm),
0097               1e-4);
0098   EXPECT_NEAR(silicon.Ar(), 28.03f, tol);
0099   EXPECT_NEAR(silicon.Z(), 14.f, tol);
0100   EXPECT_NEAR(silicon.molar_density(),
0101               static_cast<scalar>(2.32 / silicon_tml<double>().Ar() *
0102                                   unit<double>::mol / unit<double>::cm3),
0103               tol);
0104 }
0105
0106 GTEST_TEST(detray_material, mixture) {
0107   // Check if material property doesn't change after mixing with other
0108   // material of 0 ratio
0109   constexpr oxygen_gas<scalar> pure_oxygen;
0110   constexpr mixture<scalar, oxygen_gas<scalar>,
0111                     aluminium<scalar, std::ratio<0, 1>>>
0112       oxygen_mix;
0113
0114   static_assert(concepts::material_params<decltype(oxygen_mix)>);
0115   static_assert(concepts::homogeneous_material<decltype(oxygen_mix)>);
0116   static_assert(!concepts::surface_material<decltype(oxygen_mix)>);
0117   static_assert(concepts::volume_material<decltype(oxygen_mix)>);
0118   static_assert(!concepts::material_map<decltype(oxygen_mix)>);
0119
0120   EXPECT_NEAR(pure_oxygen.X0(), oxygen_mix.X0(), 0.02f);
0121   EXPECT_NEAR(pure_oxygen.L0(), oxygen_mix.L0(), tol);
0122   EXPECT_NEAR(pure_oxygen.Ar(), oxygen_mix.Ar(), tol);
0123   EXPECT_NEAR(pure_oxygen.Z(), oxygen_mix.Z(), tol);
0124   EXPECT_NEAR(pure_oxygen.mass_density(), oxygen_mix.mass_density(), tol);
0125   EXPECT_NEAR(pure_oxygen.molar_density(), oxygen_mix.molar_density(), tol);
0126
0127   // Air mixture check
0128   constexpr mixture<scalar, carbon_gas<scalar, std::ratio<0, 100>>,
0129                     nitrogen_gas<scalar, std::ratio<76, 100>>,
0130                     oxygen_gas<scalar, std::ratio<23, 100>>,
0131                     argon_gas<scalar, std::ratio<1, 100>>>
0132       air_mix;
0133   constexpr air<scalar> pure_air;
0134
0135   static_assert(concepts::material_params<decltype(pure_air)>);
0136   static_assert(concepts::homogeneous_material<decltype(pure_air)>);
0137   static_assert(!concepts::surface_material<decltype(pure_air)>);
0138   static_assert(concepts::volume_material<decltype(pure_air)>);
0139   static_assert(!concepts::material_map<decltype(pure_air)>);
0140
0141   EXPECT_TRUE(std::abs(air_mix.X0() - pure_air.X0()) / pure_air.X0() < 0.01f);
0142   EXPECT_TRUE(std::abs(air_mix.L0() - pure_air.L0()) / pure_air.L0() < 0.01f);
0143   EXPECT_TRUE(std::abs(air_mix.Ar() - pure_air.Ar()) / pure_air.Ar() < 0.01f);
0144   EXPECT_TRUE(std::abs(air_mix.Z() - pure_air.Z()) / pure_air.Z() < 0.01f);
0145   EXPECT_TRUE(std::abs(air_mix.mass_density() - pure_air.mass_density()) /
0146                   pure_air.mass_density() <
0147               0.01f);
0148
0149   // Vector check
0150   material_slab<scalar> slab1(air_mix, 5.5f);
0151   material_slab<scalar> slab2(pure_air, 2.3f);
0152   material_slab<scalar> slab3(oxygen_gas<scalar>(), 2.f);
0153
0154   std::vector<material_slab<scalar>> slab_vec;
0155
0156   slab_vec.push_back(slab1);
0157   slab_vec.push_back(slab2);
0158   slab_vec.push_back(slab3);
0159
0160   EXPECT_NEAR(slab_vec[0].thickness_in_X0(),
0161               slab1.thickness() / slab1.get_material().X0(), tol);
0162   EXPECT_NEAR(slab_vec[1].thickness_in_X0(),
0163               slab2.thickness() / slab2.get_material().X0(), tol);
0164   EXPECT_NEAR(slab_vec[2].thickness_in_X0(),
0165               slab3.thickness() / slab3.get_material().X0(), tol);
0166 }
0167
0168 GTEST_TEST(detray_material, mixture2) {
0169   constexpr mixture<scalar, silicon_with_ded<scalar, std::ratio<1, 3>>,
0170                     cesium_iodide_with_ded<scalar, std::ratio<2, 3>>>
0171       mix;
0172
0173   static_assert(concepts::material_params<decltype(mix)>);
0174   static_assert(concepts::homogeneous_material<decltype(mix)>);
0175   static_assert(!concepts::surface_material<decltype(mix)>);
0176   static_assert(concepts::volume_material<decltype(mix)>);
0177   static_assert(!concepts::material_map<decltype(mix)>);
0178
0179   // For the moment, the mixture is not supposed to have density effect data
0180   // in any case
0181   EXPECT_EQ(mix.has_density_effect_data(), false);
0182   EXPECT_EQ(mix.Ar(), silicon<scalar>().Ar() * 1.f / 3.f +
0183                           cesium_iodide<scalar>().Ar() * 2.f / 3.f);
0184   EXPECT_EQ(mix.Z(), silicon<scalar>().Z() * 1.f / 3.f +
0185                          cesium_iodide<scalar>().Z() * 2.f / 3.f);
0186 }
0187
0188 // This tests the material slab functionalities
0189 GTEST_TEST(detray_material, material_slab) {
0190   constexpr material_slab<scalar> slab(oxygen_gas<scalar>(),
0191                                        2.f * unit<scalar>::mm);
0192
0193   static_assert(concepts::material_slab<decltype(slab)>);
0194   static_assert(concepts::homogeneous_material<decltype(slab)>);
0195   static_assert(concepts::surface_material<decltype(slab)>);
0196   static_assert(!concepts::volume_material<decltype(slab)>);
0197   static_assert(!concepts::material_map<decltype(slab)>);
0198
0199   const scalar cos_inc_ang{0.3f};
0200
0201   EXPECT_NEAR(slab.path_segment(cos_inc_ang), 2.f * unit<scalar>::mm / 0.3f,
0202               tol);
0203   EXPECT_NEAR(slab.path_segment_in_X0(cos_inc_ang),
0204               slab.path_segment(cos_inc_ang) / slab.get_material().X0(), tol);
0205   EXPECT_NEAR(slab.path_segment_in_L0(cos_inc_ang),
0206               slab.path_segment(cos_inc_ang) / slab.get_material().L0(), tol);
0207 }
0208
0209 // This tests the material rod functionalities
0210 GTEST_TEST(detray_material, material_rod) {
0211   // Rod with 1 mm radius
0212   constexpr material_rod<scalar> rod(oxygen_gas<scalar>(),
0213                                      1.f * unit<scalar>::mm);
0214
0215   static_assert(concepts::material_rod<decltype(rod)>);
0216   static_assert(concepts::homogeneous_material<decltype(rod)>);
0217   static_assert(concepts::surface_material<decltype(rod)>);
0218   static_assert(!concepts::volume_material<decltype(rod)>);
0219   static_assert(!concepts::material_map<decltype(rod)>);
0220
0221   // tf3 with Identity rotation and no translation
0222   const vector3 x{1.f, 0.f, 0.f};
0223   const vector3 z{0.f, 0.f, 1.f};
0224   const vector3 t{0.f, 0.f, 0.f};
0225   const transform3 tf{t, vector::normalize(z), vector::normalize(x)};
0226
0227   // Create a track
0228   const point3 pos{-1.f / 6.f, -10.f, 0.f};
0229   const vector3 dir{0.f, 1.f, 3.f};
0230   const free_track_parameters<test_algebra> trk(pos, 0.f, dir, -1.f);
0231
0232   // Infinite wire with 1 mm radial cell size
0233   const mask<line_circular, test_algebra> ln{
0234       0u, 1.f * unit<scalar>::mm, std::numeric_limits<scalar>::infinity()};
0235
0236   auto is = ray_intersector<line_circular, test_algebra, true>{}(
0237       detail::ray<test_algebra>(trk), surface_descriptor<>{}, ln, tf);
0238
0239   const scalar cos_inc_ang{std::abs(vector::dot(
0240       line2D<test_algebra>::normal(tf, is.local()), vector::normalize(dir)))};
0241   const scalar approach{is.local()[0]};
0242
0243   EXPECT_NEAR(rod.path_segment(cos_inc_ang, approach),
0244               2.f * std::sqrt(10.f - 10.f / 36.f), 1e-5f);
0245   EXPECT_NEAR(rod.path_segment_in_X0(cos_inc_ang, approach),
0246               rod.path_segment(cos_inc_ang, approach) / rod.get_material().X0(),
0247               tol);
0248   EXPECT_NEAR(rod.path_segment_in_L0(cos_inc_ang, approach),
0249               rod.path_segment(cos_inc_ang, approach) / rod.get_material().L0(),
0250               tol);
0251 }