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