feat: geo model detector elements, converters (acts-project#3256)

Add converters for various shape types to the GeoModel plugin, as well as a factory to read surfaces from a database. Surfaces can be written via a python script to *.obj files for visualization.

Built on top of acts-project#3213, add some missing things and restructures the code a bit.

Co-authored-by: Andreas Salzburger <[email protected]>

Core/include/Acts/Surfaces/detail/AnnulusBoundsHelper.hpp

@@ -0,0 +1,40 @@
+// 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 <memory>
+#include <tuple>
+#include <vector>
+
+namespace Acts {
+
+class AnnulusBounds;
+
+namespace detail::AnnulusBoundsHelper {
+
+/// @brief The factory function to create an annulus bounds
+///
+/// @param transform the transform of final surface object
+/// @param rMin minimal radius in disc frame
+/// @param rMax maximal radius in disc frame
+/// @param vertices the vertices of the cutout trapezoid
+///
+/// @note - for the time being, the vertices follow ROOT::TGeo convention
+/// i.e. [0 - 3] vertices defince clockwise (sic!) the plane at -z
+/// and [4 - 7] vertices define clockwise (sic!) the plane at +z
+///
+/// @return AnnulusBounds
+std::tuple<std::shared_ptr<AnnulusBounds>, Transform3> create(
+    const Transform3& transform, ActsScalar rMin, ActsScalar rMax,
+    std::vector<Vector2> vertices);
+
+}  // namespace detail::AnnulusBoundsHelper
+}  // namespace Acts

Core/src/Surfaces/CMakeLists.txt

@@ -23,8 +23,9 @@ target_sources(
     SurfaceArray.cpp
     SurfaceError.cpp
     TrapezoidBounds.cpp
-    detail/AlignmentHelper.cpp
     VerticesHelper.cpp
     RegularSurface.cpp
     CurvilinearSurface.cpp
+    detail/AlignmentHelper.cpp
+    detail/AnnulusBoundsHelper.cpp
 )

Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp

@@ -0,0 +1,65 @@
+// 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/Surfaces/detail/AnnulusBoundsHelper.hpp"
+
+#include "Acts/Surfaces/AnnulusBounds.hpp"
+#include "Acts/Utilities/VectorHelpers.hpp"
+
+#include <iostream>
+
+std::tuple<std::shared_ptr<Acts::AnnulusBounds>, Acts::Transform3>
+Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform,
+                                          ActsScalar rMin, ActsScalar rMax,
+                                          std::vector<Vector2> vertices) {
+  using Line2D = Eigen::Hyperplane<double, 2>;
+
+  // Construct the bound lines
+  std::vector<std::pair<Vector2, Vector2>> boundLines;
+  for (std::size_t i = 0; i < vertices.size(); ++i) {
+    Vector2 a = vertices.at(i);
+    Vector2 b = vertices.at((i + 1) % vertices.size());
+    Vector2 ab = b - a;
+    double phi = VectorHelpers::phi(ab);
+
+    if (std::abs(phi) > 3 * M_PI / 4. || std::abs(phi) < M_PI / 4.) {
+      if (a.norm() < b.norm()) {
+        boundLines.push_back(std::make_pair(a, b));
+      } else {
+        boundLines.push_back(std::make_pair(b, a));
+      }
+    }
+  }
+
+  if (boundLines.size() != 2) {
+    throw std::logic_error(
+        "Input DiscPoly bounds type does not have sensible edges.");
+  }
+
+  Line2D lA = Line2D::Through(boundLines[0].first, boundLines[0].second);
+  Line2D lB = Line2D::Through(boundLines[1].first, boundLines[1].second);
+  Vector2 ix = lA.intersection(lB);
+
+  const Eigen::Translation3d originTranslation(ix.x(), ix.y(), 0.);
+  const Vector2 originShift = -ix;
+
+  // Update transform by prepending the origin shift translation
+  Transform3 boundsTransform = transform * originTranslation;
+  // Transform phi line point to new origin and get phi
+  double phi1 = VectorHelpers::phi(boundLines[0].second - boundLines[0].first);
+  double phi2 = VectorHelpers::phi(boundLines[1].second - boundLines[1].first);
+  double phiMax = std::max(phi1, phi2);
+  double phiMin = std::min(phi1, phi2);
+  double phiShift = 0.;
+
+  // Create the bounds
+  auto annulusBounds = std::make_shared<AnnulusBounds>(
+      rMin, rMax, phiMin, phiMax, originShift, phiShift);
+
+  return std::make_tuple(annulusBounds, boundsTransform);
+}

Examples/Python/src/GeoModel.cpp

@@ -6,9 +6,19 @@
 // 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/.
 
+// Must be on top to avoid some conflict between forward declare and typedef
+// clang-format off
+#include "GeoModelRead/ReadGeoModel.h"
+// clang-format on
+
+#include "Acts/Plugins/GeoModel/GeoModelConverters.hpp"
+#include "Acts/Plugins/GeoModel/GeoModelDetectorElement.hpp"
+#include "Acts/Plugins/GeoModel/GeoModelDetectorSurfaceFactory.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelReader.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelTree.hpp"
+#include "Acts/Plugins/GeoModel/IGeoShapeConverter.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 
 #include <string>
 
@@ -29,5 +39,97 @@ void addGeoModel(Context& ctx) {
   py::class_<Acts::GeoModelTree>(gm, "GeoModelTree").def(py::init<>());
 
   gm.def("readFromDb", &Acts::GeoModelReader::readFromDb);
+
+  py::class_<Acts::GeoModelDetectorElement,
+             std::shared_ptr<Acts::GeoModelDetectorElement>>(
+      gm, "GeoModelDetectorElement");
+
+  // Shape converters
+  {
+    py::class_<Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::IGeoShapeConverter>>(gm,
+                                                          "IGeoShapeConverter");
+
+    py::class_<Acts::GeoBoxConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoBoxConverter>>(gm, "GeoBoxConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface", &Acts::GeoBoxConverter::toSensitiveSurface)
+        .def("toPassiveSurface", &Acts::GeoBoxConverter::toPassiveSurface);
+
+    py::class_<Acts::GeoTrdConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoTrdConverter>>(gm, "GeoTrdConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface", &Acts::GeoTrdConverter::toSensitiveSurface)
+        .def("toPassiveSurface", &Acts::GeoTrdConverter::toPassiveSurface);
+
+    py::class_<Acts::GeoTubeConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoTubeConverter>>(gm, "GeoTubeConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface", &Acts::GeoTubeConverter::toSensitiveSurface)
+        .def("toPassiveSurface", &Acts::GeoTubeConverter::toPassiveSurface);
+
+    py::class_<Acts::GeoUnionDoubleTrdConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoUnionDoubleTrdConverter>>(
+        gm, "GeoUnionDoubleTrdConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface",
+             &Acts::GeoUnionDoubleTrdConverter::toSensitiveSurface)
+        .def("toPassiveSurface",
+             &Acts::GeoUnionDoubleTrdConverter::toPassiveSurface);
+
+    py::class_<Acts::GeoIntersectionAnnulusConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoIntersectionAnnulusConverter>>(
+        gm, "GeoIntersectionAnnulusConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface",
+             &Acts::GeoIntersectionAnnulusConverter::toSensitiveSurface)
+        .def("toPassiveSurface",
+             &Acts::GeoIntersectionAnnulusConverter::toPassiveSurface);
+
+    py::class_<Acts::GeoShiftConverter, Acts::IGeoShapeConverter,
+               std::shared_ptr<Acts::GeoShiftConverter>>(gm,
+                                                         "GeoShiftConverter")
+        .def(py::init<>())
+        .def("toSensitiveSurface", &Acts::GeoShiftConverter::toSensitiveSurface)
+        .def("toPassiveSurface", &Acts::GeoShiftConverter::toPassiveSurface);
+  }
+
+  // Surface factory
+  {
+    auto f =
+        py::class_<Acts::GeoModelDetectorSurfaceFactory,
+                   std::shared_ptr<Acts::GeoModelDetectorSurfaceFactory>>(
+            gm, "GeoModelDetectorSurfaceFactory")
+            .def(py::init(
+                [](const Acts::GeoModelDetectorSurfaceFactory::Config& cfg,
+                   Acts::Logging::Level level) {
+                  return std::make_shared<Acts::GeoModelDetectorSurfaceFactory>(
+                      cfg, Acts::getDefaultLogger(
+                               "GeoModelDetectorSurfaceFactory", level));
+                }))
+            .def("construct", &Acts::GeoModelDetectorSurfaceFactory::construct);
+
+    py::class_<Acts::GeoModelDetectorSurfaceFactory::Config>(f, "Config")
+        .def(py::init<>())
+        .def_readwrite(
+            "shapeConverters",
+            &Acts::GeoModelDetectorSurfaceFactory::Config::shapeConverters)
+        .def_readwrite("nameList",
+                       &Acts::GeoModelDetectorSurfaceFactory::Config::nameList)
+        .def_readwrite(
+            "materialList",
+            &Acts::GeoModelDetectorSurfaceFactory::Config::materialList);
+
+    py::class_<Acts::GeoModelDetectorSurfaceFactory::Cache>(f, "Cache")
+        .def(py::init<>())
+        .def_readwrite(
+            "sensitiveSurfaces",
+            &Acts::GeoModelDetectorSurfaceFactory::Cache::sensitiveSurfaces);
+
+    py::class_<Acts::GeoModelDetectorSurfaceFactory::Options>(f, "Options")
+        .def(py::init<>())
+        .def_readwrite("queries",
+                       &Acts::GeoModelDetectorSurfaceFactory::Options::queries);
+  }
 }
 }  // namespace Acts::Python

Examples/Python/src/Obj.cpp

@@ -34,13 +34,16 @@ void addObj(Context& ctx) {
     /// @param surfaces is the collection of surfaces
     /// @param viewContext is the geometry context
     /// @param viewRgb is the color of the surfaces
+    /// @param viewSegements is the number of segments to approximate a full circle
     /// @param fileName is the path to the output file
     ///
     mex.def("writeSurfacesObj",
             [](const std::vector<std::shared_ptr<Surface>>& surfaces,
                const GeometryContext& viewContext,
-               const std::array<int, 3>& viewRgb, const std::string& fileName) {
+               const std::array<int, 3>& viewRgb, unsigned int viewSegements,
+               const std::string& fileName) {
               Acts::ViewConfig sConfig = Acts::ViewConfig{viewRgb};
+              sConfig.nSegments = viewSegements;
               Acts::GeometryView3D view3D;
               Acts::ObjVisualization3D obj;
 

Examples/Scripts/Python/geomodel.py

@@ -0,0 +1,69 @@
+import acts
+import argparse
+from acts import logging, GeometryContext
+from acts import geomodel as gm
+from acts import examples
+
+
+def main():
+    p = argparse.ArgumentParser()
+
+    p.add_argument("-i", "--input", type=str, default="", help="Input SQL file")
+
+    p.add_argument(
+        "-q",
+        "--queries",
+        type=str,
+        nargs="+",
+        default="GeoModelXML",
+        help="List of Queries for Published full phys volumes",
+    )
+
+    p.add_argument(
+        "--output-obj",
+        help="Write the surfaces to OBJ files",
+        action="store_true",
+        default=False,
+    )
+
+    args = p.parse_args()
+
+    gContext = acts.GeometryContext()
+
+    # Read the geometry model from the database
+    gmTree = acts.geomodel.readFromDb(args.input)
+
+    gmFactoryConfig = gm.GeoModelDetectorSurfaceFactory.Config()
+    gmFactoryConfig.shapeConverters = [
+        gm.GeoBoxConverter(),
+        gm.GeoTrdConverter(),
+        gm.GeoIntersectionAnnulusConverter(),
+        gm.GeoShiftConverter(),
+        gm.GeoUnionDoubleTrdConverter(),
+    ]
+    gmFactoryConfig.nameList = []
+    gmFactoryConfig.materialList = [
+        "std::Silicon",
+    ]
+
+    gmFactory = gm.GeoModelDetectorSurfaceFactory(gmFactoryConfig, logging.VERBOSE)
+    # The options
+    gmFactoryOptions = gm.GeoModelDetectorSurfaceFactory.Options()
+    gmFactoryOptions.queries = args.queries
+    # The Cache & construct call
+    gmFactoryCache = gm.GeoModelDetectorSurfaceFactory.Cache()
+    gmFactory.construct(gmFactoryCache, gContext, gmTree, gmFactoryOptions)
+
+    # Output the surface to an OBJ file
+    segments = 720
+    if args.output_obj:
+        ssurfaces = [ss[1] for ss in gmFactoryCache.sensitiveSurfaces]
+        acts.examples.writeSurfacesObj(
+            ssurfaces, gContext, [75, 220, 100], segments, "geomodel.obj"
+        )
+
+    return
+
+
+if "__main__" == __name__:
+    main()

Plugins/GeoModel/CMakeLists.txt

@@ -2,7 +2,16 @@ include(FetchContent)
 
 add_library(
   ActsPluginGeoModel SHARED
-  src/GeoModelReader.cpp)
+  src/GeoModelConversionError.cpp
+  src/GeoModelReader.cpp
+  src/GeoModelDetectorElement.cpp
+  src/GeoModelDetectorSurfaceFactory.cpp
+  src/detail/GeoBoxConverter.cpp
+  src/detail/GeoTrdConverter.cpp
+  src/detail/GeoTubeConverter.cpp
+  src/detail//GeoShiftConverter.cpp
+  src/detail//GeoIntersectionAnnulusConverter.cpp
+  src/detail/GeoUnionDoubleTrdConverter.cpp)
 target_include_directories(
     ActsPluginGeoModel
   PUBLIC

Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelConversionError.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 <system_error>
+#include <type_traits>
+
+class GeoFullPhysVol;
+
+namespace Acts {
+
+enum class GeoModelConversionError {
+  // ensure all values are non-zero
+  WrongShapeForConverter = 1,
+  InvalidShapeParameters,
+  UnkownShape,
+  MissingLogicalVolume
+};
+
+std::error_code make_error_code(Acts::GeoModelConversionError e);
+
+}  // namespace Acts
+
+namespace std {
+// register with STL
+template <>
+struct is_error_code_enum<Acts::GeoModelConversionError> : std::true_type {};
+}  // namespace std