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   data.extra_id = surface.geometryId().extra();
0077   // center position
0078   auto center = surface.center(geoCtx);
0079   data.cx = center.x() / Acts::UnitConstants::mm;
0080   data.cy = center.y() / Acts::UnitConstants::mm;
0081   data.cz = center.z() / Acts::UnitConstants::mm;
0082   // rotation matrix components are unit-less
0083   auto transform = surface.transform(geoCtx);
0084   data.rot_xu = transform(0, 0);
0085   data.rot_xv = transform(0, 1);
0086   data.rot_xw = transform(0, 2);
0087   data.rot_yu = transform(1, 0);
0088   data.rot_yv = transform(1, 1);
0089   data.rot_yw = transform(1, 2);
0090   data.rot_zu = transform(2, 0);
0091   data.rot_zv = transform(2, 1);
0092   data.rot_zw = transform(2, 2);
0093
0094   std::array<float*, 7> dataBoundParameters = {
0095       &data.bound_param0, &data.bound_param1, &data.bound_param2,
0096       &data.bound_param3, &data.bound_param4, &data.bound_param5,
0097       &data.bound_param6};
0098
0099   const auto& bounds = surface.bounds();
0100   data.bounds_type = static_cast<int>(bounds.type());
0101   auto boundValues = bounds.values();
0102
0103   if (boundValues.size() > dataBoundParameters.size()) {
0104     throw std::invalid_argument(
0105         "Bound types with too many parameters. Should never happen.");
0106   }
0107
0108   for (std::size_t ipar = 0; ipar < boundValues.size(); ++ipar) {
0109     (*dataBoundParameters[ipar]) = boundValues[ipar];
0110   }
0111
0112   if (surface.associatedDetectorElement() != nullptr) {
0113     data.module_t = surface.associatedDetectorElement()->thickness() /
0114                     Acts::UnitConstants::mm;
0115   }
0116 }
0117
0118 /// Write a single surface.
0119 void writeSurface(SurfaceWriter& sfWriter, const Acts::Surface& surface,
0120                   const Acts::GeometryContext& geoCtx) {
0121   SurfaceData data;
0122   fillSurfaceData(data, surface, geoCtx);
0123   sfWriter.append(data);
0124 }
0125
0126 /// Helper method for layer volume writing
0127 ///
0128 /// @param lv the layer volume to be written
0129 /// @param transform the layer transform
0130 /// @param representingBoundValues [in,out] the bound values
0131 /// @param last is the last layer
0132 void writeCylinderLayerVolume(LayerVolumeWriter& lvWriter,
0133                               const Acts::Layer& lv,
0134                               const Acts::Transform3& transform,
0135                               std::vector<double>& representingBoundValues,
0136                               std::vector<double>& volumeBoundValues,
0137                               std::vector<double>& lastBoundValues, bool last) {
0138   // The layer volume to be written
0139   LayerVolumeData lvDims;
0140   lvDims.geometry_id = lv.geometryId().value();
0141   lvDims.volume_id = lv.geometryId().volume();
0142   lvDims.layer_id = lv.geometryId().layer();
0143   bool isCylinderLayer = (lv.surfaceRepresentation().bounds().type() ==
0144                           Acts::SurfaceBounds::eCylinder);
0145
0146   auto lTranslation = transform.translation();
0147   // Change volume Bound values to r min, r max, z min, z max, phi min,
0148   // phi max
0149   representingBoundValues = {
0150       representingBoundValues[0],
0151       representingBoundValues[1],
0152       lTranslation.z() - representingBoundValues[2],
0153       lTranslation.z() + representingBoundValues[2],
0154       representingBoundValues[4] - representingBoundValues[3],
0155       representingBoundValues[4] + representingBoundValues[3]};
0156
0157   // Synchronize
0158   lvDims.min_v0 =
0159       isCylinderLayer ? representingBoundValues[0] : volumeBoundValues[0];
0160   lvDims.max_v0 =
0161       isCylinderLayer ? representingBoundValues[1] : volumeBoundValues[1];
0162   lvDims.min_v1 =
0163       isCylinderLayer ? volumeBoundValues[2] : representingBoundValues[2];
0164   lvDims.max_v1 =
0165       isCylinderLayer ? volumeBoundValues[3] : representingBoundValues[3];
0166   lvDims.min_v2 = representingBoundValues[4];
0167   lvDims.max_v2 = representingBoundValues[5];
0168
0169   // Write the prior navigation layer
0170   LayerVolumeData nlvDims;
0171   nlvDims.geometry_id = lv.geometryId().value();
0172   nlvDims.volume_id = lv.geometryId().volume();
0173   nlvDims.layer_id = lv.geometryId().layer() - 1;
0174   if (isCylinderLayer) {
0175     nlvDims.min_v0 = lastBoundValues[0];
0176     nlvDims.max_v0 = representingBoundValues[0];
0177     nlvDims.min_v1 = lastBoundValues[2];
0178     nlvDims.max_v1 = lastBoundValues[3];
0179     // Reset the r min boundary
0180     lastBoundValues[0] = representingBoundValues[1];
0181   } else {
0182     nlvDims.min_v0 = lastBoundValues[0];
0183     nlvDims.max_v0 = lastBoundValues[1];
0184     nlvDims.min_v1 = lastBoundValues[2];
0185     nlvDims.max_v1 = representingBoundValues[2];
0186     // Reset the r min boundary
0187     lastBoundValues[2] = representingBoundValues[3];
0188   }
0189   nlvDims.min_v2 = representingBoundValues[4];
0190   nlvDims.max_v2 = representingBoundValues[5];
0191   lvWriter.append(nlvDims);
0192   // Write the volume dimensions for the sensitive layer
0193   lvWriter.append(lvDims);
0194
0195   // Write last if needed
0196   if (last) {
0197     // Write the last navigation layer volume
0198     LayerVolumeData llvDims;
0199     llvDims.geometry_id = lv.geometryId().value();
0200     llvDims.volume_id = lv.geometryId().volume();
0201     llvDims.layer_id = lv.geometryId().layer() + 1;
0202     llvDims.min_v0 = lastBoundValues[0];
0203     llvDims.max_v0 = lastBoundValues[1];
0204     llvDims.min_v1 = lastBoundValues[2];
0205     llvDims.max_v1 = lastBoundValues[3];
0206     llvDims.min_v2 = representingBoundValues[4];
0207     llvDims.max_v2 = representingBoundValues[5];
0208     // Close up volume
0209     lvWriter.append(llvDims);
0210   }
0211 }
0212
0213 /// Write a single surface.
0214 void writeBoundarySurface(SurfaceWriter& writer,
0215                           const BoundarySurface& bsurface,
0216                           const Acts::GeometryContext& geoCtx) {
0217   SurfaceData data;
0218   fillSurfaceData(data, bsurface.surfaceRepresentation(), geoCtx);
0219   writer.append(data);
0220 }
0221
0222 /// Write all child surfaces and descend into confined volumes.
0223 void writeVolume(SurfaceWriter& sfWriter, SurfaceGridWriter& sfGridWriter,
0224                  LayerVolumeWriter& lvWriter,
0225                  const Acts::TrackingVolume& volume, bool writeSensitive,
0226                  bool writeBoundary, bool writeSurfaceGrid,
0227                  bool writeLayerVolume, const Acts::GeometryContext& geoCtx) {
0228   // process all layers that are directly stored within this volume
0229   if (volume.confinedLayers() != nullptr) {
0230     const auto& vTransform = volume.transform();
0231
0232     // Get the values of the volume boundaries
0233     std::vector<double> volumeBoundValues = volume.volumeBounds().values();
0234     std::vector<double> lastBoundValues;
0235
0236     if (volume.volumeBounds().type() == Acts::VolumeBounds::eCylinder) {
0237       auto vTranslation = vTransform.translation();
0238       // values to r min, r max, z min, z max, phi min, phi max
0239       volumeBoundValues = {
0240           volumeBoundValues[0],
0241           volumeBoundValues[1],
0242           vTranslation.z() - volumeBoundValues[2],
0243           vTranslation.z() + volumeBoundValues[2],
0244           volumeBoundValues[4] - volumeBoundValues[3],
0245           volumeBoundValues[4] + volumeBoundValues[3],
0246       };
0247       lastBoundValues = volumeBoundValues;
0248     }
0249
0250     unsigned int layerIdx = 0;
0251     const auto& layers = volume.confinedLayers()->arrayObjects();
0252
0253     // If we only have three layers, then the volume is the layer volume
0254     // so let's write it - this case will be excluded afterwards
0255     if (layers.size() == 3 && writeLayerVolume) {
0256       auto slayer = layers[1];
0257       LayerVolumeData plvDims;
0258       plvDims.geometry_id = slayer->geometryId().value();
0259       plvDims.volume_id = slayer->geometryId().volume();
0260       plvDims.layer_id = slayer->geometryId().layer();
0261       plvDims.min_v0 = volumeBoundValues[0];
0262       plvDims.max_v0 = volumeBoundValues[1];
0263       plvDims.min_v1 = volumeBoundValues[2];
0264       plvDims.max_v1 = volumeBoundValues[3];
0265       lvWriter.append(plvDims);
0266     }
0267
0268     // Now loop over the layer and write them
0269     for (const auto& layer : layers) {
0270       // We skip over navigation layers for layer volume writing
0271       // they will be written with the sensitive/passive parts for
0272       // synchronization
0273       if (layer->layerType() == Acts::navigation) {
0274         ++layerIdx;
0275         // For a correct layer volume setup, we need the navigation layers
0276         if (writeLayerVolume) {
0277           writeSurface(sfWriter, layer->surfaceRepresentation(), geoCtx);
0278         }
0279         continue;
0280
0281       } else {
0282         // Get the representing volume
0283         const auto* rVolume = layer->representingVolume();
0284
0285         // Write the layer volume, exclude single layer volumes (written above)
0286         if (rVolume != nullptr && writeLayerVolume && layers.size() > 3) {
0287           // Get the values of the representing volume
0288           std::vector<double> representingBoundValues =
0289               rVolume->volumeBounds().values();
0290           if (rVolume->volumeBounds().type() == Acts::VolumeBounds::eCylinder) {
0291             bool last = (layerIdx + 2 ==
0292                          volume.confinedLayers()->arrayObjects().size());
0293             writeCylinderLayerVolume(lvWriter, *layer, rVolume->transform(),
0294                                      representingBoundValues, volumeBoundValues,
0295                                      lastBoundValues, last);
0296           }
0297         }
0298
0299         // Surface has sub surfaces
0300         if (layer->surfaceArray() != nullptr) {
0301           auto sfArray = layer->surfaceArray();
0302
0303           // Write the surface grid itself if configured
0304           if (writeSurfaceGrid) {
0305             SurfaceGridData sfGrid;
0306             sfGrid.geometry_id = layer->geometryId().value();
0307             sfGrid.volume_id = layer->geometryId().volume();
0308             sfGrid.layer_id = layer->geometryId().layer();
0309
0310             // Draw the grid itself
0311             auto binning = sfArray->binningValues();
0312             auto axes = sfArray->getAxes();
0313             if (!binning.empty() && binning.size() == 2 && axes.size() == 2) {
0314               auto loc0Values = axes[0]->getBinEdges();
0315               sfGrid.nbins_loc0 = loc0Values.size() - 1;
0316               sfGrid.type_loc0 = static_cast<int>(binning[0]);
0317               sfGrid.min_loc0 = loc0Values[0];
0318               sfGrid.max_loc0 = loc0Values[loc0Values.size() - 1];
0319
0320               auto loc1Values = axes[1]->getBinEdges();
0321               sfGrid.nbins_loc1 = loc1Values.size() - 1;
0322               sfGrid.type_loc1 = static_cast<int>(binning[1]);
0323               sfGrid.min_loc1 = loc1Values[0];
0324               sfGrid.max_loc1 = loc1Values[loc1Values.size() - 1];
0325             }
0326             sfGridWriter.append(sfGrid);
0327           }
0328
0329           // Write the sensitive surface if configured
0330           if (writeSensitive) {
0331             for (auto surface : sfArray->surfaces()) {
0332               if (surface != nullptr) {
0333                 writeSurface(sfWriter, *surface, geoCtx);
0334               }
0335             }
0336           }
0337         } else {
0338           // Write the passive surface
0339           writeSurface(sfWriter, layer->surfaceRepresentation(), geoCtx);
0340         }
0341       }
0342       ++layerIdx;
0343     }  // end of layer loop
0344
0345     // This is a navigation volume, write the boundaries
0346     if (writeBoundary) {
0347       for (const auto& bsurface : volume.boundarySurfaces()) {
0348         writeBoundarySurface(sfWriter, *bsurface, geoCtx);
0349       }
0350     }
0351   }
0352   // step down into hierarchy to process all child volumnes
0353   if (volume.confinedVolumes()) {
0354     for (const auto& confined : volume.confinedVolumes()->arrayObjects()) {
0355       writeVolume(sfWriter, sfGridWriter, lvWriter, *confined, writeSensitive,
0356                   writeBoundary, writeSurfaceGrid, 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 }