Merge branch 'main' into positron-brems

CMakeLists.txt

@@ -6,7 +6,7 @@
 
 # Set up the project.
 cmake_minimum_required(VERSION 3.11)
-project(detray VERSION 0.81.0 LANGUAGES CXX)
+project(detray VERSION 0.83.0 LANGUAGES CXX)
 
 # Set up the used C++ standard(s).
 set(CMAKE_CXX_STANDARD 20 CACHE STRING "The (host) C++ standard to use")

core/include/detray/core/detail/container_buffers.hpp

@@ -157,11 +157,10 @@ struct dmulti_buffer : public detail::dbase_buffer {
 
 /// @brief Get the buffer representation of a vecmem vector - non-const
 template <class T>
-auto get_buffer(
-    const dvector_view<T>& vec_view, vecmem::memory_resource& mr,
-    vecmem::copy& cpy, detray::copy cpy_type = detray::copy::sync,
-    vecmem::data::buffer_type buff_type = vecmem::data::buffer_type::fixed_size
-    /*, stream*/) {
+auto get_buffer(const dvector_view<T>& vec_view, vecmem::memory_resource& mr,
+                vecmem::copy& cpy, detray::copy cpy_type = detray::copy::sync,
+                vecmem::data::buffer_type buff_type =
+                    vecmem::data::buffer_type::fixed_size) {
 
     // In case the view references a const object, return a non-const buffer
     using ret_buffer_t = dvector_buffer<std::remove_cv_t<T>>;
@@ -171,7 +170,7 @@ auto get_buffer(
     // TODO: Move this to detray copy util, which bundles vecmem copy object and
     // stream handle and gets this switch case right automatically
     if (cpy_type == detray::copy::async) {
-        cpy(vec_view, buff /*, stream*/);
+        cpy(vec_view, buff)->ignore();
     } else {
         cpy(vec_view, buff)->wait();
     }

core/include/detray/definitions/detail/math.hpp

@@ -1,6 +1,6 @@
 /** Detray library, part of the ACTS project (R&D line)
  *
- * (c) 2023 CERN for the benefit of the ACTS project
+ * (c) 2023-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -9,7 +9,7 @@
 
 // SYCL include(s).
 #if defined(CL_SYCL_LANGUAGE_VERSION) || defined(SYCL_LANGUAGE_VERSION)
-#include <CL/sycl.hpp>
+#include <sycl/sycl.hpp>
 #endif
 
 // System include(s).
@@ -19,7 +19,7 @@ namespace detray {
 
 /// Namespace to pick up math functions from
 #if defined(CL_SYCL_LANGUAGE_VERSION) || defined(SYCL_LANGUAGE_VERSION)
-namespace math = cl::sycl;
+namespace math = ::sycl;
 #elif IS_SOA
 
 namespace math {

core/include/detray/propagator/detail/jacobian_cartesian.hpp

@@ -80,9 +80,9 @@ struct jacobian<cartesian2D<algebra_t>> {
         const matrix_type<3, 2> bound_pos_to_free_pos_derivative =
             matrix_operator().template block<3, 2>(frame, 0u, 0u);
 
-        matrix_operator().template set_block(bound_to_free_jacobian,
-                                             bound_pos_to_free_pos_derivative,
-                                             e_free_pos0, e_bound_loc0);
+        matrix_operator().set_block(bound_to_free_jacobian,
+                                    bound_pos_to_free_pos_derivative,
+                                    e_free_pos0, e_bound_loc0);
     }
 
     DETRAY_HOST_DEVICE
@@ -98,9 +98,9 @@ struct jacobian<cartesian2D<algebra_t>> {
         const matrix_type<2, 3> free_pos_to_bound_pos_derivative =
             matrix_operator().template block<2, 3>(frameT, 0, 0);
 
-        matrix_operator().template set_block(free_to_bound_jacobian,
-                                             free_pos_to_bound_pos_derivative,
-                                             e_bound_loc0, e_free_pos0);
+        matrix_operator().set_block(free_to_bound_jacobian,
+                                    free_pos_to_bound_pos_derivative,
+                                    e_bound_loc0, e_free_pos0);
     }
 
     DETRAY_HOST_DEVICE

core/include/detray/propagator/detail/jacobian_cylindrical.hpp

@@ -102,9 +102,9 @@ struct jacobian<cylindrical2D<algebra_t>> {
         const auto bound_pos_to_free_pos_derivative =
             matrix_operator().template block<3, 2>(frame, 0u, 0u);
 
-        matrix_operator().template set_block(bound_to_free_jacobian,
-                                             bound_pos_to_free_pos_derivative,
-                                             e_free_pos0, e_bound_loc0);
+        matrix_operator().set_block(bound_to_free_jacobian,
+                                    bound_pos_to_free_pos_derivative,
+                                    e_free_pos0, e_bound_loc0);
     }
 
     DETRAY_HOST_DEVICE
@@ -120,9 +120,9 @@ struct jacobian<cylindrical2D<algebra_t>> {
         const auto free_pos_to_bound_pos_derivative =
             matrix_operator().template block<2, 3>(frameT, 0u, 0u);
 
-        matrix_operator().template set_block(free_to_bound_jacobian,
-                                             free_pos_to_bound_pos_derivative,
-                                             e_bound_loc0, e_free_pos0);
+        matrix_operator().set_block(free_to_bound_jacobian,
+                                    free_pos_to_bound_pos_derivative,
+                                    e_bound_loc0, e_free_pos0);
     }
 
     DETRAY_HOST_DEVICE

core/include/detray/propagator/detail/jacobian_line.hpp

@@ -120,9 +120,9 @@ struct jacobian<line2D<algebra_t>> {
         const auto bound_pos_to_free_pos_derivative =
             matrix_operator().template block<3, 2>(frame, 0u, 0u);
 
-        matrix_operator().template set_block(bound_to_free_jacobian,
-                                             bound_pos_to_free_pos_derivative,
-                                             e_free_pos0, e_bound_loc0);
+        matrix_operator().set_block(bound_to_free_jacobian,
+                                    bound_pos_to_free_pos_derivative,
+                                    e_free_pos0, e_bound_loc0);
     }
 
     DETRAY_HOST_DEVICE
@@ -138,9 +138,9 @@ struct jacobian<line2D<algebra_t>> {
         const auto free_pos_to_bound_pos_derivative =
             matrix_operator().template block<2, 3>(frameT, 0u, 0u);
 
-        matrix_operator().template set_block(free_to_bound_jacobian,
-                                             free_pos_to_bound_pos_derivative,
-                                             e_bound_loc0, e_free_pos0);
+        matrix_operator().set_block(free_to_bound_jacobian,
+                                    free_pos_to_bound_pos_derivative,
+                                    e_bound_loc0, e_free_pos0);
     }
 
     DETRAY_HOST_DEVICE

core/include/detray/propagator/detail/jacobian_polar.hpp

@@ -106,9 +106,9 @@ struct jacobian<polar2D<algebra_t>> {
         matrix_operator().template set_block<3, 1>(
             bound_pos_to_free_pos_derivative, col1, e_free_pos0, e_bound_loc1);
 
-        matrix_operator().template set_block(bound_to_free_jacobian,
-                                             bound_pos_to_free_pos_derivative,
-                                             e_free_pos0, e_bound_loc0);
+        matrix_operator().set_block(bound_to_free_jacobian,
+                                    bound_pos_to_free_pos_derivative,
+                                    e_free_pos0, e_bound_loc0);
     }
 
     DETRAY_HOST_DEVICE
@@ -149,9 +149,9 @@ struct jacobian<polar2D<algebra_t>> {
         matrix_operator().template set_block<1, 3>(
             free_pos_to_bound_pos_derivative, row1, e_bound_loc1, e_free_pos0);
 
-        matrix_operator().template set_block(free_to_bound_jacobian,
-                                             free_pos_to_bound_pos_derivative,
-                                             e_bound_loc0, e_free_pos0);
+        matrix_operator().set_block(free_to_bound_jacobian,
+                                    free_pos_to_bound_pos_derivative,
+                                    e_bound_loc0, e_free_pos0);
     }
 
     DETRAY_HOST_DEVICE

extern/vecmem/CMakeLists.txt

@@ -18,7 +18,7 @@ message(STATUS "Building VecMem as part of the Detray project")
 
 # Declare where to get VecMem from.
 set(DETRAY_VECMEM_SOURCE
-    "URL;https://github.com/acts-project/vecmem/archive/refs/tags/v1.8.0.tar.gz;URL_MD5;afddf52d9568964f25062e1c887246b7"
+    "URL;https://github.com/acts-project/vecmem/archive/refs/tags/v1.13.0.tar.gz;URL_MD5;02fe327552c21779f501c224b8c95e87"
     CACHE STRING
     "Source for VecMem, when built as part of this project"
 )

io/include/detray/io/common/geometry_reader.hpp

@@ -117,7 +117,7 @@ class geometry_reader {
         }
 
         // @TODO: Implement voume finder IO
-        det_builder.template set_volume_finder();
+        det_builder.set_volume_finder();
     }
 
     /// @returns a surface transform from its io payload @param trf_data

tests/benchmarks/cpu/benchmark_propagator.cpp

@@ -76,20 +76,50 @@ auto toy_cfg =
     toy_det_config{}.n_brl_layers(4u).n_edc_layers(7u).do_check(false);
 
 void fill_tracks(vecmem::vector<free_track_parameters<algebra_t>> &tracks,
-                 const std::size_t theta_steps, const std::size_t phi_steps) {
-    // Set momentum of tracks
-    const scalar mom_mag{10.f * unit<scalar>::GeV};
+                 const std::size_t n_tracks, bool do_sort = true) {
+    using scalar_t = dscalar<algebra_t>;
+    using uniform_gen_t =
+        detail::random_numbers<scalar_t,
+                               std::uniform_real_distribution<scalar_t>>;
+    using trk_generator_t =
+        random_track_generator<free_track_parameters<algebra_t>, uniform_gen_t>;
+
+    trk_generator_t::configuration trk_gen_cfg{};
+    trk_gen_cfg.seed(42u);
+    trk_gen_cfg.n_tracks(n_tracks);
+    trk_gen_cfg.randomize_charge(true);
+    trk_gen_cfg.phi_range(-constant<scalar_t>::pi, constant<scalar_t>::pi);
+    trk_gen_cfg.eta_range(-3.f, 3.f);
+    trk_gen_cfg.mom_range(1.f * unit<scalar_t>::GeV,
+                          100.f * unit<scalar_t>::GeV);
+    trk_gen_cfg.origin({0.f, 0.f, 0.f});
+    trk_gen_cfg.origin_stddev({0.f * unit<scalar_t>::mm,
+                               0.f * unit<scalar_t>::mm,
+                               0.f * unit<scalar_t>::mm});
 
     // Iterate through uniformly distributed momentum directions
-    for (auto traj : uniform_track_generator<free_track_parameters<algebra_t>>(
-             phi_steps, theta_steps, mom_mag)) {
+    for (auto traj : trk_generator_t{trk_gen_cfg}) {
         tracks.push_back(traj);
     }
+
+    if (do_sort) {
+        // Sort by theta angle
+        const auto traj_comp = [](const auto &lhs, const auto &rhs) {
+            constexpr auto pi_2{constant<scalar_t>::pi_2};
+            return math::fabs(pi_2 - getter::theta(lhs.dir())) <
+                   math::fabs(pi_2 - getter::theta(rhs.dir()));
+        };
+
+        std::ranges::sort(tracks, traj_comp);
+    }
 }
 
 template <propagate_option opt>
 static void BM_PROPAGATOR_CPU(benchmark::State &state) {
 
+    std::size_t n_tracks{static_cast<std::size_t>(state.range(0)) *
+                         static_cast<std::size_t>(state.range(0))};
+
     // Create the toy geometry and bfield
     auto [det, names] = build_toy_detector(host_mr, toy_cfg);
     test::vector3 B{0.f, 0.f, 2.f * unit<scalar>::T};
@@ -109,8 +139,7 @@ static void BM_PROPAGATOR_CPU(benchmark::State &state) {
 
         // Get tracks
         vecmem::vector<free_track_parameters<algebra_t>> tracks(&host_mr);
-        fill_tracks(tracks, static_cast<std::size_t>(state.range(0)),
-                    static_cast<std::size_t>(state.range(0)));
+        fill_tracks(tracks, n_tracks);
 
         total_tracks += tracks.size();
 

tests/benchmarks/cuda/benchmark_propagator_cuda.cpp

@@ -35,30 +35,61 @@ auto toy_cfg =
     toy_det_config{}.n_brl_layers(4u).n_edc_layers(7u).do_check(false);
 
 void fill_tracks(vecmem::vector<free_track_parameters<algebra_t>> &tracks,
-                 const std::size_t theta_steps, const std::size_t phi_steps) {
-    // Set momentum of tracks
-    const scalar mom_mag{10.f * unit<scalar>::GeV};
+                 const std::size_t n_tracks, bool do_sort = true) {
+    using scalar_t = dscalar<algebra_t>;
+    using uniform_gen_t =
+        detail::random_numbers<scalar_t,
+                               std::uniform_real_distribution<scalar_t>>;
+    using trk_generator_t =
+        random_track_generator<free_track_parameters<algebra_t>, uniform_gen_t>;
+
+    trk_generator_t::configuration trk_gen_cfg{};
+    trk_gen_cfg.seed(42u);
+    trk_gen_cfg.n_tracks(n_tracks);
+    trk_gen_cfg.randomize_charge(true);
+    trk_gen_cfg.phi_range(-constant<scalar_t>::pi, constant<scalar_t>::pi);
+    trk_gen_cfg.eta_range(-3.f, 3.f);
+    trk_gen_cfg.mom_range(1.f * unit<scalar_t>::GeV,
+                          100.f * unit<scalar_t>::GeV);
+    trk_gen_cfg.origin({0.f, 0.f, 0.f});
+    trk_gen_cfg.origin_stddev({0.f * unit<scalar_t>::mm,
+                               0.f * unit<scalar_t>::mm,
+                               0.f * unit<scalar_t>::mm});
 
     // Iterate through uniformly distributed momentum directions
-    for (auto traj : uniform_track_generator<free_track_parameters<algebra_t>>(
-             phi_steps, theta_steps, mom_mag)) {
+    for (auto traj : trk_generator_t{trk_gen_cfg}) {
         tracks.push_back(traj);
     }
+
+    if (do_sort) {
+        // Sort by theta angle
+        const auto traj_comp = [](const auto &lhs, const auto &rhs) {
+            constexpr auto pi_2{constant<scalar_t>::pi_2};
+            return math::fabs(pi_2 - getter::theta(lhs.dir())) <
+                   math::fabs(pi_2 - getter::theta(rhs.dir()));
+        };
+
+        std::ranges::sort(tracks, traj_comp);
+    }
 }
 
 template <propagate_option opt>
 static void BM_PROPAGATOR_CUDA(benchmark::State &state) {
 
+    std::size_t n_tracks{static_cast<std::size_t>(state.range(0)) *
+                         static_cast<std::size_t>(state.range(0))};
+
     // Create the toy geometry
-    auto [det, names] = build_toy_detector(bp_mng_mr, toy_cfg);
+    auto [det, names] = build_toy_detector(host_mr, toy_cfg);
     test::vector3 B{0.f, 0.f, 2.f * unit<scalar>::T};
     auto bfield = bfield::create_const_field(B);
 
-    // Get detector data
-    auto det_data = detray::get_data(det);
-
     // vecmem copy helper object
-    vecmem::cuda::copy copy;
+    vecmem::cuda::copy cuda_cpy;
+
+    // Copy detector to device
+    auto det_buff = detray::get_buffer(det, dev_mr, cuda_cpy);
+    auto det_view = detray::get_data(det_buff);
 
     std::size_t total_tracks = 0;
 
@@ -68,8 +99,7 @@ static void BM_PROPAGATOR_CUDA(benchmark::State &state) {
 
         // Get tracks
         vecmem::vector<free_track_parameters<algebra_t>> tracks(&bp_mng_mr);
-        fill_tracks(tracks, static_cast<std::size_t>(state.range(0)),
-                    static_cast<std::size_t>(state.range(0)));
+        fill_tracks(tracks, n_tracks);
 
         total_tracks += tracks.size();
 
@@ -79,7 +109,7 @@ static void BM_PROPAGATOR_CUDA(benchmark::State &state) {
         auto tracks_data = vecmem::get_data(tracks);
 
         // Run the propagator test for GPU device
-        propagator_benchmark(det_data, bfield, tracks_data, opt);
+        propagator_benchmark(det_view, bfield, tracks_data, opt);
     }
 
     state.counters["TracksPropagated"] = benchmark::Counter(

tests/include/detray/test/device/cuda/material_validation.hpp

@@ -81,13 +81,13 @@ struct run_material_validation {
         vecmem::data::vector_buffer<material_record_t> mat_records_buffer(
             static_cast<unsigned int>(tracks.size()), *dev_mr,
             vecmem::data::buffer_type::fixed_size);
-        cuda_cpy.setup(mat_records_buffer);
+        cuda_cpy.setup(mat_records_buffer)->wait();
         auto mat_records_view = vecmem::get_data(mat_records_buffer);
 
         // Buffer for the material parameters at every surface per track
         vecmem::data::jagged_vector_buffer<material_params_t> mat_steps_buffer(
             capacities, *dev_mr, host_mr, vecmem::data::buffer_type::resizable);
-        cuda_cpy.setup(mat_steps_buffer);
+        cuda_cpy.setup(mat_steps_buffer)->wait();
         auto mat_steps_view = vecmem::get_data(mat_steps_buffer);
 
         // Run the material tracing on device
@@ -96,10 +96,10 @@ struct run_material_validation {
 
         // Get the results back to the host and pass them on to be checked
         vecmem::vector<material_record_t> mat_records(host_mr);
-        cuda_cpy(mat_records_buffer, mat_records);
+        cuda_cpy(mat_records_buffer, mat_records)->wait();
 
         vecmem::jagged_vector<material_params_t> mat_steps(host_mr);
-        cuda_cpy(mat_steps_buffer, mat_steps);
+        cuda_cpy(mat_steps_buffer, mat_steps)->wait();
 
         return std::make_tuple(mat_records, mat_steps);
     }