Csv/src/CsvPlanarClusterWriter.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017 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/CsvPlanarClusterWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Digitization/DigitizationCell.hpp"
0014 #include "Acts/Digitization/DigitizationSourceLink.hpp"
0015 #include "Acts/Digitization/PlanarModuleCluster.hpp"
0016 #include "Acts/EventData/SourceLink.hpp"
0017 #include "Acts/Geometry/GeometryIdentifier.hpp"
0018 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0021 #include "ActsExamples/Utilities/GroupBy.hpp"
0022 #include "ActsExamples/Utilities/Paths.hpp"
0023 #include "ActsExamples/Utilities/Range.hpp"
0024 #include "ActsFatras/EventData/Barcode.hpp"
0025 #include "ActsFatras/EventData/Hit.hpp"
0026
0027 #include <ostream>
0028 #include <stdexcept>
0029 #include <utility>
0030 #include <vector>
0031
0032 #include <dfe/dfe_io_dsv.hpp>
0033
0034 #include "CsvOutputData.hpp"
0035
0036 ActsExamples::CsvPlanarClusterWriter::CsvPlanarClusterWriter(
0037     const ActsExamples::CsvPlanarClusterWriter::Config& config,
0038     Acts::Logging::Level level)
0039     : WriterT(config.inputClusters, "CsvPlanarClusterWriter", level),
0040       m_cfg(config) {
0041   // inputClusters 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.trackingGeometry) {
0046     throw std::invalid_argument("Missing tracking geometry");
0047   }
0048
0049   m_inputSimHits.initialize(m_cfg.inputSimHits);
0050 }
0051
0052 ActsExamples::ProcessCode ActsExamples::CsvPlanarClusterWriter::writeT(
0053     const AlgorithmContext& ctx,
0054     const ActsExamples::GeometryIdMultimap<Acts::PlanarModuleCluster>&
0055         clusters) {
0056   // retrieve simulated hits
0057   const auto& simHits = m_inputSimHits(ctx);
0058
0059   // open per-event file for all components
0060   std::string pathHits =
0061       perEventFilepath(m_cfg.outputDir, "hits.csv", ctx.eventNumber);
0062   std::string pathCells =
0063       perEventFilepath(m_cfg.outputDir, "cells.csv", ctx.eventNumber);
0064   std::string pathTruth =
0065       perEventFilepath(m_cfg.outputDir, "truth.csv", ctx.eventNumber);
0066
0067   dfe::NamedTupleCsvWriter<HitData> writerHits(pathHits, m_cfg.outputPrecision);
0068   dfe::NamedTupleCsvWriter<CellDataLegacy> writerCells(pathCells,
0069                                                        m_cfg.outputPrecision);
0070   dfe::NamedTupleCsvWriter<TruthHitData> writerTruth(pathTruth,
0071                                                      m_cfg.outputPrecision);
0072
0073   HitData hit;
0074   CellDataLegacy cell;
0075   TruthHitData truth;
0076   // will be reused as hit counter
0077   hit.hit_id = 0;
0078
0079   for (auto [moduleGeoId, moduleClusters] : groupByModule(clusters)) {
0080     const Acts::Surface* surfacePtr =
0081         m_cfg.trackingGeometry->findSurface(moduleGeoId);
0082     if (surfacePtr == nullptr) {
0083       ACTS_ERROR("Could not find surface for " << moduleGeoId);
0084       return ProcessCode::ABORT;
0085     }
0086
0087     // encoded geometry identifier. same for all hits on the module
0088     hit.geometry_id = moduleGeoId.value();
0089
0090     for (const auto& entry : moduleClusters) {
0091       const Acts::PlanarModuleCluster& cluster = entry.second;
0092       // local cluster information
0093       const auto& parameters = cluster.parameters();
0094       Acts::Vector2 localPos(parameters[0], parameters[1]);
0095       Acts::Vector3 globalFakeMom(1, 1, 1);
0096       Acts::Vector3 globalPos =
0097           surfacePtr->localToGlobal(ctx.geoContext, localPos, globalFakeMom);
0098
0099       // write global hit information
0100       hit.x = globalPos.x() / Acts::UnitConstants::mm;
0101       hit.y = globalPos.y() / Acts::UnitConstants::mm;
0102       hit.z = globalPos.z() / Acts::UnitConstants::mm;
0103       hit.t = parameters[2] / Acts::UnitConstants::mm;
0104       writerHits.append(hit);
0105
0106       // write local cell information
0107       cell.hit_id = hit.hit_id;
0108       for (auto& c : cluster.digitizationCells()) {
0109         cell.channel0 = c.channel0;
0110         cell.channel1 = c.channel1;
0111         // TODO store digital timestamp once added to the cell definition
0112         cell.timestamp = 0;
0113         cell.value = c.data;
0114         writerCells.append(cell);
0115       }
0116
0117       // write hit-particle truth association
0118       // each hit can have multiple particles, e.g. in a dense environment
0119       truth.hit_id = hit.hit_id;
0120       truth.geometry_id = hit.geometry_id;
0121       const auto& sl = cluster.sourceLink().get<Acts::DigitizationSourceLink>();
0122       for (auto idx : sl.indices()) {
0123         auto it = simHits.nth(idx);
0124         if (it == simHits.end()) {
0125           ACTS_FATAL("Simulation hit with index " << idx << " does not exist");
0126           return ProcessCode::ABORT;
0127         }
0128
0129         const auto& simHit = *it;
0130         truth.particle_id = simHit.particleId().value();
0131         // hit position
0132         truth.tx = simHit.position().x() / Acts::UnitConstants::mm;
0133         truth.ty = simHit.position().y() / Acts::UnitConstants::mm;
0134         truth.tz = simHit.position().z() / Acts::UnitConstants::mm;
0135         truth.tt = simHit.time() / Acts::UnitConstants::mm;
0136         // particle four-momentum before interaction
0137         truth.tpx = simHit.momentum4Before().x() / Acts::UnitConstants::GeV;
0138         truth.tpy = simHit.momentum4Before().y() / Acts::UnitConstants::GeV;
0139         truth.tpz = simHit.momentum4Before().z() / Acts::UnitConstants::GeV;
0140         truth.te = simHit.momentum4Before().w() / Acts::UnitConstants::GeV;
0141         // particle four-momentum change due to interaction
0142         const auto delta4 = simHit.momentum4After() - simHit.momentum4Before();
0143         truth.deltapx = delta4.x() / Acts::UnitConstants::GeV;
0144         truth.deltapy = delta4.y() / Acts::UnitConstants::GeV;
0145         truth.deltapz = delta4.z() / Acts::UnitConstants::GeV;
0146         truth.deltae = delta4.w() / Acts::UnitConstants::GeV;
0147         // TODO write hit index along the particle trajectory
0148         truth.index = simHit.index();
0149         writerTruth.append(truth);
0150       }
0151
0152       // increase hit id for next iteration
0153       hit.hit_id += 1;
0154     }
0155   }
0156
0157   return ActsExamples::ProcessCode::SUCCESS;
0158 }