cpu/material/material_interaction.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/pdg_particle.hpp"
0011 #include "detray/definitions/units.hpp"
0012 #include "detray/geometry/mask.hpp"
0013 #include "detray/geometry/shapes/unbounded.hpp"
0014 #include "detray/material/interaction.hpp"
0015 #include "detray/material/material.hpp"
0016 #include "detray/material/material_slab.hpp"
0017 #include "detray/material/predefined_materials.hpp"
0018 #include "detray/navigation/caching_navigator.hpp"
0019 #include "detray/propagator/actors.hpp"
0020 #include "detray/propagator/actors/parameter_updater.hpp"
0021 #include "detray/propagator/line_stepper.hpp"
0022 #include "detray/propagator/propagator.hpp"
0023 #include "detray/propagator/rk_stepper.hpp"
0024 #include "detray/test/common/bfield.hpp"
0025
0026 // Detray test include(s)
0027 #include "detray/test/common/build_telescope_detector.hpp"
0028 #include "detray/test/framework/types.hpp"
0029 #include "detray/test/utils/inspectors.hpp"
0030 #include "detray/test/utils/random_scatterer.hpp"
0031 #include "detray/test/utils/statistics.hpp"
0032
0033 // VecMem include(s).
0034 #include <vecmem/memory/host_memory_resource.hpp>
0035
0036 // GTest include(s).
0037 #include <gtest/gtest.h>
0038
0039 using namespace detray;
0040
0041 using test_algebra = test::algebra;
0042 using scalar = test::scalar;
0043 using covariance_t =
0044     typename bound_track_parameters<test_algebra>::covariance_type;
0045 using interactor_t = actor::pointwise_material_interactor<test_algebra>;
0046
0047 static_assert(detray::concepts::actor<interactor_t>);
0048
0049 // Test is done for muon
0050 namespace {
0051 pdg_particle ptc = muon<scalar>();
0052 }
0053
0054 // Material interaction test with telescope Geometry
0055 GTEST_TEST(detray_material, telescope_geometry_energy_loss) {
0056   vecmem::host_memory_resource host_mr;
0057
0058   // Build in x-direction from given module positions
0059   detail::ray<test_algebra> traj{{0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f}, -1.f};
0060   std::vector<scalar> positions = {0.f,   50.f,  100.f, 150.f, 200.f, 250.f,
0061                                    300.f, 350.f, 400.f, 450.f, 500.f};
0062
0063   const auto mat = silicon_tml<scalar>();
0064   constexpr scalar thickness{0.17f * unit<scalar>::cm};
0065
0066   tel_det_config<test_algebra, rectangle2D> tel_cfg{20.f * unit<scalar>::mm,
0067                                                     20.f * unit<scalar>::mm};
0068   tel_cfg.positions(positions)
0069       .pilot_track(traj)
0070       .module_material(mat)
0071       .mat_thickness(thickness);
0072
0073   const auto [det, names] =
0074       build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0075
0076   using navigator_t = caching_navigator<decltype(det)>;
0077   using stepper_t = line_stepper<test_algebra>;
0078   using pathlimit_aborter_t = actor::pathlimit_aborter<scalar>;
0079   using actor_chain_t =
0080       actor_chain<pathlimit_aborter_t,
0081                   actor::parameter_updater<test_algebra, interactor_t>>;
0082   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0083
0084   // Propagator is built from the stepper and navigator
0085   propagation::config prop_cfg{};
0086   prop_cfg.navigation.intersection.overstep_tolerance =
0087       -100.f * unit<float>::um;
0088   propagator_t p{prop_cfg};
0089
0090   constexpr scalar iniP{10.f * unit<scalar>::GeV};
0091
0092   // Bound vector
0093   bound_parameters_vector<test_algebra> bound_vector{};
0094   bound_vector.set_theta(constant<scalar>::pi_2);
0095   bound_vector.set_qop(ptc.charge() / iniP);
0096
0097   auto bound_cov = matrix::identity<covariance_t>();
0098   getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 0.f;
0099
0100   // bound track parameter at first physical plane
0101   const bound_track_parameters<test_algebra> bound_param(
0102       det.surface(0u).identifier(), bound_vector, bound_cov);
0103
0104   pathlimit_aborter_t::state aborter_state{};
0105   actor::parameter_updater_state<test_algebra> updater_state{prop_cfg,
0106                                                              bound_param};
0107   interactor_t::state interactor_state{};
0108
0109   // Create actor states tuples
0110   auto actor_states =
0111       detray::tie(aborter_state, updater_state, interactor_state);
0112
0113   propagator_t::state state(bound_param, det);
0114   state.debug(true);
0115
0116   // Propagate the entire detector
0117   ASSERT_TRUE(p.propagate(state, actor_states));
0118
0119   // new momentum
0120   const scalar newP{updater_state.bound_params().p(ptc.charge())};
0121
0122   // mass
0123   const auto mass = ptc.mass();
0124
0125   // new energy
0126   const scalar newE{std::hypot(newP, mass)};
0127
0128   // Initial energy
0129   const scalar iniE{std::hypot(iniP, mass)};
0130
0131   // New qop variance
0132   const scalar new_var_qop{
0133       getter::element(updater_state.bound_params().covariance(), e_bound_qoverp,
0134                       e_bound_qoverp)};
0135
0136   // Interaction object
0137   interaction<scalar> I;
0138
0139   // Zero incidence angle
0140   const scalar cos_inc_ang{1.f};
0141
0142   // Same material used for default telescope detector
0143   material_slab<scalar> slab(mat, thickness);
0144
0145   // Path segment in the material
0146   const scalar path_segment{slab.path_segment(cos_inc_ang)};
0147
0148   // Expected Bethe Stopping power for telescope geometry is estimated
0149   // as (number of planes * energy loss per plane assuming 1 GeV muon).
0150   // It is not perfectly precise as the track loses its energy during
0151   // propagation. However, since the energy loss << the track momentum,
0152   // the assumption is not very bad
0153   const scalar dE{
0154       I.compute_energy_loss_bethe_bloch(path_segment, slab.get_material(), ptc,
0155                                         {ptc, ptc.charge() / iniP}) *
0156       static_cast<scalar>(positions.size())};
0157
0158   // Check if the new energy after propagation is close enough to the
0159   // expected value
0160   EXPECT_NEAR(newE, iniE - dE, 1e-5f) << "dE = " << dE;
0161
0162   const scalar sigma_qop{I.compute_energy_loss_landau_sigma_QOverP(
0163       path_segment, slab.get_material(), ptc, {ptc, ptc.charge() / iniP})};
0164
0165   const scalar dvar_qop{sigma_qop * sigma_qop *
0166                         static_cast<scalar>(positions.size() - 1u)};
0167
0168   EXPECT_NEAR(new_var_qop, dvar_qop, 1e-10f);
0169
0170   // @todo: Validate the backward direction case as well?
0171 }
0172
0173 // Material interaction test with telescope Geometry
0174 GTEST_TEST(detray_material, telescope_geometry_scattering_angle) {
0175   vecmem::host_memory_resource host_mr;
0176
0177   // Build in x-direction from given module positions
0178   detail::ray<test_algebra> traj{{0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f}, -1.f};
0179   std::vector<scalar> positions = {0.f};
0180
0181   // To make sure that someone won't put more planes than one by accident
0182   EXPECT_EQ(positions.size(), 1u);
0183
0184   // Material
0185   const auto mat = silicon_tml<scalar>();
0186   const scalar thickness = 100.f * unit<scalar>::cm;
0187
0188   // Create telescope geometry
0189   tel_det_config<test_algebra, rectangle2D> tel_cfg{2000.f * unit<scalar>::mm,
0190                                                     2000.f * unit<scalar>::mm};
0191   tel_cfg.positions(positions)
0192       .pilot_track(traj)
0193       .module_material(mat)
0194       .mat_thickness(thickness);
0195
0196   const auto [det, names] =
0197       build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0198
0199   using navigator_t = caching_navigator<decltype(det)>;
0200   using stepper_t = line_stepper<test_algebra>;
0201   using simulator_t = actor::random_scatterer<test_algebra>;
0202   using pathlimit_aborter_t = actor::pathlimit_aborter<scalar>;
0203   using actor_chain_t =
0204       actor_chain<pathlimit_aborter_t,
0205                   actor::parameter_updater<test_algebra, simulator_t>>;
0206   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0207
0208   // Propagator is built from the stepper and navigator
0209   propagation::config prop_cfg{};
0210   prop_cfg.navigation.intersection.overstep_tolerance =
0211       -100.f * unit<float>::um;
0212   propagator_t p{prop_cfg};
0213
0214   constexpr scalar q{-1.f};
0215   constexpr scalar iniP{10.f * unit<scalar>::GeV};
0216
0217   // Initial track parameters directing x-axis
0218   bound_parameters_vector<test_algebra> bound_vector{};
0219   bound_vector.set_theta(constant<scalar>::pi_2);
0220   bound_vector.set_qop(q / iniP);
0221
0222   auto bound_cov = matrix::identity<covariance_t>();
0223   getter::element(bound_cov, e_bound_phi, e_bound_phi) = 0.f;
0224   getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.f;
0225
0226   // bound track parameter
0227   const bound_track_parameters<test_algebra> bound_param(
0228       det.surface(0u).identifier(), bound_vector, bound_cov);
0229
0230   std::size_t n_samples{100000u};
0231   std::vector<scalar> phis;
0232   std::vector<scalar> thetas;
0233
0234   for (std::size_t i = 0u; i < n_samples; i++) {
0235     pathlimit_aborter_t::state aborter_state{};
0236     // Seed = sample id
0237     actor::parameter_updater_state<test_algebra> updater_state{prop_cfg,
0238                                                                bound_param};
0239     simulator_t::state simulator_state{i};
0240     simulator_state.do_energy_loss = false;
0241
0242     // Create actor states tuples
0243     auto actor_states =
0244         detray::tie(aborter_state, updater_state, simulator_state);
0245
0246     propagator_t::state state(bound_param, det);
0247     state.debug(true);
0248
0249     // Propagate the entire detector
0250     ASSERT_TRUE(p.propagate(state, actor_states));
0251
0252     const auto& final_param = updater_state.bound_params();
0253
0254     // Updated phi and theta variance
0255     if (i == 0u) {
0256       interactor_t{}.update_angle_variance(
0257           bound_cov, traj.dir(), simulator_state.projected_scattering_angle);
0258     }
0259
0260     phis.push_back(final_param.phi());
0261     thetas.push_back(final_param.theta());
0262   }
0263   scalar phi_variance{statistics::rms(phis, bound_param.phi())};
0264   scalar theta_variance{statistics::rms(thetas, bound_param.theta())};
0265
0266   // Get the phi and theta variance
0267   scalar ref_phi_variance =
0268       getter::element(bound_cov, e_bound_phi, e_bound_phi);
0269   scalar ref_theta_variance =
0270       getter::element(bound_cov, e_bound_theta, e_bound_theta);
0271
0272   // Tolerate upto 1% difference
0273   EXPECT_NEAR((phi_variance - ref_phi_variance) / ref_phi_variance, 0.f, 1e-2f);
0274   EXPECT_NEAR((theta_variance - ref_theta_variance) / ref_theta_variance, 0.f,
0275               1e-2f);
0276
0277   // To make sure that the variances are not zero
0278   EXPECT_TRUE(ref_phi_variance > 1e-9f && ref_theta_variance > 1e-9f);
0279 }
0280
0281 // Material interaction test with telescope Geometry with volume material
0282 GTEST_TEST(detray_material, telescope_geometry_volume_material) {
0283   vecmem::host_memory_resource host_mr;
0284
0285   // Propagator types
0286   using bfield_t = bfield::const_field_t<scalar>;
0287   using stepper_t = rk_stepper<bfield_t::view_t, test_algebra>;
0288
0289   using pathlimit_aborter_t = actor::pathlimit_aborter<scalar>;
0290   using actor_chain_t = actor_chain<pathlimit_aborter_t>;
0291
0292   constexpr std::size_t cache_size{navigation::default_cache_size};
0293
0294   // Bfield setup
0295   test::vector3 B_z{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0296                     2.f * unit<scalar>::T};
0297   const bfield_t const_bfield = create_const_field<scalar>(B_z);
0298
0299   // Track setup
0300   constexpr scalar q{-1.f};
0301   constexpr scalar iniP{10.f * unit<scalar>::GeV};
0302
0303   bound_parameters_vector<test_algebra> bound_vector{};
0304   bound_vector.set_theta(constant<scalar>::pi_2);
0305   bound_vector.set_qop(q / iniP);
0306
0307   auto bound_cov = matrix::identity<covariance_t>();
0308   getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 0.f;
0309
0310   // Create actor states tuples
0311   const scalar path_limit = 100 * unit<scalar>::mm;
0312
0313   // Build in x-direction from given module positions
0314   detail::ray<test_algebra> traj{{0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f}, -1.f};
0315   std::vector<scalar> positions = {0.f, 5000.f * unit<scalar>::mm};
0316
0317   // NO material at modules
0318   const auto module_mat = vacuum<scalar>();
0319
0320   // Create telescope geometry
0321   tel_det_config<test_algebra, rectangle2D> tel_cfg{50000.f * unit<scalar>::mm,
0322                                                     50000.f * unit<scalar>::mm};
0323   tel_cfg.positions(positions).pilot_track(traj).module_material(module_mat);
0324
0325   std::vector<material<scalar>> vol_mats = {
0326       vacuum<scalar>(), isobutane<scalar>(), silicon<scalar>(),
0327       tungsten<scalar>()};
0328
0329   for (const auto& mat : vol_mats) {
0330     tel_cfg.volume_material(mat);
0331     const auto [det, names] =
0332         build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0333
0334     // bound track parameter at first physical plane
0335     const bound_track_parameters<test_algebra> bound_param(
0336         det.surface(0u).identifier(), bound_vector, bound_cov);
0337
0338     using navigator_t = caching_navigator<decltype(det), cache_size,
0339                                           navigation::print_inspector>;
0340
0341     using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0342
0343     // Propagator is built from the stepper and navigator
0344     propagation::config prop_cfg{};
0345     prop_cfg.navigation.intersection.overstep_tolerance =
0346         -100.f * unit<float>::um;
0347     propagator_t p{prop_cfg};
0348
0349     propagator_t::stepper_type::magnetic_field_type bfield_view(const_bfield);
0350     propagator_t::state state(bound_param, bfield_view, det);
0351
0352     pathlimit_aborter_t::state abrt_state{path_limit};
0353     auto actor_states = detray::tie(abrt_state);
0354
0355     p.propagate(state, actor_states);
0356
0357     const auto newP = state.stepping()().p(ptc.charge());
0358     const auto mass = ptc.mass();
0359
0360     const auto eloss_approx =
0361         interaction<scalar>().compute_energy_loss_bethe_bloch(
0362             state.stepping().path_length(), mat, ptc, {ptc, bound_param.qop()});
0363
0364     const auto iniE = std::sqrt(iniP * iniP + mass * mass);
0365     const auto newE = std::sqrt(newP * newP + mass * mass);
0366     const auto eloss = iniE - newE;
0367
0368     if (mat == vacuum<scalar>()) {
0369       ASSERT_FLOAT_EQ(float(eloss), 0.f);
0370     } else {
0371       ASSERT_TRUE(eloss > 0.f) << "mat.: " << mat << "\n"
0372                                << state.navigation().inspector().to_string();
0373     }
0374
0375     ASSERT_NEAR(eloss, eloss_approx, eloss * 0.01);
0376   }
0377 }