cpu/simulation/scattering.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/propagator/actors/pointwise_material_interactor.hpp"
0011 #include "detray/tracks/tracks.hpp"
0012
0013 // Detray test include(s)
0014 #include "detray/test/framework/types.hpp"
0015 #include "detray/test/utils/random_scatterer.hpp"
0016 #include "detray/test/utils/scattering_helper.hpp"
0017 #include "detray/test/utils/statistics.hpp"
0018
0019 // google-test include(s).
0020 #include <gtest/gtest.h>
0021
0022 // System include(s).
0023 #include <algorithm>
0024 #include <random>
0025 #include <vector>
0026
0027 using namespace detray;
0028
0029 using test_algebra = test::algebra;
0030 using scalar = test::scalar;
0031 using vector3 = test::vector3;
0032
0033 namespace {
0034 std::random_device rd{};
0035 std::mt19937_64 generator{rd()};
0036 }  // namespace
0037
0038 // Test scattering helper
0039 GTEST_TEST(detray_simulation, scattering_helper) {
0040   generator.seed(0u);
0041
0042   const vector3 dir{2.f, 2.f, 4.f};
0043   const scalar scattering_angle{0.1f};
0044
0045   // xaxis of curvilinear plane
0046   const vector3 u = unit_vectors<vector3>().make_curvilinear_unit_u(dir);
0047
0048   // normal vector of the plane defined by u and dir
0049   const vector3 w = vector::cross(dir, u);
0050
0051   std::vector<scalar> angle_diffs;
0052
0053   std::size_t n_samples{100000u};
0054   for (std::size_t i = 0u; i < n_samples; i++) {
0055     // Scatter the input direction
0056     const vector3 new_dir =
0057         scattering_helper<test_algebra>()(dir, scattering_angle, generator);
0058
0059     // Assign a plus sign if the dot product of new direction and w vector
0060     // is plus
0061     const scalar sign = vector::dot(w, new_dir) > 0 ? 1.f : -1.f;
0062
0063     // the cosine of angle between original direction and new one
0064     scalar cos_theta =
0065         vector::dot(dir, new_dir) / (vector::norm(dir) * vector::norm(new_dir));
0066
0067     // To prevent rounding error where cos_theta is out of [-1, 1]
0068     cos_theta = std::clamp(cos_theta, scalar{-1.f}, scalar{1.f});
0069
0070     // Geht the angle between original direction and new one
0071     const scalar angle = sign * std::acos(cos_theta);
0072
0073     angle_diffs.push_back(angle);
0074   }
0075
0076   EXPECT_NEAR(statistics::mean(angle_diffs), 0.f, 1e-3f);
0077
0078   // Tolerate upto 1% difference
0079   const scalar stddev = std::sqrt(statistics::rms(angle_diffs, 0.f));
0080   EXPECT_NEAR((stddev - scattering_angle) / scattering_angle, 0.f, 1e-2f);
0081 }
0082
0083 // Test angle update
0084 GTEST_TEST(detray_simulation, angle_update) {
0085   generator.seed(0u);
0086
0087   // Initial direction
0088   vector3 dir{1.f, 2.f, 3.f};
0089   dir = vector::normalize(dir);
0090
0091   const scalar phi0 = vector::phi(dir);
0092   const scalar theta0 = vector::theta(dir);
0093
0094   // Projected scattering angle (Tests will fail with relatively large angles)
0095   const scalar projected_scattering_angle{0.01f};
0096
0097   // Initial bound covariance
0098   auto bound_cov = matrix::zero<
0099       typename bound_track_parameters<test_algebra>::covariance_type>();
0100
0101   // We are comparing the bound covariances (var[phi] and var[theta]) to the
0102   // variance of samples taken by random scattering
0103
0104   // Update the bound covariance with projected scattering angle
0105   actor::pointwise_material_interactor<test_algebra>().update_angle_variance(
0106       bound_cov, dir, projected_scattering_angle);
0107
0108   // Get the samples of phi and theta after the random scattering
0109   std::vector<scalar> phis;
0110   std::vector<scalar> thetas;
0111   std::size_t n_samples{100000u};
0112   for (std::size_t i = 0u; i < n_samples; i++) {
0113     const auto new_dir = actor::random_scatterer<test_algebra>().scatter(
0114         dir, projected_scattering_angle, generator);
0115     phis.push_back(vector::phi(new_dir));
0116     thetas.push_back(vector::theta(new_dir));
0117   }
0118
0119   // Tolerate upto 1% difference
0120   const auto var_phi = getter::element(bound_cov, e_bound_phi, e_bound_phi);
0121   const auto var_theta =
0122       getter::element(bound_cov, e_bound_theta, e_bound_theta);
0123   EXPECT_NEAR((var_phi - statistics::rms(phis, phi0)) / var_phi, 0.f, 1e-2f);
0124   EXPECT_NEAR((var_theta - statistics::rms(thetas, theta0)) / var_theta, 0.f,
0125               1e-2f);
0126 }