Core/Geometry/BlueprintApiTests.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 <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/old/interface.hpp>
0011 #include <boost/test/unit_test.hpp>
0012 #include <boost/test/unit_test_suite.hpp>
0013
0014 #include "Acts/Definitions/Algebra.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/Geometry/Blueprint.hpp"
0017 #include "Acts/Geometry/ContainerBlueprintNode.hpp"
0018 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0019 #include "Acts/Geometry/GeometryContext.hpp"
0020 #include "Acts/Geometry/LayerBlueprintNode.hpp"
0021 #include "Acts/Geometry/MaterialDesignatorBlueprintNode.hpp"
0022 #include "Acts/Geometry/TrackingGeometry.hpp"
0023 #include "Acts/Geometry/TrackingVolume.hpp"
0024 #include "Acts/Geometry/VolumeAttachmentStrategy.hpp"
0025 #include "Acts/Geometry/VolumeResizeStrategy.hpp"
0026 #include "Acts/Navigation/INavigationPolicy.hpp"
0027 #include "Acts/Navigation/NavigationStream.hpp"
0028 #include "Acts/Surfaces/RectangleBounds.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "Acts/Utilities/ProtoAxis.hpp"
0031 #include "Acts/Visualization/GeometryView3D.hpp"
0032 #include "Acts/Visualization/ObjVisualization3D.hpp"
0033 #include "ActsTests/CommonHelpers/DetectorElementStub.hpp"
0034
0035 #include <fstream>
0036 #include <random>
0037 #include <vector>
0038
0039 using namespace Acts;
0040 using namespace UnitLiterals;
0041
0042 using Experimental::Blueprint;
0043 using Experimental::LayerBlueprintNode;
0044 using Experimental::MaterialDesignatorBlueprintNode;
0045
0046 namespace ActsTests {
0047
0048 auto logger = getDefaultLogger("UnitTests", Logging::DEBUG);
0049
0050 auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0051
0052 inline std::vector<std::shared_ptr<Surface>> makeFanLayer(
0053     const Transform3& base,
0054     std::vector<std::unique_ptr<SurfacePlacementBase>>& elements,
0055     double r = 300_mm, std::size_t nSensors = 8, double thickness = 0) {
0056   auto recBounds = std::make_shared<RectangleBounds>(40_mm, 60_mm);
0057
0058   double deltaPhi = 2 * std::numbers::pi / nSensors;
0059   std::vector<std::shared_ptr<Surface>> surfaces;
0060   for (std::size_t i = 0; i < nSensors; i++) {
0061     // Create a fan of sensors
0062
0063     Transform3 trf = base * AngleAxis3{deltaPhi * i, Vector3::UnitZ()} *
0064                      Translation3(Vector3::UnitX() * r);
0065
0066     if (i % 2 == 0) {
0067       trf = trf * Translation3{Vector3::UnitZ() * 5_mm};
0068     }
0069
0070     auto& element = elements.emplace_back(
0071         std::make_unique<DetectorElementStub>(trf, recBounds, thickness));
0072
0073     element->surface().assignSurfacePlacement(*element);
0074
0075     surfaces.push_back(element->surface().getSharedPtr());
0076   }
0077   return surfaces;
0078 }
0079
0080 inline std::vector<std::shared_ptr<Surface>> makeBarrelLayer(
0081     const Transform3& base,
0082     std::vector<std::unique_ptr<SurfacePlacementBase>>& elements,
0083     double r = 300_mm, std::size_t nStaves = 10, int nSensorsPerStave = 8,
0084     double thickness = 0, double hlPhi = 40_mm, double hlZ = 60_mm) {
0085   auto recBounds = std::make_shared<RectangleBounds>(hlPhi, hlZ);
0086
0087   double deltaPhi = 2 * std::numbers::pi / nStaves;
0088   std::vector<std::shared_ptr<Surface>> surfaces;
0089
0090   for (std::size_t istave = 0; istave < nStaves; istave++) {
0091     for (int isensor = -nSensorsPerStave; isensor <= nSensorsPerStave;
0092          isensor++) {
0093       double z = isensor * (2 * hlZ + 5_mm);
0094
0095       Transform3 trf = base * Translation3(Vector3::UnitZ() * z) *
0096                        AngleAxis3{deltaPhi * istave, Vector3::UnitZ()} *
0097                        Translation3(Vector3::UnitX() * r) *
0098                        AngleAxis3{10_degree, Vector3::UnitZ()} *
0099                        AngleAxis3{90_degree, Vector3::UnitY()} *
0100                        AngleAxis3{90_degree, Vector3::UnitZ()};
0101       auto& element = elements.emplace_back(
0102           std::make_unique<DetectorElementStub>(trf, recBounds, thickness));
0103       element->surface().assignSurfacePlacement(*element);
0104       surfaces.push_back(element->surface().getSharedPtr());
0105     }
0106   }
0107
0108   return surfaces;
0109 }
0110
0111 }  // namespace ActsTests
0112
0113 using namespace ActsTests;
0114
0115 BOOST_AUTO_TEST_SUITE(GeometrySuite);
0116
0117 void pseudoNavigation(const TrackingGeometry& trackingGeometry,
0118                       Vector3 position, const Vector3& direction,
0119                       std::ostream& csv, std::size_t run,
0120                       std::size_t substepsPerCm, const Logger& logger) {
0121   ACTS_VERBOSE("start navigation " << run);
0122   ACTS_VERBOSE("dir: " << direction.transpose());
0123   ACTS_VERBOSE(direction.norm());
0124
0125   std::mt19937 rng{static_cast<unsigned int>(run)};
0126   std::uniform_real_distribution<> dist{0.01, 0.99};
0127
0128   const auto* volume = trackingGeometry.lowestTrackingVolume(gctx, position);
0129   BOOST_REQUIRE_NE(volume, nullptr);
0130   ACTS_VERBOSE(volume->volumeName());
0131
0132   NavigationStream main;
0133   const TrackingVolume* currentVolume = volume;
0134
0135   csv << run << "," << position[0] << "," << position[1] << "," << position[2];
0136   csv << "," << volume->geometryId().volume();
0137   csv << "," << volume->geometryId().boundary();
0138   csv << "," << volume->geometryId().sensitive();
0139   csv << std::endl;
0140
0141   ACTS_VERBOSE("start pseudo navigation");
0142
0143   for (std::size_t i = 0; i < 100; i++) {
0144     main = NavigationStream{};
0145     AppendOnlyNavigationStream stream{main};
0146
0147     NavigationArguments navArgs{.position = position, .direction = direction};
0148     NavigationPolicyStateManager stateManager;
0149     currentVolume->navigationPolicy()->createState(gctx, navArgs, stateManager,
0150                                                    logger);
0151     auto policyState = stateManager.currentState();
0152     currentVolume->initializeNavigationCandidates(gctx, navArgs, policyState,
0153                                                   stream, logger);
0154
0155     ACTS_VERBOSE(main.candidates().size() << " candidates");
0156
0157     for (const auto& candidate : main.candidates()) {
0158       ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0159       ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0160     }
0161
0162     ACTS_VERBOSE("initializing candidates");
0163     main.initialize(gctx, {position, direction}, BoundaryTolerance::None());
0164
0165     ACTS_VERBOSE(main.candidates().size() << " candidates remaining");
0166
0167     for (const auto& candidate : main.candidates()) {
0168       ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0169       ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0170     }
0171
0172     if (main.currentCandidate().surface().isOnSurface(gctx, position,
0173                                                       direction)) {
0174       ACTS_VERBOSE("Already on surface at initialization, skipping candidate");
0175
0176       auto id = main.currentCandidate().surface().geometryId();
0177       csv << run << "," << position[0] << "," << position[1] << ","
0178           << position[2];
0179       csv << "," << id.volume();
0180       csv << "," << id.boundary();
0181       csv << "," << id.sensitive();
0182       csv << std::endl;
0183       if (!main.switchToNextCandidate()) {
0184         ACTS_WARNING("candidates exhausted unexpectedly");
0185         break;
0186       }
0187     }
0188
0189     auto writeIntersection = [&](const Vector3& pos, const Surface& surface) {
0190       csv << run << "," << pos[0] << "," << pos[1] << "," << pos[2];
0191       csv << "," << surface.geometryId().volume();
0192       csv << "," << surface.geometryId().boundary();
0193       csv << "," << surface.geometryId().sensitive();
0194       csv << std::endl;
0195     };
0196
0197     bool terminated = false;
0198     while (main.remainingCandidates() > 0) {
0199       const auto& candidate = main.currentCandidate();
0200
0201       ACTS_VERBOSE(candidate.position().transpose());
0202
0203       ACTS_VERBOSE("moving to position: " << position.transpose() << " (r="
0204                                           << VectorHelpers::perp(position)
0205                                           << ")");
0206
0207       Vector3 delta = candidate.position() - position;
0208
0209       std::size_t substeps =
0210           std::max(1l, std::lround(delta.norm() / 10_cm * substepsPerCm));
0211
0212       for (std::size_t j = 0; j < substeps; j++) {
0213         // position += delta / (substeps + 1);
0214         Vector3 subpos = position + dist(rng) * delta;
0215         csv << run << "," << subpos[0] << "," << subpos[1] << "," << subpos[2];
0216         csv << "," << currentVolume->geometryId().volume();
0217         csv << ",0,0";  // zero boundary and sensitive ids
0218         csv << std::endl;
0219       }
0220
0221       position = candidate.position();
0222       ACTS_VERBOSE("                 -> "
0223                    << position.transpose()
0224                    << " (r=" << VectorHelpers::perp(position) << ")");
0225
0226       writeIntersection(position, candidate.surface());
0227
0228       if (candidate.isPortalTarget()) {
0229         ACTS_VERBOSE("On portal: " << candidate.surface().toStream(gctx));
0230         currentVolume =
0231             candidate.portal().resolveVolume(gctx, position, direction).value();
0232
0233         if (currentVolume == nullptr) {
0234           ACTS_VERBOSE("switched to nullptr -> we're done");
0235           terminated = true;
0236         }
0237         break;
0238
0239       } else {
0240         ACTS_VERBOSE("Not on portal");
0241       }
0242
0243       main.switchToNextCandidate();
0244     }
0245
0246     if (terminated) {
0247       ACTS_VERBOSE("Terminate pseudo navigation");
0248       break;
0249     }
0250
0251     ACTS_VERBOSE("switched to " << currentVolume->volumeName());
0252
0253     ACTS_VERBOSE("-----");
0254   }
0255 }
0256
0257 BOOST_AUTO_TEST_CASE(NodeApiTestContainers) {
0258   // Transform3 base{AngleAxis3{30_degree, Vector3{1, 0, 0}}};
0259   Transform3 base{Transform3::Identity()};
0260
0261   std::vector<std::unique_ptr<SurfacePlacementBase>> detectorElements;
0262   auto makeFan = [&](const Transform3& layerBase, auto&&..., double r,
0263                      std::size_t nSensors, double thickness) {
0264     return makeFanLayer(layerBase, detectorElements, r, nSensors, thickness);
0265   };
0266
0267   Blueprint::Config cfg;
0268   cfg.envelope[AxisDirection::AxisZ] = {20_mm, 20_mm};
0269   cfg.envelope[AxisDirection::AxisR] = {0_mm, 20_mm};
0270   auto root = std::make_unique<Blueprint>(cfg);
0271
0272   root->addMaterial("GlobalMaterial", [&](MaterialDesignatorBlueprintNode&
0273                                               mat) {
0274     using enum AxisDirection;
0275     using enum AxisBoundaryType;
0276     using enum CylinderVolumeBounds::Face;
0277
0278     // Configure cylinder faces with proper binning
0279     mat.configureFace(OuterCylinder, {AxisRPhi, Bound, 20}, {AxisZ, Bound, 20});
0280     mat.configureFace(NegativeDisc, {AxisR, Bound, 15}, {AxisPhi, Bound, 25});
0281     mat.configureFace(PositiveDisc, {AxisR, Bound, 15}, {AxisPhi, Bound, 25});
0282
0283     mat.addCylinderContainer("Detector", AxisDirection::AxisR, [&](auto& det) {
0284       det.addCylinderContainer("Pixel", AxisDirection::AxisZ, [&](auto& cyl) {
0285         cyl.setAttachmentStrategy(VolumeAttachmentStrategy::Gap)
0286             .setResizeStrategy(VolumeResizeStrategy::Gap);
0287
0288         cyl.addCylinderContainer(
0289             "PixelNegativeEndcap", AxisDirection::AxisZ, [&](auto& ec) {
0290               ec.setAttachmentStrategy(VolumeAttachmentStrategy::Gap);
0291
0292               auto makeLayer = [&](const Transform3& trf, auto& layer) {
0293                 std::vector<std::shared_ptr<Surface>> surfaces;
0294                 auto layerSurfaces = makeFan(trf, 300_mm, 10, 2_mm);
0295                 std::copy(layerSurfaces.begin(), layerSurfaces.end(),
0296                           std::back_inserter(surfaces));
0297                 layerSurfaces = makeFan(trf, 500_mm, 16, 2_mm);
0298                 std::copy(layerSurfaces.begin(), layerSurfaces.end(),
0299                           std::back_inserter(surfaces));
0300
0301                 layer.setSurfaces(surfaces)
0302                     .setLayerType(LayerBlueprintNode::LayerType::Disc)
0303                     .setEnvelope(ExtentEnvelope{{
0304                         .z = {5_mm, 5_mm},
0305                         .r = {10_mm, 20_mm},
0306                     }})
0307                     .setTransform(base);
0308               };
0309
0310               ec.addLayer("PixelNeg1", [&](auto& layer) {
0311                 makeLayer(base * Translation3{Vector3{0, 0, -700_mm}}, layer);
0312               });
0313
0314               ec.addLayer("PixelNeg2", [&](auto& layer) {
0315                 makeLayer(base * Translation3{Vector3{0, 0, -500_mm}}, layer);
0316               });
0317             });
0318
0319         cyl.addCylinderContainer(
0320             "PixelBarrel", AxisDirection::AxisR, [&](auto& brl) {
0321               brl.setAttachmentStrategy(VolumeAttachmentStrategy::Gap)
0322                   .setResizeStrategy(VolumeResizeStrategy::Gap);
0323
0324               auto makeLayer = [&](const std::string& name, double r,
0325                                    std::size_t nStaves, int nSensorsPerStave) {
0326                 brl.addLayer(name, [&](auto& layer) {
0327                   std::vector<std::shared_ptr<Surface>> surfaces =
0328                       makeBarrelLayer(base, detectorElements, r, nStaves,
0329                                       nSensorsPerStave, 2.5_mm, 10_mm, 20_mm);
0330
0331                   layer.setSurfaces(surfaces)
0332                       .setLayerType(LayerBlueprintNode::LayerType::Cylinder)
0333                       .setEnvelope(ExtentEnvelope{{
0334                           .z = {5_mm, 5_mm},
0335                           .r = {1_mm, 1_mm},
0336                       }})
0337                       .setTransform(base);
0338                 });
0339               };
0340
0341               makeLayer("PixelLayer0", 30_mm, 18, 5);
0342               makeLayer("PixelLayer1", 90_mm, 30, 6);
0343
0344               brl.addStaticVolume(base,
0345                                   std::make_shared<CylinderVolumeBounds>(
0346                                       100_mm, 110_mm, 250_mm),
0347                                   "PixelSupport");
0348
0349               makeLayer("PixelLayer2", 150_mm, 40, 7);
0350               makeLayer("PixelLayer3", 250_mm, 70, 8);
0351             });
0352
0353         auto& ec =
0354             cyl.addCylinderContainer("PixelPosWrapper", AxisDirection::AxisR);
0355         ec.setResizeStrategy(VolumeResizeStrategy::Gap);
0356         ec.addStaticVolume(std::make_unique<TrackingVolume>(
0357             base * Translation3{Vector3{0, 0, 600_mm}},
0358             std::make_shared<CylinderVolumeBounds>(150_mm, 390_mm, 200_mm),
0359             "PixelPositiveEndcap"));
0360       });
0361
0362       det.addStaticVolume(
0363           base, std::make_shared<CylinderVolumeBounds>(0_mm, 23_mm, 1000_mm),
0364           "BeamPipe");
0365     });
0366   });
0367
0368   std::ofstream dot{"api_test_container.dot"};
0369   root->graphviz(dot);
0370
0371   auto trackingGeometry = root->construct({}, gctx, *logger);
0372
0373   BOOST_REQUIRE(trackingGeometry);
0374   BOOST_CHECK(trackingGeometry->geometryVersion() ==
0375               TrackingGeometry::GeometryVersion::Gen3);
0376
0377   trackingGeometry->visitVolumes([&](const TrackingVolume* volume) {
0378     std::cout << volume->volumeName() << std::endl;
0379     std::cout << " -> id: " << volume->geometryId() << std::endl;
0380     std::cout << " -> " << volume->portals().size() << " portals" << std::endl;
0381   });
0382
0383   ObjVisualization3D vis;
0384
0385   trackingGeometry->visualize(vis, gctx, {}, {});
0386
0387   vis.write("api_test_container.obj");
0388
0389   Vector3 position = Vector3::Zero();
0390   std::ofstream csv{"api_test_container.csv"};
0391   csv << "x,y,z,volume,boundary,sensitive" << std::endl;
0392
0393   std::mt19937 rnd{42};
0394
0395   std::uniform_real_distribution<> dist{-1, 1};
0396
0397   double etaWidth = 3;
0398   double thetaMin = 2 * std::atan(std::exp(-etaWidth));
0399   double thetaMax = 2 * std::atan(std::exp(etaWidth));
0400   std::uniform_real_distribution<> thetaDist{thetaMin, thetaMax};
0401
0402   using namespace UnitLiterals;
0403
0404   for (std::size_t i = 0; i < 5000; i++) {
0405     double theta = thetaDist(rnd);
0406     double phi = 2 * std::numbers::pi * dist(rnd);
0407
0408     Vector3 direction;
0409     direction[0] = std::sin(theta) * std::cos(phi);
0410     direction[1] = std::sin(theta) * std::sin(phi);
0411     direction[2] = std::cos(theta);
0412
0413     pseudoNavigation(*trackingGeometry, position, direction, csv, i, 2,
0414                      *logger->clone(std::nullopt, Logging::DEBUG));
0415   }
0416 }
0417
0418 BOOST_AUTO_TEST_CASE(NodeApiTestCuboid) {
0419   Transform3 base{Transform3::Identity()};
0420
0421   Blueprint::Config cfg;
0422   cfg.envelope[AxisDirection::AxisZ] = {20_mm, 20_mm};
0423   cfg.envelope[AxisDirection::AxisR] = {0_mm, 20_mm};
0424   auto root = std::make_unique<Blueprint>(cfg);
0425
0426   root->addMaterial("GlobalMaterial", [&](MaterialDesignatorBlueprintNode&
0427                                               mat) {
0428     using enum AxisDirection;
0429     using enum AxisBoundaryType;
0430     using enum CuboidVolumeBounds::Face;
0431
0432     // Configure valid axis combinations for each face type
0433     mat.configureFace(NegativeXFace, {AxisX, Bound, 20}, {AxisY, Bound, 20});
0434     mat.configureFace(PositiveXFace, {AxisX, Bound, 20}, {AxisY, Bound, 20});
0435     mat.configureFace(NegativeYFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0436     mat.configureFace(PositiveYFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0437     mat.configureFace(NegativeZFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0438     mat.configureFace(PositiveZFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0439
0440     mat.addStaticVolume(
0441         base, std::make_shared<CuboidVolumeBounds>(100_mm, 100_mm, 100_mm),
0442         "TestVolume");
0443   });
0444
0445   auto trackingGeometry = root->construct({}, gctx, *logger);
0446   BOOST_REQUIRE(trackingGeometry);
0447   BOOST_CHECK(trackingGeometry->geometryVersion() ==
0448               TrackingGeometry::GeometryVersion::Gen3);
0449 }
0450
0451 BOOST_AUTO_TEST_SUITE_END();