cpu/propagator/guided_navigator.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 #include "detray/definitions/units.hpp"
0010 #include "detray/geometry/mask.hpp"
0011 #include "detray/geometry/shapes/unbounded.hpp"
0012 #include "detray/navigation/caching_navigator.hpp"
0013 #include "detray/navigation/policies.hpp"
0014 #include "detray/propagator/actor_chain.hpp"
0015 #include "detray/propagator/actors/aborters.hpp"
0016 #include "detray/propagator/propagator.hpp"
0017 #include "detray/propagator/rk_stepper.hpp"
0018 #include "detray/tracks/tracks.hpp"
0019
0020 // Detray test include(s)
0021 #include "detray/test/common/bfield.hpp"
0022 #include "detray/test/common/build_telescope_detector.hpp"
0023 #include "detray/test/framework/types.hpp"
0024 #include "detray/test/utils/inspectors.hpp"
0025
0026 // vecmem include(s)
0027 #include <vecmem/memory/host_memory_resource.hpp>
0028
0029 // GTest include(s)
0030 #include <gtest/gtest.h>
0031
0032 // This tests the construction and general methods of the navigator
0033 GTEST_TEST(detray_navigation, guided_navigator) {
0034   using namespace detray;
0035   using namespace navigation;
0036
0037   using test_algebra = test::algebra;
0038   using scalar = test::scalar;
0039   using point3 = test::point3;
0040   using vector3 = test::vector3;
0041
0042   vecmem::host_memory_resource host_mr;
0043
0044   // Module positions along z-axis
0045   const std::vector<scalar> positions = {0.f,  10.f, 20.f, 30.f, 40.f, 50.f,
0046                                          60.f, 70.f, 80.f, 90.f, 100.f};
0047
0048   // Build telescope detector with unbounded rectangles
0049   tel_det_config<test_algebra, unbounded<rectangle2D>> tel_cfg{
0050       20.f * unit<scalar>::mm, 20.f * unit<scalar>::mm};
0051   tel_cfg.positions(positions).envelope(0.2f * unit<scalar>::mm);
0052
0053   const auto [telescope_det, names] =
0054       build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0055
0056   // Inspectors are optional, of course
0057   using detector_t = decltype(telescope_det);
0058   using intersection_t =
0059       intersection2D<typename detector_t::surface_type, test_algebra,
0060                      intersection::contains_pos>;
0061   using object_tracer_t =
0062       object_tracer<intersection_t, dvector, status::e_on_portal,
0063                     status::e_on_object>;
0064   using inspector_t = aggregate_inspector<object_tracer_t, print_inspector>;
0065   using b_field_t = bfield::const_field_t<scalar>;
0066   using runge_kutta_stepper =
0067       rk_stepper<b_field_t::view_t, test_algebra, unconstrained_step<scalar>,
0068                  guided_navigation>;
0069   using guided_navigator =
0070       caching_navigator<detector_t, navigation::default_cache_size, inspector_t,
0071                         intersection_t>;
0072   using pathlimit_aborter_t = actor::pathlimit_aborter<scalar>;
0073   using actor_chain_t = actor_chain<pathlimit_aborter_t>;
0074   using propagator_t =
0075       propagator<runge_kutta_stepper, guided_navigator, actor_chain_t>;
0076
0077   // track must point into the direction of the telescope
0078   const point3 pos{0.f, 0.f, 0.f};
0079   const vector3 mom{0.f, 0.f, 1.f};
0080   free_track_parameters<test_algebra> track(pos, 0.f, mom, -1.f);
0081   const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0082                   1.f * unit<scalar>::T};
0083   const b_field_t b_field = create_const_field<scalar>(B);
0084
0085   // Actors
0086   pathlimit_aborter_t::state pathlimit{200.f * unit<scalar>::cm};
0087
0088   // Propagator
0089   propagation::config prop_cfg{};
0090   propagator_t p{prop_cfg};
0091   propagator_t::stepper_type::magnetic_field_type bfield_view(b_field);
0092   propagator_t::state guided_state(track, bfield_view, telescope_det);
0093
0094   // Propagate
0095   p.propagate(guided_state, detray::tie(pathlimit));
0096
0097   auto &nav_state = guided_state.navigation();
0098   auto &debug_printer = nav_state.inspector().template get<print_inspector>();
0099   auto &obj_tracer = nav_state.inspector().template get<object_tracer_t>();
0100
0101   // Check that navigator exited
0102   ASSERT_TRUE(nav_state.finished()) << debug_printer.to_string();
0103
0104   // Sequence of surface ids we expect to see
0105   const std::vector<unsigned int> sf_sequence = {0u, 1u, 2u, 3u, 4u,  5u,
0106                                                  6u, 7u, 8u, 9u, 10u, 11u};
0107   // Check the surfaces that have been visited by the navigation
0108   EXPECT_EQ(obj_tracer.object_trace.size(), sf_sequence.size())
0109       << debug_printer.to_string();
0110   for (std::size_t i = 0u; i < sf_sequence.size(); ++i) {
0111     const auto &candidate = obj_tracer[i].intersection;
0112     auto geo_id = geometry::identifier{};
0113     geo_id.set_volume(0u).set_index(sf_sequence[i]);
0114     // The first transform in the detector belongs to the volume
0115     geo_id.set_transform(sf_sequence[i] + 1);
0116     geo_id.set_id((i == 11u) ? surface_id::e_portal : surface_id::e_sensitive);
0117     EXPECT_TRUE(candidate.surface().identifier() == geo_id)
0118         << "error at intersection on surface:\n"
0119         << "Expected: " << geo_id
0120         << "\nFound: " << candidate.surface().identifier();
0121   }
0122 }