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_fakeRatePlotTool(m_cfg.fakeRatePlotToolConfig, lvl),
0036       m_duplicationPlotTool(m_cfg.duplicationPlotToolConfig, lvl),
0037       m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig, lvl) {
0038   // tracks collection name is already checked by base ctor
0039   if (m_cfg.inputParticles.empty()) {
0040     throw std::invalid_argument("Missing particles input collection");
0041   }
0042   if (m_cfg.inputTrackParticleMatching.empty()) {
0043     throw std::invalid_argument("Missing input track particles matching");
0044   }
0045   if (m_cfg.inputParticleTrackMatching.empty()) {
0046     throw std::invalid_argument("Missing input particle track matching");
0047   }
0048   if (m_cfg.filePath.empty()) {
0049     throw std::invalid_argument("Missing output filename");
0050   }
0051
0052   m_inputParticles.initialize(m_cfg.inputParticles);
0053   m_inputTrackParticleMatching.initialize(m_cfg.inputTrackParticleMatching);
0054   m_inputParticleTrackMatching.initialize(m_cfg.inputParticleTrackMatching);
0055
0056   // the output file can not be given externally since TFile accesses to the
0057   // same file from multiple threads are unsafe.
0058   // must always be opened internally
0059   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0060   if (m_outputFile == nullptr) {
0061     throw std::invalid_argument("Could not open '" + m_cfg.filePath + "'");
0062   }
0063
0064   if (m_cfg.writeMatchingDetails) {
0065     m_matchingTree = new TTree("matchingdetails", "matchingdetails");
0066
0067     m_matchingTree->Branch("event_nr", &m_treeEventNr);
0068     m_matchingTree->Branch("particle_id", &m_treeParticleId);
0069     m_matchingTree->Branch("matched", &m_treeIsMatched);
0070   }
0071
0072   // initialize the plot tools
0073   m_effPlotTool.book(m_effPlotCache);
0074   m_fakeRatePlotTool.book(m_fakeRatePlotCache);
0075   m_duplicationPlotTool.book(m_duplicationPlotCache);
0076   m_trackSummaryPlotTool.book(m_trackSummaryPlotCache);
0077 }
0078
0079 TrackFinderPerformanceWriter::~TrackFinderPerformanceWriter() {
0080   m_effPlotTool.clear(m_effPlotCache);
0081   m_fakeRatePlotTool.clear(m_fakeRatePlotCache);
0082   m_duplicationPlotTool.clear(m_duplicationPlotCache);
0083   m_trackSummaryPlotTool.clear(m_trackSummaryPlotCache);
0084   if (m_outputFile != nullptr) {
0085     m_outputFile->Close();
0086   }
0087 }
0088
0089 ProcessCode TrackFinderPerformanceWriter::finalize() {
0090   float eff_tracks = static_cast<float>(m_nTotalMatchedTracks) / m_nTotalTracks;
0091   float fakeRate_tracks =
0092       static_cast<float>(m_nTotalFakeTracks) / m_nTotalTracks;
0093   float duplicationRate_tracks =
0094       static_cast<float>(m_nTotalDuplicateTracks) / m_nTotalTracks;
0095
0096   float eff_particle =
0097       static_cast<float>(m_nTotalMatchedParticles) / m_nTotalParticles;
0098   float fakeRate_particle =
0099       static_cast<float>(m_nTotalFakeParticles) / m_nTotalParticles;
0100   float duplicationRate_particle =
0101       static_cast<float>(m_nTotalDuplicateParticles) / m_nTotalParticles;
0102
0103   ACTS_DEBUG("nTotalTracks                = " << m_nTotalTracks);
0104   ACTS_DEBUG("nTotalMatchedTracks         = " << m_nTotalMatchedTracks);
0105   ACTS_DEBUG("nTotalDuplicateTracks       = " << m_nTotalDuplicateTracks);
0106   ACTS_DEBUG("nTotalFakeTracks            = " << m_nTotalFakeTracks);
0107
0108   ACTS_INFO(
0109       "Efficiency with tracks (nMatchedTracks/ nAllTracks) = " << eff_tracks);
0110   ACTS_INFO(
0111       "Fake rate with tracks (nFakeTracks/nAllTracks) = " << fakeRate_tracks);
0112   ACTS_INFO("Duplicate rate with tracks (nDuplicateTracks/nAllTracks) = "
0113             << duplicationRate_tracks);
0114   ACTS_INFO("Efficiency with particles (nMatchedParticles/nTrueParticles) = "
0115             << eff_particle);
0116   ACTS_INFO("Fake rate with particles (nFakeParticles/nTrueParticles) = "
0117             << fakeRate_particle);
0118   ACTS_INFO(
0119       "Duplicate rate with particles (nDuplicateParticles/nTrueParticles) = "
0120       << duplicationRate_particle);
0121
0122   auto writeFloat = [&](float f, const char* name) {
0123     TVectorF v(1);
0124     v[0] = f;
0125     m_outputFile->WriteObject(&v, name);
0126   };
0127
0128   if (m_outputFile != nullptr) {
0129     m_outputFile->cd();
0130     m_effPlotTool.write(m_effPlotCache);
0131     m_fakeRatePlotTool.write(m_fakeRatePlotCache);
0132     m_duplicationPlotTool.write(m_duplicationPlotCache);
0133     m_trackSummaryPlotTool.write(m_trackSummaryPlotCache);
0134     writeFloat(eff_tracks, "eff_tracks");
0135     writeFloat(fakeRate_tracks, "fakerate_tracks");
0136     writeFloat(duplicationRate_tracks, "duplicaterate_tracks");
0137     writeFloat(eff_particle, "eff_particles");
0138     writeFloat(fakeRate_particle, "fakerate_particles");
0139     writeFloat(duplicationRate_particle, "duplicaterate_particles");
0140
0141     if (m_matchingTree != nullptr) {
0142       m_matchingTree->Write();
0143     }
0144
0145     ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
0146   }
0147   return ProcessCode::SUCCESS;
0148 }
0149
0150 ProcessCode TrackFinderPerformanceWriter::writeT(
0151     const AlgorithmContext& ctx, const ConstTrackContainer& tracks) {
0152   // The number of majority particle hits and fitted track parameters
0153   using Acts::VectorHelpers::perp;
0154
0155   // Read truth input collections
0156   const auto& particles = m_inputParticles(ctx);
0157   const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx);
0158   const auto& particleTrackMatching = m_inputParticleTrackMatching(ctx);
0159
0160   // Exclusive access to the tree while writing
0161   std::lock_guard<std::mutex> lock(m_writeMutex);
0162
0163   // Vector of input features for neural network classification
0164   std::vector<float> inputFeatures(3);
0165
0166   for (const auto& track : tracks) {
0167     // Counting number of total trajectories
0168     m_nTotalTracks++;
0169
0170     // Check if the reco track has fitted track parameters
0171     if (!track.hasReferenceSurface()) {
0172       ACTS_WARNING("No fitted track parameters for track, index = "
0173                    << track.index() << " tip index = " << track.tipIndex());
0174       continue;
0175     }
0176
0177     Acts::BoundTrackParameters fittedParameters =
0178         track.createParametersAtReference();
0179
0180     // Fill the trajectory summary info
0181     m_trackSummaryPlotTool.fill(m_trackSummaryPlotCache, fittedParameters,
0182                                 track.nTrackStates(), track.nMeasurements(),
0183                                 track.nOutliers(), track.nHoles(),
0184                                 track.nSharedHits());
0185
0186     // Get the truth matching information
0187     auto imatched = trackParticleMatching.find(track.index());
0188     if (imatched == trackParticleMatching.end()) {
0189       ACTS_DEBUG("No truth matching information for this track, index = "
0190                  << track.index() << " tip index = " << track.tipIndex());
0191       continue;
0192     }
0193
0194     const auto& particleMatch = imatched->second;
0195
0196     if (particleMatch.classification == TrackMatchClassification::Fake) {
0197       m_nTotalFakeTracks++;
0198     }
0199
0200     if (particleMatch.classification == TrackMatchClassification::Duplicate) {
0201       m_nTotalDuplicateTracks++;
0202     }
0203
0204     // Fill fake rate plots
0205     m_fakeRatePlotTool.fill(
0206         m_fakeRatePlotCache, fittedParameters,
0207         particleMatch.classification == TrackMatchClassification::Fake);
0208
0209     // Fill the duplication rate
0210     m_duplicationPlotTool.fill(
0211         m_duplicationPlotCache, fittedParameters,
0212         particleMatch.classification == TrackMatchClassification::Duplicate);
0213   }
0214
0215   // Loop over all truth particles for efficiency plots and reco details.
0216   for (const auto& particle : particles) {
0217     auto particleId = particle.particleId();
0218
0219     // Investigate the truth-matched tracks
0220     std::size_t nMatchedTracks = 0;
0221     std::size_t nFakeTracks = 0;
0222     bool isReconstructed = false;
0223     if (auto imatched = particleTrackMatching.find(particleId);
0224         imatched != particleTrackMatching.end()) {
0225       nMatchedTracks = (imatched->second.track.has_value() ? 1 : 0) +
0226                        imatched->second.duplicates;
0227
0228       // Add number for total matched tracks here
0229       m_nTotalMatchedTracks += nMatchedTracks;
0230       m_nTotalMatchedParticles += 1;
0231
0232       // Check if the particle has more than one matched track for the duplicate
0233       // rate
0234       if (nMatchedTracks > 1) {
0235         m_nTotalDuplicateParticles += 1;
0236       }
0237       isReconstructed = imatched->second.track.has_value();
0238
0239       nFakeTracks = imatched->second.fakes;
0240       if (nFakeTracks > 0) {
0241         m_nTotalFakeParticles += 1;
0242       }
0243     }
0244
0245     // Loop over all the other truth particle and find the distance to the
0246     // closest one
0247     double minDeltaR = -1;
0248     for (const auto& closeParticle : particles) {
0249       if (closeParticle.particleId() == particleId) {
0250         continue;
0251       }
0252       double p_phi = phi(particle.direction());
0253       double p_eta = eta(particle.direction());
0254       double c_phi = phi(closeParticle.direction());
0255       double c_eta = eta(closeParticle.direction());
0256       double distance = sqrt(pow(p_phi - c_phi, 2) + pow(p_eta - c_eta, 2));
0257       if (minDeltaR == -1 || distance < minDeltaR) {
0258         minDeltaR = distance;
0259       }
0260     }
0261
0262     // Fill efficiency plots
0263     m_effPlotTool.fill(m_effPlotCache, particle.initial(), minDeltaR,
0264                        isReconstructed);
0265     // Fill number of duplicated tracks for this particle
0266     m_duplicationPlotTool.fill(m_duplicationPlotCache, particle.initial(),
0267                                nMatchedTracks - 1);
0268
0269     // Fill number of reconstructed/truth-matched/fake tracks for this particle
0270     m_fakeRatePlotTool.fill(m_fakeRatePlotCache, particle.initial(),
0271                             nMatchedTracks, nFakeTracks);
0272
0273     m_nTotalParticles += 1;
0274   }
0275
0276   // Write additional stuff to TTree
0277   if (m_cfg.writeMatchingDetails && m_matchingTree != nullptr) {
0278     for (const auto& particle : particles) {
0279       auto particleId = particle.particleId();
0280
0281       m_treeEventNr = ctx.eventNumber;
0282       m_treeParticleId = particleId.value();
0283
0284       m_treeIsMatched = false;
0285       if (auto imatched = particleTrackMatching.find(particleId);
0286           imatched != particleTrackMatching.end()) {
0287         m_treeIsMatched = imatched->second.track.has_value();
0288       }
0289
0290       m_matchingTree->Fill();
0291     }
0292   }
0293
0294   return ProcessCode::SUCCESS;
0295 }
0296
0297 }  // namespace ActsExamples