feat: Improvements to Cuboid Volume Builder (#1166)

Adds a number of changes to the CuboidVolumeBuilder to make it more robust.

Core/include/Acts/Geometry/CuboidVolumeBuilder.hpp

@@ -16,6 +16,7 @@
 #include <functional>
 #include <iosfwd>
 #include <memory>
+#include <optional>
 #include <string>
 #include <utility>
 #include <vector>
@@ -28,7 +29,7 @@ class RectangleBounds;
 class ISurfaceMaterial;
 class IVolumeMaterial;
 class DetectorElementBase;
-class PlaneSurface;
+class Surface;
 class Layer;
 
 /// @brief This class builds a box detector with a configurable amount of
@@ -59,18 +60,20 @@ class CuboidVolumeBuilder : public ITrackingVolumeBuilder {
   };
 
   /// @brief This struct stores the data for the construction of a PlaneLayer
-  /// that has a single PlaneSurface encapsulated
   struct LayerConfig {
     // Configuration of the surface
-    SurfaceConfig surfaceCfg;
+    std::vector<SurfaceConfig> surfaceCfg;
     // Encapsulated surface
-    std::shared_ptr<const PlaneSurface> surface = nullptr;
+    std::vector<std::shared_ptr<const Surface>> surfaces;
     // Boolean flag if layer is active
     bool active = false;
     // Bins in Y direction
     size_t binsY = 1;
     // Bins in Z direction
     size_t binsZ = 1;
+    // Envelope in X (along layer normal)
+    std::pair<double, double> envelopeX{0, 0};
+    std::optional<RotationMatrix3> rotation{std::nullopt};
   };
 
   /// @brief This struct stores the data for the construction of a cuboid
@@ -124,8 +127,8 @@ class CuboidVolumeBuilder : public ITrackingVolumeBuilder {
   /// @param [in] cfg Configuration of the surface
   ///
   /// @return Pointer to the created surface
-  std::shared_ptr<const PlaneSurface> buildSurface(
-      const GeometryContext& gctx, const SurfaceConfig& cfg) const;
+  std::shared_ptr<const Surface> buildSurface(const GeometryContext& gctx,
+                                              const SurfaceConfig& cfg) const;
 
   /// @brief This function creates a layer with a surface encaspulated with a
   /// given configuration. The surface gets a detector element attached if the

Core/src/Geometry/CuboidVolumeBuilder.cpp

@@ -33,10 +33,8 @@
 #include "Acts/Utilities/Logger.hpp"
 
 #include <limits>
-#include <optional>
 
-std::shared_ptr<const Acts::PlaneSurface>
-Acts::CuboidVolumeBuilder::buildSurface(
+std::shared_ptr<const Acts::Surface> Acts::CuboidVolumeBuilder::buildSurface(
     const GeometryContext& /*gctx*/,
     const CuboidVolumeBuilder::SurfaceConfig& cfg) const {
   std::shared_ptr<PlaneSurface> surface;
@@ -60,33 +58,73 @@ Acts::CuboidVolumeBuilder::buildSurface(
 std::shared_ptr<const Acts::Layer> Acts::CuboidVolumeBuilder::buildLayer(
     const GeometryContext& gctx,
     Acts::CuboidVolumeBuilder::LayerConfig& cfg) const {
+  if (cfg.surfaces.empty() && cfg.surfaceCfg.empty()) {
+    throw std::runtime_error{
+        "Neither surfaces nor config to build surfaces was provided. Cannot "
+        "proceed"};
+  }
+
   // Build the surface
-  if (cfg.surface == nullptr) {
-    cfg.surface = buildSurface(gctx, cfg.surfaceCfg);
+  if (cfg.surfaces.empty()) {
+    for (const auto& sCfg : cfg.surfaceCfg) {
+      cfg.surfaces.push_back(buildSurface(gctx, sCfg));
+    }
   }
   // Build transformation centered at the surface position
-  Transform3 trafo(Transform3::Identity() * cfg.surfaceCfg.rotation);
-  trafo.translation() = cfg.surfaceCfg.position;
+  Vector3 centroid{0., 0., 0.};
+
+  for (const auto& surface : cfg.surfaces) {
+    centroid += surface->transform(gctx).translation();
+  }
+
+  centroid /= cfg.surfaces.size();
+
+  // In the case the layer configuration doesn't define the rotation of the
+  // layer use the orientation of the first surface to define the layer rotation
+  // in space.
+  Transform3 trafo = Transform3::Identity();
+  trafo.translation() = centroid;
+  if (cfg.rotation) {
+    trafo.linear() = *cfg.rotation;
+  } else {
+    trafo.linear() = cfg.surfaces.front()->transform(gctx).rotation();
+  }
 
   LayerCreator::Config lCfg;
   lCfg.surfaceArrayCreator = std::make_shared<const SurfaceArrayCreator>();
   LayerCreator layerCreator(lCfg);
-  return layerCreator.planeLayer(gctx, {cfg.surface}, cfg.binsY, cfg.binsZ,
-                                 BinningValue::binX, std::nullopt, trafo);
+  ProtoLayer pl{gctx, cfg.surfaces};
+  pl.envelope[binX] = cfg.envelopeX;
+  return layerCreator.planeLayer(gctx, cfg.surfaces, cfg.binsY, cfg.binsZ,
+                                 BinningValue::binX, pl, trafo);
 }
 
 std::pair<double, double> Acts::CuboidVolumeBuilder::binningRange(
-    const GeometryContext& /*gctx*/,
+    const GeometryContext& gctx,
     const Acts::CuboidVolumeBuilder::VolumeConfig& cfg) const {
   using namespace UnitLiterals;
   // Construct return value
   std::pair<double, double> minMax = std::make_pair(
       std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());
+
+  // Compute the min volume boundaries for computing the binning start
+  // See
+  // https://acts.readthedocs.io/en/latest/core/geometry.html#geometry-building
+  // !! IMPORTANT !! The volume is assumed to be already rotated into the
+  // telescope geometry
+  Vector3 minVolumeBoundaries = cfg.position - 0.5 * cfg.length;
+  Vector3 maxVolumeBoundaries = cfg.position + 0.5 * cfg.length;
+
+  // Compute first the min-max from the layers
+
   for (const auto& layercfg : cfg.layerCfg) {
-    auto surfacePosMin = layercfg.surfaceCfg.position.x() -
-                         layercfg.surfaceCfg.thickness / 2. - 1._um;
-    auto surfacePosMax = layercfg.surfaceCfg.position.x() +
-                         layercfg.surfaceCfg.thickness / 2. + 1._um;
+    // recreating the protolayer for each layer => slow, but only few sensors
+    ProtoLayer pl{gctx, layercfg.surfaces};
+    pl.envelope[binX] = layercfg.envelopeX;
+
+    double surfacePosMin = pl.min(binX);
+    double surfacePosMax = pl.max(binX);
+
     // Test if new extreme is found and set it
     if (surfacePosMin < minMax.first) {
       minMax.first = surfacePosMin;
@@ -95,6 +133,11 @@ std::pair<double, double> Acts::CuboidVolumeBuilder::binningRange(
       minMax.second = surfacePosMax;
     }
   }
+
+  // Use the volume boundaries as limits for the binning
+  minMax.first = std::min(minMax.first, minVolumeBoundaries(binX));
+  minMax.second = std::max(minMax.second, maxVolumeBoundaries(binX));
+
   return minMax;
 }
 
@@ -124,7 +167,7 @@ std::shared_ptr<Acts::TrackingVolume> Acts::CuboidVolumeBuilder::buildVolume(
         RectangleBounds(cfg.length.y() * 0.5, cfg.length.z() * 0.5));
 
     LayerConfig lCfg;
-    lCfg.surfaceCfg = sCfg;
+    lCfg.surfaceCfg = {sCfg};
 
     cfg.layerCfg.push_back(lCfg);
   }
@@ -183,6 +226,13 @@ Acts::MutableTrackingVolumePtr Acts::CuboidVolumeBuilder::trackingVolume(
     volumes.push_back(buildVolume(gctx, volCfg));
   }
 
+  // Sort the volumes vectors according to the center location, otherwise the
+  // binning boundaries will fail
+  std::sort(volumes.begin(), volumes.end(),
+            [](const TrackingVolumePtr& lhs, const TrackingVolumePtr& rhs) {
+              return lhs->center().x() < rhs->center().x();
+            });
+
   // Glue volumes
   for (unsigned int i = 0; i < volumes.size() - 1; i++) {
     volumes[i + 1]->glueTrackingVolume(

Tests/UnitTests/Core/Geometry/CuboidVolumeBuilderTests.cpp

@@ -102,7 +102,7 @@ BOOST_AUTO_TEST_CASE(CuboidVolumeBuilderTest) {
 
   // Test that 4 surfaces can be built
   for (const auto& cfg : surfaceConfig) {
-    std::shared_ptr<const PlaneSurface> pSur = cvb.buildSurface(tgContext, cfg);
+    std::shared_ptr<const Surface> pSur = cvb.buildSurface(tgContext, cfg);
     BOOST_CHECK_NE(pSur, nullptr);
     CHECK_CLOSE_ABS(pSur->center(tgContext), cfg.position, 1e-9);
     BOOST_CHECK_NE(pSur->surfaceMaterial(), nullptr);
@@ -114,7 +114,7 @@ BOOST_AUTO_TEST_CASE(CuboidVolumeBuilderTest) {
   std::vector<CuboidVolumeBuilder::LayerConfig> layerConfig;
   for (auto& sCfg : surfaceConfig) {
     CuboidVolumeBuilder::LayerConfig cfg;
-    cfg.surfaceCfg = sCfg;
+    cfg.surfaceCfg = {sCfg};
     layerConfig.push_back(cfg);
   }
 
@@ -125,13 +125,13 @@ BOOST_AUTO_TEST_CASE(CuboidVolumeBuilderTest) {
   for (auto& cfg : layerConfig) {
     LayerPtr layer = cvb.buildLayer(tgContext, cfg);
     BOOST_CHECK_NE(layer, nullptr);
-    BOOST_CHECK_NE(cfg.surface, nullptr);
+    BOOST_CHECK(!cfg.surfaces.empty());
     BOOST_CHECK_EQUAL(layer->surfaceArray()->surfaces().size(), 1u);
     BOOST_CHECK_EQUAL(layer->layerType(), LayerType::active);
   }
 
   for (auto& cfg : layerConfig) {
-    cfg.surface = nullptr;
+    cfg.surfaces = {};
   }
 
   // Build volume configuration
@@ -164,7 +164,7 @@ BOOST_AUTO_TEST_CASE(CuboidVolumeBuilderTest) {
 
   volumeConfig.layers.clear();
   for (auto& lay : volumeConfig.layerCfg) {
-    lay.surface = nullptr;
+    lay.surfaces = {};
     lay.active = true;
   }
   trVol = cvb.buildVolume(tgContext, volumeConfig);
@@ -218,7 +218,7 @@ BOOST_AUTO_TEST_CASE(CuboidVolumeBuilderTest) {
   std::vector<CuboidVolumeBuilder::LayerConfig> layerConfig2;
   for (auto& sCfg : surfaceConfig2) {
     CuboidVolumeBuilder::LayerConfig cfg;
-    cfg.surfaceCfg = sCfg;
+    cfg.surfaceCfg = {sCfg};
     layerConfig2.push_back(cfg);
   }
   CuboidVolumeBuilder::VolumeConfig volumeConfig2;

Tests/UnitTests/Core/Propagator/StepperTests.cpp

@@ -983,7 +983,7 @@ BOOST_AUTO_TEST_CASE(step_extension_trackercalomdt_test) {
       new HomogeneousSurfaceMaterial(matProp));
   sConf1.thickness = 1._mm;
   CuboidVolumeBuilder::LayerConfig lConf1;
-  lConf1.surfaceCfg = sConf1;
+  lConf1.surfaceCfg = {sConf1};
 
   CuboidVolumeBuilder::SurfaceConfig sConf2;
   sConf2.position = Vector3(0.6_m, 0., 0.);
@@ -996,7 +996,7 @@ BOOST_AUTO_TEST_CASE(step_extension_trackercalomdt_test) {
       new HomogeneousSurfaceMaterial(matProp));
   sConf2.thickness = 1._mm;
   CuboidVolumeBuilder::LayerConfig lConf2;
-  lConf2.surfaceCfg = sConf2;
+  lConf2.surfaceCfg = {sConf2};
 
   CuboidVolumeBuilder::VolumeConfig muConf1;
   muConf1.position = {2.3_m, 0., 0.};

Tests/UnitTests/Core/Visualization/EventDataView3DBase.hpp

@@ -160,7 +160,7 @@ static inline std::string testMultiTrajectory(IVisualization3D& helper) {
           return new Test::DetectorElementStub(trans, bounds, thickness);
         };
     CuboidVolumeBuilder::LayerConfig lConf;
-    lConf.surfaceCfg = sConf;
+    lConf.surfaceCfg = {sConf};
     lConfs.push_back(lConf);
   }