Core/src/GeometryGen3.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Definitions/Units.hpp"
0010 #include "Acts/Geometry/Blueprint.hpp"
0011 #include "Acts/Geometry/BlueprintNode.hpp"
0012 #include "Acts/Geometry/ContainerBlueprintNode.hpp"
0013 #include "Acts/Geometry/GeometryIdentifierBlueprintNode.hpp"
0014 #include "Acts/Geometry/LayerBlueprintNode.hpp"
0015 #include "Acts/Geometry/MaterialDesignatorBlueprintNode.hpp"
0016 #include "Acts/Geometry/PortalLinkBase.hpp"
0017 #include "Acts/Geometry/StaticBlueprintNode.hpp"
0018 #include "Acts/Geometry/VolumeAttachmentStrategy.hpp"
0019 #include "Acts/Geometry/VolumeResizeStrategy.hpp"
0020 #include "Acts/Navigation/NavigationStream.hpp"
0021 #include "Acts/Utilities/AxisDefinitions.hpp"
0022 #include "Acts/Utilities/Logger.hpp"
0023 #include "ActsPython/Utilities/Helpers.hpp"
0024 #include "ActsPython/Utilities/Macros.hpp"
0025
0026 #include <fstream>
0027 #include <random>
0028 #include <stdexcept>
0029 #include <utility>
0030
0031 #include <pybind11/functional.h>
0032 #include <pybind11/pybind11.h>
0033 #include <pybind11/pytypes.h>
0034 #include <pybind11/stl.h>
0035 #include <pybind11/stl/filesystem.h>
0036
0037 namespace py = pybind11;
0038 using namespace pybind11::literals;
0039
0040 using namespace Acts;
0041
0042 namespace ActsPython {
0043 namespace {
0044 using std::uniform_real_distribution;
0045
0046 // This is temporary!
0047 void pseudoNavigation(const TrackingGeometry& trackingGeometry,
0048                       const GeometryContext& gctx, std::filesystem::path& path,
0049                       std::size_t runs, std::size_t substepsPerCm,
0050                       std::pair<double, double> etaRange,
0051                       Logging::Level logLevel) {
0052   using namespace Acts;
0053   using namespace UnitLiterals;
0054
0055   ACTS_LOCAL_LOGGER(getDefaultLogger("pseudoNavigation", logLevel));
0056
0057   std::ofstream csv{path};
0058   csv << "x,y,z,volume,boundary,sensitive,material" << std::endl;
0059
0060   std::mt19937 rnd{42};
0061
0062   std::uniform_real_distribution<> dist{-1, 1};
0063   std::uniform_real_distribution<> subStepDist{0.01, 0.99};
0064
0065   double thetaMin = 2 * std::atan(std::exp(-etaRange.first));
0066   double thetaMax = 2 * std::atan(std::exp(-etaRange.second));
0067   std::uniform_real_distribution<> thetaDist{thetaMin, thetaMax};
0068
0069   using namespace UnitLiterals;
0070
0071   for (std::size_t run = 0; run < runs; run++) {
0072     Vector3 position = Vector3::Zero();
0073
0074     double theta = thetaDist(rnd);
0075     double phi = 2 * std::numbers::pi * dist(rnd);
0076
0077     Vector3 direction;
0078     direction[0] = std::sin(theta) * std::cos(phi);
0079     direction[1] = std::sin(theta) * std::sin(phi);
0080     direction[2] = std::cos(theta);
0081
0082     ACTS_VERBOSE("start navigation " << run);
0083     ACTS_VERBOSE("pos: " << position.transpose());
0084     ACTS_VERBOSE("dir: " << direction.transpose());
0085     ACTS_VERBOSE(direction.norm());
0086
0087     std::mt19937 rng{static_cast<unsigned int>(run)};
0088
0089     const auto* volume = trackingGeometry.lowestTrackingVolume(gctx, position);
0090     assert(volume != nullptr);
0091     ACTS_VERBOSE(volume->volumeName());
0092
0093     NavigationStream main;
0094     const TrackingVolume* currentVolume = volume;
0095
0096     csv << run << "," << position[0] << "," << position[1] << ","
0097         << position[2];
0098     csv << "," << volume->geometryId().volume();
0099     csv << "," << volume->geometryId().boundary();
0100     csv << "," << volume->geometryId().sensitive();
0101     csv << "," << 0;
0102     csv << std::endl;
0103
0104     ACTS_VERBOSE("start pseudo navigation");
0105
0106     auto writeIntersection = [&](const Vector3& pos, const Surface& surface) {
0107       csv << run << "," << pos[0] << "," << pos[1] << "," << pos[2];
0108       csv << "," << surface.geometryId().volume();
0109       csv << "," << surface.geometryId().boundary();
0110       csv << "," << surface.geometryId().sensitive();
0111       csv << "," << (surface.surfaceMaterial() != nullptr ? 1 : 0);
0112       csv << std::endl;
0113     };
0114
0115     for (std::size_t i = 0; i < 100; i++) {
0116       assert(currentVolume != nullptr);
0117       main = NavigationStream{};
0118
0119       AppendOnlyNavigationStream navStream{main};
0120       currentVolume->initializeNavigationCandidates(
0121           gctx, {.position = position, .direction = direction}, navStream,
0122           logger());
0123
0124       ACTS_VERBOSE(main.candidates().size() << " candidates");
0125
0126       for (const auto& candidate : main.candidates()) {
0127         ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0128         ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0129       }
0130
0131       ACTS_VERBOSE("initializing candidates");
0132       main.initialize(gctx, {position, direction}, BoundaryTolerance::None());
0133
0134       ACTS_VERBOSE(main.candidates().size() << " candidates remaining");
0135
0136       for (const auto& candidate : main.candidates()) {
0137         ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0138         ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0139       }
0140
0141       if (main.currentCandidate().surface().isOnSurface(gctx, position,
0142                                                         direction)) {
0143         ACTS_VERBOSE(
0144             "Already on surface at initialization, skipping candidate");
0145
0146         writeIntersection(position, main.currentCandidate().surface());
0147
0148         if (!main.switchToNextCandidate()) {
0149           ACTS_WARNING("candidates exhausted unexpectedly");
0150           break;
0151         }
0152       }
0153
0154       bool terminated = false;
0155       while (main.remainingCandidates() > 0) {
0156         const auto& candidate = main.currentCandidate();
0157
0158         ACTS_VERBOSE(candidate.position().transpose());
0159
0160         ACTS_VERBOSE("moving to position: " << position.transpose() << " (r="
0161                                             << VectorHelpers::perp(position)
0162                                             << ")");
0163
0164         Vector3 delta = candidate.position() - position;
0165
0166         std::size_t substeps =
0167             std::max(1l, std::lround(delta.norm() / 10_cm * substepsPerCm));
0168
0169         for (std::size_t j = 0; j < substeps; j++) {
0170           Vector3 subpos = position + subStepDist(rng) * delta;
0171           csv << run << "," << subpos[0] << "," << subpos[1] << ","
0172               << subpos[2];
0173           csv << "," << currentVolume->geometryId().volume();
0174           csv << ",0,0,0";  // zero boundary and sensitive ids
0175           csv << std::endl;
0176         }
0177
0178         position = candidate.position();
0179         ACTS_VERBOSE("                 -> "
0180                      << position.transpose()
0181                      << " (r=" << VectorHelpers::perp(position) << ")");
0182
0183         writeIntersection(position, candidate.surface());
0184
0185         if (candidate.isPortalTarget()) {
0186           ACTS_VERBOSE("On portal: " << candidate.surface().toStream(gctx));
0187           currentVolume = candidate.portal()
0188                               .resolveVolume(gctx, position, direction)
0189                               .value();
0190
0191           if (currentVolume == nullptr) {
0192             ACTS_VERBOSE("switched to nullptr -> we're done");
0193             terminated = true;
0194           }
0195           break;
0196
0197         } else {
0198           ACTS_VERBOSE("Not on portal");
0199         }
0200
0201         main.switchToNextCandidate();
0202       }
0203
0204       if (terminated) {
0205         ACTS_VERBOSE("Terminate pseudo navigation");
0206         break;
0207       }
0208
0209       ACTS_VERBOSE("switched to " << currentVolume->volumeName());
0210
0211       ACTS_VERBOSE("-----");
0212     }
0213   }
0214 }
0215
0216 }  // namespace
0217
0218 /// This adds the geometry building bindings for the Gen3 geometry
0219 /// @param m the module to add the bindings to
0220 void addGeometryGen3(py::module_& m) {
0221   using Experimental::Blueprint;
0222   using Experimental::BlueprintNode;
0223   using Experimental::BlueprintOptions;
0224   using Experimental::CuboidContainerBlueprintNode;
0225   using Experimental::CylinderContainerBlueprintNode;
0226   using Experimental::GeometryIdentifierBlueprintNode;
0227   using Experimental::LayerBlueprintNode;
0228   using Experimental::MaterialDesignatorBlueprintNode;
0229   using Experimental::StaticBlueprintNode;
0230
0231   py::class_<Portal>(m, "Portal");
0232
0233   auto blueprintNode =
0234       py::class_<BlueprintNode, std::shared_ptr<BlueprintNode>>(
0235           m, "BlueprintNode");
0236
0237   auto rootNode =
0238       py::class_<Blueprint, BlueprintNode, std::shared_ptr<Blueprint>>(
0239           m, "Blueprint");
0240
0241   rootNode
0242       .def(py::init<const Blueprint::Config&>())
0243       // Return value needs to be shared pointer because python otherwise
0244       // can't manage the lifetime
0245       .def(
0246           "construct",
0247           [](Blueprint& self, const BlueprintOptions& options,
0248              const GeometryContext& gctx,
0249              Logging::Level level) -> std::shared_ptr<TrackingGeometry> {
0250             return self.construct(options, gctx,
0251                                   *getDefaultLogger("Blueprint", level));
0252           },
0253           py::arg("options"), py::arg("gctx"),
0254           py::arg("level") = Logging::INFO);
0255
0256   {
0257     auto c = py::class_<Blueprint::Config>(rootNode, "Config").def(py::init());
0258     ACTS_PYTHON_STRUCT(c, envelope);
0259   }
0260
0261   auto addContextManagerProtocol = []<typename class_>(class_& cls) {
0262     using type = typename class_::type;
0263     cls.def("__enter__", [](type& self) -> type& { return self; })
0264         .def("__exit__", [](type& /*self*/, const py::object& /*exc_type*/,
0265                             const py::object& /*exc_value*/,
0266                             const py::object& /*traceback*/) {
0267           // No action needed on exit
0268         });
0269   };
0270
0271   auto addNodeMethods = [&blueprintNode](
0272                             std::initializer_list<std::string> names,
0273                             auto&& callable, auto&&... args) {
0274     for (const auto& name : names) {
0275       blueprintNode.def(name.c_str(), callable, args...);
0276     }
0277   };
0278
0279   blueprintNode
0280       .def("__str__",
0281            [](const BlueprintNode& self) {
0282              std::stringstream ss;
0283              ss << self;
0284              return ss.str();
0285            })
0286       .def("addChild", &BlueprintNode::addChild)
0287       .def_property_readonly("children",
0288                              py::overload_cast<>(&BlueprintNode::children))
0289       .def("clearChildren", &BlueprintNode::clearChildren)
0290       .def_property_readonly("name", &BlueprintNode::name)
0291       .def_property_readonly("depth", &BlueprintNode::depth)
0292       .def("graphviz", [](BlueprintNode& self, const py::object& fh) {
0293         std::stringstream ss;
0294         self.graphviz(ss);
0295         fh.attr("write")(ss.str());
0296       });
0297
0298   py::class_<BlueprintOptions>(m, "BlueprintOptions")
0299       .def(py::init<>())
0300       .def_readwrite("defaultNavigationPolicyFactory",
0301                      &BlueprintOptions::defaultNavigationPolicyFactory);
0302
0303   py::class_<BlueprintNode::MutableChildRange>(blueprintNode,
0304                                                "MutableChildRange")
0305       .def(
0306           "__iter__",
0307           [](BlueprintNode::MutableChildRange& self) {
0308             return py::make_iterator(self.begin(), self.end());
0309           },
0310           py::keep_alive<0, 1>())
0311       .def(
0312           "__getitem__",
0313           [](BlueprintNode::MutableChildRange& self, int i) -> BlueprintNode& {
0314             if (i < 0) {
0315               i += self.size();
0316             }
0317             return self.at(i);
0318           },
0319           py::return_value_policy::reference_internal)
0320       .def("__len__", [](const BlueprintNode::MutableChildRange& self) {
0321         return self.size();
0322       });
0323
0324   auto staticNode =
0325       py::class_<StaticBlueprintNode, BlueprintNode,
0326                  std::shared_ptr<StaticBlueprintNode>>(m, "StaticBlueprintNode")
0327           .def(py::init([](const Transform3& transform,
0328                            const std::shared_ptr<VolumeBounds>& bounds,
0329                            const std::string& name) {
0330                  return std::make_shared<StaticBlueprintNode>(
0331                      std::make_unique<TrackingVolume>(transform, bounds, name));
0332                }),
0333                py::arg("transform"), py::arg("bounds"),
0334                py::arg("name") = "undefined")
0335           .def_property("navigationPolicyFactory",
0336                         &StaticBlueprintNode::navigationPolicyFactory,
0337                         &StaticBlueprintNode::setNavigationPolicyFactory);
0338
0339   addContextManagerProtocol(staticNode);
0340
0341   addNodeMethods(
0342       {"StaticVolume", "addStaticVolume"},
0343       [](BlueprintNode& self, const Transform3& transform,
0344          const std::shared_ptr<VolumeBounds>& bounds, const std::string& name) {
0345         auto node = std::make_shared<StaticBlueprintNode>(
0346             std::make_unique<TrackingVolume>(transform, bounds, name));
0347         self.addChild(node);
0348         return node;
0349       },
0350       py::arg("transform"), py::arg("bounds"), py::arg("name") = "undefined");
0351
0352   auto cylNode =
0353       py::class_<CylinderContainerBlueprintNode, BlueprintNode,
0354                  std::shared_ptr<CylinderContainerBlueprintNode>>(
0355           m, "CylinderContainerBlueprintNode")
0356           .def(py::init<const std::string&, AxisDirection,
0357                         VolumeAttachmentStrategy, VolumeResizeStrategy>(),
0358                py::arg("name"), py::arg("direction"),
0359                py::arg("attachmentStrategy") = VolumeAttachmentStrategy::Gap,
0360                py::arg("resizeStrategy") = VolumeResizeStrategy::Gap)
0361           .def_property("attachmentStrategy",
0362                         &CylinderContainerBlueprintNode::attachmentStrategy,
0363                         &CylinderContainerBlueprintNode::setAttachmentStrategy)
0364           .def_property("resizeStrategies",
0365                         &CylinderContainerBlueprintNode::resizeStrategies,
0366                         [](CylinderContainerBlueprintNode& self,
0367                            std::pair<VolumeResizeStrategy, VolumeResizeStrategy>
0368                                strategies) {
0369                           self.setResizeStrategies(strategies.first,
0370                                                    strategies.second);
0371                         })
0372           .def_property("direction", &CylinderContainerBlueprintNode::direction,
0373                         &CylinderContainerBlueprintNode::setDirection);
0374
0375   addContextManagerProtocol(cylNode);
0376
0377   addNodeMethods(
0378       {"CylinderContainer", "addCylinderContainer"},
0379       [](BlueprintNode& self, const std::string& name,
0380          AxisDirection direction) {
0381         auto cylinder =
0382             std::make_shared<CylinderContainerBlueprintNode>(name, direction);
0383         self.addChild(cylinder);
0384         return cylinder;
0385       },
0386       py::arg("name"), py::arg("direction"));
0387
0388   auto boxNode =
0389       py::class_<CuboidContainerBlueprintNode, BlueprintNode,
0390                  std::shared_ptr<CuboidContainerBlueprintNode>>(
0391           m, "CuboidContainerBlueprintNode")
0392           .def(py::init<const std::string&, AxisDirection,
0393                         VolumeAttachmentStrategy, VolumeResizeStrategy>(),
0394                py::arg("name"), py::arg("direction"),
0395                py::arg("attachmentStrategy") = VolumeAttachmentStrategy::Gap,
0396                py::arg("resizeStrategy") = VolumeResizeStrategy::Gap)
0397           .def_property("attachmentStrategy",
0398                         &CuboidContainerBlueprintNode::attachmentStrategy,
0399                         &CuboidContainerBlueprintNode::setAttachmentStrategy)
0400           .def_property("resizeStrategies",
0401                         &CuboidContainerBlueprintNode::resizeStrategies,
0402                         &CuboidContainerBlueprintNode::setResizeStrategies)
0403           .def_property("direction", &CuboidContainerBlueprintNode::direction,
0404                         &CuboidContainerBlueprintNode::setDirection);
0405
0406   addContextManagerProtocol(boxNode);
0407
0408   addNodeMethods(
0409       {"CuboidContainer", "addCuboidContainer"},
0410       [](BlueprintNode& self, const std::string& name,
0411          AxisDirection direction) {
0412         auto cylinder =
0413             std::make_shared<CuboidContainerBlueprintNode>(name, direction);
0414         self.addChild(cylinder);
0415         return cylinder;
0416       },
0417       py::arg("name"), py::arg("direction"));
0418
0419   auto matNode = py::class_<MaterialDesignatorBlueprintNode, BlueprintNode,
0420                             std::shared_ptr<MaterialDesignatorBlueprintNode>>(
0421                      m, "MaterialDesignatorBlueprintNode")
0422                      .def(py::init<const std::string&>(), "name"_a)
0423                      .def("configureFace",
0424                           py::overload_cast<CylinderVolumeBounds::Face,
0425                                             const DirectedProtoAxis&,
0426                                             const DirectedProtoAxis&>(
0427                               &MaterialDesignatorBlueprintNode::configureFace),
0428                           "face"_a, "loc0"_a, "loc1"_a)
0429                      .def("configureFace",
0430                           py::overload_cast<CuboidVolumeBounds::Face,
0431                                             const DirectedProtoAxis&,
0432                                             const DirectedProtoAxis&>(
0433                               &MaterialDesignatorBlueprintNode::configureFace),
0434                           "face"_a, "loc0"_a, "loc1"_a);
0435
0436   addContextManagerProtocol(matNode);
0437
0438   addNodeMethods(
0439       {"Material", "addMaterial"},
0440       [](BlueprintNode& self, const std::string& name) {
0441         auto child = std::make_shared<MaterialDesignatorBlueprintNode>(name);
0442         self.addChild(child);
0443         return child;
0444       },
0445       "name"_a);
0446
0447   auto layerNode =
0448       py::class_<LayerBlueprintNode, StaticBlueprintNode,
0449                  std::shared_ptr<LayerBlueprintNode>>(m, "LayerBlueprintNode")
0450           .def(py::init<const std::string&>(), py::arg("name"))
0451           .def_property_readonly("name", &LayerBlueprintNode::name)
0452           .def_property("surfaces", &LayerBlueprintNode::surfaces,
0453                         &LayerBlueprintNode::setSurfaces)
0454           .def_property("transform", &LayerBlueprintNode::transform,
0455                         &LayerBlueprintNode::setTransform)
0456           .def_property("envelope", &LayerBlueprintNode::envelope,
0457                         &LayerBlueprintNode::setEnvelope)
0458           .def_property("layerType", &LayerBlueprintNode::layerType,
0459                         &LayerBlueprintNode::setLayerType)
0460           .def_property("navigationPolicyFactory",
0461                         &LayerBlueprintNode::navigationPolicyFactory,
0462                         &LayerBlueprintNode::setNavigationPolicyFactory);
0463
0464   py::enum_<LayerBlueprintNode::LayerType>(layerNode, "LayerType")
0465       .value("Cylinder", LayerBlueprintNode::LayerType::Cylinder)
0466       .value("Disc", LayerBlueprintNode::LayerType::Disc)
0467       .value("Plane", LayerBlueprintNode::LayerType::Plane);
0468
0469   addContextManagerProtocol(layerNode);
0470
0471   addNodeMethods(
0472       {"Layer", "addLayer"},
0473       [](BlueprintNode& self, const std::string& name) {
0474         auto child = std::make_shared<LayerBlueprintNode>(name);
0475         self.addChild(child);
0476         return child;
0477       },
0478       py::arg("name"));
0479
0480   auto geoIdNode =
0481       py::class_<GeometryIdentifierBlueprintNode, BlueprintNode,
0482                  std::shared_ptr<GeometryIdentifierBlueprintNode>>(
0483           m, "GeometryIdentifierBlueprintNode")
0484           .def(py::init<>())
0485           .def("setLayerIdTo", &GeometryIdentifierBlueprintNode::setLayerIdTo,
0486                py::arg("value"))
0487           .def("incrementLayerIds",
0488                &GeometryIdentifierBlueprintNode::incrementLayerIds,
0489                py::arg("start") = 0)
0490           .def("setAllVolumeIdsTo",
0491                &GeometryIdentifierBlueprintNode::setAllVolumeIdsTo,
0492                py::arg("value"))
0493           // Need to do some massaging to avoid copy issues
0494           .def(
0495               "sortBy",
0496               [](GeometryIdentifierBlueprintNode& self,
0497                  const py::function& func) -> GeometryIdentifierBlueprintNode& {
0498                 if (func.is_none()) {
0499                   throw std::invalid_argument(
0500                       "sortBy requires a comparison function");
0501                 }
0502                 return self.sortBy(
0503                     [func](const TrackingVolume& a, const TrackingVolume& b) {
0504                       return func(&a, &b).cast<bool>();
0505                     });
0506               },
0507               py::arg("compare"));
0508
0509   auto geoIdFactory = [](BlueprintNode& self) {
0510     auto child = std::make_shared<GeometryIdentifierBlueprintNode>();
0511     self.addChild(child);
0512     return child;
0513   };
0514
0515   addNodeMethods({"GeometryIdentifier", "withGeometryIdentifier"},
0516                  geoIdFactory);
0517   addContextManagerProtocol(geoIdNode);
0518
0519   // TEMPORARY
0520   m.def("pseudoNavigation", &pseudoNavigation, "trackingGeometry"_a, "gctx"_a,
0521         "path"_a, "runs"_a, "substepsPerCm"_a = 2,
0522         "etaRange"_a = std::pair{-4.5, 4.5}, "logLevel"_a = Logging::INFO);
0523 }
0524
0525 }  // namespace ActsPython