Csv/src/CsvMeasurementWriter.cpp

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