src/utils/event_data.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2024-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 "traccc/utils/event_data.hpp"
0010
0011 #include "traccc/edm/measurement_collection.hpp"
0012 #include "traccc/io/csv/make_cell_reader.hpp"
0013 #include "traccc/io/csv/make_hit_reader.hpp"
0014 #include "traccc/io/csv/make_measurement_edm.hpp"
0015 #include "traccc/io/csv/make_measurement_hit_id_reader.hpp"
0016 #include "traccc/io/csv/make_measurement_reader.hpp"
0017 #include "traccc/io/csv/make_particle_reader.hpp"
0018 #include "traccc/io/read_cells.hpp"
0019 #include "traccc/io/read_digitization_config.hpp"
0020 #include "traccc/io/utils.hpp"
0021 #include "traccc/utils/particle.hpp"
0022
0023 // System include(s).
0024 #include <algorithm>
0025 #include <filesystem>
0026
0027 namespace traccc {
0028
0029 event_data::event_data(const std::string& event_dir, const std::size_t event_id,
0030                        vecmem::memory_resource& resource,
0031                        bool use_acts_geom_source, const host_detector* det,
0032                        data_format format, bool include_silicon_cells)
0033     : m_measurements(resource),
0034       m_event_dir(event_dir),
0035       m_event_id(event_id),
0036       m_mr(resource) {
0037   // Currently, we only support csv type for event data
0038   assert(format == data_format::csv);
0039   if (format == data_format::csv) {
0040     setup_csv(use_acts_geom_source, det, include_silicon_cells);
0041   }
0042 }
0043
0044 void event_data::setup_csv(bool use_acts_geom_source, const host_detector* det,
0045                            bool include_silicon_cells) {
0046   /********************
0047    *  Read Csv files  *
0048    ********************/
0049
0050   // CSV IO EDM containers
0051   std::vector<io::csv::cell> csv_cells;
0052   std::vector<io::csv::hit> csv_hits;
0053   std::vector<io::csv::measurement> csv_measurements;
0054   std::vector<io::csv::measurement_hit_id> csv_meas_hit_ids;
0055   std::vector<io::csv::particle> csv_particles;
0056
0057   if (include_silicon_cells) {
0058     // Read the cells from the relevant event file
0059     std::string io_cells_file =
0060         io::get_absolute_path((std::filesystem::path(m_event_dir) /
0061                                std::filesystem::path(io::get_event_filename(
0062                                    m_event_id, "-cells.csv")))
0063                                   .native());
0064
0065     auto creader = io::csv::make_cell_reader(io_cells_file);
0066     io::csv::cell iocell;
0067
0068     while (creader.read(iocell)) {
0069       csv_cells.push_back(iocell);
0070     }
0071   }
0072
0073   // Read the hits from the relevant event file
0074   std::string io_hits_file = io::get_absolute_path(
0075       (std::filesystem::path(m_event_dir) /
0076        std::filesystem::path(io::get_event_filename(m_event_id, "-hits.csv")))
0077           .native());
0078
0079   auto hreader = io::csv::make_hit_reader(io_hits_file);
0080   {
0081     io::csv::hit iohit;
0082     while (hreader.read(iohit)) {
0083       csv_hits.push_back(iohit);
0084     }
0085   }
0086
0087   // Read the measurements from the relevant event file
0088   std::string io_measurements_file =
0089       io::get_absolute_path((std::filesystem::path(m_event_dir) /
0090                              std::filesystem::path(io::get_event_filename(
0091                                  m_event_id, "-measurements.csv")))
0092                                 .native());
0093
0094   auto mreader = io::csv::make_measurement_reader(io_measurements_file);
0095   {
0096     io::csv::measurement iomeas;
0097     while (mreader.read(iomeas)) {
0098       csv_measurements.push_back(iomeas);
0099     }
0100   }
0101
0102   // Read the measurement hit id from the relevant event file
0103   std::string io_measurement_hit_id_file =
0104       io::get_absolute_path((std::filesystem::path(m_event_dir) /
0105                              std::filesystem::path(io::get_event_filename(
0106                                  m_event_id, "-measurement-simhit-map.csv")))
0107                                 .native());
0108
0109   auto mhid_reader =
0110       io::csv::make_measurement_hit_id_reader(io_measurement_hit_id_file);
0111   {
0112     io::csv::measurement_hit_id mh_id;
0113     while (mhid_reader.read(mh_id)) {
0114       csv_meas_hit_ids.push_back(mh_id);
0115     }
0116   }
0117
0118   // Read the particles from the relevant event file
0119   std::string io_particles_file =
0120       io::get_absolute_path((std::filesystem::path(m_event_dir) /
0121                              std::filesystem::path(io::get_event_filename(
0122                                  m_event_id, "-particles_initial.csv")))
0123                                 .native());
0124
0125   auto preader = io::csv::make_particle_reader(io_particles_file);
0126   {
0127     io::csv::particle ioptc;
0128     while (preader.read(ioptc)) {
0129       csv_particles.push_back(ioptc);
0130     }
0131   }
0132
0133   /********************
0134    * Make geom_id map *
0135    ********************/
0136
0137   // For Acts data, build a map of acts->detray geometry IDs
0138   std::map<geometry_id, geometry_id> acts_to_detray_id;
0139   if (use_acts_geom_source) {
0140     host_detector_visitor<detector_type_list>(
0141         *det, [&acts_to_detray_id]<typename detector_traits_t>(
0142                   const typename detector_traits_t::host& d) {
0143           for (const auto& surface_desc : d.surfaces()) {
0144             acts_to_detray_id[surface_desc.source] =
0145                 surface_desc.identifier().value();
0146           }
0147         });
0148   }
0149
0150   /********************
0151    * Make EDM maps    *
0152    ********************/
0153
0154   bool csv_measurements_have_time = false;
0155
0156   for (const auto& iomeas : csv_measurements) {
0157     if (std::fabs(iomeas.time) != 0.f) {
0158       csv_measurements_have_time = true;
0159       break;
0160     }
0161   }
0162
0163   // TODO: Put some proper logging here in a future PR
0164   if (csv_measurements_have_time) {
0165     std::cout << "Using measurement time" << std::endl;
0166   } else {
0167     std::cout << "Using hit time" << std::endl;
0168   }
0169
0170   // Measurement map
0171   for (const auto& iomeas : csv_measurements) {
0172     // Construct the measurement object.
0173     m_measurements.resize(m_measurements.size() + 1u);
0174     auto meas = m_measurements.at(m_measurements.size() - 1u);
0175     if (use_acts_geom_source) {
0176       traccc::io::csv::make_measurement_edm(iomeas, meas, &acts_to_detray_id);
0177     } else {
0178       traccc::io::csv::make_measurement_edm(iomeas, meas, nullptr);
0179     }
0180
0181     if (!csv_measurements_have_time) {
0182       const auto hid = csv_meas_hit_ids.at(iomeas.measurement_id).hit_id;
0183       const auto& iohit = csv_hits.at(hid);
0184
0185       meas.time() = iohit.tt;
0186     }
0187
0188     if (iomeas.measurement_id <=
0189         std::numeric_limits<measurement_id_type>::max()) {
0190       m_measurement_map[static_cast<measurement_id_type>(
0191           iomeas.measurement_id)] = meas;
0192     } else {
0193       throw std::runtime_error("Measurement ID exceeds the bound");
0194     }
0195   }
0196
0197   // Particle map
0198   for (const auto& ioptc : csv_particles) {
0199     point3 pos{ioptc.vx, ioptc.vy, ioptc.vz};
0200     vector3 mom{ioptc.px, ioptc.py, ioptc.pz};
0201
0202     m_particle_map[ioptc.particle_id] =
0203         traccc::particle{ioptc.particle_id, ioptc.particle_type,
0204                          ioptc.process,     pos,
0205                          ioptc.vt,          mom,
0206                          ioptc.m,           ioptc.q};
0207   }
0208
0209   /************************************
0210    * Make measurement to particle map *
0211    ************************************/
0212
0213   // When including silicon cells
0214   if (include_silicon_cells) {
0215     std::map<measurement_proxy, std::vector<io::csv::cell>> meas_to_cluster_map;
0216
0217     for (const auto& iocell : csv_cells) {
0218       // Fill the measurement_to_cluster_map
0219       auto meas_id = iocell.measurement_id;
0220       auto hid = csv_meas_hit_ids[meas_id].hit_id;
0221       const auto& iohit = csv_hits[hid];
0222
0223       measurement_proxy meas =
0224           m_measurement_map.at(static_cast<measurement_id_type>(meas_id));
0225       meas_to_cluster_map[meas].push_back(iocell);
0226
0227       const auto& ptc = m_particle_map.at(iohit.particle_id);
0228       m_cell_to_particle_map[iocell] = ptc;
0229     }
0230
0231     // Fill the meas_to_particle_map
0232     for (auto const& [ms, cluster] : meas_to_cluster_map) {
0233       for (const auto& cell : cluster) {
0234         const auto& ptc = m_cell_to_particle_map.at(cell);
0235         m_meas_to_ptc_map[ms][ptc]++;
0236       }
0237     }
0238   }
0239
0240   for (const auto& iomeas : csv_measurements) {
0241     // Hit index
0242     const auto hid = csv_meas_hit_ids.at(iomeas.measurement_id).hit_id;
0243
0244     // Make spacepoint
0245     const auto& iohit = csv_hits.at(hid);
0246     point3 global_pos{iohit.tx, iohit.ty, iohit.tz};
0247     point3 global_mom{iohit.tpx, iohit.tpy, iohit.tpz};
0248
0249     // Make particle
0250     const auto& ptc = m_particle_map.at(iohit.particle_id);
0251
0252     // Construct the measurement object.
0253     measurement_proxy meas = m_measurement_map.at(
0254         static_cast<measurement_id_type>(iomeas.measurement_id));
0255
0256     // Fill measurement to truth global position and momentum map
0257     m_meas_to_param_map[meas] = std::make_pair(global_pos, global_mom);
0258
0259     // Fill particle to measurement map
0260     auto& meas_vec = m_ptc_to_meas_map[ptc];
0261
0262     meas_vec.insert(std::upper_bound(meas_vec.begin(), meas_vec.end(), meas,
0263                                      [](const measurement_proxy& val,
0264                                         const measurement_proxy& old) {
0265                                        return old.time() > val.time();
0266                                      }),
0267                     meas);
0268
0269     if (!include_silicon_cells) {
0270       auto insert_return = m_meas_to_ptc_map.insert({meas, {}});
0271       if (insert_return.second == false) {
0272         // TODO: Put some logging here when that's ready
0273       } else {
0274         (*(insert_return.first)).second[ptc] = 1u;
0275       }
0276     }
0277   }
0278 }
0279
0280 void event_data::fill_cca_result(
0281     const edm::silicon_cell_collection::host& cells,
0282     const edm::silicon_cluster_collection::host& cca_clusters,
0283     const edm::measurement_collection::host& cca_measurements,
0284     const detector_conditions_description::host& det_cond) {
0285   const std::size_t n_cca_clusters = cca_measurements.size();
0286
0287   std::map<measurement_proxy, std::vector<io::csv::cell>>
0288       found_meas_to_cluster_map;
0289
0290   for (std::size_t i = 0; i < n_cca_clusters; i++) {
0291     const auto meas = cca_measurements.at(i);
0292     const auto cluster = cca_clusters[i];
0293
0294     std::vector<io::csv::cell> iocells;
0295     for (const unsigned int cell_idx : cluster.cell_indices()) {
0296       const auto cell = cells.at(cell_idx);
0297       io::csv::cell iocell{det_cond.acts_geometry_id().at(cell.module_index()),
0298                            0u,
0299                            cell.channel0(),
0300                            cell.channel1(),
0301                            static_cast<float>(cell.time()),
0302                            static_cast<float>(cell.activation())};
0303
0304       iocells.push_back(iocell);
0305     }
0306     found_meas_to_cluster_map[meas] = iocells;
0307   }
0308
0309   // Construct a mapping of geometry IDs to cells and the associated
0310   // particles, in order to get much faster lookups later.
0311   std::map<geometry_id, std::vector<std::reference_wrapper<const std::decay_t<
0312                             decltype(m_cell_to_particle_map)>::value_type>>>
0313       geo_id_to_cell_to_particle_map;
0314
0315   for (auto const& v : m_cell_to_particle_map) {
0316     geo_id_to_cell_to_particle_map[v.first.geometry_id].emplace_back(v);
0317   }
0318
0319   for (auto const& [ms, cluster] : found_meas_to_cluster_map) {
0320     std::map<uint64_t, std::size_t> meas_counts;
0321
0322     // Cells from CCL
0323     for (const auto& cell1 : cluster) {
0324       // Cells from truth [cell, particle] map
0325       for (auto const& v : geo_id_to_cell_to_particle_map[cell1.geometry_id]) {
0326         const auto& cell2 = v.get().first;
0327         const auto& ptc = v.get().second;
0328         assert(cell1.geometry_id == cell2.geometry_id);
0329         // Increase the particle number if the cell is the same
0330         if (cell1.channel0 == cell2.channel0 &&
0331             cell1.channel1 == cell2.channel1) {
0332           m_found_meas_to_ptc_map[ms][ptc]++;
0333           meas_counts[cell2.measurement_id]++;
0334         }
0335       }
0336     }
0337
0338     // Find most contributing measurement and its corresponding hit
0339     using pair_type = decltype(meas_counts)::value_type;
0340     auto pr = std::max_element(std::begin(meas_counts), std::end(meas_counts),
0341                                [](const pair_type& p1, const pair_type& p2) {
0342                                  return p1.second < p2.second;
0343                                });
0344
0345     const auto& meas_id = pr->first;
0346
0347     if (meas_id <= std::numeric_limits<measurement_id_type>::max()) {
0348       m_found_meas_to_param_map[ms] = m_meas_to_param_map
0349           [m_measurement_map[static_cast<measurement_id_type>(meas_id)]];
0350
0351     } else {
0352       throw std::runtime_error("Measurement ID exceeds the bound");
0353     }
0354   }
0355 }
0356 }  // namespace traccc