Root/src/RootTrackFinderPerformanceWriter.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/RootTrackFinderPerformanceWriter.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 RootTrackFinderPerformanceWriter::RootTrackFinderPerformanceWriter(
0031     RootTrackFinderPerformanceWriter::Config cfg, Acts::Logging::Level lvl)
0032     : WriterT(cfg.inputTracks, "RootTrackFinderPerformanceWriter", 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_vertex_primary",
0074                            &m_treeParticleVertexPrimary);
0075     m_matchingTree->Branch("particle_id_vertex_secondary",
0076                            &m_treeParticleVertexSecondary);
0077     m_matchingTree->Branch("particle_id_particle", &m_treeParticleParticle);
0078     m_matchingTree->Branch("particle_id_generation", &m_treeParticleGeneration);
0079     m_matchingTree->Branch("particle_id_sub_particle",
0080                            &m_treeParticleSubParticle);
0081     m_matchingTree->Branch("matched", &m_treeIsMatched);
0082   }
0083
0084   // initialize the plot tools
0085   m_effPlotTool.book(m_effPlotCache);
0086   m_fakePlotTool.book(m_fakePlotCache);
0087   m_duplicationPlotTool.book(m_duplicationPlotCache);
0088   m_trackSummaryPlotTool.book(m_trackSummaryPlotCache);
0089   for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0090     m_trackSummaryPlotTool.book(m_subDetectorSummaryCaches[key], key);
0091   }
0092   m_trackQualityPlotTool.book(m_trackQualityPlotCache);
0093 }
0094
0095 RootTrackFinderPerformanceWriter::~RootTrackFinderPerformanceWriter() {
0096   m_effPlotTool.clear(m_effPlotCache);
0097   m_fakePlotTool.clear(m_fakePlotCache);
0098   m_duplicationPlotTool.clear(m_duplicationPlotCache);
0099   m_trackSummaryPlotTool.clear(m_trackSummaryPlotCache);
0100   for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0101     m_trackSummaryPlotTool.clear(m_subDetectorSummaryCaches.at(key));
0102   }
0103   m_trackQualityPlotTool.clear(m_trackQualityPlotCache);
0104   if (m_outputFile != nullptr) {
0105     m_outputFile->Close();
0106   }
0107 }
0108
0109 ProcessCode RootTrackFinderPerformanceWriter::finalize() {
0110   float eff_tracks = static_cast<float>(m_nTotalMatchedTracks) / m_nTotalTracks;
0111   float fakeRatio_tracks =
0112       static_cast<float>(m_nTotalFakeTracks) / m_nTotalTracks;
0113   float duplicationRatio_tracks =
0114       static_cast<float>(m_nTotalDuplicateTracks) / m_nTotalTracks;
0115
0116   float eff_particle =
0117       static_cast<float>(m_nTotalMatchedParticles) / m_nTotalParticles;
0118   float fakeRatio_particle =
0119       static_cast<float>(m_nTotalFakeParticles) / m_nTotalParticles;
0120   float duplicationRatio_particle =
0121       static_cast<float>(m_nTotalDuplicateParticles) / m_nTotalParticles;
0122
0123   ACTS_DEBUG("nTotalTracks                = " << m_nTotalTracks);
0124   ACTS_DEBUG("nTotalMatchedTracks         = " << m_nTotalMatchedTracks);
0125   ACTS_DEBUG("nTotalDuplicateTracks       = " << m_nTotalDuplicateTracks);
0126   ACTS_DEBUG("nTotalFakeTracks            = " << m_nTotalFakeTracks);
0127
0128   ACTS_INFO(
0129       "Efficiency with tracks (nMatchedTracks/ nAllTracks) = " << eff_tracks);
0130   ACTS_INFO(
0131       "Fake ratio with tracks (nFakeTracks/nAllTracks) = " << fakeRatio_tracks);
0132   ACTS_INFO("Duplicate ratio with tracks (nDuplicateTracks/nAllTracks) = "
0133             << duplicationRatio_tracks);
0134   ACTS_INFO("Efficiency with particles (nMatchedParticles/nTrueParticles) = "
0135             << eff_particle);
0136   ACTS_INFO("Fake ratio with particles (nFakeParticles/nTrueParticles) = "
0137             << fakeRatio_particle);
0138   ACTS_INFO(
0139       "Duplicate ratio with particles (nDuplicateParticles/nTrueParticles) = "
0140       << duplicationRatio_particle);
0141
0142   auto writeFloat = [&](float f, const char* name) {
0143     TVectorF v(1);
0144     v[0] = f;
0145     m_outputFile->WriteObject(&v, name);
0146   };
0147
0148   if (m_outputFile != nullptr) {
0149     m_outputFile->cd();
0150     m_effPlotTool.write(m_effPlotCache);
0151     m_fakePlotTool.write(m_fakePlotCache);
0152     m_duplicationPlotTool.write(m_duplicationPlotCache);
0153     m_trackSummaryPlotTool.write(m_trackSummaryPlotCache);
0154     for (const auto& [key, _] : m_cfg.subDetectorTrackSummaryVolumes) {
0155       m_trackSummaryPlotTool.write(m_subDetectorSummaryCaches.at(key));
0156     }
0157     m_trackQualityPlotTool.write(m_trackQualityPlotCache);
0158
0159     writeFloat(eff_tracks, "eff_tracks");
0160     writeFloat(fakeRatio_tracks, "fakeratio_tracks");
0161     writeFloat(duplicationRatio_tracks, "duplicateratio_tracks");
0162     writeFloat(eff_particle, "eff_particles");
0163     writeFloat(fakeRatio_particle, "fakeratio_particles");
0164     writeFloat(duplicationRatio_particle, "duplicateratio_particles");
0165
0166     if (m_matchingTree != nullptr) {
0167       m_matchingTree->Write();
0168     }
0169
0170     ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
0171   }
0172   return ProcessCode::SUCCESS;
0173 }
0174
0175 ProcessCode RootTrackFinderPerformanceWriter::writeT(
0176     const AlgorithmContext& ctx, const ConstTrackContainer& tracks) {
0177   // The number of majority particle hits and fitted track parameters
0178   using Acts::VectorHelpers::perp;
0179
0180   // Read truth input collections
0181   const auto& particles = m_inputParticles(ctx);
0182   const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx);
0183   const auto& particleTrackMatching = m_inputParticleTrackMatching(ctx);
0184   const auto& particleMeasurementsMap = m_inputParticleMeasurementsMap(ctx);
0185
0186   // Exclusive access to the tree while writing
0187   std::lock_guard<std::mutex> lock(m_writeMutex);
0188
0189   // Vector of input features for neural network classification
0190   std::vector<float> inputFeatures(3);
0191
0192   ACTS_DEBUG("Collect information from " << tracks.size() << " tracks");
0193   std::size_t unmatched = 0, missingRefSurface = 0;
0194   for (const auto& track : tracks) {
0195     // Counting number of total trajectories
0196     m_nTotalTracks++;
0197
0198     // Check if the reco track has fitted track parameters
0199     if (!track.hasReferenceSurface()) {
0200       ACTS_VERBOSE("No fitted track parameters for track, index = "
0201                    << track.index() << " tip index = " << track.tipIndex());
0202       missingRefSurface++;
0203       continue;
0204     }
0205
0206     Acts::BoundTrackParameters fittedParameters =
0207         track.createParametersAtReference();
0208
0209     // Fill the trajectory summary info
0210     m_trackSummaryPlotTool.fill(m_trackSummaryPlotCache, fittedParameters,
0211                                 track.nTrackStates(), track.nMeasurements(),
0212                                 track.nOutliers(), track.nHoles(),
0213                                 track.nSharedHits());
0214
0215     // Potentially fill other track summary caches for the given volumes
0216     for (const auto& [key, volumes] : m_cfg.subDetectorTrackSummaryVolumes) {
0217       ACTS_VERBOSE("Fill track summary stats for subset " << key);
0218       std::size_t nTrackStates{}, nMeasurements{}, nOutliers{}, nHoles{},
0219           nSharedHits{};
0220       for (auto state : track.trackStatesReversed()) {
0221         if (!state.hasReferenceSurface() ||
0222             !volumes.contains(state.referenceSurface().geometryId().volume())) {
0223           continue;
0224         }
0225
0226         nTrackStates++;
0227         nMeasurements +=
0228             static_cast<std::size_t>(state.typeFlags().isMeasurement());
0229         nOutliers += static_cast<std::size_t>(state.typeFlags().isOutlier());
0230         nHoles += static_cast<std::size_t>(state.typeFlags().isHole());
0231         nSharedHits +=
0232             static_cast<std::size_t>(state.typeFlags().isSharedHit());
0233       }
0234       m_trackSummaryPlotTool.fill(m_subDetectorSummaryCaches.at(key),
0235                                   fittedParameters, nTrackStates, nMeasurements,
0236                                   nOutliers, nHoles, nSharedHits);
0237     }
0238
0239     // Get the truth matching information
0240     auto imatched = trackParticleMatching.find(track.index());
0241     if (imatched == trackParticleMatching.end()) {
0242       ACTS_DEBUG("No truth matching information for this track, index = "
0243                  << track.index() << " tip index = " << track.tipIndex());
0244       unmatched++;
0245       continue;
0246     }
0247
0248     const auto& particleMatch = imatched->second;
0249
0250     if (particleMatch.classification == TrackMatchClassification::Fake) {
0251       m_nTotalFakeTracks++;
0252     }
0253
0254     if (particleMatch.classification == TrackMatchClassification::Duplicate) {
0255       m_nTotalDuplicateTracks++;
0256     }
0257
0258     // Fill fake ratio plots
0259     m_fakePlotTool.fill(
0260         m_fakePlotCache, fittedParameters,
0261         particleMatch.classification == TrackMatchClassification::Fake);
0262
0263     // Fill the duplication ratio
0264     m_duplicationPlotTool.fill(
0265         m_duplicationPlotCache, fittedParameters,
0266         particleMatch.classification == TrackMatchClassification::Duplicate);
0267
0268     if (particleMatch.particle.has_value() &&
0269         particleMeasurementsMap.contains(particleMatch.particle.value())) {
0270       const auto measurements =
0271           particleMeasurementsMap.equal_range(particleMatch.particle.value());
0272
0273       std::size_t nTrackMeasurements =
0274           track.nMeasurements() + track.nOutliers();
0275       std::size_t nMatchedHits =
0276           particleMatch.contributingParticles.front().hitCount;
0277       std::size_t nParticleHits =
0278           std::distance(measurements.first, measurements.second);
0279
0280       double completeness = static_cast<double>(nMatchedHits) / nParticleHits;
0281       double purity = static_cast<double>(nMatchedHits) / nTrackMeasurements;
0282
0283       // Fill the track quality plots
0284       m_trackQualityPlotTool.fill(m_trackQualityPlotCache, fittedParameters,
0285                                   completeness, purity);
0286     }
0287   }
0288
0289   if (unmatched > 0) {
0290     ACTS_DEBUG("No matching information found for " << unmatched << " tracks");
0291   }
0292   if (missingRefSurface > 0) {
0293     ACTS_DEBUG("Reference surface was missing for " << missingRefSurface
0294                                                     << " tracks");
0295   }
0296
0297   // Loop over all truth particles for efficiency plots and reco details.
0298   for (const auto& particle : particles) {
0299     auto particleId = particle.particleId();
0300
0301     // Investigate the truth-matched tracks
0302     std::size_t nMatchedTracks = 0;
0303     std::size_t nFakeTracks = 0;
0304     bool isReconstructed = false;
0305     if (auto imatched = particleTrackMatching.find(particleId);
0306         imatched != particleTrackMatching.end()) {
0307       isReconstructed = imatched->second.track.has_value();
0308       nMatchedTracks = (isReconstructed ? 1 : 0) + imatched->second.duplicates;
0309
0310       // Add number for total matched tracks here
0311       m_nTotalMatchedTracks += nMatchedTracks;
0312       m_nTotalMatchedParticles += isReconstructed ? 1 : 0;
0313
0314       // Check if the particle has more than one matched track for the duplicate
0315       // rate/ratio
0316       if (nMatchedTracks > 1) {
0317         m_nTotalDuplicateParticles += 1;
0318       }
0319
0320       nFakeTracks = imatched->second.fakes;
0321       if (nFakeTracks > 0) {
0322         m_nTotalFakeParticles += 1;
0323       }
0324     }
0325
0326     // Loop over all the other truth particle and find the distance to the
0327     // closest one
0328     double minDeltaR = -1;
0329     for (const auto& closeParticle : particles) {
0330       if (closeParticle.particleId() == particleId) {
0331         continue;
0332       }
0333       double distance = Acts::VectorHelpers::deltaR(particle.direction(),
0334                                                     closeParticle.direction());
0335       if (minDeltaR == -1 || distance < minDeltaR) {
0336         minDeltaR = distance;
0337       }
0338     }
0339
0340     // Fill efficiency plots
0341     m_effPlotTool.fill(ctx.geoContext, m_effPlotCache, particle.initialState(),
0342                        minDeltaR, isReconstructed);
0343     // Fill duplication plots
0344     m_duplicationPlotTool.fill(m_duplicationPlotCache, particle.initialState(),
0345                                nMatchedTracks);
0346
0347     // Fill number of reconstructed/truth-matched/fake tracks for this particle
0348     m_fakePlotTool.fill(m_fakePlotCache, particle.initialState(),
0349                         nMatchedTracks, nFakeTracks);
0350
0351     m_nTotalParticles += 1;
0352   }
0353
0354   // Write additional stuff to TTree
0355   if (m_cfg.writeMatchingDetails && m_matchingTree != nullptr) {
0356     for (const auto& particle : particles) {
0357       auto particleId = particle.particleId();
0358
0359       m_treeEventNr = ctx.eventNumber;
0360       m_treeParticleVertexPrimary = particleId.vertexPrimary();
0361       m_treeParticleVertexSecondary = particleId.vertexSecondary();
0362       m_treeParticleParticle = particleId.particle();
0363       m_treeParticleGeneration = particleId.generation();
0364       m_treeParticleSubParticle = particleId.subParticle();
0365
0366       m_treeIsMatched = false;
0367       if (auto imatched = particleTrackMatching.find(particleId);
0368           imatched != particleTrackMatching.end()) {
0369         m_treeIsMatched = imatched->second.track.has_value();
0370       }
0371
0372       m_matchingTree->Fill();
0373     }
0374   }
0375
0376   return ProcessCode::SUCCESS;
0377 }
0378
0379 }  // namespace ActsExamples