src/efficiency/seeding_performance_writer.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 "traccc/efficiency/seeding_performance_writer.hpp"
0010
0011 #include "duplication_plot_tool.hpp"
0012 #include "eff_plot_tool.hpp"
0013 #include "track_classification.hpp"
0014
0015 // ROOT include(s).
0016 #ifdef TRACCC_HAVE_ROOT
0017 #include <TFile.h>
0018 #endif  // TRACCC_HAVE_ROOT
0019
0020 // System include(s).
0021 #include <iostream>
0022 #include <memory>
0023 #include <stdexcept>
0024
0025 namespace traccc {
0026 namespace details {
0027
0028 struct seeding_performance_writer_data {
0029   /// Constructor
0030   seeding_performance_writer_data(const seeding_performance_writer::config& cfg)
0031       : m_eff_plot_tool({cfg.var_binning}),
0032         m_duplication_plot_tool({cfg.var_binning}) {}
0033
0034   /// Plot tool for efficiency
0035   eff_plot_tool m_eff_plot_tool;
0036   eff_plot_tool::eff_plot_cache m_eff_plot_cache;
0037
0038   /// Plot tool for duplication rate
0039   duplication_plot_tool m_duplication_plot_tool;
0040   duplication_plot_tool::duplication_plot_cache m_duplication_plot_cache;
0041
0042   std::map<event_data::measurement_proxy, std::map<particle, std::size_t>>
0043       m_measurement_particle_map;
0044   std::map<std::uint64_t, particle> m_particle_map;
0045
0046 };  // struct seeding_performance_writer_data
0047
0048 }  // namespace details
0049
0050 seeding_performance_writer::seeding_performance_writer(
0051     const config& cfg, std::unique_ptr<const traccc::Logger> logger)
0052     : messaging(std::move(logger)),
0053       m_cfg(cfg),
0054       m_data(std::make_unique<details::seeding_performance_writer_data>(cfg)) {
0055   m_data->m_eff_plot_tool.book("seeding", m_data->m_eff_plot_cache);
0056   m_data->m_duplication_plot_tool.book("seeding",
0057                                        m_data->m_duplication_plot_cache);
0058 }
0059
0060 seeding_performance_writer::~seeding_performance_writer() {}
0061
0062 void seeding_performance_writer::write(
0063     const edm::seed_collection::const_view& seeds_view,
0064     const edm::spacepoint_collection::const_view& spacepoints_view,
0065     const edm::measurement_collection::const_view& measurements_view,
0066     const event_data& evt_data) {
0067   std::map<particle_id, std::size_t> match_counter;
0068
0069   // Create the device collections.
0070   const edm::seed_collection::const_device seeds(seeds_view);
0071   const edm::spacepoint_collection::const_device spacepoints(spacepoints_view);
0072   const edm::measurement_collection::const_device measurements(
0073       measurements_view);
0074
0075   const std::size_t n_seeds = seeds.size();
0076
0077   std::size_t total_fake_seeds = 0;
0078
0079   // Iterate over the seeds.
0080   for (edm::seed_collection::const_device::size_type i = 0u; i < seeds.size();
0081        ++i) {
0082     const auto sd = seeds.at(i);
0083
0084     std::vector<particle_hit_count> particle_hit_counts;
0085
0086     if (sd.bottom_index() >= spacepoints.size() ||
0087         sd.middle_index() >= spacepoints.size() ||
0088         sd.top_index() >= spacepoints.size()) {
0089       TRACCC_WARNING("Seed with index "
0090                      << i << " (out of " << seeds.size()
0091                      << ") has invalid spacepoint (B: " << sd.bottom_index()
0092                      << ", M: " << sd.middle_index()
0093                      << ", T: " << sd.top_index()
0094                      << ", #SP: " << spacepoints.size() << ")");
0095       continue;
0096     }
0097
0098     if (spacepoints.at(sd.bottom_index()).measurement_index_1() >=
0099             measurements.size() ||
0100         spacepoints.at(sd.middle_index()).measurement_index_1() >=
0101             measurements.size() ||
0102         spacepoints.at(sd.top_index()).measurement_index_1() >=
0103             measurements.size()) {
0104       TRACCC_WARNING(
0105           "Seed with index "
0106           << i << " (out of " << seeds.size()
0107           << ") has invalid measurement (B: "
0108           << spacepoints.at(sd.bottom_index()).measurement_index_1()
0109           << ", M: " << spacepoints.at(sd.middle_index()).measurement_index_1()
0110           << ", T: " << spacepoints.at(sd.top_index()).measurement_index_1()
0111           << ", #SP: " << spacepoints.size() << ")");
0112       continue;
0113     }
0114
0115     // Get the measurements for this seed.
0116     assert(spacepoints.at(sd.bottom_index()).measurement_index_2() ==
0117            edm::spacepoint_collection::host::INVALID_MEASUREMENT_INDEX);
0118     assert(spacepoints.at(sd.middle_index()).measurement_index_2() ==
0119            edm::spacepoint_collection::host::INVALID_MEASUREMENT_INDEX);
0120     assert(spacepoints.at(sd.top_index()).measurement_index_2() ==
0121            edm::spacepoint_collection::host::INVALID_MEASUREMENT_INDEX);
0122     std::array<event_data::measurement_proxy, 3> seed_measurements{
0123         measurements.at(
0124             spacepoints.at(sd.bottom_index()).measurement_index_1()),
0125         measurements.at(
0126             spacepoints.at(sd.middle_index()).measurement_index_1()),
0127         measurements.at(spacepoints.at(sd.top_index()).measurement_index_1())};
0128
0129     if (!evt_data.m_found_meas_to_ptc_map.empty()) {
0130       particle_hit_counts = identify_contributing_particles(
0131           seed_measurements, evt_data.m_found_meas_to_ptc_map);
0132     } else {
0133       particle_hit_counts = identify_contributing_particles(
0134           seed_measurements, evt_data.m_meas_to_ptc_map);
0135     }
0136
0137     // Consider it being matched if hit counts is larger than the half
0138     // of the number of measurements
0139     assert(seed_measurements.size() > 0u);
0140     if (static_cast<double>(particle_hit_counts.at(0).hit_counts) /
0141             static_cast<double>(seed_measurements.size()) >
0142         m_cfg.seed_truth_config.matching_ratio) {
0143       auto pid = particle_hit_counts.at(0).ptc.particle_id;
0144       match_counter[pid]++;
0145     } else {
0146       total_fake_seeds++;
0147     }
0148   }
0149
0150   std::size_t total_truth_particles = 0;
0151   std::size_t total_matched_truth_particles = 0;
0152   std::size_t total_duplicate_seeds = 0;
0153
0154   for (auto const& [pid, ptc] : evt_data.m_particle_map) {
0155     std::size_t num_measurements = 0;
0156
0157     // Find the number of measurements belonging to this track
0158     if (auto it = evt_data.m_ptc_to_meas_map.find(ptc);
0159         it != evt_data.m_ptc_to_meas_map.end()) {
0160       num_measurements = it->second.size();
0161     } else {
0162       continue;
0163     }
0164
0165     // Count only charged particles which satisfiy pT_cut
0166     if (ptc.charge == 0 ||
0167         vector::perp(ptc.momentum) < m_cfg.truth_config.pT_min ||
0168         ptc.vertex[2] < m_cfg.truth_config.z_min ||
0169         ptc.vertex[2] > m_cfg.truth_config.z_max ||
0170         vector::perp(ptc.vertex) > m_cfg.truth_config.r_max ||
0171         std::abs(vector::eta(ptc.momentum)) > m_cfg.truth_config.eta_max ||
0172         (m_cfg.truth_config.process_id >= 0 &&
0173          m_cfg.truth_config.process_id != ptc.process) ||
0174         num_measurements < m_cfg.truth_config.min_track_candidates) {
0175       continue;
0176     }
0177
0178     total_truth_particles++;
0179
0180     bool is_matched = false;
0181     std::size_t n_matched_seeds_for_particle = 0;
0182     auto it = match_counter.find(pid);
0183     if (it != match_counter.end()) {
0184       is_matched = true;
0185       n_matched_seeds_for_particle = it->second;
0186       total_matched_truth_particles++;
0187       total_duplicate_seeds += n_matched_seeds_for_particle - 1;
0188     }
0189
0190     m_data->m_eff_plot_tool.fill(m_data->m_eff_plot_cache, ptc, is_matched);
0191     m_data->m_duplication_plot_tool.fill(m_data->m_duplication_plot_cache, ptc,
0192                                          n_matched_seeds_for_particle - 1);
0193   }
0194
0195   TRACCC_INFO("Total number of truth particles was " << total_truth_particles);
0196   TRACCC_INFO("Total number of found seeds was " << n_seeds);
0197   TRACCC_INFO("Total number of seed-matched particles was "
0198               << total_matched_truth_particles);
0199   TRACCC_INFO("Total number of duplicate seeds was " << total_duplicate_seeds);
0200   TRACCC_INFO("Total number of fake seeds was " << total_fake_seeds);
0201   TRACCC_INFO("Total seed efficiency was "
0202               << (100. * static_cast<double>(total_matched_truth_particles) /
0203                   static_cast<double>(total_truth_particles))
0204               << "%");
0205   TRACCC_INFO("Total seed duplicate rate was "
0206               << (static_cast<double>(total_duplicate_seeds) /
0207                   static_cast<double>(total_matched_truth_particles)));
0208   TRACCC_INFO("Total seed fake rate was "
0209               << (static_cast<double>(total_fake_seeds) /
0210                   static_cast<double>(total_truth_particles)));
0211 }
0212
0213 void seeding_performance_writer::finalize() {
0214 #ifdef TRACCC_HAVE_ROOT
0215   // Open the output file.
0216   std::unique_ptr<TFile> ofile(
0217       TFile::Open(m_cfg.file_path.c_str(), m_cfg.file_mode.c_str()));
0218   if ((!ofile) || ofile->IsZombie()) {
0219     throw std::runtime_error("Could not open output file \"" + m_cfg.file_path +
0220                              "\" in mode \"" + m_cfg.file_mode + "\"");
0221   }
0222   ofile->cd();
0223 #else
0224   std::cout << "ROOT file \"" << m_cfg.file_path << "\" is NOT created"
0225             << std::endl;
0226 #endif  // TRACCC_HAVE_ROOT
0227
0228   m_data->m_eff_plot_tool.write(m_data->m_eff_plot_cache);
0229   m_data->m_duplication_plot_tool.write(m_data->m_duplication_plot_cache);
0230 }
0231
0232 }  // namespace traccc