Root/src/RootSpacePointPerformanceWriter.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/RootSpacePointPerformanceWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "ActsExamples/EventData/AverageSimHits.hpp"
0013 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0014 #include "ActsExamples/EventData/SimParticle.hpp"
0015 #include "ActsExamples/EventData/SpacePoint.hpp"
0016 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0017 #include "ActsExamples/Utilities/StripModulePairing.hpp"
0018 #include "ActsPlugins/Root/HistogramConverter.hpp"
0019
0020 #include <ios>
0021 #include <ranges>
0022 #include <stdexcept>
0023
0024 #include <TEfficiency.h>
0025 #include <TFile.h>
0026 #include <TH1.h>
0027 #include <TTree.h>
0028
0029 namespace ActsExamples {
0030
0031 namespace {
0032
0033 template <typename Begin, typename End>
0034 auto toRange(const std::pair<Begin, End>& pair) {
0035   return std::ranges::subrange(pair.first, pair.second);
0036 }
0037
0038 template <typename Range1, typename Range2>
0039 bool hasCommon(Range1&& a, Range2&& b) {
0040   for (const auto& elementA : a) {
0041     for (const auto& elementB : b) {
0042       if (elementA == elementB) {
0043         return true;
0044       }
0045     }
0046   }
0047   return false;
0048 }
0049
0050 }  // namespace
0051
0052 RootSpacePointPerformanceWriter::RootSpacePointPerformanceWriter(
0053     const Config& config, Acts::Logging::Level level)
0054     : WriterT(config.inputSpacePoints, "RootSpacePointPerformanceWriter",
0055               level),
0056       m_cfg(config) {
0057   // Input container for space points is already checked by base constructor
0058   if (m_cfg.inputParticles.empty()) {
0059     throw std::invalid_argument("Missing particles input collection");
0060   }
0061   if (m_cfg.inputMeasurements.empty()) {
0062     throw std::invalid_argument("Missing measurements input collection");
0063   }
0064   if (m_cfg.inputSimHits.empty()) {
0065     throw std::invalid_argument("Missing sim hits input collection");
0066   }
0067   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0068     throw std::invalid_argument(
0069         "Missing measurement-to-hits map input collection");
0070   }
0071   if (m_cfg.inputMeasurementParticlesMap.empty()) {
0072     throw std::invalid_argument(
0073         "Missing measurement-to-particles map input collection");
0074   }
0075   if (m_cfg.trackingGeometry == nullptr) {
0076     throw std::invalid_argument("Missing tracking geometry");
0077   }
0078
0079   m_inputParticles.initialize(m_cfg.inputParticles);
0080   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0081   m_inputSimHits.initialize(m_cfg.inputSimHits);
0082   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0083   m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
0084
0085   m_stripModulePairMap = pairStripModules(
0086       *m_cfg.trackingGeometry, m_cfg.stripGeometrySelection, logger());
0087
0088   // Setup ROOT File
0089   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0090   if (m_outputFile == nullptr) {
0091     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0092   }
0093
0094   m_outputFile->cd();
0095
0096   m_fakeVsZ.emplace("fake_vs_z", "Space point fake ratio vs z",
0097                     std::array{m_cfg.zAxis});
0098   m_fakeVsR.emplace("fake_vs_r", "Space point fake ratio vs r",
0099                     std::array{m_cfg.rAxis});
0100   m_fakeVsEta.emplace("fake_vs_eta", "Space point fake ratio vs eta",
0101                       std::array{m_cfg.etaAxis});
0102   m_fakeVsPhi.emplace("fake_vs_phi", "Space point fake ratio vs phi",
0103                       std::array{m_cfg.phiAxis});
0104
0105   m_effVsZ.emplace("efficiency_vs_z", "Space point efficiency vs z",
0106                    std::array{m_cfg.zAxis});
0107   m_effVsR.emplace("efficiency_vs_r", "Space point efficiency vs r",
0108                    std::array{m_cfg.rAxis});
0109   m_effVsEta.emplace("efficiency_vs_eta", "Space point efficiency vs eta",
0110                      std::array{m_cfg.etaAxis});
0111   m_effVsPhi.emplace("efficiency_vs_phi", "Space point efficiency vs phi",
0112                      std::array{m_cfg.phiAxis});
0113 }
0114
0115 RootSpacePointPerformanceWriter::~RootSpacePointPerformanceWriter() {
0116   if (m_outputFile != nullptr) {
0117     m_outputFile->Close();
0118   }
0119 }
0120
0121 ProcessCode RootSpacePointPerformanceWriter::finalize() {
0122   // Close the file if it's yours
0123   m_outputFile->cd();
0124
0125   ActsPlugins::toRoot(*m_fakeVsZ)->Write();
0126   ActsPlugins::toRoot(*m_fakeVsR)->Write();
0127   ActsPlugins::toRoot(*m_fakeVsEta)->Write();
0128   ActsPlugins::toRoot(*m_fakeVsPhi)->Write();
0129
0130   ActsPlugins::toRoot(*m_effVsZ)->Write();
0131   ActsPlugins::toRoot(*m_effVsR)->Write();
0132   ActsPlugins::toRoot(*m_effVsEta)->Write();
0133   ActsPlugins::toRoot(*m_effVsPhi)->Write();
0134
0135   m_outputFile->Close();
0136
0137   return ProcessCode::SUCCESS;
0138 }
0139
0140 ProcessCode RootSpacePointPerformanceWriter::writeT(
0141     const AlgorithmContext& ctx, const SpacePointContainer& spacePoints) {
0142   const auto& particles = m_inputParticles(ctx);
0143   const auto& measurements = m_inputMeasurements(ctx);
0144   const auto& simHits = m_inputSimHits(ctx);
0145   const auto& measurementSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0146   const auto& measurementParticlesMap = m_inputMeasurementParticlesMap(ctx);
0147
0148   const auto transformSecond =
0149       std::views::transform([](auto& p) -> auto& { return p.second; });
0150   const auto filterParticles = std::views::filter(
0151       [&](const SimBarcode& particle) { return particles.contains(particle); });
0152
0153   std::lock_guard<std::mutex> lock(m_writeMutex);
0154
0155   std::map<std::pair<Index, Index>, SpacePointIndex> trueSpacePoints;
0156
0157   for (const auto& spacePoint : spacePoints) {
0158     const Acts::Vector3 position(spacePoint.x(), spacePoint.y(),
0159                                  spacePoint.z());
0160     const double z = position.z();
0161     const double r = Acts::VectorHelpers::perp(position);
0162     const double phi = Acts::VectorHelpers::phi(position);
0163     const double eta = Acts::VectorHelpers::eta(position);
0164
0165     const auto& sourceLinks = spacePoint.sourceLinks();
0166
0167     bool isFake = false;
0168
0169     if (sourceLinks.size() == 2) {
0170       const auto measurementIndex1 =
0171           sourceLinks[0].get<IndexSourceLink>().index();
0172       const auto measurementIndex2 =
0173           sourceLinks[1].get<IndexSourceLink>().index();
0174
0175       auto contributingParticles1 =
0176           toRange(measurementParticlesMap.equal_range(measurementIndex1)) |
0177           transformSecond | filterParticles;
0178       auto contributingParticles2 =
0179           toRange(measurementParticlesMap.equal_range(measurementIndex2)) |
0180           transformSecond | filterParticles;
0181
0182       isFake = !hasCommon(contributingParticles1, contributingParticles2);
0183
0184       if (!isFake) {
0185         trueSpacePoints.emplace(
0186             std::pair<Index, Index>{measurementIndex1, measurementIndex2},
0187             spacePoint.index());
0188       }
0189     }
0190
0191     m_fakeVsZ->fill({z}, isFake);
0192     m_fakeVsR->fill({r}, isFake);
0193     m_fakeVsEta->fill({eta}, isFake);
0194     m_fakeVsPhi->fill({phi}, isFake);
0195   }
0196
0197   for (const auto [module1, module2] : m_stripModulePairMap) {
0198     const Acts::Surface* surface1 =
0199         m_cfg.trackingGeometry->findSurface(module1);
0200     const Acts::Surface* surface2 =
0201         m_cfg.trackingGeometry->findSurface(module2);
0202
0203     if (surface1 == nullptr || surface2 == nullptr) {
0204       ACTS_WARNING("Could not find surfaces for modules " << module1 << " and "
0205                                                           << module2);
0206       continue;
0207     }
0208
0209     const auto sourceLinks1 =
0210         measurements.orderedIndices().equal_range(module1);
0211     const auto sourceLinks2 =
0212         measurements.orderedIndices().equal_range(module2);
0213
0214     for (const auto& sourceLink1 : toRange(sourceLinks1)) {
0215       const auto measurementIndex1 = sourceLink1.index();
0216       auto contributingParticles1 =
0217           toRange(measurementParticlesMap.equal_range(measurementIndex1)) |
0218           transformSecond | filterParticles;
0219
0220       for (const auto& sourceLink2 : toRange(sourceLinks2)) {
0221         const auto measurementIndex2 = sourceLink2.index();
0222         auto contributingParticles2 =
0223             toRange(measurementParticlesMap.equal_range(measurementIndex2)) |
0224             transformSecond | filterParticles;
0225
0226         if (!hasCommon(contributingParticles1, contributingParticles2)) {
0227           continue;
0228         }
0229
0230         // Find the contributing simulated hits
0231         const auto hitRange1 =
0232             makeRange(measurementSimHitsMap.equal_range(measurementIndex1));
0233         // Use average truth in the case of multiple contributing sim hits
0234         const auto [local1, position1, dir1] = averageSimHits(
0235             ctx.geoContext, *surface1, simHits, hitRange1, logger());
0236
0237         const double z = position1.z();
0238         const double r = Acts::VectorHelpers::perp(position1);
0239         const double phi = Acts::VectorHelpers::phi(position1);
0240         const double eta = Acts::VectorHelpers::eta(position1.head<3>());
0241
0242         const bool isReconstructed =
0243             trueSpacePoints.contains(std::pair<Index, Index>{
0244                 measurementIndex1, measurementIndex2}) ||
0245             trueSpacePoints.contains(
0246                 std::pair<Index, Index>{measurementIndex2, measurementIndex1});
0247
0248         m_effVsZ->fill({z}, isReconstructed);
0249         m_effVsR->fill({r}, isReconstructed);
0250         m_effVsEta->fill({eta}, isReconstructed);
0251         m_effVsPhi->fill({phi}, isReconstructed);
0252       }
0253     }
0254   }
0255
0256   return ProcessCode::SUCCESS;
0257 }
0258
0259 }  // namespace ActsExamples