src/efficiency/finding_performance_writer.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2023-2025 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/finding_performance_writer.hpp"
0010
0011 #include "../resolution/stat_plot_tool.hpp"
0012 #include "duplication_plot_tool.hpp"
0013 #include "eff_plot_tool.hpp"
0014 #include "fake_tracks_plot_tool.hpp"
0015 #include "traccc/edm/track_fit_outcome.hpp"
0016 #include "traccc/utils/logging.hpp"
0017 #include "track_classification.hpp"
0018
0019 // ROOT include(s).
0020 #ifdef TRACCC_HAVE_ROOT
0021 #include <TFile.h>
0022 #endif  // TRACCC_HAVE_ROOT
0023
0024 // System include(s).
0025 #include <iostream>
0026 #include <memory>
0027 #include <stdexcept>
0028
0029 namespace traccc {
0030 namespace details {
0031
0032 struct finding_performance_writer_data {
0033   /// Constructor
0034   finding_performance_writer_data(const finding_performance_writer::config& cfg)
0035       : m_eff_plot_tool({cfg.var_binning}),
0036         m_duplication_plot_tool({cfg.var_binning}),
0037         m_fake_tracks_plot_tool({cfg.var_binning}),
0038         m_stat_plot_tool(cfg.stat_config) {}
0039
0040   /// Plot tool for efficiency
0041   eff_plot_tool m_eff_plot_tool;
0042   eff_plot_tool::eff_plot_cache m_eff_plot_cache;
0043
0044   /// Plot tool for duplication rate
0045   duplication_plot_tool m_duplication_plot_tool;
0046   duplication_plot_tool::duplication_plot_cache m_duplication_plot_cache;
0047
0048   // Plot tool for fake tracks monitoring
0049   fake_tracks_plot_tool m_fake_tracks_plot_tool;
0050   fake_tracks_plot_tool::fake_tracks_plot_cache m_fake_tracks_plot_cache;
0051
0052   std::map<event_data::measurement_proxy, std::map<particle, std::size_t>>
0053       m_measurement_particle_map;
0054   std::map<std::uint64_t, particle> m_particle_map;
0055
0056   /// Plot tool for statistics
0057   stat_plot_tool m_stat_plot_tool;
0058   stat_plot_tool::stat_plot_cache m_stat_plot_cache;
0059
0060 };  // struct finding_performance_writer_data
0061
0062 }  // namespace details
0063
0064 finding_performance_writer::finding_performance_writer(
0065     const config& cfg, std::unique_ptr<const traccc::Logger> logger)
0066     : messaging(std::move(logger)),
0067       m_cfg(cfg),
0068       m_data(std::make_unique<details::finding_performance_writer_data>(cfg)) {
0069   m_data->m_eff_plot_tool.book(m_cfg.algorithm_name, m_data->m_eff_plot_cache);
0070   m_data->m_duplication_plot_tool.book(m_cfg.algorithm_name,
0071                                        m_data->m_duplication_plot_cache);
0072   m_data->m_fake_tracks_plot_tool.book(m_cfg.algorithm_name,
0073                                        m_data->m_fake_tracks_plot_cache);
0074   m_data->m_stat_plot_tool.book(m_data->m_stat_plot_cache);
0075 }
0076
0077 finding_performance_writer::~finding_performance_writer() {}
0078
0079 namespace {
0080
0081 /**
0082  * @brief For ambiguity resolution only. Associates each reconstructed track
0083  * with its measurements.
0084  *
0085  * @param track_view the track candidates found by the finding algorithm.
0086  * @return std::vector<std::vector<measurement>> Associates each track index
0087  * with its corresponding measurements.
0088  */
0089 std::vector<std::vector<event_data::measurement_proxy>> prepare_data(
0090     const edm::track_container<default_algebra>::const_view& track_view,
0091     bool require_fit = false) {
0092   std::vector<std::vector<event_data::measurement_proxy>> result;
0093
0094   // Set up the input containers.
0095   const edm::track_container<default_algebra>::const_device tracks(track_view);
0096
0097   // Iterate over the tracks.
0098   const unsigned int n_tracks = tracks.tracks.size();
0099   result.reserve(n_tracks);
0100
0101   for (unsigned int i = 0; i < n_tracks; i++) {
0102     if (require_fit &&
0103         tracks.tracks.at(i).fit_outcome() != track_fit_outcome::SUCCESS) {
0104       continue;
0105     }
0106     std::vector<event_data::measurement_proxy> result_measurements;
0107     for (const edm::track_constituent_link& link :
0108          tracks.tracks.constituent_links().at(i)) {
0109       if (link.type == edm::track_constituent_link::measurement) {
0110         result_measurements.push_back(tracks.measurements.at(link.index));
0111       } else if (link.type == edm::track_constituent_link::track_state) {
0112         result_measurements.push_back(tracks.measurements.at(
0113             tracks.states.at(link.index).measurement_index()));
0114       }
0115     }
0116     result.push_back(std::move(result_measurements));
0117   }
0118   return result;
0119 }
0120
0121 }  // namespace
0122
0123 void finding_performance_writer::write_common(
0124     const std::vector<std::vector<event_data::measurement_proxy>>& tracks,
0125     const event_data& evt_data) {
0126   // Associates truth particle_ids with the number of tracks made entirely of
0127   // some (or all) of its hits.
0128   std::map<particle_id, std::size_t> match_counter;
0129
0130   // Associates truth particle_ids with the number of tracks sharing hits from
0131   // more than one truth particle.
0132   std::map<particle_id, std::size_t> fake_counter;
0133
0134   // Iterate over the tracks.
0135   const std::size_t n_tracks = tracks.size();
0136
0137   std::size_t total_fake_tracks = 0;
0138
0139   for (std::size_t i = 0; i < n_tracks; i++) {
0140     const std::vector<event_data::measurement_proxy>& found_measurements =
0141         tracks[i];
0142
0143     // Check which particle matches this seed.
0144     // Input :
0145     //    - the list of measurements for this track
0146     //    - the truth particles map
0147     // Output :
0148     //    - a list of particles, having for each of them a particle_id and
0149     //      a count value.
0150     // If there is only a single truth particle contributing to this track,
0151     // then increment the match_counter for this truth particle id.
0152     // If there are at least two particles contributing to the hit list of
0153     // this track, increment the fake_counter for each truth particle.
0154     std::vector<particle_hit_count> particle_hit_counts;
0155
0156     if (!evt_data.m_found_meas_to_ptc_map.empty()) {
0157       particle_hit_counts = identify_contributing_particles(
0158           found_measurements, evt_data.m_found_meas_to_ptc_map);
0159     } else {
0160       particle_hit_counts = identify_contributing_particles(
0161           found_measurements, evt_data.m_meas_to_ptc_map);
0162     }
0163
0164     if (particle_hit_counts.empty()) {
0165       TRACCC_VERBOSE("No contributing particles found for track: " << n_tracks);
0166       continue;
0167     }
0168
0169     const auto major_ptc = particle_hit_counts.at(0).ptc;
0170     const auto n_major_hits = particle_hit_counts.at(0).hit_counts;
0171
0172     // Truth measureemnt from the particle
0173     const std::vector<event_data::measurement_proxy> truth_measurements =
0174         evt_data.m_ptc_to_meas_map.at(major_ptc);
0175
0176     // Consider it being matched if hit counts is larger than the half
0177     // of the number of measurements
0178     assert(found_measurements.size() > 0u);
0179     assert(truth_measurements.size() > 0u);
0180
0181     const double purity = static_cast<double>(n_major_hits) /
0182                           static_cast<double>(found_measurements.size());
0183     const double completeness = static_cast<double>(n_major_hits) /
0184                                 static_cast<double>(truth_measurements.size());
0185
0186     const bool reco_matched = purity >= m_cfg.track_truth_config.matching_ratio;
0187     const bool truth_matched =
0188         completeness >= m_cfg.track_truth_config.matching_ratio;
0189
0190     m_data->m_stat_plot_tool.fill(m_data->m_stat_plot_cache, purity,
0191                                   completeness);
0192
0193     if ((!m_cfg.track_truth_config.double_matching && reco_matched) ||
0194         (m_cfg.track_truth_config.double_matching && reco_matched &&
0195          truth_matched)) {
0196       const auto pid = major_ptc.particle_id;
0197       match_counter[pid]++;
0198     } else {
0199       for (particle_hit_count const& phc : particle_hit_counts) {
0200         const auto pid = phc.ptc.particle_id;
0201         fake_counter[pid]++;
0202       }
0203       total_fake_tracks++;
0204     }
0205   }
0206
0207   std::size_t total_truth_particles = 0;
0208   std::size_t total_matched_truth_particles = 0;
0209   std::size_t total_duplicate_tracks = 0;
0210
0211   // For each truth particle...
0212   for (auto const& [pid, ptc] : evt_data.m_particle_map) {
0213     auto ptc_particle =
0214         detail::particle_from_pdg_number<scalar>(ptc.particle_type);
0215     if (ptc_particle.pdg_num() == 0) {
0216       // TODO: Add some debug logging here.
0217       continue;
0218     }
0219
0220     // Find the number of measurements belonging to this track
0221     std::size_t num_measurements = 0;
0222     if (auto it = evt_data.m_ptc_to_meas_map.find(ptc);
0223         it != evt_data.m_ptc_to_meas_map.end()) {
0224       num_measurements = it->second.size();
0225     } else {
0226       continue;
0227     }
0228
0229     // Count only charged particles which satisfy pT_cut and vertex cut
0230     if (ptc.charge == 0 ||
0231         vector::perp(ptc.momentum) < m_cfg.truth_config.pT_min ||
0232         ptc.vertex[2] < m_cfg.truth_config.z_min ||
0233         ptc.vertex[2] > m_cfg.truth_config.z_max ||
0234         vector::perp(ptc.vertex) > m_cfg.truth_config.r_max ||
0235         std::abs(vector::eta(ptc.momentum)) > m_cfg.truth_config.eta_max ||
0236         (m_cfg.truth_config.process_id >= 0 &&
0237          m_cfg.truth_config.process_id != ptc.process) ||
0238         num_measurements < m_cfg.truth_config.min_track_candidates) {
0239       continue;
0240     }
0241
0242     total_truth_particles++;
0243
0244     // Finds how many tracks were made solely by hits from the current truth
0245     // particle
0246     bool is_matched = false;
0247     std::size_t n_matched_seeds_for_particle = 0;
0248     auto it = match_counter.find(pid);
0249     if (it != match_counter.end()) {
0250       is_matched = true;
0251       total_matched_truth_particles++;
0252       n_matched_seeds_for_particle = it->second;
0253       total_duplicate_tracks += n_matched_seeds_for_particle - 1;
0254     } else {
0255       TRACCC_DEBUG("Not matched: " << pid);
0256     }
0257
0258     // Finds how many (fake) tracks were made with at least one hit from the
0259     // current truth particle
0260     std::size_t fake_count = 0;
0261     auto itf = fake_counter.find(pid);
0262     if (itf != fake_counter.end()) {
0263       fake_count = itf->second;
0264     }
0265
0266     m_data->m_eff_plot_tool.fill(m_data->m_eff_plot_cache, ptc, is_matched);
0267     m_data->m_duplication_plot_tool.fill(m_data->m_duplication_plot_cache, ptc,
0268                                          n_matched_seeds_for_particle - 1);
0269     m_data->m_fake_tracks_plot_tool.fill(m_data->m_fake_tracks_plot_cache, ptc,
0270                                          fake_count);
0271   }
0272
0273   TRACCC_INFO("Total number of truth particles was " << total_truth_particles);
0274   TRACCC_INFO("Total number of found tracks was " << n_tracks);
0275   TRACCC_INFO("Total number of track-matched particles was "
0276               << total_matched_truth_particles);
0277   TRACCC_INFO("Total number of duplicated tracks was "
0278               << total_duplicate_tracks);
0279   TRACCC_INFO("Total number of fake tracks was " << total_fake_tracks);
0280   TRACCC_INFO("Total track efficiency was "
0281               << (100. * static_cast<double>(total_matched_truth_particles) /
0282                   static_cast<double>(total_truth_particles))
0283               << "%");
0284   TRACCC_INFO("Total track duplicate rate was "
0285               << (static_cast<double>(total_duplicate_tracks) /
0286                   static_cast<double>(total_matched_truth_particles)));
0287   TRACCC_INFO("Total track fake rate was "
0288               << (static_cast<double>(total_fake_tracks) /
0289                   static_cast<double>(total_truth_particles)));
0290 }
0291
0292 /// For ambiguity resolution
0293 void finding_performance_writer::write(
0294     const edm::track_container<default_algebra>::const_view& track_view,
0295     const event_data& evt_data) {
0296   // Set up the input containers.
0297   const edm::track_container<default_algebra>::const_device tracks{track_view};
0298
0299   const unsigned int n_tracks = tracks.tracks.size();
0300   for (unsigned int i = 0; i < n_tracks; i++) {
0301     if (m_cfg.require_fit &&
0302         tracks.tracks.at(i).fit_outcome() != track_fit_outcome::SUCCESS) {
0303       continue;
0304     }
0305
0306     // Fill stat plots
0307     m_data->m_stat_plot_tool.fill(m_data->m_stat_plot_cache,
0308                                   tracks.tracks.at(i));
0309   }
0310
0311   auto prep_data = prepare_data(track_view, m_cfg.require_fit);
0312   write_common(prep_data, evt_data);
0313 }
0314
0315 void finding_performance_writer::finalize() {
0316 #ifdef TRACCC_HAVE_ROOT
0317   // Open the output file.
0318   std::unique_ptr<TFile> ofile(
0319       TFile::Open(m_cfg.file_path.c_str(), m_cfg.file_mode.c_str()));
0320   if ((!ofile) || ofile->IsZombie()) {
0321     throw std::runtime_error("Could not open output file \"" + m_cfg.file_path +
0322                              "\" in mode \"" + m_cfg.file_mode + "\"");
0323   }
0324   ofile->cd();
0325 #else
0326   std::cout << "ROOT file \"" << m_cfg.file_path << "\" is NOT created"
0327             << std::endl;
0328 #endif  // TRACCC_HAVE_ROOT
0329
0330   m_data->m_eff_plot_tool.write(m_data->m_eff_plot_cache);
0331   m_data->m_duplication_plot_tool.write(m_data->m_duplication_plot_cache);
0332   m_data->m_fake_tracks_plot_tool.write(m_data->m_fake_tracks_plot_cache);
0333   m_data->m_stat_plot_tool.write(m_data->m_stat_plot_cache);
0334 }
0335
0336 }  // namespace traccc