Root/src/TrackFinderPerformanceWriter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/Root/TrackFinderPerformanceWriter.hpp"
0010
0011 #include "Acts/EventData/TrackParameters.hpp"
0012 #include "Acts/Utilities/VectorHelpers.hpp"
0013
0014 #include <cstddef>
0015 #include <map>
0016 #include <ostream>
0017 #include <stdexcept>
0018 #include <utility>
0019
0020 #include <TFile.h>
0021 #include <TTree.h>
0022 #include <TVectorFfwd.h>
0023 #include <TVectorT.h>
0024
0025 using Acts::VectorHelpers::eta;
0026 using Acts::VectorHelpers::phi;
0027
0028 namespace ActsExamples {
0029
0030 TrackFinderPerformanceWriter::TrackFinderPerformanceWriter(
0031     TrackFinderPerformanceWriter::Config cfg, Acts::Logging::Level lvl)
0032     : WriterT(cfg.inputTracks, "TrackFinderPerformanceWriter", lvl),
0033       m_cfg(std::move(cfg)),
0034       m_effPlotTool(m_cfg.effPlotToolConfig, lvl),
0035       m_fakePlotTool(m_cfg.fakePlotToolConfig, lvl),
0036       m_duplicationPlotTool(m_cfg.duplicationPlotToolConfig, lvl),
0037       m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig, lvl),
0038       m_trackQualityPlotTool(m_cfg.trackQualityPlotToolConfig, lvl) {
0039   // tracks collection name is already checked by base ctor
0040   if (m_cfg.inputParticles.empty()) {
0041     throw std::invalid_argument("Missing particles input collection");
0042   }
0043   if (m_cfg.inputTrackParticleMatching.empty()) {
0044     throw std::invalid_argument("Missing input track particles matching");
0045   }
0046   if (m_cfg.inputParticleTrackMatching.empty()) {
0047     throw std::invalid_argument("Missing input particle track matching");
0048   }
0049   if (m_cfg.inputParticleMeasurementsMap.empty()) {
0050     throw std::invalid_argument("Missing input measurement particles map");
0051   }
0052   if (m_cfg.filePath.empty()) {
0053     throw std::invalid_argument("Missing output filename");
0054   }
0055
0056   m_inputParticles.initialize(m_cfg.inputParticles);
0057   m_inputTrackParticleMatching.initialize(m_cfg.inputTrackParticleMatching);
0058   m_inputParticleTrackMatching.initialize(m_cfg.inputParticleTrackMatching);
0059   m_inputParticleMeasurementsMap.initialize(m_cfg.inputParticleMeasurementsMap);
0060
0061   // the output file can not be given externally since TFile accesses to the
0062   // same file from multiple threads are unsafe.
0063   // must always be opened internally
0064   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0065   if (m_outputFile == nullptr) {
0066     throw std::invalid_argument("Could not open '" + m_cfg.filePath + "'");
0067   }
0068
0069   if (m_cfg.writeMatchingDetails) {
0070     m_matchingTree = new TTree("matchingdetails", "matchingdetails");
0071
0072     m_matchingTree->Branch("event_nr", &m_treeEventNr);
0073     m_matchingTree->Branch("particle_id", &m_treeParticleId);
0074     m_matchingTree->Branch("matched", &m_treeIsMatched);
0075   }
0076
0077   // initialize the plot tools
0078   m_effPlotTool.book(m_effPlotCache);
0079   m_fakePlotTool.book(m_fakePlotCache);
0080   m_duplicationPlotTool.book(m_duplicationPlotCache);
0081   m_trackSummaryPlotTool.book(m_trackSummaryPlotCache);
0082   for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0083     m_trackSummaryPlotTool.book(m_subDetectorSummaryCaches[key], key);
0084   }
0085   m_trackQualityPlotTool.book(m_trackQualityPlotCache);
0086 }
0087
0088 TrackFinderPerformanceWriter::~TrackFinderPerformanceWriter() {
0089   m_effPlotTool.clear(m_effPlotCache);
0090   m_fakePlotTool.clear(m_fakePlotCache);
0091   m_duplicationPlotTool.clear(m_duplicationPlotCache);
0092   m_trackSummaryPlotTool.clear(m_trackSummaryPlotCache);
0093   for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0094     m_trackSummaryPlotTool.clear(m_subDetectorSummaryCaches.at(key));
0095   }
0096   m_trackQualityPlotTool.clear(m_trackQualityPlotCache);
0097   if (m_outputFile != nullptr) {
0098     m_outputFile->Close();
0099   }
0100 }
0101
0102 ProcessCode TrackFinderPerformanceWriter::finalize() {
0103   float eff_tracks = static_cast<float>(m_nTotalMatchedTracks) / m_nTotalTracks;
0104   float fakeRatio_tracks =
0105       static_cast<float>(m_nTotalFakeTracks) / m_nTotalTracks;
0106   float duplicationRatio_tracks =
0107       static_cast<float>(m_nTotalDuplicateTracks) / m_nTotalTracks;
0108
0109   float eff_particle =
0110       static_cast<float>(m_nTotalMatchedParticles) / m_nTotalParticles;
0111   float fakeRatio_particle =
0112       static_cast<float>(m_nTotalFakeParticles) / m_nTotalParticles;
0113   float duplicationRatio_particle =
0114       static_cast<float>(m_nTotalDuplicateParticles) / m_nTotalParticles;
0115
0116   ACTS_DEBUG("nTotalTracks                = " << m_nTotalTracks);
0117   ACTS_DEBUG("nTotalMatchedTracks         = " << m_nTotalMatchedTracks);
0118   ACTS_DEBUG("nTotalDuplicateTracks       = " << m_nTotalDuplicateTracks);
0119   ACTS_DEBUG("nTotalFakeTracks            = " << m_nTotalFakeTracks);
0120
0121   ACTS_INFO(
0122       "Efficiency with tracks (nMatchedTracks/ nAllTracks) = " << eff_tracks);
0123   ACTS_INFO(
0124       "Fake ratio with tracks (nFakeTracks/nAllTracks) = " << fakeRatio_tracks);
0125   ACTS_INFO("Duplicate ratio with tracks (nDuplicateTracks/nAllTracks) = "
0126             << duplicationRatio_tracks);
0127   ACTS_INFO("Efficiency with particles (nMatchedParticles/nTrueParticles) = "
0128             << eff_particle);
0129   ACTS_INFO("Fake ratio with particles (nFakeParticles/nTrueParticles) = "
0130             << fakeRatio_particle);
0131   ACTS_INFO(
0132       "Duplicate ratio with particles (nDuplicateParticles/nTrueParticles) = "
0133       << duplicationRatio_particle);
0134
0135   auto writeFloat = [&](float f, const char* name) {
0136     TVectorF v(1);
0137     v[0] = f;
0138     m_outputFile->WriteObject(&v, name);
0139   };
0140
0141   if (m_outputFile != nullptr) {
0142     m_outputFile->cd();
0143     m_effPlotTool.write(m_effPlotCache);
0144     m_fakePlotTool.write(m_fakePlotCache);
0145     m_duplicationPlotTool.write(m_duplicationPlotCache);
0146     m_trackSummaryPlotTool.write(m_trackSummaryPlotCache);
0147     for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0148       m_trackSummaryPlotTool.write(m_subDetectorSummaryCaches.at(key));
0149     }
0150     m_trackQualityPlotTool.write(m_trackQualityPlotCache);
0151
0152     writeFloat(eff_tracks, "eff_tracks");
0153     writeFloat(fakeRatio_tracks, "fakeratio_tracks");
0154     writeFloat(duplicationRatio_tracks, "duplicateratio_tracks");
0155     writeFloat(eff_particle, "eff_particles");
0156     writeFloat(fakeRatio_particle, "fakeratio_particles");
0157     writeFloat(duplicationRatio_particle, "duplicateratio_particles");
0158
0159     if (m_matchingTree != nullptr) {
0160       m_matchingTree->Write();
0161     }
0162
0163     ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
0164   }
0165   return ProcessCode::SUCCESS;
0166 }
0167
0168 ProcessCode TrackFinderPerformanceWriter::writeT(
0169     const AlgorithmContext& ctx, const ConstTrackContainer& tracks) {
0170   // The number of majority particle hits and fitted track parameters
0171   using Acts::VectorHelpers::perp;
0172
0173   // Read truth input collections
0174   const auto& particles = m_inputParticles(ctx);
0175   const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx);
0176   const auto& particleTrackMatching = m_inputParticleTrackMatching(ctx);
0177   const auto& particleMeasurementsMap = m_inputParticleMeasurementsMap(ctx);
0178
0179   // Exclusive access to the tree while writing
0180   std::lock_guard<std::mutex> lock(m_writeMutex);
0181
0182   // Vector of input features for neural network classification
0183   std::vector<float> inputFeatures(3);
0184
0185   ACTS_DEBUG("Collect information from " << tracks.size() << " tracks");
0186   std::size_t unmatched = 0, missingRefSurface = 0;
0187   for (const auto& track : tracks) {
0188     // Counting number of total trajectories
0189     m_nTotalTracks++;
0190
0191     // Check if the reco track has fitted track parameters
0192     if (!track.hasReferenceSurface()) {
0193       ACTS_VERBOSE("No fitted track parameters for track, index = "
0194                    << track.index() << " tip index = " << track.tipIndex());
0195       missingRefSurface++;
0196       continue;
0197     }
0198
0199     Acts::BoundTrackParameters fittedParameters =
0200         track.createParametersAtReference();
0201
0202     // Fill the trajectory summary info
0203     m_trackSummaryPlotTool.fill(m_trackSummaryPlotCache, fittedParameters,
0204                                 track.nTrackStates(), track.nMeasurements(),
0205                                 track.nOutliers(), track.nHoles(),
0206                                 track.nSharedHits());
0207
0208     // Potentially fill other track summary caches for the given volumes
0209     for (const auto& [key, volumes] : m_cfg.subDetectorTrackSummaryVolumes) {
0210       ACTS_VERBOSE("Fill track summary stats for subset " << key);
0211       std::size_t nTrackStates{}, nMeasurements{}, nOutliers{}, nHoles{},
0212           nSharedHits{};
0213       for (auto state : track.trackStatesReversed()) {
0214         if (!state.hasReferenceSurface() ||
0215             !volumes.contains(state.referenceSurface().geometryId().volume())) {
0216           continue;
0217         }
0218
0219         nTrackStates++;
0220         nMeasurements += static_cast<std::size_t>(
0221             state.typeFlags().test(Acts::MeasurementFlag));
0222         nOutliers +=
0223             static_cast<std::size_t>(state.typeFlags().test(Acts::OutlierFlag));
0224         nHoles +=
0225             static_cast<std::size_t>(state.typeFlags().test(Acts::HoleFlag));
0226         nSharedHits += static_cast<std::size_t>(
0227             state.typeFlags().test(Acts::SharedHitFlag));
0228       }
0229       m_trackSummaryPlotTool.fill(m_subDetectorSummaryCaches.at(key),
0230                                   fittedParameters, nTrackStates, nMeasurements,
0231                                   nOutliers, nHoles, nSharedHits);
0232     }
0233
0234     // Get the truth matching information
0235     auto imatched = trackParticleMatching.find(track.index());
0236     if (imatched == trackParticleMatching.end()) {
0237       ACTS_DEBUG("No truth matching information for this track, index = "
0238                  << track.index() << " tip index = " << track.tipIndex());
0239       unmatched++;
0240       continue;
0241     }
0242
0243     const auto& particleMatch = imatched->second;
0244
0245     if (particleMatch.classification == TrackMatchClassification::Fake) {
0246       m_nTotalFakeTracks++;
0247     }
0248
0249     if (particleMatch.classification == TrackMatchClassification::Duplicate) {
0250       m_nTotalDuplicateTracks++;
0251     }
0252
0253     // Fill fake ratio plots
0254     m_fakePlotTool.fill(
0255         m_fakePlotCache, fittedParameters,
0256         particleMatch.classification == TrackMatchClassification::Fake);
0257
0258     // Fill the duplication ratio
0259     m_duplicationPlotTool.fill(
0260         m_duplicationPlotCache, fittedParameters,
0261         particleMatch.classification == TrackMatchClassification::Duplicate);
0262
0263     if (particleMatch.particle.has_value() &&
0264         particleMeasurementsMap.contains(particleMatch.particle.value())) {
0265       const auto measurements =
0266           particleMeasurementsMap.equal_range(particleMatch.particle.value());
0267
0268       std::size_t nTrackMeasurements =
0269           track.nMeasurements() + track.nOutliers();
0270       std::size_t nMatchedHits =
0271           particleMatch.contributingParticles.front().hitCount;
0272       std::size_t nParticleHits =
0273           std::distance(measurements.first, measurements.second);
0274
0275       double completeness = static_cast<double>(nMatchedHits) / nParticleHits;
0276       double purity = static_cast<double>(nMatchedHits) / nTrackMeasurements;
0277
0278       // Fill the track quality plots
0279       m_trackQualityPlotTool.fill(m_trackQualityPlotCache, fittedParameters,
0280                                   completeness, purity);
0281     }
0282   }
0283
0284   if (unmatched > 0) {
0285     ACTS_DEBUG("No matching information found for " << unmatched << " tracks");
0286   }
0287   if (missingRefSurface > 0) {
0288     ACTS_DEBUG("Reference surface was missing for " << missingRefSurface
0289                                                     << " tracks");
0290   }
0291
0292   // Loop over all truth particles for efficiency plots and reco details.
0293   for (const auto& particle : particles) {
0294     auto particleId = particle.particleId();
0295
0296     // Investigate the truth-matched tracks
0297     std::size_t nMatchedTracks = 0;
0298     std::size_t nFakeTracks = 0;
0299     bool isReconstructed = false;
0300     if (auto imatched = particleTrackMatching.find(particleId);
0301         imatched != particleTrackMatching.end()) {
0302       nMatchedTracks = (imatched->second.track.has_value() ? 1 : 0) +
0303                        imatched->second.duplicates;
0304
0305       // Add number for total matched tracks here
0306       m_nTotalMatchedTracks += nMatchedTracks;
0307       m_nTotalMatchedParticles += 1;
0308
0309       // Check if the particle has more than one matched track for the duplicate
0310       // rate/ratio
0311       if (nMatchedTracks > 1) {
0312         m_nTotalDuplicateParticles += 1;
0313       }
0314       isReconstructed = imatched->second.track.has_value();
0315
0316       nFakeTracks = imatched->second.fakes;
0317       if (nFakeTracks > 0) {
0318         m_nTotalFakeParticles += 1;
0319       }
0320     }
0321
0322     // Loop over all the other truth particle and find the distance to the
0323     // closest one
0324     double minDeltaR = -1;
0325     for (const auto& closeParticle : particles) {
0326       if (closeParticle.particleId() == particleId) {
0327         continue;
0328       }
0329       double p_phi = phi(particle.direction());
0330       double p_eta = eta(particle.direction());
0331       double c_phi = phi(closeParticle.direction());
0332       double c_eta = eta(closeParticle.direction());
0333       double distance = sqrt(pow(p_phi - c_phi, 2) + pow(p_eta - c_eta, 2));
0334       if (minDeltaR == -1 || distance < minDeltaR) {
0335         minDeltaR = distance;
0336       }
0337     }
0338
0339     // Fill efficiency plots
0340     m_effPlotTool.fill(m_effPlotCache, particle.initial(), minDeltaR,
0341                        isReconstructed);
0342     // Fill number of duplicated tracks for this particle
0343     m_duplicationPlotTool.fill(m_duplicationPlotCache, particle.initial(),
0344                                nMatchedTracks - 1);
0345
0346     // Fill number of reconstructed/truth-matched/fake tracks for this particle
0347     m_fakePlotTool.fill(m_fakePlotCache, particle.initial(), nMatchedTracks,
0348                         nFakeTracks);
0349
0350     m_nTotalParticles += 1;
0351   }
0352
0353   // Write additional stuff to TTree
0354   if (m_cfg.writeMatchingDetails && m_matchingTree != nullptr) {
0355     for (const auto& particle : particles) {
0356       auto particleId = particle.particleId();
0357
0358       m_treeEventNr = ctx.eventNumber;
0359       m_treeParticleId = particleId.value();
0360
0361       m_treeIsMatched = false;
0362       if (auto imatched = particleTrackMatching.find(particleId);
0363           imatched != particleTrackMatching.end()) {
0364         m_treeIsMatched = imatched->second.track.has_value();
0365       }
0366
0367       m_matchingTree->Fill();
0368     }
0369   }
0370
0371   return ProcessCode::SUCCESS;
0372 }
0373
0374 }  // namespace ActsExamples