Root/src/RootMeasurementPerformanceWriter.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/RootMeasurementPerformanceWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "ActsExamples/EventData/Measurement.hpp"
0013 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0014 #include "ActsPlugins/Root/HistogramConverter.hpp"
0015
0016 #include <ios>
0017 #include <ranges>
0018 #include <stdexcept>
0019
0020 #include <TEfficiency.h>
0021 #include <TFile.h>
0022 #include <TH1.h>
0023 #include <TTree.h>
0024
0025 namespace ActsExamples {
0026
0027 namespace {
0028
0029 template <typename Begin, typename End>
0030 auto toRange(const std::pair<Begin, End>& pair) {
0031   return std::ranges::subrange(pair.first, pair.second);
0032 }
0033
0034 }  // namespace
0035
0036 RootMeasurementPerformanceWriter::RootMeasurementPerformanceWriter(
0037     const Config& config, Acts::Logging::Level level)
0038     : WriterT(config.inputMeasurements, "RootMeasurementPerformanceWriter",
0039               level),
0040       m_cfg(config) {
0041   // Input container for measurements is already checked by base constructor
0042   if (m_cfg.inputSimHits.empty()) {
0043     throw std::invalid_argument("Missing simulated hits input collection");
0044   }
0045   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0046     throw std::invalid_argument(
0047         "Missing measurement-to-hits map input collection");
0048   }
0049   if (m_cfg.inputMeasurementParticlesMap.empty()) {
0050     throw std::invalid_argument(
0051         "Missing measurement-to-particles map input collection");
0052   }
0053   if (m_cfg.inputSimHitMeasurementsMap.empty()) {
0054     throw std::invalid_argument(
0055         "Missing hits-to-measurements map input collection");
0056   }
0057
0058   m_inputSimHits.initialize(m_cfg.inputSimHits);
0059   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0060   m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
0061   m_inputSimHitMeasurementsMap.initialize(m_cfg.inputSimHitMeasurementsMap);
0062
0063   // Setup ROOT File
0064   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0065   if (m_outputFile == nullptr) {
0066     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0067   }
0068
0069   m_outputFile->cd();
0070
0071   m_measurementContributingHits.emplace(
0072       "measurement_n_contributing_hits",
0073       "Measurement number of contributing simulated hits",
0074       std::array{m_cfg.countAxis});
0075   m_measurementContributingParticles.emplace(
0076       "measurement_n_contributing_particles",
0077       "Measurement number of contributing particles",
0078       std::array{m_cfg.countAxis});
0079   m_measurementPurity.emplace("measurement_purity", "Measurements purity",
0080                               std::array{m_cfg.purityAxis});
0081   m_measurementClassification.emplace(
0082       "measurement_classification", "Measurement classification",
0083       std::array{Acts::Experimental::AxisVariant(
0084           Acts::Experimental::BoostRegularAxis{4, 0, 4, "Classification"})});
0085
0086   m_hitEffVsZ.emplace("hit_efficiency_vs_z",
0087                       "Hit to measurement efficiency vs z",
0088                       std::array{m_cfg.zAxis});
0089   m_hitEffVsR.emplace("hit_efficiency_vs_r",
0090                       "Hit to measurement efficiency vs r",
0091                       std::array{m_cfg.rAxis});
0092   m_hitEffVsEta.emplace("hit_efficiency_vs_eta",
0093                         "Hit to measurement efficiency vs eta",
0094                         std::array{m_cfg.etaAxis});
0095   m_hitEffVsPhi.emplace("hit_efficiency_vs_phi",
0096                         "Hit to measurement efficiency vs phi",
0097                         std::array{m_cfg.phiAxis});
0098 }
0099
0100 RootMeasurementPerformanceWriter::~RootMeasurementPerformanceWriter() {
0101   if (m_outputFile != nullptr) {
0102     m_outputFile->Close();
0103   }
0104 }
0105
0106 ProcessCode RootMeasurementPerformanceWriter::finalize() {
0107   // Close the file if it's yours
0108   m_outputFile->cd();
0109
0110   ActsPlugins::toRoot(*m_measurementContributingHits)->Write();
0111   ActsPlugins::toRoot(*m_measurementContributingParticles)->Write();
0112   ActsPlugins::toRoot(*m_measurementPurity)->Write();
0113   ActsPlugins::toRoot(*m_measurementClassification)->Write();
0114
0115   ActsPlugins::toRoot(*m_hitEffVsZ)->Write();
0116   ActsPlugins::toRoot(*m_hitEffVsR)->Write();
0117   ActsPlugins::toRoot(*m_hitEffVsEta)->Write();
0118   ActsPlugins::toRoot(*m_hitEffVsPhi)->Write();
0119
0120   m_outputFile->Close();
0121
0122   return ProcessCode::SUCCESS;
0123 }
0124
0125 ProcessCode RootMeasurementPerformanceWriter::writeT(
0126     const AlgorithmContext& ctx, const MeasurementContainer& measurements) {
0127   const auto& simHits = m_inputSimHits(ctx);
0128   const auto& measurementSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0129   const auto& measurementParticlesMap = m_inputMeasurementParticlesMap(ctx);
0130   const auto& simHitMeasurementsMap = m_inputSimHitMeasurementsMap(ctx);
0131
0132   std::lock_guard<std::mutex> lock(m_writeMutex);
0133
0134   for (const auto& measurement : measurements) {
0135     const auto contributingSimHits =
0136         toRange(measurementSimHitsMap.equal_range(measurement.index()));
0137     const auto contributingParticles =
0138         toRange(measurementParticlesMap.equal_range(measurement.index()));
0139
0140     const double nContributingSimHits = contributingSimHits.size();
0141     const double nContributingParticles = contributingParticles.size();
0142     const double purity = 1.0 / nContributingParticles;
0143     const MeasurementClassification classification = [&]() {
0144       if (nContributingParticles == 0) {
0145         return MeasurementClassification::Fake;
0146       } else if (purity > m_cfg.matchingRatio) {
0147         return MeasurementClassification::Matched;
0148       } else {
0149         return MeasurementClassification::Merged;
0150       }
0151     }();
0152
0153     m_measurementContributingHits->fill({nContributingSimHits});
0154     m_measurementContributingParticles->fill({nContributingParticles});
0155     m_measurementPurity->fill({purity});
0156     m_measurementClassification->fill({static_cast<double>(classification)});
0157   }
0158
0159   for (const auto& simHit : simHits) {
0160     const auto relatedMeasurements =
0161         toRange(simHitMeasurementsMap.equal_range(simHit.index()));
0162     const std::size_t nRelatedMeasurements = relatedMeasurements.size();
0163     const bool efficient = nRelatedMeasurements > 0;
0164
0165     const Acts::Vector3& position = simHit.position();
0166     const double z = position.z();
0167     const double r = Acts::VectorHelpers::perp(position);
0168     const double phi = Acts::VectorHelpers::phi(position);
0169     const double eta = Acts::VectorHelpers::eta(position);
0170
0171     m_hitEffVsZ->fill({z}, efficient);
0172     m_hitEffVsR->fill({r}, efficient);
0173     m_hitEffVsEta->fill({eta}, efficient);
0174     m_hitEffVsPhi->fill({phi}, efficient);
0175   }
0176
0177   return ProcessCode::SUCCESS;
0178 }
0179
0180 }  // namespace ActsExamples