src/csv/read_particles.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2022-2026 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 // Local include(s).
0009 #include "read_particles.hpp"
0010
0011 #include "read_measurements.hpp"
0012 #include "traccc/io/csv/make_hit_reader.hpp"
0013 #include "traccc/io/csv/make_measurement_hit_id_reader.hpp"
0014 #include "traccc/io/csv/make_particle_reader.hpp"
0015
0016 // VecMem include(s).
0017 #include <vecmem/containers/jagged_vector.hpp>
0018 #include <vecmem/memory/host_memory_resource.hpp>
0019
0020 // System include(s).
0021 #include <algorithm>
0022 #include <stdexcept>
0023 #include <unordered_map>
0024
0025 namespace traccc::io::csv {
0026
0027 void read_particles(particle_collection_types::host& particles,
0028                     std::string_view filename) {
0029   // Construct the particle reader object.
0030   auto reader = make_particle_reader(filename);
0031
0032   // Read the particles from the input file.
0033   csv::particle part;
0034   while (reader.read(part)) {
0035     particles.push_back({part.particle_id, part.particle_type, part.process,
0036                          point3{part.vx, part.vy, part.vz}, part.vt,
0037                          vector3{part.px, part.py, part.pz}, part.m, part.q});
0038   }
0039 }
0040
0041 void read_particles(
0042     particle_container_types::host& particles,
0043     edm::measurement_collection::host& measurements,
0044     std::string_view particles_file, std::string_view hits_file,
0045     std::string_view measurements_file, std::string_view hit_map_file,
0046     const traccc::host_detector* detector,
0047     const traccc::detector_design_description::host* det_desc,
0048     const traccc::detector_conditions_description::host* det_cond,
0049     const bool sort_measurements) {
0050   // Memory resource used by the temporary collections.
0051   vecmem::host_memory_resource mr;
0052
0053   // Construct all necessary reader objects.
0054   auto hit_reader = make_hit_reader(hits_file);
0055   auto measurement_hit_id_reader = make_measurement_hit_id_reader(hit_map_file);
0056
0057   // Read in all particles, into a temporary collection.
0058   particle_collection_types::host temp_particles{&mr};
0059   read_particles(temp_particles, particles_file);
0060
0061   // Read in all measurements.
0062   const std::vector<measurement_id_type> new_idx_map =
0063       read_measurements(measurements, measurements_file, detector, det_desc,
0064                         det_cond, sort_measurements);
0065
0066   // Make a hit to measurement map.
0067   std::unordered_map<std::size_t, measurement_id_type> hit_to_measurement;
0068   measurement_hit_id mhid;
0069   while (measurement_hit_id_reader.read(mhid)) {
0070     if (sort_measurements) {
0071       mhid.measurement_id = new_idx_map[mhid.measurement_id];
0072     }
0073     if (hit_to_measurement.insert({mhid.hit_id, mhid.measurement_id}).second ==
0074         false) {
0075       throw std::runtime_error("Duplicate hit ID in hit->measurement map");
0076     }
0077   }
0078
0079   // Construct the indices of the measurements belonging to each particle.
0080   vecmem::jagged_vector<unsigned int> particle_measurements{
0081       temp_particles.size(), vecmem::vector<unsigned int>{&mr}, &mr};
0082   hit h;
0083   std::size_t hit_id = 0u;
0084   while (hit_reader.read(h)) {
0085     // Find the particle belonging to this hit.
0086     auto particle_it =
0087         std::find_if(temp_particles.begin(), temp_particles.end(),
0088                      [&h](const traccc::particle& p) {
0089                        return p.particle_id == h.particle_id;
0090                      });
0091     if (particle_it == temp_particles.end()) {
0092       throw std::runtime_error("Hit without corresponding particle");
0093     }
0094
0095     // Find the measurement belonging to this hit.
0096     auto hit_to_measurement_it = hit_to_measurement.find(hit_id);
0097     if (hit_to_measurement_it == hit_to_measurement.end()) {
0098       throw std::runtime_error("Hit without corresponding measurement");
0099     }
0100
0101     // Find the index of the found particle in its collection.
0102     auto particle_index = std::distance(temp_particles.begin(), particle_it);
0103
0104     // Add the measurement to the particle's collection.
0105     particle_measurements
0106         [static_cast<decltype(particle_measurements)::size_type>(
0107              particle_index)]
0108             .push_back(hit_to_measurement_it->second);
0109
0110     // Increment the hit ID.
0111     ++hit_id;
0112   }
0113
0114   // Construct the final particle collection.
0115   for (std::size_t i = 0; i < temp_particles.size(); ++i) {
0116     particles.push_back(temp_particles[i], particle_measurements[i]);
0117   }
0118 }
0119
0120 }  // namespace traccc::io::csv