Csv/src/CsvTrackingGeometryWriter.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/CsvTrackingGeometryWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/BoundarySurfaceT.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "Acts/Geometry/GeometryIdentifier.hpp"
0016 #include "Acts/Geometry/Layer.hpp"
0017 #include "Acts/Geometry/TrackingGeometry.hpp"
0018 #include "Acts/Geometry/TrackingVolume.hpp"
0019 #include "Acts/Geometry/Volume.hpp"
0020 #include "Acts/Geometry/VolumeBounds.hpp"
0021 #include "Acts/Surfaces/Surface.hpp"
0022 #include "Acts/Surfaces/SurfaceArray.hpp"
0023 #include "Acts/Surfaces/SurfaceBounds.hpp"
0024 #include "Acts/Utilities/BinnedArray.hpp"
0025 #include "Acts/Utilities/IAxis.hpp"
0026 #include "Acts/Utilities/Logger.hpp"
0027 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0028 #include "ActsExamples/Io/Csv/CsvInputOutput.hpp"
0029 #include "ActsExamples/Utilities/Paths.hpp"
0030
0031 #include <array>
0032 #include <cstddef>
0033 #include <stdexcept>
0034 #include <utility>
0035 #include <vector>
0036
0037 #include "CsvOutputData.hpp"
0038
0039 using namespace ActsExamples;
0040
0041 CsvTrackingGeometryWriter::CsvTrackingGeometryWriter(
0042     const CsvTrackingGeometryWriter::Config& config, Acts::Logging::Level level)
0043     : m_cfg(config),
0044       m_logger(Acts::getDefaultLogger("CsvTrackingGeometryWriter", level))
0045
0046 {
0047   if (!m_cfg.trackingGeometry) {
0048     throw std::invalid_argument("Missing tracking geometry");
0049   }
0050   m_world = m_cfg.trackingGeometry->highestTrackingVolume();
0051   if (m_world == nullptr) {
0052     throw std::invalid_argument("Could not identify the world volume");
0053   }
0054 }
0055
0056 std::string CsvTrackingGeometryWriter::name() const {
0057   return "CsvTrackingGeometryWriter";
0058 }
0059
0060 namespace {
0061
0062 using SurfaceWriter = ActsExamples::NamedTupleCsvWriter<SurfaceData>;
0063 using SurfaceGridWriter = ActsExamples::NamedTupleCsvWriter<SurfaceGridData>;
0064 using LayerVolumeWriter = ActsExamples::NamedTupleCsvWriter<LayerVolumeData>;
0065 using BoundarySurface = Acts::BoundarySurfaceT<Acts::TrackingVolume>;
0066
0067 /// Write a single surface.
0068 void fillSurfaceData(SurfaceData& data, const Acts::Surface& surface,
0069                      const Acts::GeometryContext& geoCtx) noexcept(false) {
0070   // encoded and partially decoded geometry identifier
0071   data.geometry_id = surface.geometryId().value();
0072   data.volume_id = surface.geometryId().volume();
0073   data.boundary_id = surface.geometryId().boundary();
0074   data.layer_id = surface.geometryId().layer();
0075   data.module_id = surface.geometryId().sensitive();
0076   // center position
0077   auto center = surface.center(geoCtx);
0078   data.cx = center.x() / Acts::UnitConstants::mm;
0079   data.cy = center.y() / Acts::UnitConstants::mm;
0080   data.cz = center.z() / Acts::UnitConstants::mm;
0081   // rotation matrix components are unit-less
0082   auto transform = surface.transform(geoCtx);
0083   data.rot_xu = transform(0, 0);
0084   data.rot_xv = transform(0, 1);
0085   data.rot_xw = transform(0, 2);
0086   data.rot_yu = transform(1, 0);
0087   data.rot_yv = transform(1, 1);
0088   data.rot_yw = transform(1, 2);
0089   data.rot_zu = transform(2, 0);
0090   data.rot_zv = transform(2, 1);
0091   data.rot_zw = transform(2, 2);
0092
0093   std::array<float*, 7> dataBoundParameters = {
0094       &data.bound_param0, &data.bound_param1, &data.bound_param2,
0095       &data.bound_param3, &data.bound_param4, &data.bound_param5,
0096       &data.bound_param6};
0097
0098   const auto& bounds = surface.bounds();
0099   data.bounds_type = static_cast<int>(bounds.type());
0100   auto boundValues = bounds.values();
0101
0102   if (boundValues.size() > dataBoundParameters.size()) {
0103     throw std::invalid_argument(
0104         "Bound types with too many parameters. Should never happen.");
0105   }
0106
0107   for (std::size_t ipar = 0; ipar < boundValues.size(); ++ipar) {
0108     (*dataBoundParameters[ipar]) = boundValues[ipar];
0109   }
0110
0111   if (surface.associatedDetectorElement() != nullptr) {
0112     data.module_t = surface.associatedDetectorElement()->thickness() /
0113                     Acts::UnitConstants::mm;
0114   }
0115 }
0116
0117 /// Write a single surface.
0118 void writeSurface(SurfaceWriter& sfWriter, const Acts::Surface& surface,
0119                   const Acts::GeometryContext& geoCtx) {
0120   SurfaceData data;
0121   fillSurfaceData(data, surface, geoCtx);
0122   sfWriter.append(data);
0123 }
0124
0125 /// Helper method for layer volume writing
0126 ///
0127 /// @param lv the layer volume to be written
0128 /// @param transform the layer transform
0129 /// @param representingBoundValues [in,out] the bound values
0130 /// @param last is the last layer
0131 void writeCylinderLayerVolume(LayerVolumeWriter& lvWriter,
0132                               const Acts::Layer& lv,
0133                               const Acts::Transform3& transform,
0134                               std::vector<double>& representingBoundValues,
0135                               std::vector<double>& volumeBoundValues,
0136                               std::vector<double>& lastBoundValues, bool last) {
0137   // The layer volume to be written
0138   LayerVolumeData lvDims;
0139   lvDims.geometry_id = lv.geometryId().value();
0140   lvDims.volume_id = lv.geometryId().volume();
0141   lvDims.layer_id = lv.geometryId().layer();
0142   bool isCylinderLayer = (lv.surfaceRepresentation().bounds().type() ==
0143                           Acts::SurfaceBounds::eCylinder);
0144
0145   auto lTranslation = transform.translation();
0146   // Change volume Bound values to r min, r max, z min, z max, phi min,
0147   // phi max
0148   representingBoundValues = {
0149       representingBoundValues[0],
0150       representingBoundValues[1],
0151       lTranslation.z() - representingBoundValues[2],
0152       lTranslation.z() + representingBoundValues[2],
0153       representingBoundValues[4] - representingBoundValues[3],
0154       representingBoundValues[4] + representingBoundValues[3]};
0155
0156   // Synchronize
0157   lvDims.min_v0 =
0158       isCylinderLayer ? representingBoundValues[0] : volumeBoundValues[0];
0159   lvDims.max_v0 =
0160       isCylinderLayer ? representingBoundValues[1] : volumeBoundValues[1];
0161   lvDims.min_v1 =
0162       isCylinderLayer ? volumeBoundValues[2] : representingBoundValues[2];
0163   lvDims.max_v1 =
0164       isCylinderLayer ? volumeBoundValues[3] : representingBoundValues[3];
0165   lvDims.min_v2 = representingBoundValues[4];
0166   lvDims.max_v2 = representingBoundValues[5];
0167
0168   // Write the prior navigation layer
0169   LayerVolumeData nlvDims;
0170   nlvDims.geometry_id = lv.geometryId().value();
0171   nlvDims.volume_id = lv.geometryId().volume();
0172   nlvDims.layer_id = lv.geometryId().layer() - 1;
0173   if (isCylinderLayer) {
0174     nlvDims.min_v0 = lastBoundValues[0];
0175     nlvDims.max_v0 = representingBoundValues[0];
0176     nlvDims.min_v1 = lastBoundValues[2];
0177     nlvDims.max_v1 = lastBoundValues[3];
0178     // Reset the r min boundary
0179     lastBoundValues[0] = representingBoundValues[1];
0180   } else {
0181     nlvDims.min_v0 = lastBoundValues[0];
0182     nlvDims.max_v0 = lastBoundValues[1];
0183     nlvDims.min_v1 = lastBoundValues[2];
0184     nlvDims.max_v1 = representingBoundValues[2];
0185     // Reset the r min boundary
0186     lastBoundValues[2] = representingBoundValues[3];
0187   }
0188   nlvDims.min_v2 = representingBoundValues[4];
0189   nlvDims.max_v2 = representingBoundValues[5];
0190   lvWriter.append(nlvDims);
0191   // Write the volume dimensions for the sensitive layer
0192   lvWriter.append(lvDims);
0193
0194   // Write last if needed
0195   if (last) {
0196     // Write the last navigation layer volume
0197     LayerVolumeData llvDims;
0198     llvDims.geometry_id = lv.geometryId().value();
0199     llvDims.volume_id = lv.geometryId().volume();
0200     llvDims.layer_id = lv.geometryId().layer() + 1;
0201     llvDims.min_v0 = lastBoundValues[0];
0202     llvDims.max_v0 = lastBoundValues[1];
0203     llvDims.min_v1 = lastBoundValues[2];
0204     llvDims.max_v1 = lastBoundValues[3];
0205     llvDims.min_v2 = representingBoundValues[4];
0206     llvDims.max_v2 = representingBoundValues[5];
0207     // Close up volume
0208     lvWriter.append(llvDims);
0209   }
0210 }
0211
0212 /// Write a single surface.
0213 void writeBoundarySurface(SurfaceWriter& writer,
0214                           const BoundarySurface& bsurface,
0215                           const Acts::GeometryContext& geoCtx) {
0216   SurfaceData data;
0217   fillSurfaceData(data, bsurface.surfaceRepresentation(), geoCtx);
0218   writer.append(data);
0219 }
0220
0221 /// Write all child surfaces and descend into confined volumes.
0222 void writeVolume(SurfaceWriter& sfWriter, SurfaceGridWriter& sfGridWriter,
0223                  LayerVolumeWriter& lvWriter,
0224                  const Acts::TrackingVolume& volume, bool writeSensitive,
0225                  bool writeBoundary, bool writeSurfaceGrid,
0226                  bool writeLayerVolume, const Acts::GeometryContext& geoCtx) {
0227   // process all layers that are directly stored within this volume
0228   if (volume.confinedLayers() != nullptr) {
0229     const auto& vTransform = volume.transform();
0230
0231     // Get the values of the volume boundaries
0232     std::vector<double> volumeBoundValues = volume.volumeBounds().values();
0233     std::vector<double> lastBoundValues;
0234
0235     if (volume.volumeBounds().type() == Acts::VolumeBounds::eCylinder) {
0236       auto vTranslation = vTransform.translation();
0237       // values to r min, r max, z min, z max, phi min, phi max
0238       volumeBoundValues = {
0239           volumeBoundValues[0],
0240           volumeBoundValues[1],
0241           vTranslation.z() - volumeBoundValues[2],
0242           vTranslation.z() + volumeBoundValues[2],
0243           volumeBoundValues[4] - volumeBoundValues[3],
0244           volumeBoundValues[4] + volumeBoundValues[3],
0245       };
0246       lastBoundValues = volumeBoundValues;
0247     }
0248
0249     unsigned int layerIdx = 0;
0250     const auto& layers = volume.confinedLayers()->arrayObjects();
0251
0252     // If we only have three layers, then the volume is the layer volume
0253     // so let's write it - this case will be excluded afterwards
0254     if (layers.size() == 3 && writeLayerVolume) {
0255       auto slayer = layers[1];
0256       LayerVolumeData plvDims;
0257       plvDims.geometry_id = slayer->geometryId().value();
0258       plvDims.volume_id = slayer->geometryId().volume();
0259       plvDims.layer_id = slayer->geometryId().layer();
0260       plvDims.min_v0 = volumeBoundValues[0];
0261       plvDims.max_v0 = volumeBoundValues[1];
0262       plvDims.min_v1 = volumeBoundValues[2];
0263       plvDims.max_v1 = volumeBoundValues[3];
0264       lvWriter.append(plvDims);
0265     }
0266
0267     // Now loop over the layer and write them
0268     for (const auto& layer : layers) {
0269       // We skip over navigation layers for layer volume writing
0270       // they will be written with the sensitive/passive parts for
0271       // synchronization
0272       if (layer->layerType() == Acts::navigation) {
0273         ++layerIdx;
0274         // For a correct layer volume setup, we need the navigation layers
0275         if (writeLayerVolume) {
0276           writeSurface(sfWriter, layer->surfaceRepresentation(), geoCtx);
0277         }
0278         continue;
0279
0280       } else {
0281         // Get the representing volume
0282         const auto* rVolume = layer->representingVolume();
0283
0284         // Write the layer volume, exclude single layer volumes (written above)
0285         if (rVolume != nullptr && writeLayerVolume && layers.size() > 3) {
0286           // Get the values of the representing volume
0287           std::vector<double> representingBoundValues =
0288               rVolume->volumeBounds().values();
0289           if (rVolume->volumeBounds().type() == Acts::VolumeBounds::eCylinder) {
0290             bool last = (layerIdx + 2 ==
0291                          volume.confinedLayers()->arrayObjects().size());
0292             writeCylinderLayerVolume(lvWriter, *layer, rVolume->transform(),
0293                                      representingBoundValues, volumeBoundValues,
0294                                      lastBoundValues, last);
0295           }
0296         }
0297
0298         // Surface has sub surfaces
0299         if (layer->surfaceArray() != nullptr) {
0300           auto sfArray = layer->surfaceArray();
0301
0302           // Write the surface grid itself if configured
0303           if (writeSurfaceGrid) {
0304             SurfaceGridData sfGrid;
0305             sfGrid.geometry_id = layer->geometryId().value();
0306             sfGrid.volume_id = layer->geometryId().volume();
0307             sfGrid.layer_id = layer->geometryId().layer();
0308
0309             // Draw the grid itself
0310             auto binning = sfArray->binningValues();
0311             auto axes = sfArray->getAxes();
0312             if (!binning.empty() && binning.size() == 2 && axes.size() == 2) {
0313               auto loc0Values = axes[0]->getBinEdges();
0314               sfGrid.nbins_loc0 = loc0Values.size() - 1;
0315               sfGrid.type_loc0 = static_cast<int>(binning[0]);
0316               sfGrid.min_loc0 = loc0Values[0];
0317               sfGrid.max_loc0 = loc0Values[loc0Values.size() - 1];
0318
0319               auto loc1Values = axes[1]->getBinEdges();
0320               sfGrid.nbins_loc1 = loc1Values.size() - 1;
0321               sfGrid.type_loc1 = static_cast<int>(binning[1]);
0322               sfGrid.min_loc1 = loc1Values[0];
0323               sfGrid.max_loc1 = loc1Values[loc1Values.size() - 1];
0324             }
0325             sfGridWriter.append(sfGrid);
0326           }
0327
0328           // Write the sensitive surface if configured
0329           if (writeSensitive) {
0330             for (auto surface : sfArray->surfaces()) {
0331               if (surface != nullptr) {
0332                 writeSurface(sfWriter, *surface, geoCtx);
0333               }
0334             }
0335           }
0336         } else {
0337           // Write the passive surface
0338           writeSurface(sfWriter, layer->surfaceRepresentation(), geoCtx);
0339         }
0340       }
0341       ++layerIdx;
0342     }  // end of layer loop
0343
0344     // This is a navigation volume, write the boundaries
0345     if (writeBoundary) {
0346       for (const auto& bsurface : volume.boundarySurfaces()) {
0347         writeBoundarySurface(sfWriter, *bsurface, geoCtx);
0348       }
0349     }
0350   }
0351   // step down into hierarchy to process all child volumnes
0352   if (volume.confinedVolumes()) {
0353     for (const auto& confined : volume.confinedVolumes()->arrayObjects()) {
0354       writeVolume(sfWriter, sfGridWriter, lvWriter, *confined.get(),
0355                   writeSensitive, writeBoundary, writeSurfaceGrid,
0356                   writeLayerVolume, geoCtx);
0357     }
0358   }
0359 }
0360 }  // namespace
0361
0362 ProcessCode CsvTrackingGeometryWriter::write(const AlgorithmContext& ctx) {
0363   if (!m_cfg.writePerEvent) {
0364     return ProcessCode::SUCCESS;
0365   }
0366
0367   SurfaceWriter sfWriter(
0368       perEventFilepath(m_cfg.outputDir, "detectors.csv", ctx.eventNumber),
0369       m_cfg.outputPrecision);
0370
0371   SurfaceGridWriter sfGridWriter(
0372       perEventFilepath(m_cfg.outputDir, "surface-grids.csv", ctx.eventNumber),
0373       m_cfg.outputPrecision);
0374
0375   LayerVolumeWriter lvWriter(
0376       perEventFilepath(m_cfg.outputDir, "layer-volumes.csv", ctx.eventNumber),
0377       m_cfg.outputPrecision);
0378
0379   writeVolume(sfWriter, sfGridWriter, lvWriter, *m_world, m_cfg.writeSensitive,
0380               m_cfg.writeBoundary, m_cfg.writeSurfaceGrid,
0381               m_cfg.writeLayerVolume, ctx.geoContext);
0382   return ProcessCode::SUCCESS;
0383 }
0384
0385 ProcessCode CsvTrackingGeometryWriter::finalize() {
0386   SurfaceWriter sfWriter(joinPaths(m_cfg.outputDir, "detectors.csv"),
0387                          m_cfg.outputPrecision);
0388   SurfaceGridWriter sfGridWriter(
0389       joinPaths(m_cfg.outputDir, "surface-grids.csv"), m_cfg.outputPrecision);
0390
0391   LayerVolumeWriter lvWriter(joinPaths(m_cfg.outputDir, "layer-volumes.csv"),
0392                              m_cfg.outputPrecision);
0393
0394   writeVolume(sfWriter, sfGridWriter, lvWriter, *m_world, m_cfg.writeSensitive,
0395               m_cfg.writeBoundary, m_cfg.writeSurfaceGrid,
0396               m_cfg.writeLayerVolume, Acts::GeometryContext());
0397   return ProcessCode::SUCCESS;
0398 }