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