Merge pull request #1090 from krasznaa/SeedingThroughput-main-20250724

Seeding Throughput, main branch (2025.07.24.)

Author: krasznaa

examples/options/include/traccc/options/throughput.hpp

@@ -23,6 +23,14 @@ class throughput : public interface {
     /// @name Options
     /// @{
 
+    /// "Reconstruction stage" to run
+    enum class stage {
+        seeding,  ///< Run until the end of seeding
+        full      ///< Run the full chain of reconstruction
+    };
+    /// The reconstruction stage to run
+    stage reco_stage = stage::full;
+
     /// The number of events to process during the job
     std::size_t processed_events = 100;
     /// The number of events to run "cold", i.e. run without accounting for
@@ -43,6 +51,13 @@ class throughput : public interface {
     throughput();
 
     std::unique_ptr<configuration_printable> as_printable() const override;
+
+    /// Read/process the command line options
+    ///
+    /// @param vm The command line options to interpret/read
+    ///
+    void read(const boost::program_options::variables_map& vm) override;
+
 };  // class throughput
 
 }  // namespace traccc::opts

examples/options/src/throughput.cpp

@@ -18,8 +18,16 @@ namespace traccc::opts {
 /// Convenience namespace shorthand
 namespace po = boost::program_options;
 
+/// Type alias for the reconstruction stage enumeration
+using stage_type = std::string;
+/// Name of the reconstruction stage option
+static const char* stage_option = "reco-stage";
+
 throughput::throughput() : interface("Throughput Measurement Options") {
 
+    m_desc.add_options()(
+        stage_option, po::value<stage_type>()->default_value("full"),
+        "Reconstruction stage to run (\"seeding\" or \"full\")");
     m_desc.add_options()(
         "processed-events",
         po::value(&processed_events)->default_value(processed_events),
@@ -40,9 +48,38 @@ throughput::throughput() : interface("Throughput Measurement Options") {
         "File where result logs will be printed (in append mode).");
 }
 
+void throughput::read(const po::variables_map& vm) {
+
+    // Decode the input data format.
+    if (vm.count(stage_option)) {
+        const std::string stage_string = vm[stage_option].as<stage_type>();
+        if (stage_string == "full") {
+            reco_stage = stage::full;
+        } else if (stage_string == "seeding") {
+            reco_stage = stage::seeding;
+        } else {
+            throw std::invalid_argument("Unknown reconstruction stage");
+        }
+    }
+}
+
 std::unique_ptr<configuration_printable> throughput::as_printable() const {
     auto cat = std::make_unique<configuration_category>(m_description);
 
+    std::string reco_stage_string;
+    switch (reco_stage) {
+        case stage::seeding:
+            reco_stage_string = "seeding";
+            break;
+        case stage::full:
+            reco_stage_string = "full";
+            break;
+        default:
+            reco_stage_string = "unknown";
+            break;
+    }
+    cat->add_child(std::make_unique<configuration_kv_pair>(
+        "Reconstruction stage", reco_stage_string));
     cat->add_child(std::make_unique<configuration_kv_pair>(
         "Cold run events", std::to_string(cold_run_events)));
     cat->add_child(std::make_unique<configuration_kv_pair>(

examples/run/alpaka/full_chain_algorithm.cpp

@@ -200,4 +200,47 @@ full_chain_algorithm::output_type full_chain_algorithm::operator()(
     }
 }
 
+bound_track_parameters_collection_types::host full_chain_algorithm::seeding(
+    const edm::silicon_cell_collection::host& cells) const {
+
+    // Create device copy of input collections
+    edm::silicon_cell_collection::buffer cells_buffer(
+        static_cast<unsigned int>(cells.size()), *m_cached_device_mr);
+    m_vecmem_objects.async_copy()(::vecmem::get_data(cells), cells_buffer)
+        ->ignore();
+
+    // Run the clusterization (asynchronously).
+    const clusterization_algorithm::output_type measurements =
+        m_clusterization(cells_buffer, m_device_det_descr);
+    m_measurement_sorting(measurements);
+
+    // If we have a Detray detector, run the track finding and fitting.
+    if (m_detector != nullptr) {
+
+        // Run the seed-finding (asynchronously).
+        const spacepoint_formation_algorithm::output_type spacepoints =
+            m_spacepoint_formation(m_device_detector_view, measurements);
+        const track_params_estimation::output_type track_params =
+            m_track_parameter_estimation(measurements, spacepoints,
+                                         m_seeding(spacepoints), m_field_vec);
+
+        // Copy a limited amount of result data back to the host.
+        bound_track_parameters_collection_types::host result{&m_host_mr};
+        m_vecmem_objects.async_copy()(track_params, result)->wait();
+        return result;
+
+    }
+    // If not, copy the track parameters back to the host, and return a dummy
+    // object.
+    else {
+
+        // Copy the measurements back to the host.
+        measurement_collection_types::host measurements_host(&m_host_mr);
+        m_vecmem_objects.async_copy()(measurements, measurements_host)->wait();
+
+        // Return an empty object.
+        return {};
+    }
+}
+
 }  // namespace traccc::alpaka

examples/run/alpaka/full_chain_algorithm.hpp

@@ -20,6 +20,7 @@
 #include "traccc/bfield/magnetic_field.hpp"
 #include "traccc/clusterization/clustering_config.hpp"
 #include "traccc/edm/silicon_cell_collection.hpp"
+#include "traccc/edm/track_parameters.hpp"
 #include "traccc/edm/track_state.hpp"
 #include "traccc/fitting/kalman_filter/kalman_fitter.hpp"
 #include "traccc/geometry/detector.hpp"
@@ -108,6 +109,14 @@ class full_chain_algorithm
     output_type operator()(
         const edm::silicon_cell_collection::host& cells) const override;
 
+    /// Reconstruct track seeds in the entire detector
+    ///
+    /// @param cells The cells for every detector module in the event
+    /// @return The track seeds reconstructed
+    ///
+    bound_track_parameters_collection_types::host seeding(
+        const edm::silicon_cell_collection::host& cells) const;
+
     private:
     /// Alpaka Queue
     traccc::alpaka::queue m_queue;

examples/run/common/throughput_mt.ipp

@@ -54,6 +54,7 @@
 #include <cstdlib>
 #include <ctime>
 #include <fstream>
+#include <functional>
 #include <iostream>
 #include <memory>
 #include <vector>
@@ -187,6 +188,27 @@ int throughput_mt(std::string_view description, int argc, char* argv[],
              logger().clone()});
     }
 
+    // Set up a lambda that calls the correct function on the algorithms.
+    std::function<std::size_t(int, const edm::silicon_cell_collection::host&)>
+        process_event;
+    if (throughput_opts.reco_stage == opts::throughput::stage::seeding) {
+        process_event = [&](int thread,
+                            const edm::silicon_cell_collection::host& cells)
+            -> std::size_t {
+            return algs.at(static_cast<std::size_t>(thread))
+                .seeding(cells)
+                .size();
+        };
+    } else if (throughput_opts.reco_stage == opts::throughput::stage::full) {
+        process_event = [&](int thread,
+                            const edm::silicon_cell_collection::host& cells)
+            -> std::size_t {
+            return algs.at(static_cast<std::size_t>(thread))(cells).size();
+        };
+    } else {
+        throw std::invalid_argument("Unknown reconstruction stage");
+    }
+
     // Set up the TBB arena and thread group. From here on out TBB is only
     // allowed to use the specified number of threads.
     tbb::global_control global_thread_limit(
@@ -233,10 +255,9 @@ int throughput_mt(std::string_view description, int argc, char* argv[],
             // Launch the processing of the event.
             arena.execute([&, event]() {
                 group.run([&, event]() {
-                    rec_track_params.fetch_add(algs.at(static_cast<std::size_t>(
-                        tbb::this_task_arena::current_thread_index()))(
-                                                       input[event])
-                                                   .size());
+                    rec_track_params.fetch_add(process_event(
+                        tbb::this_task_arena::current_thread_index(),
+                        input[event]));
                     progress_bar.tick();
                 });
             });

examples/run/common/throughput_st.ipp

@@ -45,6 +45,7 @@
 // System include(s).
 #include <cstdlib>
 #include <ctime>
+#include <functional>
 #include <iostream>
 #include <memory>
 
@@ -151,6 +152,19 @@ int throughput_st(std::string_view description, int argc, char* argv[],
         std::srand(throughput_opts.random_seed);
     }
 
+    // Set up a lambda that calls the correct function on the algorithm.
+    std::function<std::size_t(const edm::silicon_cell_collection::host&)>
+        process_event;
+    if (throughput_opts.reco_stage == opts::throughput::stage::seeding) {
+        process_event = [&](const edm::silicon_cell_collection::host& cells)
+            -> std::size_t { return alg->seeding(cells).size(); };
+    } else if (throughput_opts.reco_stage == opts::throughput::stage::full) {
+        process_event = [&](const edm::silicon_cell_collection::host& cells)
+            -> std::size_t { return (*alg)(cells).size(); };
+    } else {
+        throw std::invalid_argument("Unknown reconstruction stage");
+    }
+
     // Dummy count uses output of tp algorithm to ensure the compiler
     // optimisations don't skip any step
     std::size_t rec_track_params = 0;
@@ -180,7 +194,7 @@ int throughput_st(std::string_view description, int argc, char* argv[],
                 input_opts.events;
 
             // Process one event.
-            rec_track_params += (*alg)(input[event]).size();
+            rec_track_params += process_event(input[event]);
             progress_bar.tick();
         }
     }

examples/run/cpu/full_chain_algorithm.cpp

@@ -90,4 +90,41 @@ full_chain_algorithm::output_type full_chain_algorithm::operator()(
     }
 }
 
+bound_track_parameters_collection_types::host full_chain_algorithm::seeding(
+    const edm::silicon_cell_collection::host& cells) const {
+
+    // Create a data object for the detector description.
+    const silicon_detector_description::const_data det_descr_data =
+        vecmem::get_data(m_det_descr.get());
+
+    // Run the clusterization.
+    auto cells_data = vecmem::get_data(cells);
+    const clustering_algorithm::output_type measurements =
+        m_clusterization(cells_data, det_descr_data);
+
+    // If we have a Detray detector, run the seeding track finding and fitting.
+    if (m_detector != nullptr) {
+
+        // Run the seed-finding.
+        const measurement_collection_types::const_view measurements_view =
+            vecmem::get_data(measurements);
+        const spacepoint_formation_algorithm::output_type spacepoints =
+            m_spacepoint_formation(*m_detector, measurements_view);
+        const edm::spacepoint_collection::const_data spacepoints_data =
+            vecmem::get_data(spacepoints);
+        const host::seeding_algorithm::output_type seeds =
+            m_seeding(spacepoints_data);
+        const edm::seed_collection::const_data seeds_data =
+            vecmem::get_data(seeds);
+        return m_track_parameter_estimation(measurements_view, spacepoints_data,
+                                            seeds_data, m_field_vec);
+    }
+    // If not, just return an empty object.
+    else {
+
+        // Return an empty object.
+        return {};
+    }
+}
+
 }  // namespace traccc

examples/run/cpu/full_chain_algorithm.hpp

@@ -11,6 +11,7 @@
 #include "traccc/bfield/magnetic_field.hpp"
 #include "traccc/clusterization/clusterization_algorithm.hpp"
 #include "traccc/edm/silicon_cell_collection.hpp"
+#include "traccc/edm/track_parameters.hpp"
 #include "traccc/edm/track_state.hpp"
 #include "traccc/finding/combinatorial_kalman_filter_algorithm.hpp"
 #include "traccc/fitting/kalman_fitting_algorithm.hpp"
@@ -87,6 +88,14 @@ class full_chain_algorithm : public algorithm<track_state_container_types::host(
     output_type operator()(
         const edm::silicon_cell_collection::host& cells) const override;
 
+    /// Reconstruct track seeds in the entire detector
+    ///
+    /// @param cells The cells for every detector module in the event
+    /// @return The track seeds reconstructed
+    ///
+    bound_track_parameters_collection_types::host seeding(
+        const edm::silicon_cell_collection::host& cells) const;
+
     private:
     /// Vecmem copy object
     std::unique_ptr<vecmem::copy> m_copy;

examples/run/cuda/full_chain_algorithm.cpp

@@ -206,4 +206,46 @@ full_chain_algorithm::output_type full_chain_algorithm::operator()(
     }
 }
 
+bound_track_parameters_collection_types::host full_chain_algorithm::seeding(
+    const edm::silicon_cell_collection::host& cells) const {
+
+    // Create device copy of input collections
+    edm::silicon_cell_collection::buffer cells_buffer(
+        static_cast<unsigned int>(cells.size()), *m_cached_device_mr);
+    m_copy(vecmem::get_data(cells), cells_buffer)->ignore();
+
+    // Run the clusterization (asynchronously).
+    const measurement_collection_types::buffer measurements =
+        m_clusterization(cells_buffer, m_device_det_descr);
+    m_measurement_sorting(measurements);
+
+    // If we have a Detray detector, run the seeding, track finding and fitting.
+    if (m_detector != nullptr) {
+
+        // Run the seed-finding (asynchronously).
+        const spacepoint_formation_algorithm::output_type spacepoints =
+            m_spacepoint_formation(m_device_detector_view, measurements);
+        const track_params_estimation::output_type track_params =
+            m_track_parameter_estimation(measurements, spacepoints,
+                                         m_seeding(spacepoints), m_field_vec);
+
+        // Copy a limited amount of result data back to the host.
+        bound_track_parameters_collection_types::host result{&m_host_mr};
+        m_copy(track_params, result)->wait();
+        return result;
+
+    }
+    // If not, copy the measurements back to the host, and return a dummy
+    // object.
+    else {
+
+        // Copy the measurements back to the host.
+        measurement_collection_types::host measurements_host(&m_host_mr);
+        m_copy(measurements, measurements_host)->wait();
+
+        // Return an empty object.
+        return {};
+    }
+}
+
 }  // namespace traccc::cuda

examples/run/cuda/full_chain_algorithm.hpp

@@ -18,6 +18,7 @@
 #include "traccc/cuda/seeding/track_params_estimation.hpp"
 #include "traccc/cuda/utils/stream.hpp"
 #include "traccc/edm/silicon_cell_collection.hpp"
+#include "traccc/edm/track_parameters.hpp"
 #include "traccc/edm/track_state.hpp"
 #include "traccc/geometry/detector.hpp"
 #include "traccc/geometry/silicon_detector_description.hpp"
@@ -109,6 +110,14 @@ class full_chain_algorithm
     output_type operator()(
         const edm::silicon_cell_collection::host& cells) const override;
 
+    /// Reconstruct track seeds in the entire detector
+    ///
+    /// @param cells The cells for every detector module in the event
+    /// @return The track seeds reconstructed
+    ///
+    bound_track_parameters_collection_types::host seeding(
+        const edm::silicon_cell_collection::host& cells) const;
+
     private:
     /// Host memory resource
     vecmem::memory_resource& m_host_mr;