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 += static_cast<std::size_t>(
0228             state.typeFlags().test(Acts::MeasurementFlag));
0229         nOutliers +=
0230             static_cast<std::size_t>(state.typeFlags().test(Acts::OutlierFlag));
0231         nHoles +=
0232             static_cast<std::size_t>(state.typeFlags().test(Acts::HoleFlag));
0233         nSharedHits += static_cast<std::size_t>(
0234             state.typeFlags().test(Acts::SharedHitFlag));
0235       }
0236       m_trackSummaryPlotTool.fill(m_subDetectorSummaryCaches.at(key),
0237                                   fittedParameters, nTrackStates, nMeasurements,
0238                                   nOutliers, nHoles, nSharedHits);
0239     }
0240
0241     // Get the truth matching information
0242     auto imatched = trackParticleMatching.find(track.index());
0243     if (imatched == trackParticleMatching.end()) {
0244       ACTS_DEBUG("No truth matching information for this track, index = "
0245                  << track.index() << " tip index = " << track.tipIndex());
0246       unmatched++;
0247       continue;
0248     }
0249
0250     const auto& particleMatch = imatched->second;
0251
0252     if (particleMatch.classification == TrackMatchClassification::Fake) {
0253       m_nTotalFakeTracks++;
0254     }
0255
0256     if (particleMatch.classification == TrackMatchClassification::Duplicate) {
0257       m_nTotalDuplicateTracks++;
0258     }
0259
0260     // Fill fake ratio plots
0261     m_fakePlotTool.fill(
0262         m_fakePlotCache, fittedParameters,
0263         particleMatch.classification == TrackMatchClassification::Fake);
0264
0265     // Fill the duplication ratio
0266     m_duplicationPlotTool.fill(
0267         m_duplicationPlotCache, fittedParameters,
0268         particleMatch.classification == TrackMatchClassification::Duplicate);
0269
0270     if (particleMatch.particle.has_value() &&
0271         particleMeasurementsMap.contains(particleMatch.particle.value())) {
0272       const auto measurements =
0273           particleMeasurementsMap.equal_range(particleMatch.particle.value());
0274
0275       std::size_t nTrackMeasurements =
0276           track.nMeasurements() + track.nOutliers();
0277       std::size_t nMatchedHits =
0278           particleMatch.contributingParticles.front().hitCount;
0279       std::size_t nParticleHits =
0280           std::distance(measurements.first, measurements.second);
0281
0282       double completeness = static_cast<double>(nMatchedHits) / nParticleHits;
0283       double purity = static_cast<double>(nMatchedHits) / nTrackMeasurements;
0284
0285       // Fill the track quality plots
0286       m_trackQualityPlotTool.fill(m_trackQualityPlotCache, fittedParameters,
0287                                   completeness, purity);
0288     }
0289   }
0290
0291   if (unmatched > 0) {
0292     ACTS_DEBUG("No matching information found for " << unmatched << " tracks");
0293   }
0294   if (missingRefSurface > 0) {
0295     ACTS_DEBUG("Reference surface was missing for " << missingRefSurface
0296                                                     << " tracks");
0297   }
0298
0299   // Loop over all truth particles for efficiency plots and reco details.
0300   for (const auto& particle : particles) {
0301     auto particleId = particle.particleId();
0302
0303     // Investigate the truth-matched tracks
0304     std::size_t nMatchedTracks = 0;
0305     std::size_t nFakeTracks = 0;
0306     bool isReconstructed = false;
0307     if (auto imatched = particleTrackMatching.find(particleId);
0308         imatched != particleTrackMatching.end()) {
0309       nMatchedTracks = (imatched->second.track.has_value() ? 1 : 0) +
0310                        imatched->second.duplicates;
0311
0312       // Add number for total matched tracks here
0313       m_nTotalMatchedTracks += nMatchedTracks;
0314       m_nTotalMatchedParticles += imatched->second.track.has_value() ? 1 : 0;
0315
0316       // Check if the particle has more than one matched track for the duplicate
0317       // rate/ratio
0318       if (nMatchedTracks > 1) {
0319         m_nTotalDuplicateParticles += 1;
0320       }
0321       isReconstructed = imatched->second.track.has_value();
0322
0323       nFakeTracks = imatched->second.fakes;
0324       if (nFakeTracks > 0) {
0325         m_nTotalFakeParticles += 1;
0326       }
0327     }
0328
0329     // Loop over all the other truth particle and find the distance to the
0330     // closest one
0331     double minDeltaR = -1;
0332     for (const auto& closeParticle : particles) {
0333       if (closeParticle.particleId() == particleId) {
0334         continue;
0335       }
0336       double distance = Acts::VectorHelpers::deltaR(particle.direction(),
0337                                                     closeParticle.direction());
0338       if (minDeltaR == -1 || distance < minDeltaR) {
0339         minDeltaR = distance;
0340       }
0341     }
0342
0343     // Fill efficiency plots
0344     m_effPlotTool.fill(ctx.geoContext, m_effPlotCache, particle.initialState(),
0345                        minDeltaR, isReconstructed);
0346     // Fill number of duplicated tracks for this particle
0347     m_duplicationPlotTool.fill(m_duplicationPlotCache, particle.initialState(),
0348                                nMatchedTracks - 1);
0349
0350     // Fill number of reconstructed/truth-matched/fake tracks for this particle
0351     m_fakePlotTool.fill(m_fakePlotCache, particle.initialState(),
0352                         nMatchedTracks, nFakeTracks);
0353
0354     m_nTotalParticles += 1;
0355   }
0356
0357   // Write additional stuff to TTree
0358   if (m_cfg.writeMatchingDetails && m_matchingTree != nullptr) {
0359     for (const auto& particle : particles) {
0360       auto particleId = particle.particleId();
0361
0362       m_treeEventNr = ctx.eventNumber;
0363       m_treeParticleVertexPrimary = particleId.vertexPrimary();
0364       m_treeParticleVertexSecondary = particleId.vertexSecondary();
0365       m_treeParticleParticle = particleId.particle();
0366       m_treeParticleGeneration = particleId.generation();
0367       m_treeParticleSubParticle = particleId.subParticle();
0368
0369       m_treeIsMatched = false;
0370       if (auto imatched = particleTrackMatching.find(particleId);
0371           imatched != particleTrackMatching.end()) {
0372         m_treeIsMatched = imatched->second.track.has_value();
0373       }
0374
0375       m_matchingTree->Fill();
0376     }
0377   }
0378
0379   return ProcessCode::SUCCESS;
0380 }
0381
0382 }  // namespace ActsExamples