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