feat: GeoModel conversion for TGC, sTGC and MicroMega

Plugins/GeoModel/CMakeLists.txt

@@ -13,6 +13,7 @@ add_library(
     src/detail/GeoBoxConverter.cpp
     src/detail/GeoTrdConverter.cpp
     src/detail/GeoTubeConverter.cpp
+    src/detail/GeoPolygonConverter.cpp
     src/detail/GeoShiftConverter.cpp
     src/detail/GeoIntersectionAnnulusConverter.cpp
     src/detail/GeoModelBinningHelper.cpp

Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelConverters.hpp

@@ -15,6 +15,7 @@
 #include "Acts/Plugins/GeoModel/detail/GenericGeoShapeConverter.hpp"
 #include "Acts/Plugins/GeoModel/detail/GeoBoxConverter.hpp"
 #include "Acts/Plugins/GeoModel/detail/GeoIntersectionAnnulusConverter.hpp"
+#include "Acts/Plugins/GeoModel/detail/GeoPolygonConverter.hpp"
 #include "Acts/Plugins/GeoModel/detail/GeoShiftConverter.hpp"
 #include "Acts/Plugins/GeoModel/detail/GeoTrdConverter.hpp"
 #include "Acts/Plugins/GeoModel/detail/GeoTubeConverter.hpp"
@@ -37,6 +38,9 @@ namespace Acts {
 using GeoBoxConverter =
     detail::GenericGeoShapeConverter<GeoBox, detail::GeoBoxConverter>;
 
+using GeoPolygonConverter =
+    detail::GenericGeoShapeConverter<GeoSimplePolygonBrep,
+                                     detail::GeoPolygonConverter>;
 /// @brief A dedicated converter for GeoInterseciton that describe annulus bounds
 ///
 /// This is very much tailored to the AnnulusBounds class
@@ -83,6 +87,8 @@ inline std::shared_ptr<const IGeoShapeConverter> geoShapesConverters(
            std::make_shared<GeoShiftConverter>()},
           {GeoTrd::getClassTypeID(), std::make_shared<GeoTrdConverter>()},
           {GeoTube::getClassTypeID(), std::make_shared<GeoTubeConverter>()},
+          {GeoSimplePolygonBrep::getClassTypeID(),
+           std::make_shared<GeoPolygonConverter>()},
           {GeoShapeUnion::getClassTypeID(),
            std::make_shared<GeoUnionDoubleTrdConverter>()}};
   auto itr = converters.find(geoShapeId);

Plugins/GeoModel/include/Acts/Plugins/GeoModel/detail/GeoPolygonConverter.hpp

@@ -0,0 +1,34 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Plugins/GeoModel/GeoModelDetectorElement.hpp"
+#include "Acts/Utilities/Result.hpp"
+
+#include <memory>
+#include <tuple>
+
+#include <GeoModelKernel/GeoSimplePolygonBrep.h>
+
+namespace Acts::detail {
+struct GeoPolygonConverter {
+  /// @brief Convert a GeoBox to a detector element and surface
+  ///
+  /// @param geoFPV The full physical volume to convert (contains shape)
+  /// @param geoBox The GeoBox to convert
+  /// @param absTransform from the GeoPhysVol
+  /// @param bool sensitive
+  ///
+  /// @return The detector element and surface
+  Result<GeoModelSensitiveSurface> operator()(
+      const PVConstLink& geoPV, const GeoSimplePolygonBrep& geoPolygon,
+      const Transform3& absTransform, bool sensitive) const;
+};
+}  // namespace Acts::detail

Plugins/GeoModel/src/GeoModelDetectorObjectFactory.cpp

@@ -139,13 +139,14 @@ Acts::GeoModelDetectorObjectFactory::findAllSubVolumes(const PVConstLink &vol) {
   }
   return sensitives;
 }
+
 bool Acts::GeoModelDetectorObjectFactory::convertBox(std::string name) {
   auto convB = std::any_of(
       m_cfg.convertBox.begin(), m_cfg.convertBox.end(),
       [&](const auto &n) { return name.find(n) != std::string::npos; });
-  // return (name.find(m_cfg.convertBox) != std::string::npos);
   return convB;
 }
+
 void Acts::GeoModelDetectorObjectFactory::convertFpv(
     const std::string &name, GeoFullPhysVol *fpv, Cache &cache,
     const GeometryContext &gctx) {

Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp

@@ -65,7 +65,6 @@ std::shared_ptr<Experimental::DetectorVolume> convertVolume(
         portalGenerator, context, name, newTransform, bounds,
         Experimental::tryAllPortalsAndSurfaces());
   } else if (shape.typeID() == GeoBox::getClassTypeID()) {
-    // TODO do the surfaces
     const GeoBox* box = dynamic_cast<const GeoBox*>(&shape);
     std::shared_ptr<CuboidVolumeBounds> bounds =
         std::make_shared<CuboidVolumeBounds>(box->getXHalfLength(),

Plugins/GeoModel/src/detail/GeoPolygonConverter.cpp

@@ -0,0 +1,125 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Plugins/GeoModel/detail/GeoPolygonConverter.hpp"
+
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Plugins/GeoModel/GeoModelConversionError.hpp"
+#include "Acts/Surfaces/DiamondBounds.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Surfaces/TrapezoidBounds.hpp"
+
+#include <GeoModelKernel/GeoBox.h>
+#include <GeoModelKernel/GeoFullPhysVol.h>
+#include <GeoModelKernel/GeoLogVol.h>
+#include <GeoModelKernel/GeoShape.h>
+#include <GeoModelKernel/Units.h>
+
+Acts::Result<Acts::GeoModelSensitiveSurface>
+Acts::detail::GeoPolygonConverter::operator()(
+    const PVConstLink& geoPV, const GeoSimplePolygonBrep& polygon,
+    const Transform3& absTransform, bool sensitive) const {
+  /// auto-calculate the unit length conversion
+  static constexpr ActsScalar unitLength =
+      Acts::UnitConstants::mm / GeoModelKernelUnits::millimeter;
+
+  // Create the surface transform
+  Transform3 transform = Transform3::Identity();
+  transform.translation() = unitLength * absTransform.translation();
+  auto rotation = absTransform.rotation();
+  // Get the half lengths
+  int nVertices = polygon.getNVertices();
+  std::vector<std::vector<double>> vertices;
+
+  for (int i = 0; i < nVertices; i++) {
+    vertices.push_back({polygon.getXVertex(i), polygon.getYVertex(i)});
+  }
+  // sort based on the y-coordinate
+  std::sort(vertices.begin(), vertices.end(),
+            [](const std::vector<double>& a, const std::vector<double>& b) {
+              return a[1] < b[1];
+            });
+  if (nVertices == 4) {
+    double hlxnegy = fabs(vertices[0][0] - vertices[1][0]) / 2;
+    double hlxposy = fabs(vertices[2][0] - vertices[3][0]) / 2;
+    double hly = fabs(vertices[0][1] - vertices[3][1]) / 2;
+    std::vector<ActsScalar> halfLengths = {hlxnegy, hlxposy, hly};
+
+    // Create the surface
+    Vector3 colX = rotation.col(0);
+    Vector3 colY = rotation.col(1);
+    Vector3 colZ = rotation.col(2);
+    rotation.col(0) = colX;
+    rotation.col(1) = colY;
+    rotation.col(2) = colZ;
+    transform.linear() = rotation;
+    // Create the surface bounds
+    ActsScalar halfXnegY = unitLength * halfLengths[0];
+    ActsScalar halfXposY = unitLength * halfLengths[1];
+    ActsScalar halfY = unitLength * halfLengths[2];
+    auto trapBounds =
+        std::make_shared<Acts::TrapezoidBounds>(halfXnegY, halfXposY, halfY);
+    if (!sensitive) {
+      auto surface = Surface::makeShared<PlaneSurface>(transform, trapBounds);
+      return std::make_tuple(nullptr, surface);
+    }
+    // Create the element and the surface
+    auto detectorElement =
+        GeoModelDetectorElement::createDetectorElement<PlaneSurface>(
+            geoPV, trapBounds, transform, unitLength * polygon.getDZ());
+    auto surface = detectorElement->surface().getSharedPtr();
+    // Return the detector element and surface
+    return std::make_tuple(detectorElement, surface);
+  } else if (nVertices == 6) {
+    double hlxnegy = fabs(vertices[0][0] - vertices[1][0]) / 2;
+    double hlxzeroy = fabs(vertices[2][0] - vertices[3][0]) / 2;
+    double hlxposy = fabs(vertices[4][0] - vertices[5][0]) / 2;
+    double hly = fabs(vertices[0][1] - vertices[4][1]) / 2;
+    std::vector<ActsScalar> halfLengths = {hlxnegy, hlxzeroy, hlxposy, hly,
+                                           hly};
+
+    // Create the surface
+
+    Vector3 colX = rotation.col(0);
+    Vector3 colY = rotation.col(1);
+    Vector3 colZ = rotation.col(2);
+    rotation.col(0) = colX;
+    rotation.col(1) = colY;
+    rotation.col(2) = colZ;
+    transform.linear() = rotation;
+
+    // Create the surface bounds
+    ActsScalar halfXnegY = unitLength * halfLengths[0];
+    ActsScalar halfXzeroY = unitLength * halfLengths[1];
+    ActsScalar halfXposY = unitLength * halfLengths[2];
+    ActsScalar halfYnegX = unitLength * halfLengths[3];
+    ActsScalar halfYposX = unitLength * halfLengths[4];
+
+    auto diamondBounds = std::make_shared<Acts::DiamondBounds>(
+        halfXnegY, halfXzeroY, halfXposY, halfYnegX, halfYposX);
+    if (!sensitive) {
+      auto surface =
+          Surface::makeShared<PlaneSurface>(transform, diamondBounds);
+      return std::make_tuple(nullptr, surface);
+    }
+    // Create the element and the surface
+    auto detectorElement =
+        GeoModelDetectorElement::createDetectorElement<PlaneSurface>(
+            geoPV, diamondBounds, transform, unitLength * polygon.getDZ());
+    auto surface = detectorElement->surface().getSharedPtr();
+    // Return the detector element and surface
+    return std::make_tuple(detectorElement, surface);
+  } else {
+    throw std::runtime_error("GeoSimplePolygonBrep with " +
+                             std::to_string(nVertices) +
+                             " can not be converted");
+  }
+}

Plugins/GeoModel/src/detail/GeoTrdConverter.cpp

@@ -36,8 +36,8 @@ Acts::detail::GeoTrdConverter::operator()(const PVConstLink& geoPV,
   Transform3 transform = Transform3::Identity();
   transform.translation() = unitLength * absTransform.translation();
 
-  // GeoTrd is defined that halfZ needs to map onto surface halfY
-  // Create the surface
+  // GeoTrd coordinates: x is the extrusion direction, y is orthogonal to the
+  // symmetry axis and z is along the symmetry axis
   ActsScalar halfX1 = geoTrd.getXHalfLength1();
   ActsScalar halfX2 = geoTrd.getXHalfLength2();
   ActsScalar halfY1 = geoTrd.getYHalfLength1();

Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp

@@ -10,10 +10,14 @@
 
 #include "Acts/Plugins/GeoModel/GeoModelDetectorObjectFactory.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelReader.hpp"
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/DiamondBounds.hpp"
 #include "Acts/Surfaces/LineBounds.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Surfaces/StrawSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Surfaces/TrapezoidBounds.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
 #include <typeinfo>
@@ -22,38 +26,69 @@
 #include <GeoModelKernel/GeoFullPhysVol.h>
 #include <GeoModelKernel/GeoLogVol.h>
 #include <GeoModelKernel/GeoMaterial.h>
+#include <GeoModelKernel/GeoSimplePolygonBrep.h>
+#include <GeoModelKernel/GeoTrd.h>
 #include <GeoModelKernel/GeoTube.h>
 
 BOOST_AUTO_TEST_SUITE(GeoModelDetObj)
 
 BOOST_AUTO_TEST_CASE(GeoModelDetectorObjectFactory) {
   // define materials
-  auto material = new GeoMaterial("Material", 1.0);
-  auto al = new GeoMaterial("Aluminium", 1.0);
+  GeoIntrusivePtr<GeoMaterial> material(new GeoMaterial("Material", 1.0));
+  GeoIntrusivePtr<GeoMaterial> al(new GeoMaterial("Aluminium", 1.0));
 
   // define dimensions
-  double gmhlx = 100, gmhly = 200, gmhlz = 2;
-  double gmrmin = 5, gmrmax = 5, gmhlzt = 100;
-  double gmhlxs = 100, gmhlys = 200, gmhlzs = 2;
+  double gmBoxHlx = 100, gmBoxHly = 200, gmBoxHlz = 2;
+  double gmTubeRmin = 5, gmTubeRmax = 6, gmTubeHlz = 100;
+  double gmRsurfHlx = 100, gmRsurfHly = 200, gmRsurfHlz = 2;
+  double gmPolyZ = 2;
+  std::vector<double> trapXVerts = {-103, 103, 183, -183};
+  std::vector<double> trapYVerts = {-50, -50, 50, 50};
+  std::vector<double> polyXVerts = {-60, 60, 153, 123, -123, -153};
+  std::vector<double> polyYVerts = {-50, -50, 0, 50, 50, 0};
 
   // create shapes
-  auto boxXY = new GeoBox(gmhlx, gmhly, gmhlz);
-  auto tube = new GeoTube(gmrmin, gmrmax, gmhlzt);
-  auto ssurface = new GeoBox(gmhlxs, gmhlys, gmhlzs);
+  GeoIntrusivePtr<GeoBox> boxXY(new GeoBox(gmBoxHlx, gmBoxHly, gmBoxHlz));
+  GeoIntrusivePtr<GeoTube> tube(new GeoTube(gmTubeRmin, gmTubeRmax, gmTubeHlz));
+  GeoIntrusivePtr<GeoBox> ssurface(
+      new GeoBox(gmRsurfHlx, gmRsurfHly, gmRsurfHlz));
+  double halfX1 = fabs(trapXVerts[0] - trapXVerts[1]) / 2;
+  double halfX2 = fabs(trapXVerts[2] - trapXVerts[3]) / 2;
+  double halfY1 = fabs(trapYVerts[0] - trapYVerts[2]) / 2;
+  GeoIntrusivePtr<GeoTrd> trd(new GeoTrd(1, 1, halfX1, halfX2, halfY1));
+  GeoIntrusivePtr<GeoSimplePolygonBrep> trap(new GeoSimplePolygonBrep(gmPolyZ));
+  for (long unsigned int i = 0; i < trapXVerts.size(); i++) {
+    trap->addVertex(trapXVerts[i], trapYVerts[i]);
+  }
+  GeoIntrusivePtr<GeoSimplePolygonBrep> poly(new GeoSimplePolygonBrep(gmPolyZ));
+  for (long unsigned int i = 0; i < polyXVerts.size(); i++) {
+    poly->addVertex(polyXVerts[i], polyYVerts[i]);
+  }
 
   // create logvols
-  auto logXY = new GeoLogVol("LogVolumeXY", boxXY, material);
-  auto logTube = new GeoLogVol("LogTube", tube, al);
-  auto logSurface = new GeoLogVol("LogSurface", ssurface, al);
+  GeoIntrusivePtr<GeoLogVol> logXY(
+      new GeoLogVol("LogVolumeXY", boxXY, material.get()));
+  GeoIntrusivePtr<GeoLogVol> logTube(new GeoLogVol("LogTube", tube, al));
+  GeoIntrusivePtr<GeoLogVol> logSurface(
+      new GeoLogVol("LogSurface", ssurface, al));
+  GeoIntrusivePtr<GeoLogVol> logTrap(new GeoLogVol("LogTrap", trap, al));
+  GeoIntrusivePtr<GeoLogVol> logTrd(new GeoLogVol("LogTrd", trd, al));
+  GeoIntrusivePtr<GeoLogVol> logPoly(new GeoLogVol("LogPoly", poly, al));
 
   // create physvols
-  auto fphysXY = new GeoFullPhysVol(logXY);
-  auto physTube = new GeoFullPhysVol(logTube);
-  auto physSurface = new GeoFullPhysVol(logSurface);
+  GeoIntrusivePtr<GeoFullPhysVol> fphysXY(new GeoFullPhysVol(logXY));
+  GeoIntrusivePtr<GeoFullPhysVol> physTube(new GeoFullPhysVol(logTube));
+  GeoIntrusivePtr<GeoFullPhysVol> physSurface(new GeoFullPhysVol(logSurface));
+  GeoIntrusivePtr<GeoFullPhysVol> physTrap(new GeoFullPhysVol(logTrap));
+  GeoIntrusivePtr<GeoFullPhysVol> physTrd(new GeoFullPhysVol(logTrd));
+  GeoIntrusivePtr<GeoFullPhysVol> physPoly(new GeoFullPhysVol(logPoly));
 
   // build hierarchy
   fphysXY->add(physTube);
   fphysXY->add(physSurface);
+  fphysXY->add(physTrap);
+  fphysXY->add(physTrd);
+  fphysXY->add(physPoly);
 
   PVConstLink physVol{fphysXY};
   auto rBounds = std::make_shared<Acts::RectangleBounds>(100, 200);
@@ -73,29 +108,50 @@ BOOST_AUTO_TEST_CASE(GeoModelDetectorObjectFactory) {
   // checking the dimension of the converted bounding boxes
   for (const auto& box : gmCache.boundingBoxes) {
     const Acts::VolumeBounds& bounds = box->volumeBounds();
-    BOOST_CHECK(gmhlx == bounds.values()[0]);
-    BOOST_CHECK(gmhly == bounds.values()[1]);
-    BOOST_CHECK(gmhlz == bounds.values()[2]);
+    BOOST_CHECK(gmBoxHlx == bounds.values()[0]);
+    BOOST_CHECK(gmBoxHly == bounds.values()[1]);
+    BOOST_CHECK(gmBoxHlz == bounds.values()[2]);
     std::vector<const Acts::Surface*> surfaces = box->surfaces();
 
     for (auto surface : surfaces) {
       const Acts::SurfaceBounds& sbounds = surface->bounds();
-      // Straw check outer radius and length without trf
+      // check straws
       if (surface->type() == Acts::Surface::SurfaceType::Straw) {
-        BOOST_CHECK(sbounds.values()[0] == gmrmax);
-        BOOST_CHECK(sbounds.values()[1] == gmhlzt);
+        const auto* lineBounds =
+            dynamic_cast<const Acts::LineBounds*>(&sbounds);
+        BOOST_CHECK(gmTubeRmax == lineBounds->get(Acts::LineBounds::eR));
+        BOOST_CHECK(gmTubeHlz ==
+                    lineBounds->get(Acts::LineBounds::eHalfLengthZ));
       }
-
-      // plane Surface check corner position without trf
-      if (surface->type() == Acts::Surface::SurfaceType::Plane) {
+      // rectangle Surface check corner position without trf
+      if (sbounds.type() == Acts::SurfaceBounds::eRectangle) {
         double csxmin = sbounds.values()[0];
         double csymin = sbounds.values()[1];
         double csxmax = sbounds.values()[2];
         double csymax = sbounds.values()[3];
-        BOOST_CHECK(gmhlxs == -csxmin);
-        BOOST_CHECK(gmhlys == -csymin);
-        BOOST_CHECK(gmhlxs == csxmax);
-        BOOST_CHECK(gmhlys == csymax);
+        BOOST_CHECK(gmRsurfHlx == -csxmin);
+        BOOST_CHECK(gmRsurfHly == -csymin);
+        BOOST_CHECK(gmRsurfHlx == csxmax);
+        BOOST_CHECK(gmRsurfHly == csymax);
+      }
+      // trap Surface without trf
+      if (sbounds.type() == Acts::SurfaceBounds::eTrapezoid) {
+        const auto* trapBounds =
+            dynamic_cast<const Acts::TrapezoidBounds*>(&sbounds);
+        std::vector<Acts::Vector2> trapVerts = trapBounds->vertices();
+        for (long unsigned int i = 0; i < trapVerts.size(); i++) {
+          BOOST_CHECK(trapVerts[i][0] == trapXVerts[i]);
+          BOOST_CHECK(trapVerts[i][1] == trapYVerts[i]);
+        }
+      }
+      if (sbounds.type() == Acts::SurfaceBounds::eDiamond) {
+        const auto* polyBounds =
+            dynamic_cast<const Acts::DiamondBounds*>(&sbounds);
+        std::vector<Acts::Vector2> polyVerts = polyBounds->vertices();
+        for (long unsigned int i = 0; i < polyVerts.size(); i++) {
+          BOOST_CHECK(polyVerts[i][0] == polyXVerts[i]);
+          BOOST_CHECK(polyVerts[i][1] == polyYVerts[i]);
+        }
       }
     }
   }