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