Root/src/RootTrackFitterPerformanceWriter.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/RootTrackFitterPerformanceWriter.hpp"
0010
0011 #include "Acts/EventData/TrackParameters.hpp"
0012 #include "Acts/Utilities/Helpers.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014 #include "ActsExamples/EventData/Track.hpp"
0015 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0016 #include "ActsExamples/Validation/TrackClassification.hpp"
0017 #include "ActsFatras/EventData/Barcode.hpp"
0018
0019 #include <cstddef>
0020 #include <ostream>
0021 #include <stdexcept>
0022 #include <utility>
0023 #include <vector>
0024
0025 #include <TFile.h>
0026
0027 using Acts::VectorHelpers::eta;
0028 using Acts::VectorHelpers::phi;
0029
0030 ActsExamples::RootTrackFitterPerformanceWriter::
0031     RootTrackFitterPerformanceWriter(
0032         ActsExamples::RootTrackFitterPerformanceWriter::Config config,
0033         Acts::Logging::Level level)
0034     : WriterT(config.inputTracks, "RootTrackFitterPerformanceWriter", level),
0035       m_cfg(std::move(config)),
0036       m_resPlotTool(m_cfg.resPlotToolConfig, level),
0037       m_effPlotTool(m_cfg.effPlotToolConfig, level),
0038       m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig, level) {
0039   // trajectories 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.filePath.empty()) {
0047     throw std::invalid_argument("Missing output filename");
0048   }
0049
0050   m_inputParticles.initialize(m_cfg.inputParticles);
0051   m_inputTrackParticleMatching.initialize(m_cfg.inputTrackParticleMatching);
0052
0053   // the output file can not be given externally since TFile accesses to the
0054   // same file from multiple threads are unsafe.
0055   // must always be opened internally
0056   auto path = m_cfg.filePath;
0057   m_outputFile = TFile::Open(path.c_str(), "RECREATE");
0058   if (m_outputFile == nullptr) {
0059     throw std::invalid_argument("Could not open '" + path + "'");
0060   }
0061
0062   // initialize the residual and efficiency plots tool
0063   m_resPlotTool.book(m_resPlotCache);
0064   m_effPlotTool.book(m_effPlotCache);
0065   m_trackSummaryPlotTool.book(m_trackSummaryPlotCache);
0066 }
0067
0068 ActsExamples::RootTrackFitterPerformanceWriter::
0069     ~RootTrackFitterPerformanceWriter() {
0070   m_resPlotTool.clear(m_resPlotCache);
0071   m_effPlotTool.clear(m_effPlotCache);
0072   m_trackSummaryPlotTool.clear(m_trackSummaryPlotCache);
0073
0074   if (m_outputFile != nullptr) {
0075     m_outputFile->Close();
0076   }
0077 }
0078
0079 ActsExamples::ProcessCode
0080 ActsExamples::RootTrackFitterPerformanceWriter::finalize() {
0081   // fill residual and pull details into additional hists
0082   m_resPlotTool.refinement(m_resPlotCache);
0083
0084   if (m_outputFile != nullptr) {
0085     m_outputFile->cd();
0086     m_resPlotTool.write(m_resPlotCache);
0087     m_effPlotTool.write(m_effPlotCache);
0088     m_trackSummaryPlotTool.write(m_trackSummaryPlotCache);
0089
0090     ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
0091   }
0092   return ProcessCode::SUCCESS;
0093 }
0094
0095 ActsExamples::ProcessCode
0096 ActsExamples::RootTrackFitterPerformanceWriter::writeT(
0097     const AlgorithmContext& ctx, const ConstTrackContainer& tracks) {
0098   // Read truth input collections
0099   const auto& particles = m_inputParticles(ctx);
0100   const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx);
0101
0102   // Truth particles with corresponding reconstructed tracks
0103   std::vector<ActsFatras::Barcode> reconParticleIds;
0104   reconParticleIds.reserve(particles.size());
0105   // For each particle within a track, how many hits did it contribute
0106   std::vector<ParticleHitCount> particleHitCounts;
0107
0108   // Exclusive access to the tree while writing
0109   std::lock_guard<std::mutex> lock(m_writeMutex);
0110
0111   // Loop over all tracks
0112   for (const auto& track : tracks) {
0113     // Select reco track with fitted parameters
0114     if (!track.hasReferenceSurface()) {
0115       ACTS_WARNING("No fitted track parameters.");
0116       continue;
0117     }
0118     Acts::BoundTrackParameters fittedParameters =
0119         track.createParametersAtReference();
0120
0121     // Get the truth-matched particle
0122     auto imatched = trackParticleMatching.find(track.index());
0123     if (imatched == trackParticleMatching.end()) {
0124       ACTS_DEBUG("No truth particle associated with this track, index = "
0125                  << track.index() << " tip index = " << track.tipIndex());
0126       continue;
0127     }
0128     const auto& particleMatch = imatched->second;
0129
0130     if (!particleMatch.particle.has_value()) {
0131       ACTS_DEBUG("No truth particle associated with this track.");
0132       continue;
0133     }
0134
0135     // Get the barcode of the majority truth particle
0136     SimBarcode majorityParticleId = particleMatch.particle.value();
0137
0138     // Find the truth particle via the barcode
0139     auto ip = particles.find(majorityParticleId);
0140     if (ip == particles.end()) {
0141       ACTS_DEBUG("Majority particle not found in the particles collection.");
0142       continue;
0143     }
0144
0145     // Record this majority particle ID of this trajectory
0146     reconParticleIds.push_back(ip->particleId());
0147     // Fill the residual plots
0148     m_resPlotTool.fill(m_resPlotCache, ctx.geoContext, ip->initialState(),
0149                        fittedParameters);
0150     // Fill the trajectory summary info
0151     m_trackSummaryPlotTool.fill(m_trackSummaryPlotCache, fittedParameters,
0152                                 track.nTrackStates(), track.nMeasurements(),
0153                                 track.nOutliers(), track.nHoles(),
0154                                 track.nSharedHits());
0155   }
0156
0157   // Fill the efficiency, defined as the ratio between number of tracks with
0158   // fitted parameter and total truth tracks (assumes one truth particle has
0159   // one truth track)
0160   for (const auto& particle : particles) {
0161     bool isReconstructed = false;
0162     // Check if the particle has been reconstructed
0163     if (Acts::rangeContainsValue(reconParticleIds, particle.particleId())) {
0164       isReconstructed = true;
0165     }
0166     // Loop over all the other truth particle and find the distance to the
0167     // closest one
0168     double minDeltaR = -1;
0169     for (const auto& closeParticle : particles) {
0170       if (closeParticle.particleId() == particle.particleId()) {
0171         continue;
0172       }
0173       double p_phi = phi(particle.direction());
0174       double p_eta = eta(particle.direction());
0175       double c_phi = phi(closeParticle.direction());
0176       double c_eta = eta(closeParticle.direction());
0177       double distance = sqrt(pow(p_phi - c_phi, 2) + pow(p_eta - c_eta, 2));
0178       if (minDeltaR == -1 || distance < minDeltaR) {
0179         minDeltaR = distance;
0180       }
0181     }
0182     m_effPlotTool.fill(ctx.geoContext, m_effPlotCache, particle.initialState(),
0183                        minDeltaR, isReconstructed);
0184   }
0185
0186   return ProcessCode::SUCCESS;
0187 }