feat: establishing Detector from GeoModel with db blueprint (Gen2) (a…

…cts-project#3289)

This PR allows to built cylindrical detectors from augmented GeoModel sqlite databases - using the Gen2 `Acts::Detector` infrastructure.

This will eventually be replaced by the Gen3 `Blueprint`, but will allow in the meantime to exercise the navigation validation and material mapping tests.

Core/include/Acts/Detector/Blueprint.hpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Detector/ProtoBinning.hpp"
+#include "Acts/Geometry/Extent.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Utilities/BinningData.hpp"
 #include "Acts/Utilities/StringHelpers.hpp"
@@ -48,15 +49,17 @@ struct Node final {
   /// @param bv the boundary values
   /// @param bss the binning values
   /// @param cs the children of the node
+  /// @param e the estimated extent of the node (optional)
   Node(const std::string& n, const Transform3& t, VolumeBounds::BoundsType bt,
        const std::vector<ActsScalar>& bv, const std::vector<BinningValue>& bss,
-       std::vector<std::unique_ptr<Node>> cs = {})
+       std::vector<std::unique_ptr<Node>> cs = {}, const Extent& e = Extent())
       : name(n),
         transform(t),
         boundsType(bt),
         boundaryValues(bv),
         children(std::move(cs)),
-        binning(bss) {
+        binning(bss),
+        extent(e) {
     for_each(children.begin(), children.end(),
              [this](std::unique_ptr<Node>& c) { c->parent = this; });
   }
@@ -68,22 +71,25 @@ struct Node final {
   /// @param bt the boundary type
   /// @param bv the boundary values
   /// @param isb the internal structure builder (optional)
+  /// @param e the estimated extent of the node (optional)
   Node(const std::string& n, const Transform3& t, VolumeBounds::BoundsType bt,
        const std::vector<ActsScalar>& bv,
-       std::shared_ptr<const IInternalStructureBuilder> isb = nullptr)
+       std::shared_ptr<const IInternalStructureBuilder> isb = nullptr,
+       const Extent& e = Extent())
       : name(n),
         transform(t),
         boundsType(bt),
         boundaryValues(bv),
-        internalsBuilder(std::move(isb)) {}
+        internalsBuilder(std::move(isb)),
+        extent(e) {}
 
   /// Name identification of this node
   std::string name = "";
   /// Transform definition of this node
   Transform3 transform = Transform3::Identity();
-  /// Boundary definition of this node
-  VolumeBounds::BoundsType boundsType = VolumeBounds::eOther;
   /// The boundary type
+  VolumeBounds::BoundsType boundsType = VolumeBounds::eOther;
+  /// The associated values
   std::vector<ActsScalar> boundaryValues = {};
   /// Parent node - nullptr for root only
   const Node* parent = nullptr;
@@ -106,6 +112,9 @@ struct Node final {
   /// Internal structure builder - for leaf nodes
   std::shared_ptr<const IInternalStructureBuilder> internalsBuilder = nullptr;
 
+  /// An optional extent object
+  Extent extent = Extent();
+
   /// @brief Check if it is a leaf node
   bool isLeaf() const { return children.empty(); }
 

Core/include/Acts/Geometry/Extent.hpp

@@ -111,6 +111,18 @@ class Extent {
   /// @param max the maximum parameter
   void set(BinningValue bValue, ActsScalar min, ActsScalar max);
 
+  /// Set a min value for a dedicated binning value
+  ///
+  /// @param bValue the binning identification
+  /// @param min the minimum parameter
+  void setMin(BinningValue bValue, ActsScalar min);
+
+  /// Set a max value for a dedicated binning value
+  ///
+  /// @param bValue the binning identification
+  /// @param max the maximum parameter
+  void setMax(BinningValue bValue, ActsScalar max);
+
   /// (re-)Set the envelope
   ///
   /// @param envelope new envelope to be set

Core/src/Geometry/Extent.cpp

@@ -100,6 +100,20 @@ void Acts::Extent::set(BinningValue bValue, ActsScalar min, ActsScalar max) {
   m_constrains.set(bValue);
 }
 
+void Acts::Extent::setMin(BinningValue bValue, ActsScalar min) {
+  ActsScalar minval = min;
+  if (bValue == binR && minval < 0.) {
+    minval = 0.;
+  }
+  m_range[bValue].setMin(0u, minval);
+  m_constrains.set(bValue);
+}
+
+void Acts::Extent::setMax(BinningValue bValue, ActsScalar max) {
+  m_range[bValue].setMax(0u, max);
+  m_constrains.set(bValue);
+}
+
 void Acts::Extent::setEnvelope(const ExtentEnvelope& envelope) {
   m_envelope = envelope;
 }

Examples/Python/src/GeoModel.cpp

@@ -11,6 +11,8 @@
 #include "GeoModelRead/ReadGeoModel.h"
 // clang-format on
 
+#include "Acts/Detector/CylindricalContainerBuilder.hpp"
+#include "Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelConverters.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelDetectorElement.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelDetectorSurfaceFactory.hpp"
@@ -131,5 +133,49 @@ void addGeoModel(Context& ctx) {
         .def_readwrite("queries",
                        &Acts::GeoModelDetectorSurfaceFactory::Options::queries);
   }
+
+  {
+    py::class_<Acts::GeoModelBlueprintCreater::Blueprint,
+               std::shared_ptr<Acts::GeoModelBlueprintCreater::Blueprint>>(
+        gm, "Blueprint")
+        .def("convertToBuilder",
+             [](Acts::GeoModelBlueprintCreater::Blueprint& self,
+                Acts::Logging::Level level) {
+               // It's a container builder
+               return std::make_shared<
+                   Acts::Experimental::CylindricalContainerBuilder>(self.node(),
+                                                                    level);
+             });
+
+    auto bpc =
+        py::class_<Acts::GeoModelBlueprintCreater,
+                   std::shared_ptr<Acts::GeoModelBlueprintCreater>>(
+            gm, "GeoModelBlueprintCreater")
+            .def(py::init([](const Acts::GeoModelBlueprintCreater::Config& cfg,
+                             Acts::Logging::Level level) {
+              return std::make_shared<Acts::GeoModelBlueprintCreater>(
+                  cfg,
+                  Acts::getDefaultLogger("GeoModelBlueprintCreater", level));
+            }))
+            .def("create", &Acts::GeoModelBlueprintCreater::create);
+
+    py::class_<Acts::GeoModelBlueprintCreater::Config>(bpc, "Config")
+        .def(py::init<>())
+        .def_readwrite(
+            "detectorSurfaces",
+            &Acts::GeoModelBlueprintCreater::Config::detectorSurfaces)
+        .def_readwrite("kdtBinning",
+                       &Acts::GeoModelBlueprintCreater::Config::kdtBinning);
+
+    py::class_<Acts::GeoModelBlueprintCreater::Options>(bpc, "Options")
+        .def(py::init<>())
+        .def_readwrite("topEntry",
+                       &Acts::GeoModelBlueprintCreater::Options::topEntry)
+        .def_readwrite(
+            "topBoundsOverride",
+            &Acts::GeoModelBlueprintCreater::Options::topBoundsOverride)
+        .def_readwrite("table",
+                       &Acts::GeoModelBlueprintCreater::Options::table);
+  }
 }
 }  // namespace Acts::Python

Examples/Scripts/MaterialMapping/material_comparison.py

@@ -90,7 +90,7 @@ def TH1D_from_TProf(tprof):
         "--eta",
         type=float,
         nargs=2,
-        default=[-3.0, 3.0],
+        default=[-4.0, 4.0],
         help="Eta range for the plotting",
     )
     p.add_argument("--eta-bins", type=int, default=60, help="Eta bins for the plotting")

Examples/Scripts/Python/geomodel.py

@@ -1,6 +1,12 @@
 import acts
 import argparse
-from acts import logging, GeometryContext
+from acts import (
+    logging,
+    GeometryContext,
+    CylindricalContainerBuilder,
+    DetectorBuilder,
+    GeometryIdGenerator,
+)
 from acts import geomodel as gm
 from acts import examples
 
@@ -10,6 +16,10 @@ def main():
 
     p.add_argument("-i", "--input", type=str, default="", help="Input SQL file")
 
+    p.add_argument(
+        "-o", "--output", type=str, default="GeoModel", help="Output file(s) base name"
+    )
+
     p.add_argument(
         "-q",
         "--queries",
@@ -19,16 +29,70 @@ def main():
         help="List of Queries for Published full phys volumes",
     )
 
+    p.add_argument(
+        "--table-name",
+        type=str,
+        default="ActsBlueprint",
+        help="Name of the blueprint table",
+    )
+
+    p.add_argument(
+        "-t",
+        "--top-node",
+        type=str,
+        default="",
+        help="Name of the top node in the blueprint tree",
+    )
+
+    p.add_argument(
+        "-b",
+        "--top-node-bounds",
+        type=str,
+        default="",
+        help="Table entry string overriding the top node bounds",
+    )
+
+    p.add_argument(
+        "-m", "--map", type=str, default="", help="Input file for the material map"
+    )
+
+    p.add_argument(
+        "--output-svg",
+        help="Write the surfaces to SVG files",
+        action="store_true",
+        default=False,
+    )
+
+    p.add_argument(
+        "--output-internals-svg",
+        help="Write the internal navigation to SVG files",
+        action="store_true",
+        default=False,
+    )
+
     p.add_argument(
         "--output-obj",
         help="Write the surfaces to OBJ files",
         action="store_true",
         default=False,
     )
 
+    p.add_argument(
+        "--output-json",
+        help="Write the surfaces to OBJ files",
+        action="store_true",
+        default=False,
+    )
+
     args = p.parse_args()
 
     gContext = acts.GeometryContext()
+    logLevel = logging.INFO
+
+    materialDecorator = None
+    if args.map != "":
+        print("Loading a material decorator from file:", args.map)
+        materialDecorator = acts.IMaterialDecorator.fromFile(args.map)
 
     # Read the geometry model from the database
     gmTree = acts.geomodel.readFromDb(args.input)
@@ -38,29 +102,105 @@ def main():
         gm.GeoBoxConverter(),
         gm.GeoTrdConverter(),
         gm.GeoIntersectionAnnulusConverter(),
-        gm.GeoShiftConverter(),
-        gm.GeoUnionDoubleTrdConverter(),
     ]
-    gmFactoryConfig.nameList = []
-    gmFactoryConfig.materialList = [
-        "std::Silicon",
-    ]
-
-    gmFactory = gm.GeoModelDetectorSurfaceFactory(gmFactoryConfig, logging.VERBOSE)
+    gmFactory = gm.GeoModelDetectorSurfaceFactory(gmFactoryConfig, logLevel)
     # The options
     gmFactoryOptions = gm.GeoModelDetectorSurfaceFactory.Options()
     gmFactoryOptions.queries = args.queries
     # The Cache & construct call
     gmFactoryCache = gm.GeoModelDetectorSurfaceFactory.Cache()
     gmFactory.construct(gmFactoryCache, gContext, gmTree, gmFactoryOptions)
 
+    # All surfaces from GeoModel
+    gmSurfaces = [ss[1] for ss in gmFactoryCache.sensitiveSurfaces]
+
+    # Construct the building hierarchy
+    gmBlueprintConfig = gm.GeoModelBlueprintCreater.Config()
+    gmBlueprintConfig.detectorSurfaces = gmSurfaces
+    gmBlueprintConfig.kdtBinning = [acts.Binning.z, acts.Binning.r]
+
+    gmBlueprintOptions = gm.GeoModelBlueprintCreater.Options()
+    gmBlueprintOptions.table = args.table_name
+    gmBlueprintOptions.topEntry = args.top_node
+    if len(args.top_node_bounds) > 0:
+        gmBlueprintOptions.topBoundsOverride = args.top_node_bounds
+
+    gmBlueprintCreater = gm.GeoModelBlueprintCreater(gmBlueprintConfig, logLevel)
+    gmBlueprint = gmBlueprintCreater.create(gContext, gmTree, gmBlueprintOptions)
+
+    gmCylindricalBuilder = gmBlueprint.convertToBuilder(logLevel)
+
+    # Top level geo id generator
+    gmGeoIdConfig = GeometryIdGenerator.Config()
+    gmGeoIdGenerator = GeometryIdGenerator(
+        gmGeoIdConfig, "GeoModelGeoIdGenerator", logLevel
+    )
+
+    # Create the detector builder
+    gmDetectorConfig = DetectorBuilder.Config()
+    gmDetectorConfig.name = args.top_node + "_DetectorBuilder"
+    gmDetectorConfig.builder = gmCylindricalBuilder
+    gmDetectorConfig.geoIdGenerator = gmGeoIdGenerator
+    gmDetectorConfig.materialDecorator = materialDecorator
+    gmDetectorConfig.auxiliary = (
+        "GeoModel based Acts::Detector from '" + args.input + "'"
+    )
+
+    gmDetectorBuilder = DetectorBuilder(gmDetectorConfig, args.top_node, logLevel)
+    detector = gmDetectorBuilder.construct(gContext)
+
+    materialSurfaces = detector.extractMaterialSurfaces()
+    print("Found ", len(materialSurfaces), " material surfaces")
+
+    # Output the detector to SVG
+    if args.output_svg:
+        surfaceStyle = acts.svg.Style()
+        surfaceStyle.fillColor = [5, 150, 245]
+        surfaceStyle.fillOpacity = 0.5
+
+        surfaceOptions = acts.svg.SurfaceOptions()
+        surfaceOptions.style = surfaceStyle
+
+        viewRange = acts.Extent([])
+        volumeOptions = acts.svg.DetectorVolumeOptions()
+        volumeOptions.surfaceOptions = surfaceOptions
+
+        xyRange = acts.Extent([[acts.Binning.z, [-50, 50]]])
+        zrRange = acts.Extent([[acts.Binning.phi, [-0.8, 0.8]]])
+
+        acts.svg.viewDetector(
+            gContext,
+            detector,
+            args.top_node,
+            [[ivol, volumeOptions] for ivol in range(detector.numberVolumes())],
+            [
+                ["xy", ["sensitives", "portals"], xyRange],
+                ["zr", ["", "", "materials"], zrRange],
+            ],
+            args.output + "_detector",
+        )
+
+        # Output the internal navigation to SVG
+        if args.output_internals_svg:
+            for vol in detector.volumes():
+                acts.svg.viewInternalNavigation(
+                    gContext, vol, [66, 111, 245, 245, 203, 66, 0.8], "/;:"
+                )
+
     # Output the surface to an OBJ file
-    segments = 720
     if args.output_obj:
+        segments = 720
         ssurfaces = [ss[1] for ss in gmFactoryCache.sensitiveSurfaces]
         acts.examples.writeSurfacesObj(
-            ssurfaces, gContext, [75, 220, 100], segments, "geomodel.obj"
+            ssurfaces,
+            gContext,
+            [75, 220, 100],
+            segments,
+            args.output + "_sensitives.obj",
         )
+    # Output to a JSON file
+    if args.output_json:
+        acts.examples.writeDetectorToJsonDetray(gContext, detector, args.output)
 
     return