Root/src/SeedingPerformanceWriter.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/SeedingPerformanceWriter.hpp"
0010
0011 #include "Acts/Utilities/MultiIndex.hpp"
0012 #include "Acts/Utilities/VectorHelpers.hpp"
0013 #include "ActsExamples/Utilities/EventDataTransforms.hpp"
0014 #include "ActsExamples/Validation/TrackClassification.hpp"
0015 #include "ActsFatras/EventData/Barcode.hpp"
0016
0017 #include <cstddef>
0018 #include <ostream>
0019 #include <stdexcept>
0020 #include <unordered_map>
0021 #include <utility>
0022 #include <vector>
0023
0024 #include <TFile.h>
0025 #include <TVectorT.h>
0026
0027 using Acts::VectorHelpers::eta;
0028 using Acts::VectorHelpers::phi;
0029
0030 ActsExamples::SeedingPerformanceWriter::SeedingPerformanceWriter(
0031     ActsExamples::SeedingPerformanceWriter::Config config,
0032     Acts::Logging::Level level)
0033     : WriterT(config.inputSeeds, "SeedingPerformanceWriter", level),
0034       m_cfg(std::move(config)),
0035       m_effPlotTool(m_cfg.effPlotToolConfig, level),
0036       m_duplicationPlotTool(m_cfg.duplicationPlotToolConfig, level) {
0037   if (m_cfg.inputMeasurementParticlesMap.empty()) {
0038     throw std::invalid_argument("Missing hit-particles map input collection");
0039   }
0040   if (m_cfg.inputParticles.empty()) {
0041     throw std::invalid_argument("Missing input particles collection");
0042   }
0043   if (m_cfg.filePath.empty()) {
0044     throw std::invalid_argument("Missing output filename");
0045   }
0046
0047   m_inputParticles.initialize(m_cfg.inputParticles);
0048   m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
0049
0050   // the output file can not be given externally since TFile accesses to the
0051   // same file from multiple threads are unsafe.
0052   // must always be opened internally
0053   auto path = m_cfg.filePath;
0054   m_outputFile = TFile::Open(path.c_str(), m_cfg.fileMode.c_str());
0055   if (m_outputFile == nullptr) {
0056     throw std::invalid_argument("Could not open '" + path + "'");
0057   }
0058   // initialize the plot tools
0059   m_effPlotTool.book(m_effPlotCache);
0060   m_duplicationPlotTool.book(m_duplicationPlotCache);
0061 }
0062
0063 ActsExamples::SeedingPerformanceWriter::~SeedingPerformanceWriter() {
0064   m_effPlotTool.clear(m_effPlotCache);
0065   m_duplicationPlotTool.clear(m_duplicationPlotCache);
0066   if (m_outputFile != nullptr) {
0067     m_outputFile->Close();
0068   }
0069 }
0070
0071 ActsExamples::ProcessCode ActsExamples::SeedingPerformanceWriter::finalize() {
0072   float eff = static_cast<float>(m_nTotalMatchedParticles) / m_nTotalParticles;
0073   float fakeRate =
0074       static_cast<float>(m_nTotalSeeds - m_nTotalMatchedSeeds) / m_nTotalSeeds;
0075   float duplicationRate = static_cast<float>(m_nTotalDuplicatedParticles) /
0076                           m_nTotalMatchedParticles;
0077   float aveNDuplicatedSeeds =
0078       static_cast<float>(m_nTotalMatchedSeeds - m_nTotalMatchedParticles) /
0079       m_nTotalMatchedParticles;
0080   float totalSeedPurity =
0081       static_cast<float>(m_nTotalMatchedSeeds) / m_nTotalSeeds;
0082   ACTS_DEBUG("nTotalSeeds               = " << m_nTotalSeeds);
0083   ACTS_DEBUG("nTotalMatchedSeeds        = " << m_nTotalMatchedSeeds);
0084   ACTS_DEBUG("nTotalParticles           = " << m_nTotalParticles);
0085   ACTS_DEBUG("nTotalMatchedParticles    = " << m_nTotalMatchedParticles);
0086   ACTS_DEBUG("nTotalDuplicatedParticles = " << m_nTotalDuplicatedParticles);
0087
0088   ACTS_INFO("Efficiency (nMatchedParticles / nAllParticles) = " << eff);
0089   ACTS_INFO("Fake rate (nUnMatchedSeeds / nAllSeeds) = " << fakeRate);
0090   ACTS_INFO("Total seed purity (nTotalMatchedSeeds / m_nTotalSeeds) = "
0091             << totalSeedPurity);
0092   ACTS_INFO(
0093       "Duplication rate (nDuplicatedMatchedParticles / nMatchedParticles) = "
0094       << duplicationRate);
0095   ACTS_INFO(
0096       "Average number of duplicated seeds ((nMatchedSeeds - nMatchedParticles) "
0097       "/ nMatchedParticles) = "
0098       << aveNDuplicatedSeeds);
0099
0100   auto writeFloat = [&](float f, const char* name) {
0101     TVectorF v(1);
0102     v[0] = f;
0103     m_outputFile->WriteObject(&v, name);
0104   };
0105
0106   if (m_outputFile != nullptr) {
0107     m_outputFile->cd();
0108     m_effPlotTool.write(m_effPlotCache);
0109     m_duplicationPlotTool.write(m_duplicationPlotCache);
0110     writeFloat(eff, "eff_seeds");
0111     writeFloat(fakeRate, "fakerate_seeds");
0112     writeFloat(duplicationRate, "duplicaterate_seeds");
0113     writeFloat(totalSeedPurity, "purity_seeds");
0114     ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'");
0115   }
0116   return ProcessCode::SUCCESS;
0117 }
0118
0119 ActsExamples::ProcessCode ActsExamples::SeedingPerformanceWriter::writeT(
0120     const AlgorithmContext& ctx, const SimSeedContainer& seeds) {
0121   // Read truth information collections
0122   const auto& particles = m_inputParticles(ctx);
0123   const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx);
0124
0125   // Exclusive access to the tree while writing
0126   std::lock_guard<std::mutex> lock(m_writeMutex);
0127
0128   std::size_t nSeeds = seeds.size();
0129   std::size_t nMatchedSeeds = 0;
0130   // Map from particles to how many times they were successfully found by a seed
0131   std::unordered_map<ActsFatras::Barcode, std::size_t> truthCount;
0132
0133   for (std::size_t itrack = 0; itrack < seeds.size(); ++itrack) {
0134     const auto& seed = seeds[itrack];
0135     const auto track = seedToPrototrack(seed);
0136     std::vector<ParticleHitCount> particleHitCounts;
0137     identifyContributingParticles(hitParticlesMap, track, particleHitCounts);
0138     // All hits matched to the same particle
0139     if (particleHitCounts.size() == 1) {
0140       auto prt = particleHitCounts.at(0);
0141       auto it = truthCount.try_emplace(prt.particleId, 0u).first;
0142       it->second += 1;
0143       nMatchedSeeds++;
0144     }
0145   }
0146
0147   int nMatchedParticles = 0;
0148   int nDuplicatedParticles = 0;
0149   // Fill the efficiency and fake rate plots
0150   for (const auto& particle : particles) {
0151     const auto it1 = truthCount.find(particle.particleId());
0152     bool isMatched = false;
0153     int nMatchedSeedsForParticle = 0;
0154     if (it1 != truthCount.end()) {
0155       isMatched = true;
0156       nMatchedParticles++;
0157       nMatchedSeedsForParticle = it1->second;
0158       if (nMatchedSeedsForParticle > 1) {
0159         nDuplicatedParticles++;
0160       }
0161     }
0162     // Loop over all the other truth particle and find the distance to the
0163     // closest one
0164     double minDeltaR = -1;
0165     for (const auto& closeParticle : particles) {
0166       if (closeParticle.particleId() == particle.particleId()) {
0167         continue;
0168       }
0169       double p_phi = phi(particle.direction());
0170       double p_eta = eta(particle.direction());
0171       double c_phi = phi(closeParticle.direction());
0172       double c_eta = eta(closeParticle.direction());
0173       double distance = sqrt(pow(p_phi - c_phi, 2) + pow(p_eta - c_eta, 2));
0174       if (minDeltaR == -1 || distance < minDeltaR) {
0175         minDeltaR = distance;
0176       }
0177     }
0178     m_effPlotTool.fill(m_effPlotCache, particle.initial(), minDeltaR,
0179                        isMatched);
0180     m_duplicationPlotTool.fill(m_duplicationPlotCache, particle.initial(),
0181                                nMatchedSeedsForParticle - 1);
0182   }
0183   ACTS_DEBUG("Number of seeds: " << nSeeds);
0184   m_nTotalSeeds += nSeeds;
0185   m_nTotalMatchedSeeds += nMatchedSeeds;
0186   m_nTotalParticles += particles.size();
0187   m_nTotalMatchedParticles += nMatchedParticles;
0188   m_nTotalDuplicatedParticles += nDuplicatedParticles;
0189
0190   return ProcessCode::SUCCESS;
0191 }