cpu/propagator/rk_stepper.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 // detray include(s)
0010 #include "detray/propagator/rk_stepper.hpp"
0011
0012 #include "detray/definitions/units.hpp"
0013 #include "detray/propagator/stepping_config.hpp"
0014 #include "detray/tracks/tracks.hpp"
0015 #include "detray/tracks/trajectories.hpp"
0016
0017 // Detray test include(s)
0018 #include "detray/test/common/bfield.hpp"
0019 #include "detray/test/common/track_generators.hpp"
0020 #include "detray/test/framework/types.hpp"
0021
0022 // System include(s)
0023 #include <memory>
0024
0025 // google-test include(s)
0026 #include <gtest/gtest.h>
0027
0028 using namespace detray;
0029
0030 using test_algebra = test::algebra;
0031 using scalar = test::scalar;
0032 using vector3 = test::vector3;
0033 using point3 = test::point3;
0034
0035 /// Runge-Kutta stepper
0036 template <typename bfield_t>
0037 using rk_stepper_t = rk_stepper<typename bfield_t::view_t, test_algebra>;
0038 template <typename bfield_t>
0039 using crk_stepper_t = rk_stepper<typename bfield_t::view_t, test_algebra,
0040                                  constrained_step<scalar>>;
0041
0042 namespace {
0043
0044 constexpr scalar tol{1e-3f};
0045
0046 const stepping::config step_cfg{};
0047 constexpr scalar step_size{1.f * unit<scalar>::mm};
0048 constexpr material<scalar> vol_mat{detray::cesium_iodide_with_ded<scalar>()};
0049
0050 }  // namespace
0051
0052 // This tests the base functionality of the Runge-Kutta stepper
0053 GTEST_TEST(detray_propagator, rk_stepper) {
0054   using namespace step;
0055
0056   // Constant magnetic field
0057   using bfield_t = bfield::const_field_t<scalar>;
0058
0059   vector3 B{1.f * unit<scalar>::T, 1.f * unit<scalar>::T,
0060             1.f * unit<scalar>::T};
0061   const bfield_t hom_bfield = create_const_field<scalar>(B);
0062   const bfield_t::view_t bfield_view(hom_bfield);
0063
0064   // RK stepper
0065   rk_stepper_t<bfield_t> rk_stepper;
0066   crk_stepper_t<bfield_t> crk_stepper;
0067
0068   // RK stepper configurations
0069   constexpr unsigned int rk_steps = 100u;
0070   constexpr scalar stepsize_constr{0.5f * unit<scalar>::mm};
0071
0072   // Track generator configuration
0073   const scalar p_mag{10.f * unit<scalar>::GeV};
0074   constexpr unsigned int theta_steps = 100u;
0075   constexpr unsigned int phi_steps = 100u;
0076
0077   // Iterate through uniformly distributed momentum directions
0078   for (auto track :
0079        uniform_track_generator<free_track_parameters<test_algebra>>(
0080            phi_steps, theta_steps, p_mag)) {
0081     // Generate track state used for propagation with constrained step size
0082     free_track_parameters<test_algebra> c_track(track);
0083
0084     // helix trajectory
0085     detail::helix helix(track, B);
0086
0087     // RK Stepping into forward direction
0088     rk_stepper_t<bfield_t>::state rk_state{track, bfield_view};
0089     crk_stepper_t<bfield_t>::state crk_state{c_track, bfield_view};
0090
0091     // Set step size constraint to half the nominal step size =>
0092     // crk_stepper will need twice as many steps
0093     crk_state.template set_constraint<constraint::e_user>(stepsize_constr);
0094     ASSERT_NEAR(crk_state.constraints().template size<>(),
0095                 0.5f * unit<scalar>::mm, tol);
0096
0097     // Forward stepping
0098     for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {
0099       rk_stepper.step(step_size, rk_state, step_cfg, true);
0100       crk_stepper.step(step_size, crk_state, step_cfg, true);
0101       crk_stepper.step(step_size, crk_state, step_cfg, true);
0102     }
0103
0104     // Check that both steppers arrive at the same point
0105     // Get relative error by dividing error with path length
0106     ASSERT_TRUE(rk_state.path_length() > 0.f);
0107     ASSERT_NEAR(rk_state.path_length(), crk_state.path_length(), tol);
0108     ASSERT_NEAR(vector::norm(rk_state().pos() - crk_state().pos()) /
0109                     rk_state.path_length(),
0110                 0.f, tol);
0111
0112     // Check that the stepper position lies on the truth helix
0113     const auto helix_pos = helix(rk_state.path_length());
0114     const auto forward_pos = rk_state().pos();
0115     const point3 forward_relative_error{(1.f / rk_state.path_length()) *
0116                                         (forward_pos - helix_pos)};
0117
0118     // Make sure that relative error is smaller than the tolerance
0119     EXPECT_NEAR(vector::norm(forward_relative_error), 0.f, tol);
0120
0121     // Roll the same track back to the origin
0122     // Use the same path length, since there is no overstepping
0123     const scalar path_length = rk_state.path_length();
0124     for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {
0125       rk_stepper.step(-step_size, rk_state, step_cfg, true);
0126       crk_stepper.step(-step_size, crk_state, step_cfg, true);
0127       crk_stepper.step(-step_size, crk_state, step_cfg, true);
0128     }
0129
0130     // Should arrive back at track origin, where path length is zero
0131     ASSERT_NEAR(rk_state.path_length(), 0.f, tol);
0132     ASSERT_NEAR(crk_state.path_length(), 0.f, tol);
0133
0134     const point3 backward_relative_error{1.f / (2.f * path_length) *
0135                                          (rk_state().pos())};
0136     // Make sure that relative error is smaller than the tolerance
0137     EXPECT_NEAR(vector::norm(backward_relative_error), 0.f, tol);
0138
0139     // The constrained stepper should be at the same position now
0140     ASSERT_NEAR(vector::norm(rk_state().pos() - crk_state().pos()) /
0141                     (2.f * path_length),
0142                 0.f, tol);
0143   }
0144 }
0145
0146 /// This tests the base functionality of the Runge-Kutta stepper in an
0147 /// in-homogeneous magnetic field, read from file
0148 TEST(detray_propagator, rk_stepper_inhomogeneous_bfield) {
0149   using namespace step;
0150
0151   // Read the magnetic field map
0152   using bfield_t = bfield::inhom_field_t<scalar>;
0153   bfield_t inhom_bfield = create_inhom_field<scalar>();
0154   const bfield_t::view_t bfield_view(inhom_bfield);
0155
0156   // RK stepper
0157   rk_stepper_t<bfield_t> rk_stepper;
0158   crk_stepper_t<bfield_t> crk_stepper;
0159
0160   // RK stepper configurations
0161   constexpr unsigned int rk_steps = 100u;
0162   constexpr scalar stepsize_constr{0.5f * unit<scalar>::mm};
0163
0164   // Track generator configuration
0165   const scalar p_mag{10.f * unit<scalar>::GeV};
0166   constexpr unsigned int theta_steps = 100u;
0167   constexpr unsigned int phi_steps = 100u;
0168
0169   // Iterate through uniformly distributed momentum directions
0170   for (auto track :
0171        uniform_track_generator<free_track_parameters<test_algebra>>(
0172            phi_steps, theta_steps, p_mag)) {
0173     // Generate track state used for propagation with constrained step size
0174     free_track_parameters<test_algebra> c_track(track);
0175
0176     // RK Stepping into forward direction
0177     rk_stepper_t<bfield_t>::state rk_state{track, bfield_view};
0178     crk_stepper_t<bfield_t>::state crk_state{c_track, bfield_view};
0179
0180     crk_state.template set_constraint<constraint::e_user>(stepsize_constr);
0181     ASSERT_NEAR(crk_state.constraints().template size<>(),
0182                 0.5f * unit<scalar>::mm, tol);
0183
0184     // Forward stepping
0185     for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {
0186       rk_stepper.step(step_size, rk_state, step_cfg, true);
0187       crk_stepper.step(step_size, crk_state, step_cfg, true);
0188       crk_stepper.step(step_size, crk_state, step_cfg, true);
0189     }
0190
0191     // Make sure the steppers moved
0192     const scalar path_length{rk_state.path_length()};
0193     ASSERT_TRUE(path_length > 0.f);
0194     ASSERT_TRUE(crk_state.path_length() > 0.f);
0195
0196     // Roll the same track back to the origin
0197     // Use the same path length, since there is no overstepping
0198     for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {
0199       rk_stepper.step(-step_size, rk_state, step_cfg, true);
0200       crk_stepper.step(-step_size, crk_state, step_cfg, true);
0201       crk_stepper.step(-step_size, crk_state, step_cfg, true);
0202     }
0203
0204     // Should arrive back at track origin, where path length is zero
0205     ASSERT_NEAR(rk_state.path_length(), 0.f, tol);
0206     ASSERT_NEAR(crk_state.path_length(), 0.f, tol);
0207
0208     const point3 backward_relative_error{1.f / (2.f * path_length) *
0209                                          (rk_state().pos())};
0210     // Make sure that relative error is smaller than the tolerance
0211     EXPECT_NEAR(vector::norm(backward_relative_error), 0.f, tol);
0212
0213     // The constrained stepper should be at the same position now
0214     ASSERT_NEAR(vector::norm(rk_state().pos() - crk_state().pos()) /
0215                     (2.f * path_length),
0216                 0.f, tol);
0217   }
0218 }
0219
0220 /// This tests dqop of the Runge-Kutta stepper
0221 TEST(detray_propagator, qop_derivative) {
0222   using namespace step;
0223
0224   // Constant magnetic field
0225   using bfield_t = bfield::const_field_t<scalar>;
0226
0227   vector3 B{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0228             2.f * unit<scalar>::T};
0229   const bfield_t hom_bfield = create_const_field<scalar>(B);
0230
0231   // RK stepper
0232   rk_stepper_t<bfield_t> rk_stepper;
0233   constexpr unsigned int rk_steps = 1000u;
0234
0235   // Theta phi for track generator
0236   const scalar p_mag{10.f * unit<scalar>::GeV};
0237   constexpr unsigned int theta_steps = 10u;
0238   constexpr unsigned int phi_steps = 10u;
0239
0240   const scalar ds = 1e-2f * unit<scalar>::mm;
0241
0242   // Iterate through uniformly distributed momentum directions
0243   for (auto track :
0244        uniform_track_generator<free_track_parameters<test_algebra>>(
0245            phi_steps, theta_steps, p_mag)) {
0246     // RK Stepping into forward direction
0247     rk_stepper_t<bfield_t>::state rk_state{track, hom_bfield};
0248
0249     for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {
0250       const scalar qop1 = rk_state().qop();
0251       const scalar d2qopdsdqop = rk_state.d2qopdsdqop(qop1, &vol_mat);
0252
0253       const scalar dqopds1 = rk_state.dqopds(qop1, &vol_mat);
0254
0255       rk_stepper.step(ds, rk_state, step_cfg, true, &vol_mat);
0256
0257       const scalar qop2 = rk_state().qop();
0258       const scalar dqopds2 = rk_state.dqopds(qop2, &vol_mat);
0259
0260       ASSERT_TRUE(qop1 > qop2);
0261       ASSERT_NEAR((qop2 - qop1) / ds, dqopds1, 1e-4);
0262       ASSERT_NEAR((dqopds2 - dqopds1) / (qop2 - qop1), d2qopdsdqop, 1e-4);
0263     }
0264   }
0265 }