Csv/src/CsvMeasurementWriter.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/Csv/CsvMeasurementWriter.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "ActsExamples/EventData/Cluster.hpp"
0015 #include "ActsExamples/EventData/Index.hpp"
0016 #include "ActsExamples/EventData/Measurement.hpp"
0017 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0018 #include "ActsExamples/Io/Csv/CsvInputOutput.hpp"
0019 #include "ActsExamples/Utilities/Paths.hpp"
0020 #include "ActsExamples/Utilities/Range.hpp"
0021
0022 #include <optional>
0023 #include <ostream>
0024 #include <stdexcept>
0025
0026 #include "CsvOutputData.hpp"
0027
0028 namespace ActsExamples {
0029
0030 CsvMeasurementWriter::CsvMeasurementWriter(const Config& config,
0031                                            Acts::Logging::Level level)
0032     : WriterT(config.inputMeasurements, "CsvMeasurementWriter", level),
0033       m_cfg(config) {
0034   // Input container for measurements is already checked by base constructor
0035   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0036     throw std::invalid_argument(
0037         "Missing hit-to-simulated-hits map input collection");
0038   }
0039
0040   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0041   m_inputClusters.maybeInitialize(m_cfg.inputClusters);
0042 }
0043
0044 CsvMeasurementWriter::~CsvMeasurementWriter() = default;
0045
0046 ProcessCode CsvMeasurementWriter::finalize() {
0047   // Write the tree
0048   return ProcessCode::SUCCESS;
0049 }
0050
0051 ProcessCode CsvMeasurementWriter::writeT(
0052     const AlgorithmContext& ctx, const MeasurementContainer& measurements) {
0053   const auto& measurementSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0054
0055   ClusterContainer clusters;
0056
0057   // Open per-event file for all components
0058   std::string pathMeasurements =
0059       perEventFilepath(m_cfg.outputDir, "measurements.csv", ctx.eventNumber);
0060   std::string pathMeasurementSimHitMap = perEventFilepath(
0061       m_cfg.outputDir, "measurement-simhit-map.csv", ctx.eventNumber);
0062
0063   BoostDescribeCsvWriter<MeasurementData> writerMeasurements(
0064       pathMeasurements, m_cfg.outputPrecision);
0065
0066   std::optional<BoostDescribeCsvWriter<CellData>> writerCells{std::nullopt};
0067   if (!m_cfg.inputClusters.empty()) {
0068     ACTS_VERBOSE(
0069         "Set up writing of clusters from collection: " << m_cfg.inputClusters);
0070     clusters = m_inputClusters(ctx);
0071     std::string pathCells =
0072         perEventFilepath(m_cfg.outputDir, "cells.csv", ctx.eventNumber);
0073     writerCells =
0074         BoostDescribeCsvWriter<CellData>{pathCells, m_cfg.outputPrecision};
0075   }
0076
0077   BoostDescribeCsvWriter<MeasurementSimHitLink> writerMeasurementSimHitMap(
0078       pathMeasurementSimHitMap, m_cfg.outputPrecision);
0079
0080   MeasurementData meas;
0081   CellData cell;
0082
0083   // Will be reused as measurement counter
0084   meas.measurement_id = 0;
0085
0086   ACTS_VERBOSE("Writing " << measurements.size()
0087                           << " measurements in this event.");
0088
0089   for (Index measIdx = 0u; measIdx < measurements.size(); ++measIdx) {
0090     const ConstVariableBoundMeasurementProxy measurement =
0091         measurements.getMeasurement(measIdx);
0092
0093     auto simHitIndices = makeRange(measurementSimHitsMap.equal_range(measIdx));
0094     for (auto [_, simHitIdx] : simHitIndices) {
0095       writerMeasurementSimHitMap.append({measIdx, simHitIdx});
0096     }
0097
0098     Acts::GeometryIdentifier geoId = measurement.geometryId();
0099     // MEASUREMENT information ------------------------------------
0100
0101     // Encoded geometry identifier. same for all hits on the module
0102     meas.geometry_id = geoId.value();
0103     meas.local_key = 0;
0104     // Create a full set of parameters
0105     auto parameters = measurement.fullParameters();
0106     meas.local0 = parameters[Acts::eBoundLoc0] / Acts::UnitConstants::mm;
0107     meas.local1 = parameters[Acts::eBoundLoc1] / Acts::UnitConstants::mm;
0108     meas.phi = parameters[Acts::eBoundPhi] / Acts::UnitConstants::rad;
0109     meas.theta = parameters[Acts::eBoundTheta] / Acts::UnitConstants::rad;
0110     meas.time = parameters[Acts::eBoundTime] / Acts::UnitConstants::mm;
0111
0112     auto covariance = measurement.fullCovariance();
0113     meas.var_local0 = covariance(Acts::eBoundLoc0, Acts::eBoundLoc0);
0114     meas.var_local1 = covariance(Acts::eBoundLoc1, Acts::eBoundLoc1);
0115     meas.var_phi = covariance(Acts::eBoundPhi, Acts::eBoundPhi);
0116     meas.var_theta = covariance(Acts::eBoundTheta, Acts::eBoundTheta);
0117     meas.var_time = covariance(Acts::eBoundTime, Acts::eBoundTime);
0118     for (unsigned int ipar = 0;
0119          ipar < static_cast<unsigned int>(Acts::eBoundSize); ++ipar) {
0120       if (measurement.contains(static_cast<Acts::BoundIndices>(ipar))) {
0121         meas.local_key = ((1 << (ipar + 1)) | meas.local_key);
0122       }
0123     }
0124
0125     if (!clusters.empty()) {
0126       const auto& cluster = clusters.at(measIdx);
0127       meas.global_x = cluster.globalPosition.x();
0128       meas.global_y = cluster.globalPosition.y();
0129       meas.global_z = cluster.globalPosition.z();
0130     } else {
0131       meas.global_x = std::numeric_limits<float>::quiet_NaN();
0132       meas.global_y = std::numeric_limits<float>::quiet_NaN();
0133       meas.global_z = std::numeric_limits<float>::quiet_NaN();
0134     }
0135
0136     writerMeasurements.append(meas);
0137
0138     // CLUSTER / channel information ------------------------------
0139     if (!clusters.empty() && writerCells) {
0140       const auto& cluster = clusters.at(measIdx);
0141       cell.geometry_id = meas.geometry_id;
0142       cell.measurement_id = meas.measurement_id;
0143       for (auto& c : cluster.channels) {
0144         cell.channel0 = c.bin[0];
0145         cell.channel1 = c.bin[1];
0146         // TODO store digital timestamp once added to the cell definition
0147         cell.timestamp = 0;
0148         cell.value = c.activation;
0149         writerCells->append(cell);
0150       }
0151     }
0152     // Increase counter
0153     meas.measurement_id += 1;
0154   }
0155   return ProcessCode::SUCCESS;
0156 }
0157
0158 }  // namespace ActsExamples