From 82de6315e790f56ed59a671462615a4c2d905040 Mon Sep 17 00:00:00 2001
From: Stephen Nicholas Swatman <stephen.nicholas.swatman@cern.ch>
Date: Tue, 25 Jun 2024 15:51:27 +0200
Subject: [PATCH 1/2] test

---
 .../traccc/seeding/doublet_finding_helper.hpp |  16 +
 device/cuda/CMakeLists.txt                    |   5 +
 .../cuda/seeding2/kernels/kd_tree_kernel.hpp  |  24 +
 .../seeding2/kernels/seed_finding_kernel.hpp  |  27 +
 .../seeding2/kernels/write_output_kernel.hpp  |  26 +
 .../traccc/cuda/seeding2/seed_finding.hpp     |  35 +
 .../cuda/seeding2/types/internal_sp.hpp       |  35 +
 .../traccc/cuda/seeding2/types/kd_tree.hpp    |  36 +
 .../traccc/cuda/seeding2/types/range3d.hpp    |  73 ++
 .../traccc/cuda/utils/device_traits.hpp       |  56 ++
 .../cuda/include/traccc/cuda/utils/sort.hpp   | 120 +++
 .../cuda/include/traccc/cuda/utils/sync.hpp   |  40 +
 .../src/seeding2/kernels/kd_tree_kernel.cu    | 833 ++++++++++++++++++
 .../seeding2/kernels/seed_finding_kernel.cu   | 706 +++++++++++++++
 .../seeding2/kernels/write_output_kernel.cu   |  82 ++
 device/cuda/src/seeding2/seed_finding2.cu     |  96 ++
 16 files changed, 2210 insertions(+)
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/types/kd_tree.hpp
 create mode 100644 device/cuda/include/traccc/cuda/seeding2/types/range3d.hpp
 create mode 100644 device/cuda/include/traccc/cuda/utils/device_traits.hpp
 create mode 100644 device/cuda/include/traccc/cuda/utils/sort.hpp
 create mode 100644 device/cuda/include/traccc/cuda/utils/sync.hpp
 create mode 100644 device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
 create mode 100644 device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
 create mode 100644 device/cuda/src/seeding2/kernels/write_output_kernel.cu
 create mode 100644 device/cuda/src/seeding2/seed_finding2.cu

diff --git a/core/include/traccc/seeding/doublet_finding_helper.hpp b/core/include/traccc/seeding/doublet_finding_helper.hpp
index 224b13da03..dfb9b96765 100644
--- a/core/include/traccc/seeding/doublet_finding_helper.hpp
+++ b/core/include/traccc/seeding/doublet_finding_helper.hpp
@@ -34,6 +34,11 @@ struct doublet_finding_helper {
         const internal_spacepoint<spacepoint>& sp2,
         const seedfinder_config& config);
 
+    static inline TRACCC_HOST_DEVICE bool isCompatible(
+        bool top, const internal_spacepoint<spacepoint>& sp1,
+        const internal_spacepoint<spacepoint>& sp2,
+        const seedfinder_config& config);
+
     /// Do the conformal transformation on doublet's coordinate
     ///
     /// @param sp1 is middle spacepoint
@@ -48,6 +53,17 @@ struct doublet_finding_helper {
                           const internal_spacepoint<spacepoint>& sp2);
 };
 
+bool doublet_finding_helper::isCompatible(
+    bool top, const internal_spacepoint<spacepoint>& sp1,
+    const internal_spacepoint<spacepoint>& sp2,
+    const seedfinder_config& config) {
+    if (top) {
+        return isCompatible<details::spacepoint_type::top>(sp1, sp2, config);
+    } else {
+        return isCompatible<details::spacepoint_type::bottom>(sp1, sp2, config);
+    }
+}
+
 template <details::spacepoint_type otherSpType>
 bool TRACCC_HOST_DEVICE
 doublet_finding_helper::isCompatible(const internal_spacepoint<spacepoint>& sp1,
diff --git a/device/cuda/CMakeLists.txt b/device/cuda/CMakeLists.txt
index 1321e547b5..0e7067222f 100644
--- a/device/cuda/CMakeLists.txt
+++ b/device/cuda/CMakeLists.txt
@@ -39,6 +39,11 @@ traccc_add_library( traccc_cuda cuda TYPE SHARED
   "src/seeding/spacepoint_binning.cu"
   "include/traccc/cuda/seeding/experimental/spacepoint_formation.hpp"
   "src/seeding/experimental/spacepoint_formation.cu"
+  # Alternate seeding.
+  "src/seeding2/seed_finding2.cu"
+  "src/seeding2/kernels/seed_finding_kernel.cu"
+  "src/seeding2/kernels/kd_tree_kernel.cu"
+  "src/seeding2/kernels/write_output_kernel.cu"
   # Clusterization
   "include/traccc/cuda/clusterization/clusterization_algorithm.hpp"
   "src/clusterization/clusterization_algorithm.cu"
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
new file mode 100644
index 0000000000..4ce1325af7
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
@@ -0,0 +1,24 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <traccc/cuda/seeding2/types/internal_sp.hpp>
+#include <traccc/cuda/seeding2/types/kd_tree.hpp>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <vector>
+
+namespace traccc::cuda {
+/**
+ * @brief Creates a k-d tree from a given set of spacepoints.
+ *
+ * @return A pair containing the k-d tree nodes as well as the number of nodes.
+ */
+std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
+    vecmem::memory_resource&, internal_sp_owning_t&&, uint32_t);
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
new file mode 100644
index 0000000000..464b4dd2af
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
@@ -0,0 +1,27 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <traccc/cuda/seeding2/types/internal_sp.hpp>
+#include <traccc/cuda/seeding2/types/kd_tree.hpp>
+#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/seeding/detail/seeding_config.hpp>
+
+namespace traccc::cuda {
+/**
+ * @brief Execute the seed finding kernel itself.
+ *
+ * @return A pair containing the list of internal seeds as well as the number
+ * of seeds.
+ */
+std::pair<vecmem::unique_alloc_ptr<alt_seed[]>, uint32_t> run_seeding(
+    seedfinder_config, seedfilter_config, vecmem::memory_resource&,
+    internal_sp_t, kd_tree_t);
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
new file mode 100644
index 0000000000..0bddf482d0
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
@@ -0,0 +1,26 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <traccc/cuda/seeding2/types/internal_sp.hpp>
+#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/seed.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/utils/memory_resource.hpp>
+
+namespace traccc::cuda {
+/**
+ * @brief Kernel to write output data back into traccc's EDM.
+ *
+ * @return A vector buffer containing the output seeds.
+ */
+alt_seed_collection_types::buffer write_output(const traccc::memory_resource &,
+                                               uint32_t, const internal_sp_t,
+                                               const alt_seed *const);
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp b/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
new file mode 100644
index 0000000000..a350a73730
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
@@ -0,0 +1,35 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/seeding/detail/seeding_config.hpp>
+#include <traccc/utils/algorithm.hpp>
+#include <traccc/utils/memory_resource.hpp>
+#include <vecmem/memory/cuda/device_memory_resource.hpp>
+
+namespace traccc::cuda {
+/**
+ * @brief Alternative seed finding algorithm, using orthogonal range search
+ * implemented through a k-d tree.
+ */
+class seed_finding2 : public algorithm<alt_seed_collection_types::buffer(
+                          const spacepoint_collection_types::const_view&)> {
+    public:
+    seed_finding2(const traccc::memory_resource& mr);
+
+    output_type operator()(
+        const spacepoint_collection_types::const_view& sps) const override;
+
+    private:
+    traccc::memory_resource m_output_mr;
+    seedfinder_config m_finder_conf;
+    seedfilter_config m_filter_conf;
+};
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp b/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
new file mode 100644
index 0000000000..94b375ff4d
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
@@ -0,0 +1,35 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <array>
+#include <cstdint>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/utils/array_wrapper.hpp>
+#include <vecmem/memory/unique_ptr.hpp>
+
+namespace traccc::cuda {
+template <template <typename> typename F>
+struct sp_t {
+    using tuple_t = std::tuple<float, float, float, float, float, unsigned int>;
+
+    F<float> x;
+    F<float> y;
+    F<float> z;
+    F<float> phi;
+    F<float> radius;
+    F<unsigned int> link;
+};
+
+using internal_sp = array_wrapper<soa, sp_t>;
+
+using internal_sp_owning_t = internal_sp::owner;
+
+using internal_sp_t = internal_sp::handle;
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/types/kd_tree.hpp b/device/cuda/include/traccc/cuda/seeding2/types/kd_tree.hpp
new file mode 100644
index 0000000000..98f103d819
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/types/kd_tree.hpp
@@ -0,0 +1,36 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <traccc/cuda/seeding2/types/range3d.hpp>
+#include <traccc/utils/array_wrapper.hpp>
+#include <vecmem/memory/unique_ptr.hpp>
+
+namespace traccc::cuda {
+enum class nodetype_e { LEAF, INTERNAL, NON_EXTANT };
+
+enum class pivot_e { Phi, R, Z };
+
+template <template <typename> typename F>
+struct node_t {
+    using tuple_t =
+        std::tuple<nodetype_e, range3d, uint32_t, uint32_t, pivot_e, float>;
+
+    F<nodetype_e> type;
+    F<range3d> range;
+    F<uint32_t> begin, end;
+    F<pivot_e> dim;
+    F<float> mid;
+};
+
+using kd_tree = array_wrapper<soa, node_t>;
+
+using kd_tree_owning_t = kd_tree::owner;
+
+using kd_tree_t = kd_tree::handle;
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/types/range3d.hpp b/device/cuda/include/traccc/cuda/seeding2/types/range3d.hpp
new file mode 100644
index 0000000000..1914b68ff5
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/seeding2/types/range3d.hpp
@@ -0,0 +1,73 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include <array>
+#include <cstdint>
+
+namespace traccc::cuda {
+/**
+ * @brief Range in three dimensions, defined by minima and maxima in the phi,
+ * r, and z dimensions.
+ */
+struct range3d {
+    float phi_min, phi_max, r_min, r_max, z_min, z_max;
+
+    __host__ __device__ __forceinline__ static range3d Infinite() {
+        range3d r;
+
+        r.phi_min = std::numeric_limits<float>::lowest();
+        r.phi_max = std::numeric_limits<float>::max();
+        r.r_min = 0.f;
+        r.r_max = std::numeric_limits<float>::max();
+        r.z_min = std::numeric_limits<float>::lowest();
+        r.z_max = std::numeric_limits<float>::max();
+
+        return r;
+    }
+
+    __host__ __device__ __forceinline__ static range3d Degenerate() {
+        range3d r;
+
+        r.phi_min = std::numeric_limits<float>::max();
+        r.phi_max = std::numeric_limits<float>::lowest();
+        r.r_min = std::numeric_limits<float>::max();
+        r.r_max = 0.f;
+        r.z_min = std::numeric_limits<float>::max();
+        r.z_max = std::numeric_limits<float>::lowest();
+
+        return r;
+    }
+
+    __host__ __device__ __forceinline__ static range3d Union(const range3d& a,
+                                                             const range3d& b) {
+        return {std::min(a.phi_min, b.phi_min), std::max(a.phi_max, b.phi_max),
+                std::min(a.r_min, b.r_min),     std::max(a.r_max, b.r_max),
+                std::min(a.z_min, b.z_min),     std::max(a.z_max, b.z_max)};
+    }
+
+    __host__ __device__ __forceinline__ bool intersects(
+        const range3d& o) const {
+        return phi_min <= o.phi_max && o.phi_min < phi_max &&
+               r_min <= o.r_max && o.r_min < r_max && z_min <= o.z_max &&
+               o.z_min < z_max;
+    }
+
+    __host__ __device__ __forceinline__ bool dominates(const range3d& o) const {
+        return phi_min <= o.phi_min && o.phi_max <= phi_max &&
+               r_min <= o.r_min && o.r_max <= r_max && z_min <= o.z_min &&
+               o.z_max <= z_max;
+    }
+
+    __host__ __device__ __forceinline__ bool contains(float phi, float r,
+                                                      float z) const {
+        return phi_min <= phi && phi <= phi_max && r_min <= r && r <= r_max &&
+               z_min <= z && z <= z_max;
+    }
+};
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/utils/device_traits.hpp b/device/cuda/include/traccc/cuda/utils/device_traits.hpp
new file mode 100644
index 0000000000..c1cc59eb1c
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/utils/device_traits.hpp
@@ -0,0 +1,56 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2021 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#if defined(__CUDACC__) && defined(__CUDA_ARCH__)
+// Set values of maximum resident threads per SM
+#if __CUDA_ARCH__ < 500
+#pragma message \
+    "Very old CUDA architecture, setting maximum resident threads per SM to 1024."
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
+#elif __CUDA_ARCH__ <= 720
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
+#elif __CUDA_ARCH__ <= 750
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
+#elif __CUDA_ARCH__ <= 800
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
+#elif __CUDA_ARCH__ <= 890
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1536u
+#elif __CUDA_ARCH__ <= 900
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
+#pragma message \
+    "Unknown CUDA architecture, setting maximum resident threads per SM to 1024."
+#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
+#endif
+
+// Set values of maximum resident blocks per SM
+#if __CUDA_ARCH__ < 500
+#pragma message \
+    "Very old CUDA architecture, setting maximum resident blocks per SM to 16."
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
+#elif __CUDA_ARCH__ <= 720
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
+#elif __CUDA_ARCH__ <= 750
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
+#elif __CUDA_ARCH__ <= 800
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
+#elif __CUDA_ARCH__ <= 870
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
+#elif __CUDA_ARCH__ <= 890
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 24u
+#elif __CUDA_ARCH__ <= 900
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
+#else
+#pragma message \
+    "Very old CUDA architecture, setting maximum resident blocks per SM to 16."
+#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
+#endif
+
+// Set values of maximum threadss per block
+#define CUDA_MAX_THREADS_PER_BLOCK 1024u
+#endif
diff --git a/device/cuda/include/traccc/cuda/utils/sort.hpp b/device/cuda/include/traccc/cuda/utils/sort.hpp
new file mode 100644
index 0000000000..5a1e0b9b15
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/utils/sort.hpp
@@ -0,0 +1,120 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2021 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+namespace traccc::cuda {
+/**
+ * @brief Swap two values of arbitrary type.
+ *
+ * @tparam T The type of values to swap.
+ *
+ * @param a The first object in the swap (will take the value of b).
+ * @param b The second object in the swap (will take the value of a).
+ */
+template <typename T>
+__device__ __forceinline__ void swap(T& a, T& b) {
+    T t = a;
+    a = b;
+    b = t;
+}
+
+/**
+ * @brief Perform a block-wide odd-even key sorting.
+ *
+ * This function performs a sorting operation across the entire block, assuming
+ * that all the threads in the block are currently active.
+ *
+ * @warning The behaviour of this function is ill-defined if any of the threads
+ * in the block have exited.
+ *
+ * @warning This method is efficient for sorting small arrays, preferable in
+ * shared memory, but given the O(n^2) worst-case performance this should not
+ * be used on larger arrays.
+ *
+ * @tparam K The type of keys to sort.
+ * @tparam C The type of the comparison function.
+ *
+ * @param keys An array of keys to sort.
+ * @param num_keys The number of keys in the array to sort.
+ * @param comparison A comparison function.
+ */
+template <typename K, typename C>
+__device__ void blockOddEvenKeySort(K* keys, uint32_t num_keys,
+                                    C&& comparison) {
+    bool sorted;
+
+    do {
+        sorted = true;
+
+        for (uint32_t j = 2 * threadIdx.x + 1; j < num_keys - 1;
+             j += 2 * blockDim.x) {
+            if (comparison(keys[j + 1], keys[j])) {
+                swap(keys[j + 1], keys[j]);
+                sorted = false;
+            }
+        }
+
+        __syncthreads();
+
+        for (uint32_t j = 2 * threadIdx.x; j < num_keys - 1;
+             j += 2 * blockDim.x) {
+            if (comparison(keys[j + 1], keys[j])) {
+                swap(keys[j + 1], keys[j]);
+                sorted = false;
+            }
+        }
+    } while (__syncthreads_or(!sorted));
+}
+
+/**
+ * @brief Perform a warp-wide odd-even key sorting.
+ *
+ * This function performs a sorting operation across a single warp, assuming
+ * that all the threads in the warp are currently active.
+ *
+ * @warning The behaviour of this function is ill-defined if any of the threads
+ * in the warp have exited.
+ *
+ * @warning This method is efficient for sorting small arrays, preferable in
+ * shared memory, but given the O(n^2) worst-case performance this should not
+ * be used on larger arrays.
+ *
+ * @tparam K The type of keys to sort.
+ * @tparam C The type of the comparison function.
+ *
+ * @param keys An array of keys to sort.
+ * @param num_keys The number of keys in the array to sort.
+ * @param comparison A comparison function.
+ */
+template <typename K, typename C>
+__device__ void warpOddEvenKeySort(K* keys, uint32_t num_keys, C&& comparison) {
+    bool sorted;
+
+    do {
+        sorted = true;
+
+        for (uint32_t j = 2 * (threadIdx.x % warpSize) + 1; j < num_keys - 1;
+             j += 2 * warpSize) {
+            if (comparison(keys[j + 1], keys[j])) {
+                swap(keys[j + 1], keys[j]);
+                sorted = false;
+            }
+        }
+
+        __syncwarp(0xFFFFFFFF);
+
+        for (uint32_t j = 2 * (threadIdx.x % warpSize); j < num_keys - 1;
+             j += 2 * warpSize) {
+            if (comparison(keys[j + 1], keys[j])) {
+                swap(keys[j + 1], keys[j]);
+                sorted = false;
+            }
+        }
+    } while (__any_sync(0xFFFFFFFF, !sorted));
+}
+}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/utils/sync.hpp b/device/cuda/include/traccc/cuda/utils/sync.hpp
new file mode 100644
index 0000000000..4e0c3c0c56
--- /dev/null
+++ b/device/cuda/include/traccc/cuda/utils/sync.hpp
@@ -0,0 +1,40 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2021-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+namespace traccc::cuda {
+/**
+ * @brief Calculate the index of the current thread in the set of threads in
+ * the warp for which a predicate evaluates to true.
+ *
+ * @param predicate The predicate to evaluate.
+ *
+ * @returns The tuple (T, I), where T is the total number of threads in the
+ * warp for which the provided predicate was true, and where I is a unique,
+ * consecutive index in the set of threads for which this was true.
+ *
+ * @example If a predicate P evaluates as `[T, F, F, T]` in a 4-lane warp, the
+ * call to `warp_indexed_ballot_sync(P)` will evaluate to the values
+ * `[(2, 0), (2, ?), (2, ?), (2, 1)]` where `?` indicates an undefined value.
+ *
+ * @warning The value of T is always well-defined, the value if I is only
+ * well-defined if the given predicate was true for the given thread.
+ *
+ * @note This function forces thread synchronization.
+ */
+__device__ __forceinline__ std::pair<uint32_t, uint32_t>
+warp_indexed_ballot_sync(bool predicate) {
+    uint32_t mask = __ballot_sync(0xFFFFFFFFu, predicate);
+
+    uint32_t tot = __popc(mask);
+    uint32_t idx =
+        __popc(mask & ~(0xFFFFFFFEu << (threadIdx.x % warpSize))) - 1u;
+
+    return {tot, idx};
+}
+}  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu b/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
new file mode 100644
index 0000000000..ff5d54c301
--- /dev/null
+++ b/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
@@ -0,0 +1,833 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include <cooperative_groups.h>
+
+#include <iostream>
+#include <limits>
+#include <memory>
+#include <traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp>
+#include <traccc/cuda/seeding2/seed_finding.hpp>
+#include <traccc/cuda/seeding2/types/kd_tree.hpp>
+#include <traccc/cuda/seeding2/types/range3d.hpp>
+#include <traccc/cuda/utils/definitions.hpp>
+#include <traccc/cuda/utils/device_traits.hpp>
+#include <traccc/cuda/utils/sort.hpp>
+#include <traccc/edm/seed.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/seeding/detail/lin_circle.hpp>
+#include <traccc/seeding/detail/seeding_config.hpp>
+#include <vecmem/memory/memory_resource.hpp>
+#include <vecmem/memory/unique_ptr.hpp>
+#include <vector>
+
+namespace traccc::cuda {
+/**
+ * @brief Bit mask for halo point indices.
+ *
+ * Our phi-axis search range wraps around, and k-d trees do not natively
+ * support this. Thankfully, we can solve this problem by simply duplicating
+ * points at the edges of the phi range. These are called phi halo points. We
+ * mark these points by setting the MSB in their index. This means that, in
+ * principle, the index is invalid, but it also saves quite a lot of space.
+ * This mask contains the exact bit that is set as the mask.
+ */
+constexpr uint32_t AUXILIARY_INDEX_MASK = (1u << 31);
+
+/**
+ * @brief Width of phi-halo range.
+ *
+ * This value defines the range in which points are copied and treated as halo
+ * points. If this width is ε, then points below -π + ε and above π - ε are
+ * halo points. This value should be at least as large as the maximum phi delta
+ * betweel spacepoints.
+ */
+// constexpr float PHI_HALO_WIDTH = 0.06f;
+
+/**
+ * @brief Get coordinate value for a halo-aware index in a given list of
+ * spacepoints.
+ *
+ * For non-halo points, this function simply returns the phi value of the
+ * point. But for halo points, it calculates the proper adjusted phi value. For
+ * points in the range -π ≤ φ ≤ -π + ε, this is φ + 2π. For points in the range
+ * π - ε ≤ φ ≤ π it is φ - 2π.
+ *
+ * @param[in] spacepoints The array of spacepoints.
+ * @param[in] i The halo-aware index in the spacepoint array.
+ *
+ * @return The ϕ value of the spacepoint with the given index.
+ */
+__device__ float get_coordinate(const internal_sp_t spacepoints, uint32_t i,
+                                pivot_e dim) {
+    bool is_real = !(i & AUXILIARY_INDEX_MASK);
+
+    if (is_real) {
+        if (dim == pivot_e::Phi) {
+            return spacepoints[i].phi;
+        } else if (dim == pivot_e::R) {
+            return spacepoints[i].radius;
+        } else {
+            return spacepoints[i].z;
+        }
+    } else {
+        uint32_t real_i = i & ~AUXILIARY_INDEX_MASK;
+
+        if (dim == pivot_e::Phi) {
+            float phi = spacepoints[real_i].phi;
+
+            if (phi >= 0) {
+                return phi - 2.f * static_cast<float>(M_PI);
+            } else {
+                return phi + 2.f * static_cast<float>(M_PI);
+            }
+        } else if (dim == pivot_e::R) {
+            return spacepoints[real_i].radius;
+        } else {
+            return spacepoints[real_i].z;
+        }
+    }
+}
+
+/**
+ * @brief Initialize the index vector used to build the k-d tree.
+ *
+ * Our k-d tree construction algorithm keeps track of spacepoints using a big
+ * index array, which is initialized here. Given n spacepoints, this has n + m
+ * elements, where m is the number of halo points.
+ *
+ * @param[in] spacepoints The list of spacepoints.
+ * @param[out] indices The output array to write to.
+ * @param[in] sps The total number of spacepoints.
+ * @param[out] extra_indices The output halo point count.
+ */
+__global__ void initialize_index_vector(internal_sp_t spacepoints,
+                                        uint32_t* __restrict__ indices,
+                                        uint32_t sps,
+                                        uint32_t* __restrict__ extra_indices) {
+    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
+
+    if (grid.thread_rank() < sps) {
+        indices[grid.thread_rank()] = grid.thread_rank();
+
+        /*
+         * If we are a halo point, add an extra index for our evil twin!
+         */
+        // if (spacepoints.phi(grid.thread_rank()) < (-M_PI + PHI_HALO_WIDTH) ||
+        //     spacepoints.phi(grid.thread_rank()) > (M_PI - PHI_HALO_WIDTH)) {
+        //     indices[sps + atomicAdd(extra_indices, 1u)] =
+        //         grid.thread_rank() | AUXILIARY_INDEX_MASK;
+        // }
+    }
+}
+
+/**
+ * @brief Kernel that initializes the k-d tree.
+ *
+ * This kernel simply sets every node in a given k-d tree to be non-extant
+ * except for the root (at index 0) which is considered an incomplete node.
+ *
+ * @param[out] tree The array of tree nodes to write to.
+ * @param[in] num_nodes The total number of tree nodes.
+ * @param[in] num_indices The total number of indices.
+ */
+__global__ void initialize_kd_tree(kd_tree_t tree, uint32_t num_nodes,
+                                   uint32_t num_indices,
+                                   uint32_t* __restrict__ work_list,
+                                   uint32_t* __restrict__ work_list_size) {
+    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
+
+    if (grid.thread_rank() == 0) {
+        /*
+         * The root node starts as an incomplete node covering all indices in
+         * the index array. Its range is initialized as an infinite range.
+         */
+        tree[0].type = nodetype_e::LEAF;
+        tree[0].begin = 0;
+        tree[0].end = num_indices;
+        tree[0].range = range3d::Infinite();
+
+        work_list[atomicAdd(work_list_size, 1u)] = 0;
+    } else if (grid.thread_rank() < num_nodes) {
+        /*
+         * All other nodes are non-extant.
+         */
+        tree[grid.thread_rank()].type = nodetype_e::NON_EXTANT;
+    }
+}
+
+/**
+ * @brief Helper function to determine the next pivot axis.
+ *
+ * k-d trees work by pivoting along a certain axis. In this implementation, we
+ * simply cycle through dimensions because that's very fast. We go in the order
+ * ϕ to r to z, and then repeat.
+ *
+ * @param[in] p The old pivot axis.
+ *
+ * @return The new pivot axis.
+ */
+__device__ pivot_e next_pivot(pivot_e p) {
+    if (p == pivot_e::Phi) {
+        return pivot_e::R;
+    } else if (p == pivot_e::R) {
+        return pivot_e::Z;
+    } else {
+        return pivot_e::Phi;
+    }
+}
+
+/**
+ * @brief Approximate mean finding algorithm for large numbers of spacepoints.
+ *
+ * When we have a lot of spacepoints to find the median for, doing so exactly
+ * is too slow. Thus, we use a sampling approach; we select a subset of points
+ * and calculate the median of those points.
+ *
+ * @param[in] spacepoints The list of spacepoints.
+ * @param[inout] indices The list of indices for _this_ computation.
+ * @param[in] num_indices The number of indices for _this_ function call.
+ * @param[in] pivot The requested axis along which to split.
+ *
+ * @return The approximate median value of the given indices in the given
+ * pivot axis.
+ */
+__device__ float get_approximate_median(internal_sp_t spacepoints,
+                                        const uint32_t* __restrict__ indices,
+                                        uint32_t num_indices, pivot_e pivot) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+
+    /*
+     * First, we declare some shared memory which will hold our samples from
+     * which to compute the median. We assume a maximum of one item per thread.
+     */
+    extern __shared__ float values[];
+
+    /*
+     * Assert that the number of indices is at least as big as the block size,
+     * otherwise this computation might not work as designed!
+     */
+    assert(num_indices >= block.size());
+
+    /*
+     * Calculate the stride between samples in the spacepoint array.
+     */
+    uint32_t delta = num_indices / block.size();
+
+    /*
+     * Calculate the index for the current thread.
+     */
+    uint32_t index = indices[(block.thread_rank() * delta) % num_indices];
+
+    /*
+     * Calculate the value of the chosen spacepoint in the requested pivot
+     * dimension and store it in the shared memory.
+     */
+    values[block.thread_rank()] = get_coordinate(spacepoints, index, pivot);
+
+    blockOddEvenKeySort(values, static_cast<uint32_t>(block.size()),
+                        std::less<float>());
+
+    /*
+     * Now that our sorting is complete, we can simply return the middle value
+     * in the range to get the median.
+     */
+    return values[(block.size() + 1) / 2];
+}
+
+/**
+ * @brief Exact mean finding algorithm for small ranges.
+ *
+ * When the array of indices is relatively small, we can afford to find the
+ * exact median. We do this by sorting the data and then returning the index
+ * (rather than the value) of the median.
+ *
+ * @warning This function has the side effect that it leaves the index array
+ * sorted!
+ *
+ * @param[in] spacepoints The array of spacepoints.
+ * @param[inout] indices The array of indices for this function call.
+ * @param[in] num_indices The number of indices for this function call.
+ * @param[in] pivot The requested pivot axis.
+ */
+__device__ uint32_t sort_and_get_median(internal_sp_t spacepoints,
+                                        const uint32_t* __restrict__ indices,
+                                        uint32_t* __restrict__ index_buffer,
+                                        uint32_t num_indices, pivot_e pivot) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+
+    for (std::size_t i = block.thread_index().x; i < num_indices;
+         i += block.size()) {
+        index_buffer[i] = indices[i];
+    }
+
+    /*
+     * Simple perform an odd-even sort to sort the range of indices by the
+     * value of the corresponding SPs in the requested pivot axis.
+     */
+    blockOddEvenKeySort(
+        index_buffer, num_indices,
+        [&spacepoints, pivot](const uint32_t& a, const uint32_t& b) {
+            return get_coordinate(spacepoints, a, pivot) <
+                   get_coordinate(spacepoints, b, pivot);
+        });
+
+    /*
+     * Remember, we are returning an index here, not a value!
+     */
+    return (num_indices + 1) / 2;
+}
+
+/**
+ * @brief Create an internal node in the k-d tree.
+ *
+ * This function constructs an internal node in the k-d tree at a given
+ * position.
+ *
+ * @warning This function has the side-effect of re-ordering the index array.
+ *
+ * @param[inout] tree The k-d tree in which to build the node.
+ * @param[in] spacepoints The spacepoint array.
+ * @param[inout] indices The indices from which to build an internal node.
+ * @param[inout] index_buffer Temporary storage for partitioning the indices.
+ * @param[in] nid The node index.
+ * @param[in] pivot The requested pivot axis.
+ */
+__device__ bool create_internal_node(kd_tree_t tree, internal_sp_t spacepoints,
+                                     const uint32_t* __restrict__ indices,
+                                     uint32_t* __restrict__ index_buffer,
+                                     int nid, pivot_e pivot) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+
+    /*
+     * Get the node and it's left and right children.
+     */
+    uint32_t lnid = 2 * nid + 1;
+    uint32_t rnid = 2 * nid + 2;
+
+    /*
+     * Start the construction of the two children, only if we are the leader
+     * thread.
+     */
+    if (block.thread_rank() == 0) {
+        tree[lnid].type = nodetype_e::LEAF;
+        tree[rnid].type = nodetype_e::LEAF;
+    }
+
+    /*
+     * Let's hope this doesn't fire, but if we have zero (or god forbid a
+     * negative number) nodes, we have a big problem.
+     */
+    assert(tree[nid].end > tree[nid].begin);
+
+    uint32_t num_indices = tree[nid].end - tree[nid].begin;
+
+    block.sync();
+
+    float mid_point_value;
+
+    /*
+     * We have two different branches here; for small nodes we take an exact
+     * approach, for big nodes we make approximations.
+     */
+    if (num_indices > block.size()) {
+        /*
+         * These counters are used for partitioning, and count how many points
+         * we have that are greater than our pivot, as well as how many are
+         * less than the pivot.
+         */
+        __shared__ uint32_t lower_idx, upper_idx;
+
+        if (block.thread_rank() == 0) {
+            lower_idx = 0;
+            upper_idx = 0;
+        }
+
+        /*
+         * Find the approximate median value of the given spacepoints in the
+         * selected dimension.
+         */
+        float mid_point =
+            get_approximate_median(spacepoints, indices, num_indices, pivot);
+
+        block.sync();
+
+        /*
+         * Now, we partition the spacepoints (or rather, their indices) into
+         * lower and upper points, counting the number of elements we have in
+         * each.
+         */
+        for (uint32_t i = block.thread_rank(); i < num_indices;
+             i += block.size()) {
+            if (get_coordinate(spacepoints, indices[i], pivot) >= mid_point) {
+                index_buffer[num_indices - (1u + atomicAdd(&upper_idx, 1u))] =
+                    indices[i];
+            } else {
+                index_buffer[atomicAdd(&lower_idx, 1u)] = indices[i];
+            }
+        }
+
+        block.sync();
+
+        /*
+         * Finally, set the beginning and end of our children, as well as the
+         * splitting point and the current node.
+         */
+        if (block.thread_rank() == 0) {
+            tree[lnid].begin = tree[nid].begin;
+            tree[lnid].end = tree[lnid].begin + lower_idx;
+
+            tree[rnid].begin = tree[lnid].end;
+            tree[rnid].end = tree[nid].end;
+
+            mid_point_value = mid_point;
+        }
+
+        block.sync();
+    } else {
+        /*
+         * For small ranges, we just get the exact median. We have a helper
+         * function for this which also sorts the indices as a side-effect, and
+         * so we do not actually need to do any partitioning.
+         */
+        uint32_t mid_point = sort_and_get_median(
+            spacepoints, indices, index_buffer, num_indices, pivot);
+
+        /*
+         * Correctly set the beginning and the of our children, as well as our
+         * own pivot value.
+         */
+        if (block.thread_rank() == 0) {
+            tree[lnid].begin = tree[nid].begin;
+            tree[lnid].end = tree[lnid].begin + mid_point;
+
+            tree[rnid].begin = tree[lnid].end;
+            tree[rnid].end = tree[nid].end;
+
+            mid_point_value =
+                get_coordinate(spacepoints, index_buffer[mid_point], pivot);
+        }
+
+        block.sync();
+    }
+
+    /*
+     * Now that the partitioning is complete, we can have our leader thread
+     * put the dots on the i's and complete this internal node.
+     */
+    if (block.thread_rank() == 0) {
+        tree[nid].type = nodetype_e::INTERNAL;
+        tree[nid].dim = pivot;
+        tree[nid].mid = mid_point_value;
+
+        tree[lnid].range = tree[nid].range;
+        tree[rnid].range = tree[nid].range;
+
+        /*
+         * We need to make sure that we update the ranges of our children!
+         */
+        if (tree[nid].dim == pivot_e::Phi) {
+            tree[lnid].range.phi_max = mid_point_value;
+            tree[rnid].range.phi_min = mid_point_value;
+        } else if (tree[nid].dim == pivot_e::R) {
+            tree[lnid].range.r_max = mid_point_value;
+            tree[rnid].range.r_min = mid_point_value;
+        } else {
+            tree[lnid].range.z_max = mid_point_value;
+            tree[rnid].range.z_min = mid_point_value;
+        }
+    }
+
+    /*
+     * This function also always succeeds.
+     */
+    return true;
+}
+
+/**
+ * @brief Main k-d tree construction kernel.
+ *
+ * A k-d tree is constructed by repeated application of this kernel. This
+ * kernel uses a one-node-per-block approach. This is not optimal for early
+ * iteration counts.
+ *
+ * @warning This kernel assumes that the tree has already been initialized
+ * through the preceding kernels.
+ *
+ * @param[inout] tree The tree array in which to build.
+ * @param[in] num_nodes The maximum size of the k-d tree.
+ * @param[in] spacepoints The spacepoints to encode in the tree.
+ * @param[inout] _indices The index array to use.
+ * @param[inout] _index_buffer Temporary buffer space for indices.
+ * @param[in] iteration The iteration count.
+ */
+__global__ void __launch_bounds__(512)
+    construct_kd_tree_big_step(kd_tree_t tree, uint32_t num_nodes,
+                               internal_sp_t spacepoints,
+                               const uint32_t* __restrict__ _indices_old,
+                               uint32_t* __restrict__ _indices_new,
+                               uint32_t iteration,
+                               uint32_t* __restrict__ work_list_current,
+                               uint32_t* __restrict__ work_list_size_current,
+                               uint32_t* __restrict__ work_list_next,
+                               uint32_t* __restrict__ work_list_size_next) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+
+    /*
+     * Calculate the current node ID. Remember that we start with 1 node
+     * indexed at [0], then 2 nodes at [1, 2], then 4 nodes at [3, 4, 5, 6],
+     * etc. Thus, the subarray for zero-indexed iteration i starts at 2^i - 1.
+     */
+    assert(blockIdx.x < *work_list_size_current);
+
+    int nid = work_list_current[blockIdx.x];
+
+    /*
+     * Note that this is not equivalent per se to the above, and that this is
+     * not an error. But it does indicate that there is no work to be done.
+     */
+    assert(tree[nid].type == nodetype_e::LEAF);
+
+    bool build_leaf = (tree[nid].end - tree[nid].begin) <= 8;
+
+    block.sync();
+
+    /*
+     * Check to see whether we could fit the range assigned to this block into
+     * a leaf node. If so, we construct one. Otherwise, we construct an
+     * internal node.
+     */
+    if (build_leaf) {
+        cooperative_groups::thread_block block =
+            cooperative_groups::this_thread_block();
+
+        /*
+         * Set the node to a leaf, and set the number of points it contains.
+         */
+        if (block.thread_rank() == 0) {
+            uint32_t begin = tree[nid].begin, end = tree[nid].end;
+
+            tree[nid].type = nodetype_e::LEAF;
+            tree[nid].begin = begin;
+            tree[nid].end = end;
+        }
+
+        block.sync();
+    } else {
+        [[maybe_unused]] bool r = false;
+
+        pivot_e pivot;
+
+        /*
+         * If we are the root node, our pivot axis is φ; otherwise, it's what-
+         * ever is next after our parent's pivot axis.
+         */
+        if (nid == 0) {
+            pivot = pivot_e::Phi;
+        } else {
+            pivot = next_pivot(tree[(nid - 1) / 2].dim);
+        }
+
+        /*
+         * Try to create an internal node with the given pivot.
+         */
+        r = create_internal_node(tree, spacepoints,
+                                 &_indices_old[tree[nid].begin],
+                                 &_indices_new[tree[nid].begin], nid, pivot);
+
+        /*
+         * Check whether the process succeeded. It should.
+         */
+        assert(r);
+
+        if (threadIdx.x == 0 && blockIdx.y == 0) {
+            uint32_t idx = atomicAdd(work_list_size_next, 2u);
+            work_list_next[idx] = 2 * nid + 1;
+            work_list_next[idx + 1] = 2 * nid + 2;
+        }
+    }
+}
+
+__global__ void __launch_bounds__(96)
+    construct_kd_tree_small_step(kd_tree_t tree, uint32_t num_nodes,
+                                 internal_sp_t spacepoints,
+                                 const uint32_t* __restrict__ _indices_old,
+                                 uint32_t* __restrict__ _indices_new,
+                                 uint32_t iteration,
+                                 uint32_t* __restrict__ work_list_current,
+                                 uint32_t* __restrict__ work_list_size_current,
+                                 uint32_t* __restrict__ work_list_next,
+                                 uint32_t* __restrict__ work_list_size_next) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+
+    /*
+     * Calculate the current node ID. Remember that we start with 1 node
+     * indexed at [0], then 2 nodes at [1, 2], then 4 nodes at [3, 4, 5, 6],
+     * etc. Thus, the subarray for zero-indexed iteration i starts at 2^i - 1.
+     */
+    assert(blockIdx.x < *work_list_size_current);
+
+    int nid = work_list_current[blockIdx.x];
+
+    /*
+     * Note that this is not equivalent per se to the above, and that this is
+     * not an error. But it does indicate that there is no work to be done.
+     */
+    assert(tree[nid].type == nodetype_e::LEAF);
+
+    bool build_leaf = (tree[nid].end - tree[nid].begin) <= 8;
+
+    block.sync();
+
+    /*
+     * Check to see whether we could fit the range assigned to this block into
+     * a leaf node. If so, we construct one. Otherwise, we construct an
+     * internal node.
+     */
+    if (build_leaf) {
+        cooperative_groups::thread_block block =
+            cooperative_groups::this_thread_block();
+
+        /*
+         * Set the node to a leaf, and set the number of points it contains.
+         */
+        if (block.thread_rank() == 0) {
+            uint32_t begin = tree[nid].begin, end = tree[nid].end;
+
+            tree[nid].type = nodetype_e::LEAF;
+            tree[nid].begin = begin;
+            tree[nid].end = end;
+        }
+
+        block.sync();
+    } else {
+        [[maybe_unused]] bool r = false;
+
+        pivot_e pivot;
+
+        /*
+         * If we are the root node, our pivot axis is φ; otherwise, it's what-
+         * ever is next after our parent's pivot axis.
+         */
+        if (nid == 0) {
+            pivot = pivot_e::Phi;
+        } else {
+            pivot = next_pivot(tree[(nid - 1) / 2].dim);
+        }
+
+        /*
+         * Try to create an internal node with the given pivot.
+         */
+        r = create_internal_node(tree, spacepoints,
+                                 &_indices_old[tree[nid].begin],
+                                 &_indices_new[tree[nid].begin], nid, pivot);
+
+        /*
+         * Check whether the process succeeded. It should.
+         */
+        assert(r);
+
+        if (threadIdx.x == 0 && blockIdx.y == 0) {
+            uint32_t idx = atomicAdd(work_list_size_next, 2u);
+            work_list_next[idx] = 2 * nid + 1;
+            work_list_next[idx + 1] = 2 * nid + 2;
+        }
+    }
+}
+
+uint32_t round_to_power_2(uint32_t i) {
+    uint32_t power = 1u;
+
+    while (power < i) {
+        power *= 2u;
+    }
+
+    return power;
+}
+
+__global__ void bake_spacepoints_kernel(uint32_t n, internal_sp_t oldsps,
+                                        internal_sp_t newsps,
+                                        uint32_t* __restrict__ indices) {
+    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
+
+    uint32_t i = grid.thread_rank();
+
+    if (i < n) {
+        uint32_t ni = indices[i] & ~AUXILIARY_INDEX_MASK;
+
+        newsps[i].x = oldsps[ni].x;
+        newsps[i].y = oldsps[ni].y;
+        newsps[i].z = oldsps[ni].z;
+        newsps[i].phi = oldsps[ni].phi;
+        newsps[i].radius = oldsps[ni].radius;
+        newsps[i].link = oldsps[ni].link;
+    }
+}
+
+std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
+    vecmem::memory_resource& mr, internal_sp_owning_t&& spacepoints,
+    uint32_t n_sp) {
+    /*
+     * Allocate space for the indices. Since each spacepoint can produce at
+     * most two indices, we allocate enough space to support this theoretical
+     * upper bound.
+     */
+    vecmem::unique_alloc_ptr<uint32_t[]> indices =
+        vecmem::make_unique_alloc<uint32_t[]>(mr, 4 * n_sp);
+
+    /*
+     * Allocate an on-device counter for the number of halo points we have,
+     * which starts out at zero.
+     */
+    vecmem::unique_alloc_ptr<uint32_t> extra_indices =
+        vecmem::make_unique_alloc<uint32_t>(mr);
+
+    CUDA_ERROR_CHECK(cudaMemset(extra_indices.get(), 0, sizeof(uint32_t)));
+
+    /*
+     * Launch the index initialization kernel, which turns n spacepoints into
+     * anywhere between n and 2n indices in the index array.
+     */
+    std::size_t threads_per_block1 = 256;
+    initialize_index_vector<<<(n_sp / threads_per_block1) +
+                                  (n_sp % threads_per_block1 == 0 ? 0 : 1),
+                              threads_per_block1>>>(spacepoints, indices.get(),
+                                                    n_sp, extra_indices.get());
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    /*
+     * Retrieve the number of halo points from the device.
+     */
+    uint32_t extra_indices_h;
+
+    CUDA_ERROR_CHECK(cudaMemcpy(&extra_indices_h, extra_indices.get(),
+                                sizeof(uint32_t), cudaMemcpyDeviceToHost));
+
+    /*
+     * The total number of indices is equal to the number of spacepoints plus
+     * the number of halo points.
+     */
+    uint32_t num_indices = n_sp + extra_indices_h;
+
+    /*
+     * In theory, the number of nodes that is needed to support N points given
+     * K points per node is 2^⌈log_2⌈N / ⌊K / 2⌋⌉ + 1⌉ - 1. In practice, we
+     * build potentially _very_ unbalanced trees. Thus, we need a bit more
+     * space.
+     */
+    uint32_t largest_layer = round_to_power_2(num_indices);
+    uint32_t num_nodes = 4 * largest_layer - 1;
+
+    /*
+     * Allocate space in which the k-d tree will live.
+     */
+    kd_tree_owning_t tree_owner(mr, num_nodes);
+
+    vecmem::unique_alloc_ptr<uint32_t[]> work_list_1 =
+        vecmem::make_unique_alloc<uint32_t[]>(mr, num_nodes);
+    vecmem::unique_alloc_ptr<uint32_t[]> work_list_2 =
+        vecmem::make_unique_alloc<uint32_t[]>(mr, num_nodes);
+    vecmem::unique_alloc_ptr<uint32_t> work_list_size_1 =
+        vecmem::make_unique_alloc<uint32_t>(mr);
+    vecmem::unique_alloc_ptr<uint32_t> work_list_size_2 =
+        vecmem::make_unique_alloc<uint32_t>(mr);
+
+    CUDA_ERROR_CHECK(cudaMemset(work_list_size_1.get(), 0, sizeof(uint32_t)));
+    CUDA_ERROR_CHECK(cudaMemset(work_list_size_2.get(), 0, sizeof(uint32_t)));
+
+    /*
+     * Run the initialization kernel for the tree itself, which sets the
+     * initial types of the nodes.
+     */
+    std::size_t threads_per_block2 = 256;
+    initialize_kd_tree<<<(num_nodes / threads_per_block2) +
+                             (num_nodes % threads_per_block2 == 0 ? 0 : 1),
+                         threads_per_block2>>>(kd_tree_t(tree_owner), num_nodes,
+                                               num_indices, work_list_1.get(),
+                                               work_list_size_1.get());
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    /*
+     * Allocate some temporary space for the index buffer, one half for "lower"
+     * points, and one for "upper" points, defined relative to some pivot
+     * point. Multiple nodes may use this buffer to do their partitions at the
+     * same time.
+     */
+    vecmem::unique_alloc_ptr<uint32_t[]> index_buffer =
+        vecmem::make_unique_alloc<uint32_t[]>(mr, num_indices);
+
+    /*
+     * Repeatedly apply the k-d tree construction kernel. The iteration width
+     * is the number of nodes that is processed in parallel, and is simply 2^i
+     * where i is the iteration. We can stop running kernels after we cover the
+     * entire last level of the tree.
+     */
+    uint32_t remaining = 1, iteration = 0;
+
+    while (remaining > 0) {
+        dim3 grid_size{remaining};
+        uint32_t threads_per_block;
+
+        if (iteration <= 4) {
+            threads_per_block = 512u;
+
+            construct_kd_tree_big_step<<<grid_size, threads_per_block,
+                                         threads_per_block * sizeof(float)>>>(
+                kd_tree_t(tree_owner), num_nodes, spacepoints, indices.get(),
+                index_buffer.get(), iteration, work_list_1.get(),
+                work_list_size_1.get(), work_list_2.get(),
+                work_list_size_2.get());
+        } else {
+            threads_per_block = 96u;
+
+            construct_kd_tree_small_step<<<grid_size, threads_per_block,
+                                           threads_per_block * sizeof(float)>>>(
+                kd_tree_t(tree_owner), num_nodes, spacepoints, indices.get(),
+                index_buffer.get(), iteration, work_list_1.get(),
+                work_list_size_1.get(), work_list_2.get(),
+                work_list_size_2.get());
+        }
+
+        CUDA_ERROR_CHECK(cudaGetLastError());
+        CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+        iteration++;
+
+        CUDA_ERROR_CHECK(cudaMemcpy(&remaining, work_list_size_2.get(),
+                                    sizeof(uint32_t), cudaMemcpyDeviceToHost));
+
+        std::swap(work_list_1, work_list_2);
+        std::swap(work_list_size_1, work_list_size_2);
+        std::swap(indices, index_buffer);
+
+        CUDA_ERROR_CHECK(
+            cudaMemset(work_list_size_2.get(), 0, sizeof(uint32_t)));
+    }
+
+    internal_sp_owning_t new_sps(mr, num_indices);
+
+    bake_spacepoints_kernel<<<
+        num_indices / 1024u + (num_indices % 1024u == 0u ? 0u : 1u), 1024u>>>(
+        num_indices, internal_sp_t(spacepoints), internal_sp_t(new_sps),
+        indices.get());
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    return {kd_tree_owning_t{std::move(tree_owner)}, num_nodes,
+            std::move(new_sps)};
+}
+}  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu b/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
new file mode 100644
index 0000000000..9e1204e7fd
--- /dev/null
+++ b/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
@@ -0,0 +1,706 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include <cooperative_groups.h>
+
+#include <iostream>
+#include <traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp>
+#include <traccc/cuda/seeding2/types/internal_sp.hpp>
+#include <traccc/cuda/seeding2/types/kd_tree.hpp>
+#include <traccc/cuda/seeding2/types/range3d.hpp>
+#include <traccc/cuda/utils/definitions.hpp>
+#include <traccc/cuda/utils/device_traits.hpp>
+#include <traccc/cuda/utils/sort.hpp>
+#include <traccc/cuda/utils/sync.hpp>
+#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/seeding/detail/lin_circle.hpp>
+#include <traccc/seeding/doublet_finding_helper.hpp>
+#include <traccc/seeding/triplet_finding_helper.hpp>
+
+#define MAX_LOWER_SP_PER_MIDDLE 100u
+#define MAX_UPPER_SP_PER_MIDDLE 100u
+#define KD_TREE_TRAVERSAL_STACK_SIZE 128u
+#define WARPS_PER_BLOCK 8u
+
+/**
+ * @brief Maximum difference in φ between two spacepoints adjacent in the same
+ * seed.
+ */
+static constexpr float MAX_DELTA_PHI = 0.025f;
+
+namespace {
+struct internal_seed {
+    uint32_t spacepoints[2];
+    float weight;
+
+    __host__ __device__ static internal_seed Invalid() {
+        internal_seed r;
+
+        r.weight = std::numeric_limits<float>::lowest();
+
+        return r;
+    }
+};
+}  // namespace
+
+namespace traccc::cuda {
+/**
+ * @brief Retrieve values from a k-d tree in a given range and write them to
+ * an output array.
+ *
+ * Performs depth-first traversal of the k-d tree, looking for candidate
+ * spacepoints.
+ *
+ * @param[in] range 3D range in which to search.
+ * @param[in] mi_idx Index of the middle spacepoint.
+ * @param[in] upper True iff we're looking for upper spacepoints.
+ * @param[out] output_arr Output array to write candidates to.
+ * @param[out] output_cnt Total number of output candidates.
+ */
+__device__ void retrieve_from_tree(
+    const seedfinder_config finder_conf, const seedfilter_config filter_conf,
+    const internal_sp_t spacepoints, const kd_tree_t tree,
+    const internal_spacepoint<spacepoint> mi, const range3d range,
+    uint32_t mi_idx, bool upper, uint32_t* __restrict__ output_arr,
+    uint32_t* __restrict__ output_cnt) {
+    __shared__ uint32_t stack[WARPS_PER_BLOCK][KD_TREE_TRAVERSAL_STACK_SIZE];
+    __shared__ uint32_t leaves[WARPS_PER_BLOCK][WARP_SIZE];
+
+    uint32_t stack_idx = 32;
+
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+    cooperative_groups::thread_block_tile<WARP_SIZE> warp =
+        cooperative_groups::tiled_partition<WARP_SIZE>(block);
+
+    uint32_t idx = WARP_SIZE - 1 + warp.thread_rank();
+    assert(tree[idx].type == nodetype_e::INTERNAL);
+    stack[warp.meta_group_rank()][warp.thread_rank()] = idx;
+
+    warp.sync();
+
+    /*
+     * The traversal runs until the stack is empty. Or until the entire kernel
+     * explodes and crashes!
+     */
+    while (stack_idx > 0) {
+        /*
+         * We perform this traversal in parallel, so multiple threads pop a
+         * value from the stack at the same time. Of course, this can only
+         * happen if the thread index fits within the current size of the
+         * stack.
+         */
+        uint32_t old_stack_idx = stack_idx;
+
+        stack_idx -= std::min(stack_idx, warp.size());
+
+        uint32_t current;
+        bool add_children = false, add_spacepoints = false;
+
+        if (warp.thread_rank() < old_stack_idx) {
+            current = stack[warp.meta_group_rank()]
+                           [old_stack_idx - (warp.thread_rank() + 1)];
+
+            /*
+             * The node must be a leaf node or an internal node, otherwise we
+             * cannot (safely) traverse it.
+             */
+            assert(tree[current].type == nodetype_e::LEAF ||
+                   tree[current].type == nodetype_e::INTERNAL);
+
+            const range3d& crange = tree[current].range;
+
+            /*
+             * Choose the appropriate code path based on whether the node is a
+             * leaf or an internal node.
+             */
+            add_spacepoints = tree[current].type == nodetype_e::LEAF ||
+                              range.dominates(crange);
+            add_children = tree[current].type == nodetype_e::INTERNAL &&
+                           range.intersects(crange);
+        }
+
+        auto [vst, vsi] = warp_indexed_ballot_sync(add_spacepoints);
+
+        if (add_spacepoints) {
+            uint32_t beg = tree[current].begin, end = tree[current].end;
+
+            for (uint32_t point_id = beg; point_id < end; point_id++) {
+                /*
+                 * Ensure that the point is actually valid using cuts which
+                 * are more specific than those encoded in the search
+                 * range.
+                 */
+                float phi = spacepoints[point_id].phi,
+                      radius = spacepoints[point_id].radius,
+                      z = spacepoints[point_id].z;
+
+                if (range.contains(phi, radius, z) && point_id != mi_idx &&
+                    doublet_finding_helper::isCompatible(
+                        upper, mi,
+                        internal_spacepoint<spacepoint>(
+                            spacepoints[point_id].x, spacepoints[point_id].y, z,
+                            radius, phi, spacepoints[point_id].link),
+                        finder_conf)) {
+                    /*
+                     * Reserve a spot in the output for this point, and
+                     * then store it if it would not exceed the output
+                     * array size.
+                     */
+                    uint32_t idx = atomicAdd(output_cnt, 1u);
+
+                    assert((upper && idx < MAX_UPPER_SP_PER_MIDDLE) ||
+                           (!upper && idx < MAX_LOWER_SP_PER_MIDDLE));
+
+                    output_arr[idx] = point_id;
+                }
+            }
+        }
+
+        auto [vct, vci] = warp_indexed_ballot_sync(add_children);
+
+        if (add_children) {
+            uint32_t idx = stack_idx + vci * 2u;
+            assert(idx + 1 < KD_TREE_TRAVERSAL_STACK_SIZE);
+
+            stack[warp.meta_group_rank()][idx] = 2 * current + 1;
+            stack[warp.meta_group_rank()][idx + 1] = 2 * current + 2;
+        }
+
+        stack_idx += 2u * vct;
+
+        /*
+         * Enforce synchronization of the block between DFS steps.
+         */
+        warp.sync();
+    }
+}
+
+/**
+ * @brief Create an internal seed from three spacepoints.
+ *
+ * @param[in] spacepoints The array of spacepoints.
+ * @param[in] lower The index of the lower spacepoint.
+ * @param[in] middle The index of the middle spacepoint.
+ * @param[in] upper The index of the upper spacepoint.
+ *
+ * @return An internal seed iff the triplet is valid, otherwise an invalid
+ * seed.
+ */
+__device__ internal_seed make_seed(
+    const seedfinder_config finder_conf, const seedfilter_config filter_conf,
+    const internal_sp_t spacepoints, const kd_tree_t tree, uint32_t lower,
+    uint32_t middle, internal_spacepoint<spacepoint> mi, uint32_t upper) {
+    const internal_spacepoint<spacepoint> lo(
+        spacepoints[lower].x, spacepoints[lower].y, spacepoints[lower].z,
+        spacepoints[lower].radius, spacepoints[lower].phi,
+        spacepoints[lower].link);
+    const internal_spacepoint<spacepoint> hi(
+        spacepoints[upper].x, spacepoints[upper].y, spacepoints[upper].z,
+        spacepoints[upper].radius, spacepoints[upper].phi,
+        spacepoints[upper].link);
+
+    /*
+     * Find the lin-circles for the bottom and top pair.
+     */
+    lin_circle lm = doublet_finding_helper::transform_coordinates<
+        details::spacepoint_type::bottom>(mi, lo);
+    lin_circle mh = doublet_finding_helper::transform_coordinates<
+        details::spacepoint_type::top>(mi, hi);
+
+    scalar iSinTheta2 = 1 + lm.cotTheta() * lm.cotTheta();
+    scalar scatteringInRegion2 = finder_conf.maxScatteringAngle2 * iSinTheta2;
+    scatteringInRegion2 *=
+        finder_conf.sigmaScattering * finder_conf.sigmaScattering;
+    scalar curvature, impact_parameter;
+
+    /*
+     * Borrow the compatibility code from the existing seed finding code.
+     */
+    if (triplet_finding_helper::isCompatible(mi, lm, mh, finder_conf,
+                                             iSinTheta2, scatteringInRegion2,
+                                             curvature, impact_parameter)) {
+        internal_seed r;
+
+        r.spacepoints[0] = lower;
+        r.spacepoints[1] = upper;
+
+        /*
+         * Add a weight based on the impact parameter to the seed.
+         */
+        r.weight = -impact_parameter * filter_conf.impactWeightFactor;
+
+        return r;
+    } else {
+        /*
+         * If the triplet is invalid, return a bogus seed.
+         */
+        return internal_seed::Invalid();
+    }
+}
+
+/**
+ * @brief Get the basic Range3D object
+ *
+ * This contains only the most basic cuts, but is designed to be refined later.
+ *
+ * @param[in] finder_conf The seed finder configuration to use.
+ *
+ * @return A 3D range object.
+ */
+__device__ range3d get_basic_range3d(const seedfinder_config finder_conf) {
+    range3d r = range3d::Infinite();
+
+    r.r_min = 0.f;
+    r.r_max = finder_conf.rMax;
+
+    r.z_min = finder_conf.zMin;
+    r.z_max = finder_conf.zMax;
+
+    return r;
+}
+
+/**
+ * @brief Return the search range for a lower spacepoint given a middle
+ * spacepoint.
+ *
+ * @param[in] finder_conf The seed finder configuration to use.
+ * @param[in] s The middle spacepoint.
+ *
+ * @return A three-dimensional search range.
+ */
+__device__ range3d get_lower_range3d(const seedfinder_config finder_conf,
+                                     const internal_spacepoint<spacepoint> s) {
+    range3d r = get_basic_range3d(finder_conf);
+
+    r.r_min = std::max(r.r_min, s.radius() - finder_conf.deltaRMax);
+    r.r_max = std::min(r.r_max, s.radius() - finder_conf.deltaRMin);
+
+    float frac_r = r.r_min / s.radius();
+
+    float z_min2 = (s.z() - finder_conf.collisionRegionMin) * frac_r +
+                   finder_conf.collisionRegionMin;
+    float z_max2 = (s.z() - finder_conf.collisionRegionMax) * frac_r +
+                   finder_conf.collisionRegionMax;
+
+    r.z_min = std::max(r.z_min, std::min(z_min2, s.z()));
+    r.z_max = std::min(r.z_max, std::max(z_max2, s.z()));
+
+    r.z_min = std::max(r.z_min, s.z() - 450.f);
+    r.z_max = std::min(r.z_max, s.z() + 450.f);
+
+    r.phi_min = std::max(r.phi_min, s.phi() - MAX_DELTA_PHI);
+    r.phi_max = std::min(r.phi_max, s.phi() + MAX_DELTA_PHI);
+
+    return r;
+}
+
+/**
+ * @brief Return the search range for an upper spacepoint given a middle
+ * spacepoint.
+ *
+ * @param[in] finder_conf The seed finder configuration to use.
+ * @param[in] s The middle spacepoint.
+ *
+ * @return A three-dimensional search range.
+ */
+__device__ range3d get_upper_range3d(const seedfinder_config finder_conf,
+                                     const internal_spacepoint<spacepoint> s) {
+    range3d r = get_basic_range3d(finder_conf);
+
+    r.r_min = std::max(r.r_min, s.radius() + finder_conf.deltaRMin);
+    r.r_max = std::min(r.r_max, s.radius() + finder_conf.deltaRMax);
+
+    float z_max2 =
+        (r.r_max / s.radius()) * (s.z() - finder_conf.collisionRegionMin) +
+        finder_conf.collisionRegionMin;
+    float z_min2 =
+        finder_conf.collisionRegionMax -
+        (r.r_max / s.radius()) * (finder_conf.collisionRegionMax - s.z());
+
+    if (s.z() > finder_conf.collisionRegionMin) {
+        r.z_max = std::min(r.z_max, z_max2);
+    } else if (s.z() < finder_conf.collisionRegionMax) {
+        r.z_min = std::max(r.z_min, z_min2);
+    }
+
+    r.z_min = std::max(
+        r.z_min, s.z() - finder_conf.cotThetaMax * (r.r_max - s.radius()));
+    r.z_max = std::min(
+        r.z_max, s.z() + finder_conf.cotThetaMax * (r.r_max - s.radius()));
+
+    r.z_min = std::max(r.z_min, s.z() - 450.f);
+    r.z_max = std::min(r.z_max, s.z() + 450.f);
+
+    r.phi_min = std::max(r.phi_min, s.phi() - MAX_DELTA_PHI);
+    r.phi_max = std::min(r.phi_max, s.phi() + MAX_DELTA_PHI);
+
+    return r;
+}
+
+/**
+ * @brief Seed finding helper function.
+ *
+ * Turns out this function actually does most of the important work.
+ *
+ * @param[inout] internal_seeds The array to write internal seeds to.
+ * @param[in] increasing_z True iff we are looking for increasing-z seeds.
+ */
+__device__ void run_helper(const seedfinder_config finder_conf,
+                           const seedfilter_config filter_conf,
+                           const internal_sp_t spacepoints,
+                           const kd_tree_t tree,
+                           const internal_spacepoint<spacepoint> sp,
+                           std::size_t sp_idx, internal_seed* internal_seeds,
+                           bool increasing_z) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+    cooperative_groups::thread_block_tile<WARP_SIZE> warp =
+        cooperative_groups::tiled_partition<WARP_SIZE>(block);
+
+    /*
+     * Allocate shared memory for the lower and upper candidates.
+     */
+    __shared__ uint32_t lower_sps[WARPS_PER_BLOCK][MAX_LOWER_SP_PER_MIDDLE];
+    __shared__ uint32_t upper_sps[WARPS_PER_BLOCK][MAX_UPPER_SP_PER_MIDDLE];
+
+    /*
+     * The leader thread initializes the candidate counts to zero.
+     */
+    __shared__ uint32_t num_lower[WARPS_PER_BLOCK], num_upper[WARPS_PER_BLOCK];
+
+    if (warp.thread_rank() == 0) {
+        num_lower[warp.meta_group_rank()] = 0;
+        num_upper[warp.meta_group_rank()] = 0;
+    }
+
+    warp.sync();
+
+    /*
+     * Retrieve the search range for the lower spacepoints, then adjust it for
+     * monotinicity in z.
+     */
+    range3d lower_range3d = get_lower_range3d(finder_conf, sp);
+
+    if (increasing_z) {
+        lower_range3d.z_max = std::min(lower_range3d.z_max, sp.z());
+    } else {
+        lower_range3d.z_min = std::max(lower_range3d.z_min, sp.z());
+    }
+
+    /*
+     * Same as above, but for the upper spacepoint. Note that some of the
+     * operations are reversed.
+     */
+    range3d upper_range3d = get_upper_range3d(finder_conf, sp);
+
+    if (increasing_z) {
+        upper_range3d.z_min = std::max(lower_range3d.z_min, sp.z());
+    } else {
+        upper_range3d.z_max = std::min(lower_range3d.z_max, sp.z());
+    }
+
+    /*
+     * If either of the search ranges is degenerate, we can never find any
+     * candidate seeds. So we can safely exit.
+     */
+    if (lower_range3d.phi_min >= lower_range3d.phi_max ||
+        lower_range3d.r_min >= lower_range3d.r_max ||
+        lower_range3d.z_min >= lower_range3d.z_max ||
+        upper_range3d.phi_min >= upper_range3d.phi_max ||
+        upper_range3d.r_min >= upper_range3d.r_max ||
+        upper_range3d.z_min >= upper_range3d.z_max) {
+        return;
+    }
+
+    /*
+     * Retrieve the lower and upper spacepoint candidates from the k-d tree.
+     */
+    retrieve_from_tree(finder_conf, filter_conf, spacepoints, tree, sp,
+                       lower_range3d, sp_idx, false,
+                       lower_sps[warp.meta_group_rank()],
+                       &num_lower[warp.meta_group_rank()]);
+
+    warp.sync();
+
+    /*
+     * If we have zero lower candidates, we might as well quit; we will never
+     * have any candidates!
+     */
+    if (num_lower[warp.meta_group_rank()] == 0) {
+        return;
+    }
+
+    retrieve_from_tree(finder_conf, filter_conf, spacepoints, tree, sp,
+                       upper_range3d, sp_idx, true,
+                       upper_sps[warp.meta_group_rank()],
+                       &num_upper[warp.meta_group_rank()]);
+
+    warp.sync();
+
+    /*
+     * Life fast and die young! (It's an early quitting condition)
+     */
+    if (num_upper[warp.meta_group_rank()] == 0) {
+        return;
+    }
+
+    /*
+     * Calculate the total number of iterations, but take into account the
+     * maximum size of the lower and upper spacepoint storage arrays.
+     */
+    uint32_t combinations =
+        std::min(num_lower[warp.meta_group_rank()], MAX_LOWER_SP_PER_MIDDLE) *
+        std::min(num_upper[warp.meta_group_rank()], MAX_UPPER_SP_PER_MIDDLE);
+
+    /*
+     * The iteration count is rounded up to a multiple of the block size to
+     * ensure that none of the threads slack off and exit prematurely.
+     */
+    uint32_t iter_count =
+        combinations + (combinations % WARP_SIZE != 0
+                            ? WARP_SIZE - (combinations % WARP_SIZE)
+                            : 0);
+
+    float min_weight = std::numeric_limits<float>::lowest();
+
+    for (uint32_t i = warp.thread_rank(); i < iter_count; i += WARP_SIZE) {
+        uint32_t lower_id = i / num_upper[warp.meta_group_rank()];
+        uint32_t upper_id = i % num_upper[warp.meta_group_rank()];
+
+        /*
+         * Consider the current combination and turn it into a seed. This seed
+         * may be bogus! But that is okay, because the sorting step below will
+         * automatically separate bogus seeds.
+         */
+        bool add_seed = false;
+        internal_seed s;
+
+        if (lower_id < std::min(num_lower[warp.meta_group_rank()],
+                                MAX_LOWER_SP_PER_MIDDLE) &&
+            upper_id < std::min(num_upper[warp.meta_group_rank()],
+                                MAX_UPPER_SP_PER_MIDDLE)) {
+            s = make_seed(finder_conf, filter_conf, spacepoints, tree,
+                          lower_sps[warp.meta_group_rank()][lower_id], sp_idx,
+                          sp, upper_sps[warp.meta_group_rank()][upper_id]);
+
+            if (s.weight > min_weight) {
+                // internal_seeds[warp.thread_rank()] = s;
+                add_seed = true;
+            }
+        }
+
+        auto [vt, vi] = warp_indexed_ballot_sync(add_seed);
+
+        if (add_seed) {
+            internal_seeds[WARP_SIZE - (1 + vi)] = s;
+        }
+
+        warp.sync();
+
+        if (vt > 0) {
+            /*
+             * Perform odd-even sorting, making sure that the seeds with the
+             * highest weight float to the top of the seed array, which is the
+             * part that will be written to the global output.
+             */
+            warpOddEvenKeySort(
+                internal_seeds, WARP_SIZE + finder_conf.maxSeedsPerSpM,
+                [](const internal_seed& a, const internal_seed& b) {
+                    return a.weight < b.weight;
+                });
+
+            min_weight = internal_seeds[WARP_SIZE].weight;
+        }
+    }
+}
+
+/**
+ * @brief Main seed finding kernel.
+ *
+ * @note This kernel assumes exactly one warp per block and operates in a
+ * one-block-per-middle-SP fashion.
+ *
+ */
+__global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
+    seed_finding_kernel(const seedfinder_config finder_conf,
+                        const seedfilter_config filter_conf,
+                        const internal_sp_t spacepoints, const kd_tree_t tree,
+                        alt_seed* output_seeds, uint32_t* output_seed_size) {
+    cooperative_groups::thread_block block =
+        cooperative_groups::this_thread_block();
+    cooperative_groups::thread_block_tile<WARP_SIZE> warp =
+        cooperative_groups::tiled_partition<WARP_SIZE>(block);
+
+    assert(warp.meta_group_size() == WARPS_PER_BLOCK);
+
+    std::size_t spacepoint_id =
+        block.group_index().x * WARPS_PER_BLOCK + warp.meta_group_rank();
+
+    /*
+     * The seed finding configuration defines a few criteria that allow us to
+     * ignore some of the middle spacepoint candidates.
+     */
+    if (spacepoint_id >= spacepoints.size() ||
+        spacepoints[spacepoint_id].phi < finder_conf.phiMin ||
+        spacepoints[spacepoint_id].phi > finder_conf.phiMax ||
+        spacepoints[spacepoint_id].radius < finder_conf.rMin ||
+        spacepoints[spacepoint_id].radius > finder_conf.rMax ||
+        spacepoints[spacepoint_id].z < finder_conf.zMin ||
+        spacepoints[spacepoint_id].z > finder_conf.zMax) {
+        return;
+    }
+
+    const internal_spacepoint<spacepoint> sp(
+        spacepoints[spacepoint_id].x, spacepoints[spacepoint_id].y,
+        spacepoints[spacepoint_id].z, spacepoints[spacepoint_id].radius,
+        spacepoints[spacepoint_id].phi, spacepoints[spacepoint_id].link);
+
+    /*
+     * The internal seeds are stored in this shared array.
+     */
+    extern __shared__ internal_seed _internal_seeds[];
+    internal_seed* internal_seeds =
+        &_internal_seeds[warp.meta_group_rank() *
+                         (WARP_SIZE + finder_conf.maxSeedsPerSpM)];
+
+    /*
+     * Initialize the internal seed array to contain only bogus seeds.
+     */
+    for (uint32_t i = warp.thread_rank();
+         i < WARP_SIZE + finder_conf.maxSeedsPerSpM; i += WARP_SIZE) {
+        internal_seeds[i].weight = std::numeric_limits<float>::lowest();
+    }
+
+    warp.sync();
+
+    /*
+     * Run the actual seed finding twice; once with monotonically increasing
+     * z values, and once with monotonically decreasing ones.
+     */
+    run_helper(finder_conf, filter_conf, spacepoints, tree, sp, spacepoint_id,
+               internal_seeds, true);
+    warp.sync();
+    run_helper(finder_conf, filter_conf, spacepoints, tree, sp, spacepoint_id,
+               internal_seeds, false);
+    warp.sync();
+
+    /*
+     * We have now completed the seeding algorithm, and we finish by writing
+     * the seeds to our output array. This task is reserved for the leader
+     * thread only, and all the other threads are idle. Thankfully, this should
+     * not take much time at all.
+     */
+    if (warp.thread_rank() == 0) {
+        /*
+         * First, we count the number of valid seeds that we will want to write
+         * to the global memory.
+         */
+        uint32_t num_valid = 0;
+
+        for (uint32_t i = 0; i < finder_conf.maxSeedsPerSpM; ++i) {
+            if (internal_seeds[WARP_SIZE + i].weight > -1000000.0f) {
+                ++num_valid;
+            }
+        }
+
+        /*
+         * Next, reserve space atomically in the output array for our new
+         * seeds.
+         */
+        uint32_t idx = atomicAdd(output_seed_size, num_valid);
+
+        /*
+         * Finally, actually write the output to global memory. Again, this is
+         * done sequentially by the leader thread, but since the number of
+         * seeds is so small it should be fine. Notice that we count down from
+         * the end of the array, because valid seeds are always sorted higher
+         * than invalid ones (if we have any).
+         */
+        for (uint32_t i = 0; i < num_valid; ++i) {
+            alt_seed& s = output_seeds[idx + i];
+            uint32_t j = WARP_SIZE + finder_conf.maxSeedsPerSpM - (i + 1);
+
+            /*
+             * Write the data back to the output.
+             */
+            s.spB_link = spacepoints[internal_seeds[j].spacepoints[0]].link;
+            s.spM_link = spacepoints[spacepoint_id].link;
+            s.spT_link = spacepoints[internal_seeds[j].spacepoints[1]].link;
+
+            s.weight = internal_seeds[i].weight;
+            s.z_vertex = 0.f;
+        }
+    }
+}
+
+/**
+ * @brief Seeding entry point for C++ code.
+ *
+ * This is basically a wrapper for the seed finding kernel.
+ *
+ * @param[in] finder_conf The seed finder configuration to use.
+ * @param[in] filter_conf The seed filter configuration to use.
+ * @param[inout] mr The memory resource for allocations.
+ * @param[in] spacepoints The spacepoint array.
+ * @param[in] tree The k-d tree to use for searching.
+ *
+ * @return A unique pointer to an array of internal seeds, and the size of that
+ * array.
+ */
+std::pair<vecmem::unique_alloc_ptr<alt_seed[]>, uint32_t> run_seeding(
+    seedfinder_config finder_conf, seedfilter_config filter_conf,
+    vecmem::memory_resource& mr, const internal_sp_t spacepoints,
+    const kd_tree_t tree) {
+    /*
+     * Allocate space for output of seeds on the device.
+     */
+    vecmem::unique_alloc_ptr<alt_seed[]> seeds_device =
+        vecmem::make_unique_alloc<alt_seed[]>(
+            mr, finder_conf.maxSeedsPerSpM * spacepoints.size());
+
+    /*
+     * Allocate space for seed count on the device.
+     */
+    vecmem::unique_alloc_ptr<uint32_t> seed_count_device =
+        vecmem::make_unique_alloc<uint32_t>(mr);
+    CUDA_ERROR_CHECK(cudaMemset(seed_count_device.get(), 0, sizeof(uint32_t)));
+
+    /*
+     * Calculate the total amount of shared memory on top of that which is
+     * statically defined in the CUDA kernels. This is equal to one seed for
+     * each block plus one seed for each output seed.
+     */
+    std::size_t total_shared_memory =
+        (WARPS_PER_BLOCK * (WARP_SIZE + finder_conf.maxSeedsPerSpM)) *
+        sizeof(internal_seed);
+    std::size_t grid_size = (spacepoints.size() / WARPS_PER_BLOCK) +
+                            (spacepoints.size() % WARPS_PER_BLOCK == 0 ? 0 : 1);
+    std::size_t block_size = WARPS_PER_BLOCK * WARP_SIZE;
+
+    /*
+     * Call the kernel and all that jazz.
+     */
+    seed_finding_kernel<<<grid_size, block_size, total_shared_memory>>>(
+        finder_conf, filter_conf, spacepoints, tree, seeds_device.get(),
+        seed_count_device.get());
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    /*
+     * Transfer the seed count back to the host and then hand it to the user.
+     */
+    uint32_t seed_count_host;
+    CUDA_ERROR_CHECK(cudaMemcpy(&seed_count_host, seed_count_device.get(),
+                                sizeof(uint32_t), cudaMemcpyDeviceToHost));
+
+    return {std::move(seeds_device), seed_count_host};
+}
+}  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/kernels/write_output_kernel.cu b/device/cuda/src/seeding2/kernels/write_output_kernel.cu
new file mode 100644
index 0000000000..d9bb4a73f2
--- /dev/null
+++ b/device/cuda/src/seeding2/kernels/write_output_kernel.cu
@@ -0,0 +1,82 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include <cooperative_groups.h>
+
+#include <cstdint>
+#include <traccc/cuda/seeding2/kernels/write_output_kernel.hpp>
+#include <traccc/cuda/seeding2/types/internal_sp.hpp>
+#include <traccc/cuda/utils/definitions.hpp>
+#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/seed.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/utils/memory_resource.hpp>
+#include <vecmem/containers/device_vector.hpp>
+#include <vecmem/containers/vector.hpp>
+#include <vecmem/utils/cuda/copy.hpp>
+
+namespace traccc::cuda {
+/**
+ * @brief Main output writing kernel.
+ *
+ * @param[in] n_seeds The number of seeds to write.
+ * @param[in] spacepoints Spacepoint array into which seeds are indexed.
+ * @param[in] in_seeds Input internal seeds to read from.
+ * @param[out] out_seeds Output seeds to write to.
+ */
+__global__ void write_output_kernel(uint32_t n_seeds,
+                                    const internal_sp_t spacepoints,
+                                    const alt_seed* const in_seeds,
+                                    alt_seed_collection_types::view out_seeds) {
+    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
+
+    /*
+     * Create a writable device vector for the seeds.
+     */
+    vecmem::device_vector<alt_seed> dvec(out_seeds);
+
+    int i = grid.thread_rank();
+
+    if (i < n_seeds) {
+        /*
+         * Get the output seed for the current thread.
+         */
+        dvec.at(i) = in_seeds[i];
+    }
+}
+
+alt_seed_collection_types::buffer write_output(
+    const traccc::memory_resource& mr, uint32_t n_seeds,
+    const internal_sp_t spacepoints, const alt_seed* const seeds) {
+    static constexpr std::size_t threads_per_block = 256;
+
+    /*
+     * Create a "vector buffer" that we can write our output to.
+     */
+    alt_seed_collection_types::buffer out(n_seeds, mr.main);
+
+    /*
+     * Fetch the current CUDA device number.
+     */
+    int device;
+    CUDA_ERROR_CHECK(cudaGetDevice(&device));
+
+    /*
+     * Execute the main output kernel.
+     */
+    write_output_kernel<<<(n_seeds / threads_per_block) +
+                              (n_seeds % threads_per_block == 0 ? 0 : 1),
+                          threads_per_block>>>(n_seeds, spacepoints, seeds,
+                                               out);
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    return out;
+}
+}  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/seed_finding2.cu b/device/cuda/src/seeding2/seed_finding2.cu
new file mode 100644
index 0000000000..0837873889
--- /dev/null
+++ b/device/cuda/src/seeding2/seed_finding2.cu
@@ -0,0 +1,96 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include <limits>
+#include <memory>
+#include <traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp>
+#include <traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp>
+#include <traccc/cuda/seeding2/kernels/write_output_kernel.hpp>
+#include <traccc/cuda/seeding2/seed_finding.hpp>
+#include <traccc/cuda/seeding2/types/kd_tree.hpp>
+#include <traccc/cuda/utils/definitions.hpp>
+#include <traccc/edm/seed.hpp>
+#include <traccc/edm/spacepoint.hpp>
+#include <traccc/seeding/detail/lin_circle.hpp>
+#include <traccc/seeding/detail/seeding_config.hpp>
+#include <vecmem/memory/memory_resource.hpp>
+#include <vecmem/memory/unique_ptr.hpp>
+#include <vecmem/utils/cuda/copy.hpp>
+#include <vector>
+
+namespace {
+/// Helper function that would produce a default seed-finder configuration
+__host__ traccc::seedfinder_config default_seedfinder2_config() {
+    traccc::seedfinder_config config;
+    traccc::seedfinder_config config_copy = config.toInternalUnits();
+    config.highland = 13.6f * std::sqrt(config_copy.radLengthPerSeed) *
+                      (1.f + 0.038f * std::log(config_copy.radLengthPerSeed));
+    float maxScatteringAngle = config.highland / config_copy.minPt;
+    config.maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle;
+    // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV
+    // and millimeter
+    config.pTPerHelixRadius = 300.f * config_copy.bFieldInZ;
+    config.minHelixDiameter2 =
+        std::pow(config_copy.minPt * 2.f / config.pTPerHelixRadius, 2.f);
+    config.pT2perRadius =
+        std::pow(config.highland / config.pTPerHelixRadius, 2.f);
+    return config;
+}
+}  // namespace
+
+namespace traccc::cuda {
+__global__ void construct_internal_sp_kernel(internal_sp_t out,
+                                             const spacepoint* in,
+                                             std::size_t n_sp) {
+    int tid = blockIdx.x * blockDim.x + threadIdx.x;
+
+    if (tid < n_sp) {
+        out[tid].x = in[tid].x();
+        out[tid].y = in[tid].y();
+        out[tid].z = in[tid].z();
+        out[tid].phi = atan2f(in[tid].y(), in[tid].x());
+        out[tid].radius = in[tid].radius();
+        out[tid].link = tid;
+    }
+}
+
+seed_finding2::seed_finding2(const traccc::memory_resource& mr)
+    : m_output_mr(mr), m_finder_conf(default_seedfinder2_config()) {}
+
+seed_finding2::output_type seed_finding2::operator()(
+    const spacepoint_collection_types::const_view& sps) const {
+    vecmem::cuda::copy copy;
+    vecmem::cuda::device_memory_resource mr;
+    std::size_t n_sp = copy.get_size(sps);
+
+    internal_sp_owning_t internal_sp_device(mr, n_sp);
+
+    static constexpr std::size_t copy_kernel_threads_per_block = 256;
+    construct_internal_sp_kernel<<<
+        (n_sp / copy_kernel_threads_per_block) +
+            (n_sp % copy_kernel_threads_per_block == 0 ? 0 : 1),
+        copy_kernel_threads_per_block>>>(internal_sp_t(internal_sp_device),
+                                         sps.ptr(), n_sp);
+
+    CUDA_ERROR_CHECK(cudaGetLastError());
+    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+
+    auto [kd_tree_device, kd_tree_size, internal_sp_device_new] =
+        create_kd_tree(mr, std::move(internal_sp_device), n_sp);
+
+    vecmem::unique_alloc_ptr<alt_seed[]> seeds_device;
+    uint32_t seed_count;
+
+    std::tie(seeds_device, seed_count) = run_seeding(
+        m_finder_conf, m_filter_conf, mr, internal_sp_t(internal_sp_device_new),
+        kd_tree_t(kd_tree_device));
+
+    return write_output(m_output_mr, seed_count,
+                        internal_sp_t(internal_sp_device_new),
+                        seeds_device.get());
+}
+}  // namespace traccc::cuda

From f3c0e1370ab6f85ee7dc490d553449fa5826661a Mon Sep 17 00:00:00 2001
From: Stephen Nicholas Swatman <stephen.nicholas.swatman@cern.ch>
Date: Tue, 25 Jun 2024 16:37:03 +0200
Subject: [PATCH 2/2] UpdatE

---
 core/include/traccc/edm/spacepoint.hpp        |   2 +
 device/cuda/CMakeLists.txt                    |   1 -
 .../cuda/seeding2/kernels/kd_tree_kernel.hpp  |   7 +-
 .../seeding2/kernels/seed_finding_kernel.hpp  |  10 +-
 .../seeding2/kernels/write_output_kernel.hpp  |  26 ---
 .../traccc/cuda/seeding2/seed_finding.hpp     |  23 ++-
 .../cuda/seeding2/types/internal_sp.hpp       |  35 ----
 .../traccc/cuda/utils/device_traits.hpp       |  56 ------
 .../src/seeding2/kernels/kd_tree_kernel.cu    | 164 ++++++++----------
 .../seeding2/kernels/seed_finding_kernel.cu   | 139 +++++++--------
 .../seeding2/kernels/write_output_kernel.cu   |  82 ---------
 device/cuda/src/seeding2/seed_finding2.cu     |  79 ++-------
 12 files changed, 182 insertions(+), 442 deletions(-)
 delete mode 100644 device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
 delete mode 100644 device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
 delete mode 100644 device/cuda/include/traccc/cuda/utils/device_traits.hpp
 delete mode 100644 device/cuda/src/seeding2/kernels/write_output_kernel.cu

diff --git a/core/include/traccc/edm/spacepoint.hpp b/core/include/traccc/edm/spacepoint.hpp
index 4814629248..add6bc1f08 100644
--- a/core/include/traccc/edm/spacepoint.hpp
+++ b/core/include/traccc/edm/spacepoint.hpp
@@ -42,6 +42,8 @@ struct spacepoint {
     scalar radius() const {
         return std::sqrt(global[0] * global[0] + global[1] * global[1]);
     }
+    TRACCC_HOST_DEVICE
+    scalar phi() const { return std::atan2(y(), x()); }
 };
 
 /// Comparison / ordering operator for spacepoints
diff --git a/device/cuda/CMakeLists.txt b/device/cuda/CMakeLists.txt
index 0e7067222f..784495698a 100644
--- a/device/cuda/CMakeLists.txt
+++ b/device/cuda/CMakeLists.txt
@@ -43,7 +43,6 @@ traccc_add_library( traccc_cuda cuda TYPE SHARED
   "src/seeding2/seed_finding2.cu"
   "src/seeding2/kernels/seed_finding_kernel.cu"
   "src/seeding2/kernels/kd_tree_kernel.cu"
-  "src/seeding2/kernels/write_output_kernel.cu"
   # Clusterization
   "include/traccc/cuda/clusterization/clusterization_algorithm.hpp"
   "src/clusterization/clusterization_algorithm.cu"
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
index 4ce1325af7..63afb2f3ac 100644
--- a/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
+++ b/device/cuda/include/traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp
@@ -7,10 +7,10 @@
 
 #pragma once
 
-#include <traccc/cuda/seeding2/types/internal_sp.hpp>
 #include <traccc/cuda/seeding2/types/kd_tree.hpp>
 #include <traccc/edm/internal_spacepoint.hpp>
 #include <traccc/edm/spacepoint.hpp>
+#include <vecmem/utils/copy.hpp>
 #include <vector>
 
 namespace traccc::cuda {
@@ -19,6 +19,7 @@ namespace traccc::cuda {
  *
  * @return A pair containing the k-d tree nodes as well as the number of nodes.
  */
-std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
-    vecmem::memory_resource&, internal_sp_owning_t&&, uint32_t);
+std::tuple<kd_tree_owning_t, uint32_t, vecmem::data::vector_buffer<std::size_t>>
+create_kd_tree(vecmem::memory_resource&, vecmem::copy& copy,
+               const spacepoint_collection_types::const_view&);
 }  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
index 464b4dd2af..a23cc37d56 100644
--- a/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
+++ b/device/cuda/include/traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp
@@ -7,12 +7,12 @@
 
 #pragma once
 
-#include <traccc/cuda/seeding2/types/internal_sp.hpp>
 #include <traccc/cuda/seeding2/types/kd_tree.hpp>
-#include <traccc/edm/alt_seed.hpp>
 #include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/seed.hpp>
 #include <traccc/edm/spacepoint.hpp>
 #include <traccc/seeding/detail/seeding_config.hpp>
+#include <vecmem/utils/copy.hpp>
 
 namespace traccc::cuda {
 /**
@@ -21,7 +21,7 @@ namespace traccc::cuda {
  * @return A pair containing the list of internal seeds as well as the number
  * of seeds.
  */
-std::pair<vecmem::unique_alloc_ptr<alt_seed[]>, uint32_t> run_seeding(
-    seedfinder_config, seedfilter_config, vecmem::memory_resource&,
-    internal_sp_t, kd_tree_t);
+seed_collection_types::buffer run_seeding(
+    seedfinder_config, seedfilter_config, vecmem::memory_resource &,
+    vecmem::copy &, const spacepoint_collection_types::const_view &, kd_tree_t);
 }  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp b/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
deleted file mode 100644
index 0bddf482d0..0000000000
--- a/device/cuda/include/traccc/cuda/seeding2/kernels/write_output_kernel.hpp
+++ /dev/null
@@ -1,26 +0,0 @@
-/** TRACCC library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-#include <traccc/cuda/seeding2/types/internal_sp.hpp>
-#include <traccc/edm/alt_seed.hpp>
-#include <traccc/edm/internal_spacepoint.hpp>
-#include <traccc/edm/seed.hpp>
-#include <traccc/edm/spacepoint.hpp>
-#include <traccc/utils/memory_resource.hpp>
-
-namespace traccc::cuda {
-/**
- * @brief Kernel to write output data back into traccc's EDM.
- *
- * @return A vector buffer containing the output seeds.
- */
-alt_seed_collection_types::buffer write_output(const traccc::memory_resource &,
-                                               uint32_t, const internal_sp_t,
-                                               const alt_seed *const);
-}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp b/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
index a350a73730..119a18188d 100644
--- a/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
+++ b/device/cuda/include/traccc/cuda/seeding2/seed_finding.hpp
@@ -7,29 +7,40 @@
 
 #pragma once
 
-#include <traccc/edm/alt_seed.hpp>
+#include <traccc/edm/seed.hpp>
 #include <traccc/edm/spacepoint.hpp>
 #include <traccc/seeding/detail/seeding_config.hpp>
 #include <traccc/utils/algorithm.hpp>
 #include <traccc/utils/memory_resource.hpp>
 #include <vecmem/memory/cuda/device_memory_resource.hpp>
+#include <vecmem/utils/copy.hpp>
+
+#include "traccc/cuda/utils/stream.hpp"
 
 namespace traccc::cuda {
 /**
  * @brief Alternative seed finding algorithm, using orthogonal range search
  * implemented through a k-d tree.
  */
-class seed_finding2 : public algorithm<alt_seed_collection_types::buffer(
+class seed_finding2 : public algorithm<seed_collection_types::buffer(
                           const spacepoint_collection_types::const_view&)> {
     public:
-    seed_finding2(const traccc::memory_resource& mr);
+    seed_finding2(const seedfinder_config& config,
+                  const seedfilter_config& filter_config,
+                  const traccc::memory_resource& mr, vecmem::copy& copy,
+                  stream& str);
 
     output_type operator()(
         const spacepoint_collection_types::const_view& sps) const override;
 
     private:
-    traccc::memory_resource m_output_mr;
-    seedfinder_config m_finder_conf;
-    seedfilter_config m_filter_conf;
+    traccc::memory_resource m_mr;
+    vecmem::copy& m_copy;
+    stream& m_stream;
+
+    seedfinder_config m_seedfinder_config;
+    seedfilter_config m_seedfilter_config;
+
+    int m_warp_size;
 };
 }  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp b/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
deleted file mode 100644
index 94b375ff4d..0000000000
--- a/device/cuda/include/traccc/cuda/seeding2/types/internal_sp.hpp
+++ /dev/null
@@ -1,35 +0,0 @@
-/** TRACCC library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-#include <array>
-#include <cstdint>
-#include <traccc/edm/internal_spacepoint.hpp>
-#include <traccc/edm/spacepoint.hpp>
-#include <traccc/utils/array_wrapper.hpp>
-#include <vecmem/memory/unique_ptr.hpp>
-
-namespace traccc::cuda {
-template <template <typename> typename F>
-struct sp_t {
-    using tuple_t = std::tuple<float, float, float, float, float, unsigned int>;
-
-    F<float> x;
-    F<float> y;
-    F<float> z;
-    F<float> phi;
-    F<float> radius;
-    F<unsigned int> link;
-};
-
-using internal_sp = array_wrapper<soa, sp_t>;
-
-using internal_sp_owning_t = internal_sp::owner;
-
-using internal_sp_t = internal_sp::handle;
-}  // namespace traccc::cuda
diff --git a/device/cuda/include/traccc/cuda/utils/device_traits.hpp b/device/cuda/include/traccc/cuda/utils/device_traits.hpp
deleted file mode 100644
index c1cc59eb1c..0000000000
--- a/device/cuda/include/traccc/cuda/utils/device_traits.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-/** TRACCC library, part of the ACTS project (R&D line)
- *
- * (c) 2021 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-#if defined(__CUDACC__) && defined(__CUDA_ARCH__)
-// Set values of maximum resident threads per SM
-#if __CUDA_ARCH__ < 500
-#pragma message \
-    "Very old CUDA architecture, setting maximum resident threads per SM to 1024."
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
-#elif __CUDA_ARCH__ <= 720
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
-#elif __CUDA_ARCH__ <= 750
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
-#elif __CUDA_ARCH__ <= 800
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
-#elif __CUDA_ARCH__ <= 890
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1536u
-#elif __CUDA_ARCH__ <= 900
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 2048u
-#pragma message \
-    "Unknown CUDA architecture, setting maximum resident threads per SM to 1024."
-#define CUDA_MAX_RESIDENT_THREADS_PER_SM 1024u
-#endif
-
-// Set values of maximum resident blocks per SM
-#if __CUDA_ARCH__ < 500
-#pragma message \
-    "Very old CUDA architecture, setting maximum resident blocks per SM to 16."
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
-#elif __CUDA_ARCH__ <= 720
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
-#elif __CUDA_ARCH__ <= 750
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
-#elif __CUDA_ARCH__ <= 800
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
-#elif __CUDA_ARCH__ <= 870
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
-#elif __CUDA_ARCH__ <= 890
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 24u
-#elif __CUDA_ARCH__ <= 900
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 32u
-#else
-#pragma message \
-    "Very old CUDA architecture, setting maximum resident blocks per SM to 16."
-#define CUDA_MAX_RESIDENT_BLOCKS_PER_SM 16u
-#endif
-
-// Set values of maximum threadss per block
-#define CUDA_MAX_THREADS_PER_BLOCK 1024u
-#endif
diff --git a/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu b/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
index ff5d54c301..e41f07cbbc 100644
--- a/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
+++ b/device/cuda/src/seeding2/kernels/kd_tree_kernel.cu
@@ -14,8 +14,6 @@
 #include <traccc/cuda/seeding2/seed_finding.hpp>
 #include <traccc/cuda/seeding2/types/kd_tree.hpp>
 #include <traccc/cuda/seeding2/types/range3d.hpp>
-#include <traccc/cuda/utils/definitions.hpp>
-#include <traccc/cuda/utils/device_traits.hpp>
 #include <traccc/cuda/utils/sort.hpp>
 #include <traccc/edm/seed.hpp>
 #include <traccc/edm/spacepoint.hpp>
@@ -25,6 +23,8 @@
 #include <vecmem/memory/unique_ptr.hpp>
 #include <vector>
 
+#include "../../utils/cuda_error_handling.hpp"
+
 namespace traccc::cuda {
 /**
  * @brief Bit mask for halo point indices.
@@ -62,23 +62,24 @@ constexpr uint32_t AUXILIARY_INDEX_MASK = (1u << 31);
  *
  * @return The ϕ value of the spacepoint with the given index.
  */
-__device__ float get_coordinate(const internal_sp_t spacepoints, uint32_t i,
-                                pivot_e dim) {
+__device__ float get_coordinate(
+    const spacepoint_collection_types::const_device& spacepoints, uint32_t i,
+    pivot_e dim) {
     bool is_real = !(i & AUXILIARY_INDEX_MASK);
 
     if (is_real) {
         if (dim == pivot_e::Phi) {
-            return spacepoints[i].phi;
+            return spacepoints[i].phi();
         } else if (dim == pivot_e::R) {
-            return spacepoints[i].radius;
+            return spacepoints[i].radius();
         } else {
-            return spacepoints[i].z;
+            return spacepoints[i].z();
         }
     } else {
         uint32_t real_i = i & ~AUXILIARY_INDEX_MASK;
 
         if (dim == pivot_e::Phi) {
-            float phi = spacepoints[real_i].phi;
+            float phi = spacepoints[real_i].phi();
 
             if (phi >= 0) {
                 return phi - 2.f * static_cast<float>(M_PI);
@@ -86,9 +87,9 @@ __device__ float get_coordinate(const internal_sp_t spacepoints, uint32_t i,
                 return phi + 2.f * static_cast<float>(M_PI);
             }
         } else if (dim == pivot_e::R) {
-            return spacepoints[real_i].radius;
+            return spacepoints[real_i].radius();
         } else {
-            return spacepoints[real_i].z;
+            return spacepoints[real_i].z();
         }
     }
 }
@@ -105,10 +106,11 @@ __device__ float get_coordinate(const internal_sp_t spacepoints, uint32_t i,
  * @param[in] sps The total number of spacepoints.
  * @param[out] extra_indices The output halo point count.
  */
-__global__ void initialize_index_vector(internal_sp_t spacepoints,
-                                        uint32_t* __restrict__ indices,
-                                        uint32_t sps,
-                                        uint32_t* __restrict__ extra_indices) {
+__global__ void initialize_index_vector(
+    const spacepoint_collection_types::const_view& _spacepoints,
+    uint32_t* __restrict__ indices, uint32_t sps,
+    uint32_t* __restrict__ extra_indices) {
+    const spacepoint_collection_types::const_device spacepoints(_spacepoints);
     cooperative_groups::grid_group grid = cooperative_groups::this_grid();
 
     if (grid.thread_rank() < sps) {
@@ -196,9 +198,9 @@ __device__ pivot_e next_pivot(pivot_e p) {
  * @return The approximate median value of the given indices in the given
  * pivot axis.
  */
-__device__ float get_approximate_median(internal_sp_t spacepoints,
-                                        const uint32_t* __restrict__ indices,
-                                        uint32_t num_indices, pivot_e pivot) {
+__device__ float get_approximate_median(
+    const spacepoint_collection_types::const_device& spacepoints,
+    const uint32_t* __restrict__ indices, uint32_t num_indices, pivot_e pivot) {
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
 
@@ -255,10 +257,10 @@ __device__ float get_approximate_median(internal_sp_t spacepoints,
  * @param[in] num_indices The number of indices for this function call.
  * @param[in] pivot The requested pivot axis.
  */
-__device__ uint32_t sort_and_get_median(internal_sp_t spacepoints,
-                                        const uint32_t* __restrict__ indices,
-                                        uint32_t* __restrict__ index_buffer,
-                                        uint32_t num_indices, pivot_e pivot) {
+__device__ uint32_t sort_and_get_median(
+    const spacepoint_collection_types::const_device& spacepoints,
+    const uint32_t* __restrict__ indices, uint32_t* __restrict__ index_buffer,
+    uint32_t num_indices, pivot_e pivot) {
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
 
@@ -299,10 +301,11 @@ __device__ uint32_t sort_and_get_median(internal_sp_t spacepoints,
  * @param[in] nid The node index.
  * @param[in] pivot The requested pivot axis.
  */
-__device__ bool create_internal_node(kd_tree_t tree, internal_sp_t spacepoints,
-                                     const uint32_t* __restrict__ indices,
-                                     uint32_t* __restrict__ index_buffer,
-                                     int nid, pivot_e pivot) {
+__device__ bool create_internal_node(
+    kd_tree_t tree,
+    const spacepoint_collection_types::const_device& spacepoints,
+    const uint32_t* __restrict__ indices, uint32_t* __restrict__ index_buffer,
+    int nid, pivot_e pivot) {
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
 
@@ -468,16 +471,16 @@ __device__ bool create_internal_node(kd_tree_t tree, internal_sp_t spacepoints,
  * @param[inout] _index_buffer Temporary buffer space for indices.
  * @param[in] iteration The iteration count.
  */
-__global__ void __launch_bounds__(512)
-    construct_kd_tree_big_step(kd_tree_t tree, uint32_t num_nodes,
-                               internal_sp_t spacepoints,
-                               const uint32_t* __restrict__ _indices_old,
-                               uint32_t* __restrict__ _indices_new,
-                               uint32_t iteration,
-                               uint32_t* __restrict__ work_list_current,
-                               uint32_t* __restrict__ work_list_size_current,
-                               uint32_t* __restrict__ work_list_next,
-                               uint32_t* __restrict__ work_list_size_next) {
+__global__ void construct_kd_tree_big_step(
+    kd_tree_t tree, uint32_t num_nodes,
+    const spacepoint_collection_types::const_view& _spacepoints,
+    const uint32_t* __restrict__ _indices_old,
+    uint32_t* __restrict__ _indices_new, uint32_t iteration,
+    uint32_t* __restrict__ work_list_current,
+    uint32_t* __restrict__ work_list_size_current,
+    uint32_t* __restrict__ work_list_next,
+    uint32_t* __restrict__ work_list_size_next) {
+    const spacepoint_collection_types::const_device spacepoints(_spacepoints);
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
 
@@ -556,16 +559,16 @@ __global__ void __launch_bounds__(512)
     }
 }
 
-__global__ void __launch_bounds__(96)
-    construct_kd_tree_small_step(kd_tree_t tree, uint32_t num_nodes,
-                                 internal_sp_t spacepoints,
-                                 const uint32_t* __restrict__ _indices_old,
-                                 uint32_t* __restrict__ _indices_new,
-                                 uint32_t iteration,
-                                 uint32_t* __restrict__ work_list_current,
-                                 uint32_t* __restrict__ work_list_size_current,
-                                 uint32_t* __restrict__ work_list_next,
-                                 uint32_t* __restrict__ work_list_size_next) {
+__global__ void construct_kd_tree_small_step(
+    kd_tree_t tree, uint32_t num_nodes,
+    const spacepoint_collection_types::const_view& _spacepoints,
+    const uint32_t* __restrict__ _indices_old,
+    uint32_t* __restrict__ _indices_new, uint32_t iteration,
+    uint32_t* __restrict__ work_list_current,
+    uint32_t* __restrict__ work_list_size_current,
+    uint32_t* __restrict__ work_list_next,
+    uint32_t* __restrict__ work_list_size_next) {
+    const spacepoint_collection_types::const_device spacepoints(_spacepoints);
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
 
@@ -654,28 +657,11 @@ uint32_t round_to_power_2(uint32_t i) {
     return power;
 }
 
-__global__ void bake_spacepoints_kernel(uint32_t n, internal_sp_t oldsps,
-                                        internal_sp_t newsps,
-                                        uint32_t* __restrict__ indices) {
-    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
-
-    uint32_t i = grid.thread_rank();
+std::tuple<kd_tree_owning_t, uint32_t, vecmem::data::vector_buffer<std::size_t>>
+create_kd_tree(vecmem::memory_resource& mr, vecmem::copy& copy,
+               const spacepoint_collection_types::const_view& spacepoints) {
+    std::size_t n_sp = copy.get_size(spacepoints);
 
-    if (i < n) {
-        uint32_t ni = indices[i] & ~AUXILIARY_INDEX_MASK;
-
-        newsps[i].x = oldsps[ni].x;
-        newsps[i].y = oldsps[ni].y;
-        newsps[i].z = oldsps[ni].z;
-        newsps[i].phi = oldsps[ni].phi;
-        newsps[i].radius = oldsps[ni].radius;
-        newsps[i].link = oldsps[ni].link;
-    }
-}
-
-std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
-    vecmem::memory_resource& mr, internal_sp_owning_t&& spacepoints,
-    uint32_t n_sp) {
     /*
      * Allocate space for the indices. Since each spacepoint can produce at
      * most two indices, we allocate enough space to support this theoretical
@@ -691,7 +677,8 @@ std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
     vecmem::unique_alloc_ptr<uint32_t> extra_indices =
         vecmem::make_unique_alloc<uint32_t>(mr);
 
-    CUDA_ERROR_CHECK(cudaMemset(extra_indices.get(), 0, sizeof(uint32_t)));
+    TRACCC_CUDA_ERROR_CHECK(
+        cudaMemset(extra_indices.get(), 0, sizeof(uint32_t)));
 
     /*
      * Launch the index initialization kernel, which turns n spacepoints into
@@ -703,16 +690,17 @@ std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
                               threads_per_block1>>>(spacepoints, indices.get(),
                                                     n_sp, extra_indices.get());
 
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+    TRACCC_CUDA_ERROR_CHECK(cudaGetLastError());
+    TRACCC_CUDA_ERROR_CHECK(cudaDeviceSynchronize());
 
     /*
      * Retrieve the number of halo points from the device.
      */
     uint32_t extra_indices_h;
 
-    CUDA_ERROR_CHECK(cudaMemcpy(&extra_indices_h, extra_indices.get(),
-                                sizeof(uint32_t), cudaMemcpyDeviceToHost));
+    TRACCC_CUDA_ERROR_CHECK(cudaMemcpy(&extra_indices_h, extra_indices.get(),
+                                       sizeof(uint32_t),
+                                       cudaMemcpyDeviceToHost));
 
     /*
      * The total number of indices is equal to the number of spacepoints plus
@@ -743,8 +731,10 @@ std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
     vecmem::unique_alloc_ptr<uint32_t> work_list_size_2 =
         vecmem::make_unique_alloc<uint32_t>(mr);
 
-    CUDA_ERROR_CHECK(cudaMemset(work_list_size_1.get(), 0, sizeof(uint32_t)));
-    CUDA_ERROR_CHECK(cudaMemset(work_list_size_2.get(), 0, sizeof(uint32_t)));
+    TRACCC_CUDA_ERROR_CHECK(
+        cudaMemset(work_list_size_1.get(), 0, sizeof(uint32_t)));
+    TRACCC_CUDA_ERROR_CHECK(
+        cudaMemset(work_list_size_2.get(), 0, sizeof(uint32_t)));
 
     /*
      * Run the initialization kernel for the tree itself, which sets the
@@ -757,8 +747,8 @@ std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
                                                num_indices, work_list_1.get(),
                                                work_list_size_1.get());
 
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+    TRACCC_CUDA_ERROR_CHECK(cudaGetLastError());
+    TRACCC_CUDA_ERROR_CHECK(cudaDeviceSynchronize());
 
     /*
      * Allocate some temporary space for the index buffer, one half for "lower"
@@ -801,33 +791,25 @@ std::tuple<kd_tree_owning_t, uint32_t, internal_sp_owning_t> create_kd_tree(
                 work_list_size_2.get());
         }
 
-        CUDA_ERROR_CHECK(cudaGetLastError());
-        CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+        TRACCC_CUDA_ERROR_CHECK(cudaGetLastError());
+        TRACCC_CUDA_ERROR_CHECK(cudaDeviceSynchronize());
 
         iteration++;
 
-        CUDA_ERROR_CHECK(cudaMemcpy(&remaining, work_list_size_2.get(),
-                                    sizeof(uint32_t), cudaMemcpyDeviceToHost));
+        TRACCC_CUDA_ERROR_CHECK(cudaMemcpy(&remaining, work_list_size_2.get(),
+                                           sizeof(uint32_t),
+                                           cudaMemcpyDeviceToHost));
 
         std::swap(work_list_1, work_list_2);
         std::swap(work_list_size_1, work_list_size_2);
         std::swap(indices, index_buffer);
 
-        CUDA_ERROR_CHECK(
+        TRACCC_CUDA_ERROR_CHECK(
             cudaMemset(work_list_size_2.get(), 0, sizeof(uint32_t)));
     }
 
-    internal_sp_owning_t new_sps(mr, num_indices);
-
-    bake_spacepoints_kernel<<<
-        num_indices / 1024u + (num_indices % 1024u == 0u ? 0u : 1u), 1024u>>>(
-        num_indices, internal_sp_t(spacepoints), internal_sp_t(new_sps),
-        indices.get());
-
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+    vecmem::data::vector_buffer<std::size_t> out;
 
-    return {kd_tree_owning_t{std::move(tree_owner)}, num_nodes,
-            std::move(new_sps)};
+    return {kd_tree_owning_t{std::move(tree_owner)}, num_nodes, std::move(out)};
 }
 }  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu b/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
index 9e1204e7fd..601cef0b41 100644
--- a/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
+++ b/device/cuda/src/seeding2/kernels/seed_finding_kernel.cu
@@ -9,24 +9,25 @@
 
 #include <iostream>
 #include <traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp>
-#include <traccc/cuda/seeding2/types/internal_sp.hpp>
 #include <traccc/cuda/seeding2/types/kd_tree.hpp>
 #include <traccc/cuda/seeding2/types/range3d.hpp>
-#include <traccc/cuda/utils/definitions.hpp>
-#include <traccc/cuda/utils/device_traits.hpp>
 #include <traccc/cuda/utils/sort.hpp>
 #include <traccc/cuda/utils/sync.hpp>
-#include <traccc/edm/alt_seed.hpp>
 #include <traccc/edm/internal_spacepoint.hpp>
+#include <traccc/edm/seed.hpp>
 #include <traccc/edm/spacepoint.hpp>
 #include <traccc/seeding/detail/lin_circle.hpp>
+#include <traccc/seeding/detail/spacepoint_type.hpp>
 #include <traccc/seeding/doublet_finding_helper.hpp>
 #include <traccc/seeding/triplet_finding_helper.hpp>
 
+#include "../../utils/cuda_error_handling.hpp"
+
 #define MAX_LOWER_SP_PER_MIDDLE 100u
 #define MAX_UPPER_SP_PER_MIDDLE 100u
 #define KD_TREE_TRAVERSAL_STACK_SIZE 128u
 #define WARPS_PER_BLOCK 8u
+#define WARP_SIZE 32u
 
 /**
  * @brief Maximum difference in φ between two spacepoints adjacent in the same
@@ -65,10 +66,10 @@ namespace traccc::cuda {
  */
 __device__ void retrieve_from_tree(
     const seedfinder_config finder_conf, const seedfilter_config filter_conf,
-    const internal_sp_t spacepoints, const kd_tree_t tree,
-    const internal_spacepoint<spacepoint> mi, const range3d range,
-    uint32_t mi_idx, bool upper, uint32_t* __restrict__ output_arr,
-    uint32_t* __restrict__ output_cnt) {
+    const spacepoint_collection_types::const_device& spacepoints,
+    const kd_tree_t tree, const internal_spacepoint<spacepoint> mi,
+    const range3d range, uint32_t mi_idx, bool upper,
+    uint32_t* __restrict__ output_arr, uint32_t* __restrict__ output_cnt) {
     __shared__ uint32_t stack[WARPS_PER_BLOCK][KD_TREE_TRAVERSAL_STACK_SIZE];
     __shared__ uint32_t leaves[WARPS_PER_BLOCK][WARP_SIZE];
 
@@ -137,16 +138,15 @@ __device__ void retrieve_from_tree(
                  * are more specific than those encoded in the search
                  * range.
                  */
-                float phi = spacepoints[point_id].phi,
-                      radius = spacepoints[point_id].radius,
-                      z = spacepoints[point_id].z;
+                float phi = spacepoints[point_id].phi(),
+                      radius = spacepoints[point_id].radius(),
+                      z = spacepoints[point_id].z();
 
                 if (range.contains(phi, radius, z) && point_id != mi_idx &&
                     doublet_finding_helper::isCompatible(
                         upper, mi,
-                        internal_spacepoint<spacepoint>(
-                            spacepoints[point_id].x, spacepoints[point_id].y, z,
-                            radius, phi, spacepoints[point_id].link),
+                        internal_spacepoint<spacepoint>(spacepoints[point_id],
+                                                        point_id, {0., 0.}),
                         finder_conf)) {
                     /*
                      * Reserve a spot in the output for this point, and
@@ -195,24 +195,21 @@ __device__ void retrieve_from_tree(
  */
 __device__ internal_seed make_seed(
     const seedfinder_config finder_conf, const seedfilter_config filter_conf,
-    const internal_sp_t spacepoints, const kd_tree_t tree, uint32_t lower,
-    uint32_t middle, internal_spacepoint<spacepoint> mi, uint32_t upper) {
-    const internal_spacepoint<spacepoint> lo(
-        spacepoints[lower].x, spacepoints[lower].y, spacepoints[lower].z,
-        spacepoints[lower].radius, spacepoints[lower].phi,
-        spacepoints[lower].link);
-    const internal_spacepoint<spacepoint> hi(
-        spacepoints[upper].x, spacepoints[upper].y, spacepoints[upper].z,
-        spacepoints[upper].radius, spacepoints[upper].phi,
-        spacepoints[upper].link);
+    const spacepoint_collection_types::const_device& spacepoints,
+    const kd_tree_t tree, uint32_t lower, uint32_t middle,
+    internal_spacepoint<spacepoint> mi, uint32_t upper) {
+    const internal_spacepoint<spacepoint> lo(spacepoints[lower], lower,
+                                             {0., 0.});
+    const internal_spacepoint<spacepoint> hi(spacepoints[upper], upper,
+                                             {0., 0.});
 
     /*
      * Find the lin-circles for the bottom and top pair.
      */
     lin_circle lm = doublet_finding_helper::transform_coordinates<
-        details::spacepoint_type::bottom>(mi, lo);
+        ::traccc::details::spacepoint_type::bottom>(mi, lo);
     lin_circle mh = doublet_finding_helper::transform_coordinates<
-        details::spacepoint_type::top>(mi, hi);
+        ::traccc::details::spacepoint_type::top>(mi, hi);
 
     scalar iSinTheta2 = 1 + lm.cotTheta() * lm.cotTheta();
     scalar scatteringInRegion2 = finder_conf.maxScatteringAngle2 * iSinTheta2;
@@ -352,13 +349,11 @@ __device__ range3d get_upper_range3d(const seedfinder_config finder_conf,
  * @param[inout] internal_seeds The array to write internal seeds to.
  * @param[in] increasing_z True iff we are looking for increasing-z seeds.
  */
-__device__ void run_helper(const seedfinder_config finder_conf,
-                           const seedfilter_config filter_conf,
-                           const internal_sp_t spacepoints,
-                           const kd_tree_t tree,
-                           const internal_spacepoint<spacepoint> sp,
-                           std::size_t sp_idx, internal_seed* internal_seeds,
-                           bool increasing_z) {
+__device__ void run_helper(
+    const seedfinder_config finder_conf, const seedfilter_config filter_conf,
+    const spacepoint_collection_types::const_device& spacepoints,
+    const kd_tree_t tree, const internal_spacepoint<spacepoint> sp,
+    std::size_t sp_idx, internal_seed* internal_seeds, bool increasing_z) {
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
     cooperative_groups::thread_block_tile<WARP_SIZE> warp =
@@ -528,11 +523,13 @@ __device__ void run_helper(const seedfinder_config finder_conf,
  * one-block-per-middle-SP fashion.
  *
  */
-__global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
-    seed_finding_kernel(const seedfinder_config finder_conf,
-                        const seedfilter_config filter_conf,
-                        const internal_sp_t spacepoints, const kd_tree_t tree,
-                        alt_seed* output_seeds, uint32_t* output_seed_size) {
+__global__ void seed_finding_kernel(
+    const seedfinder_config finder_conf, const seedfilter_config filter_conf,
+    const spacepoint_collection_types::const_view& _spacepoints,
+    const kd_tree_t tree, seed_collection_types::view _seeds) {
+    const spacepoint_collection_types::const_device spacepoints(_spacepoints);
+    seed_collection_types::device seeds(_seeds);
+
     cooperative_groups::thread_block block =
         cooperative_groups::this_thread_block();
     cooperative_groups::thread_block_tile<WARP_SIZE> warp =
@@ -548,19 +545,17 @@ __global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
      * ignore some of the middle spacepoint candidates.
      */
     if (spacepoint_id >= spacepoints.size() ||
-        spacepoints[spacepoint_id].phi < finder_conf.phiMin ||
-        spacepoints[spacepoint_id].phi > finder_conf.phiMax ||
-        spacepoints[spacepoint_id].radius < finder_conf.rMin ||
-        spacepoints[spacepoint_id].radius > finder_conf.rMax ||
-        spacepoints[spacepoint_id].z < finder_conf.zMin ||
-        spacepoints[spacepoint_id].z > finder_conf.zMax) {
+        spacepoints[spacepoint_id].phi() < finder_conf.phiMin ||
+        spacepoints[spacepoint_id].phi() > finder_conf.phiMax ||
+        spacepoints[spacepoint_id].radius() < finder_conf.rMin ||
+        spacepoints[spacepoint_id].radius() > finder_conf.rMax ||
+        spacepoints[spacepoint_id].z() < finder_conf.zMin ||
+        spacepoints[spacepoint_id].z() > finder_conf.zMax) {
         return;
     }
 
-    const internal_spacepoint<spacepoint> sp(
-        spacepoints[spacepoint_id].x, spacepoints[spacepoint_id].y,
-        spacepoints[spacepoint_id].z, spacepoints[spacepoint_id].radius,
-        spacepoints[spacepoint_id].phi, spacepoints[spacepoint_id].link);
+    const internal_spacepoint<spacepoint> sp(spacepoints[spacepoint_id],
+                                             spacepoint_id, {0., 0.});
 
     /*
      * The internal seeds are stored in this shared array.
@@ -614,7 +609,7 @@ __global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
          * Next, reserve space atomically in the output array for our new
          * seeds.
          */
-        uint32_t idx = atomicAdd(output_seed_size, num_valid);
+        uint32_t idx = seeds.bulk_append_implicit_unsafe(num_valid);
 
         /*
          * Finally, actually write the output to global memory. Again, this is
@@ -624,15 +619,15 @@ __global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
          * than invalid ones (if we have any).
          */
         for (uint32_t i = 0; i < num_valid; ++i) {
-            alt_seed& s = output_seeds[idx + i];
+            seed& s = seeds[idx + i];
             uint32_t j = WARP_SIZE + finder_conf.maxSeedsPerSpM - (i + 1);
 
             /*
              * Write the data back to the output.
              */
-            s.spB_link = spacepoints[internal_seeds[j].spacepoints[0]].link;
-            s.spM_link = spacepoints[spacepoint_id].link;
-            s.spT_link = spacepoints[internal_seeds[j].spacepoints[1]].link;
+            s.spB_link = internal_seeds[j].spacepoints[0];
+            s.spM_link = spacepoint_id;
+            s.spT_link = internal_seeds[j].spacepoints[1];
 
             s.weight = internal_seeds[i].weight;
             s.z_vertex = 0.f;
@@ -654,23 +649,29 @@ __global__ void __launch_bounds__(WARP_SIZE* WARPS_PER_BLOCK)
  * @return A unique pointer to an array of internal seeds, and the size of that
  * array.
  */
-std::pair<vecmem::unique_alloc_ptr<alt_seed[]>, uint32_t> run_seeding(
+seed_collection_types::buffer run_seeding(
     seedfinder_config finder_conf, seedfilter_config filter_conf,
-    vecmem::memory_resource& mr, const internal_sp_t spacepoints,
+    vecmem::memory_resource& mr, vecmem::copy& copy,
+    const spacepoint_collection_types::const_view& spacepoints,
     const kd_tree_t tree) {
+    std::size_t n_sp = copy.get_size(spacepoints);
+
     /*
      * Allocate space for output of seeds on the device.
      */
-    vecmem::unique_alloc_ptr<alt_seed[]> seeds_device =
-        vecmem::make_unique_alloc<alt_seed[]>(
-            mr, finder_conf.maxSeedsPerSpM * spacepoints.size());
+    std::size_t output_size = finder_conf.maxSeedsPerSpM * n_sp;
+
+    seed_collection_types::buffer seed_buffer(
+        output_size, mr, vecmem::data::buffer_type::resizable);
+    copy.setup(seed_buffer)->wait();
 
     /*
      * Allocate space for seed count on the device.
      */
     vecmem::unique_alloc_ptr<uint32_t> seed_count_device =
         vecmem::make_unique_alloc<uint32_t>(mr);
-    CUDA_ERROR_CHECK(cudaMemset(seed_count_device.get(), 0, sizeof(uint32_t)));
+    TRACCC_CUDA_ERROR_CHECK(
+        cudaMemset(seed_count_device.get(), 0, sizeof(uint32_t)));
 
     /*
      * Calculate the total amount of shared memory on top of that which is
@@ -680,27 +681,19 @@ std::pair<vecmem::unique_alloc_ptr<alt_seed[]>, uint32_t> run_seeding(
     std::size_t total_shared_memory =
         (WARPS_PER_BLOCK * (WARP_SIZE + finder_conf.maxSeedsPerSpM)) *
         sizeof(internal_seed);
-    std::size_t grid_size = (spacepoints.size() / WARPS_PER_BLOCK) +
-                            (spacepoints.size() % WARPS_PER_BLOCK == 0 ? 0 : 1);
+    std::size_t grid_size =
+        (n_sp / WARPS_PER_BLOCK) + (n_sp % WARPS_PER_BLOCK == 0 ? 0 : 1);
     std::size_t block_size = WARPS_PER_BLOCK * WARP_SIZE;
 
     /*
      * Call the kernel and all that jazz.
      */
     seed_finding_kernel<<<grid_size, block_size, total_shared_memory>>>(
-        finder_conf, filter_conf, spacepoints, tree, seeds_device.get(),
-        seed_count_device.get());
-
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
+        finder_conf, filter_conf, spacepoints, tree, seed_buffer);
 
-    /*
-     * Transfer the seed count back to the host and then hand it to the user.
-     */
-    uint32_t seed_count_host;
-    CUDA_ERROR_CHECK(cudaMemcpy(&seed_count_host, seed_count_device.get(),
-                                sizeof(uint32_t), cudaMemcpyDeviceToHost));
+    TRACCC_CUDA_ERROR_CHECK(cudaGetLastError());
+    TRACCC_CUDA_ERROR_CHECK(cudaDeviceSynchronize());
 
-    return {std::move(seeds_device), seed_count_host};
+    return seed_buffer;
 }
 }  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/kernels/write_output_kernel.cu b/device/cuda/src/seeding2/kernels/write_output_kernel.cu
deleted file mode 100644
index d9bb4a73f2..0000000000
--- a/device/cuda/src/seeding2/kernels/write_output_kernel.cu
+++ /dev/null
@@ -1,82 +0,0 @@
-/** TRACCC library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#include <cooperative_groups.h>
-
-#include <cstdint>
-#include <traccc/cuda/seeding2/kernels/write_output_kernel.hpp>
-#include <traccc/cuda/seeding2/types/internal_sp.hpp>
-#include <traccc/cuda/utils/definitions.hpp>
-#include <traccc/edm/alt_seed.hpp>
-#include <traccc/edm/internal_spacepoint.hpp>
-#include <traccc/edm/seed.hpp>
-#include <traccc/edm/spacepoint.hpp>
-#include <traccc/utils/memory_resource.hpp>
-#include <vecmem/containers/device_vector.hpp>
-#include <vecmem/containers/vector.hpp>
-#include <vecmem/utils/cuda/copy.hpp>
-
-namespace traccc::cuda {
-/**
- * @brief Main output writing kernel.
- *
- * @param[in] n_seeds The number of seeds to write.
- * @param[in] spacepoints Spacepoint array into which seeds are indexed.
- * @param[in] in_seeds Input internal seeds to read from.
- * @param[out] out_seeds Output seeds to write to.
- */
-__global__ void write_output_kernel(uint32_t n_seeds,
-                                    const internal_sp_t spacepoints,
-                                    const alt_seed* const in_seeds,
-                                    alt_seed_collection_types::view out_seeds) {
-    cooperative_groups::grid_group grid = cooperative_groups::this_grid();
-
-    /*
-     * Create a writable device vector for the seeds.
-     */
-    vecmem::device_vector<alt_seed> dvec(out_seeds);
-
-    int i = grid.thread_rank();
-
-    if (i < n_seeds) {
-        /*
-         * Get the output seed for the current thread.
-         */
-        dvec.at(i) = in_seeds[i];
-    }
-}
-
-alt_seed_collection_types::buffer write_output(
-    const traccc::memory_resource& mr, uint32_t n_seeds,
-    const internal_sp_t spacepoints, const alt_seed* const seeds) {
-    static constexpr std::size_t threads_per_block = 256;
-
-    /*
-     * Create a "vector buffer" that we can write our output to.
-     */
-    alt_seed_collection_types::buffer out(n_seeds, mr.main);
-
-    /*
-     * Fetch the current CUDA device number.
-     */
-    int device;
-    CUDA_ERROR_CHECK(cudaGetDevice(&device));
-
-    /*
-     * Execute the main output kernel.
-     */
-    write_output_kernel<<<(n_seeds / threads_per_block) +
-                              (n_seeds % threads_per_block == 0 ? 0 : 1),
-                          threads_per_block>>>(n_seeds, spacepoints, seeds,
-                                               out);
-
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
-
-    return out;
-}
-}  // namespace traccc::cuda
diff --git a/device/cuda/src/seeding2/seed_finding2.cu b/device/cuda/src/seeding2/seed_finding2.cu
index 0837873889..222660c2a6 100644
--- a/device/cuda/src/seeding2/seed_finding2.cu
+++ b/device/cuda/src/seeding2/seed_finding2.cu
@@ -9,10 +9,8 @@
 #include <memory>
 #include <traccc/cuda/seeding2/kernels/kd_tree_kernel.hpp>
 #include <traccc/cuda/seeding2/kernels/seed_finding_kernel.hpp>
-#include <traccc/cuda/seeding2/kernels/write_output_kernel.hpp>
 #include <traccc/cuda/seeding2/seed_finding.hpp>
 #include <traccc/cuda/seeding2/types/kd_tree.hpp>
-#include <traccc/cuda/utils/definitions.hpp>
 #include <traccc/edm/seed.hpp>
 #include <traccc/edm/spacepoint.hpp>
 #include <traccc/seeding/detail/lin_circle.hpp>
@@ -22,75 +20,28 @@
 #include <vecmem/utils/cuda/copy.hpp>
 #include <vector>
 
-namespace {
-/// Helper function that would produce a default seed-finder configuration
-__host__ traccc::seedfinder_config default_seedfinder2_config() {
-    traccc::seedfinder_config config;
-    traccc::seedfinder_config config_copy = config.toInternalUnits();
-    config.highland = 13.6f * std::sqrt(config_copy.radLengthPerSeed) *
-                      (1.f + 0.038f * std::log(config_copy.radLengthPerSeed));
-    float maxScatteringAngle = config.highland / config_copy.minPt;
-    config.maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle;
-    // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV
-    // and millimeter
-    config.pTPerHelixRadius = 300.f * config_copy.bFieldInZ;
-    config.minHelixDiameter2 =
-        std::pow(config_copy.minPt * 2.f / config.pTPerHelixRadius, 2.f);
-    config.pT2perRadius =
-        std::pow(config.highland / config.pTPerHelixRadius, 2.f);
-    return config;
-}
-}  // namespace
+#include "../utils/cuda_error_handling.hpp"
+#include "../utils/utils.hpp"
 
 namespace traccc::cuda {
-__global__ void construct_internal_sp_kernel(internal_sp_t out,
-                                             const spacepoint* in,
-                                             std::size_t n_sp) {
-    int tid = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (tid < n_sp) {
-        out[tid].x = in[tid].x();
-        out[tid].y = in[tid].y();
-        out[tid].z = in[tid].z();
-        out[tid].phi = atan2f(in[tid].y(), in[tid].x());
-        out[tid].radius = in[tid].radius();
-        out[tid].link = tid;
-    }
-}
-
-seed_finding2::seed_finding2(const traccc::memory_resource& mr)
-    : m_output_mr(mr), m_finder_conf(default_seedfinder2_config()) {}
+seed_finding2::seed_finding2(const seedfinder_config& config,
+                             const seedfilter_config& filter_config,
+                             const traccc::memory_resource& mr,
+                             vecmem::copy& copy, stream& str)
+    : m_seedfinder_config(config),
+      m_seedfilter_config(filter_config),
+      m_mr(mr),
+      m_copy(copy),
+      m_stream(str),
+      m_warp_size(details::get_warp_size(str.device())) {}
 
 seed_finding2::output_type seed_finding2::operator()(
     const spacepoint_collection_types::const_view& sps) const {
-    vecmem::cuda::copy copy;
-    vecmem::cuda::device_memory_resource mr;
-    std::size_t n_sp = copy.get_size(sps);
-
-    internal_sp_owning_t internal_sp_device(mr, n_sp);
-
-    static constexpr std::size_t copy_kernel_threads_per_block = 256;
-    construct_internal_sp_kernel<<<
-        (n_sp / copy_kernel_threads_per_block) +
-            (n_sp % copy_kernel_threads_per_block == 0 ? 0 : 1),
-        copy_kernel_threads_per_block>>>(internal_sp_t(internal_sp_device),
-                                         sps.ptr(), n_sp);
-
-    CUDA_ERROR_CHECK(cudaGetLastError());
-    CUDA_ERROR_CHECK(cudaDeviceSynchronize());
 
     auto [kd_tree_device, kd_tree_size, internal_sp_device_new] =
-        create_kd_tree(mr, std::move(internal_sp_device), n_sp);
-
-    vecmem::unique_alloc_ptr<alt_seed[]> seeds_device;
-    uint32_t seed_count;
-
-    std::tie(seeds_device, seed_count) = run_seeding(
-        m_finder_conf, m_filter_conf, mr, internal_sp_t(internal_sp_device_new),
-        kd_tree_t(kd_tree_device));
+        create_kd_tree(m_mr.main, m_copy, sps);
 
-    return write_output(m_output_mr, seed_count,
-                        internal_sp_t(internal_sp_device_new),
-                        seeds_device.get());
+    return run_seeding(m_seedfinder_config, m_seedfilter_config, m_mr.main,
+                       m_copy, sps, kd_tree_t(kd_tree_device));
 }
 }  // namespace traccc::cuda