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 GeometryContext gctx;
0051
0052 inline std::vector<std::shared_ptr<Surface>> makeFanLayer(
0053     const Transform3& base,
0054     std::vector<std::unique_ptr<DetectorElementBase>>& 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().assignDetectorElement(*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<DetectorElementBase>>& 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().assignDetectorElement(*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     currentVolume->initializeNavigationCandidates(
0148         {.position = position, .direction = direction}, stream, logger);
0149
0150     ACTS_VERBOSE(main.candidates().size() << " candidates");
0151
0152     for (const auto& candidate : main.candidates()) {
0153       ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0154       ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0155     }
0156
0157     ACTS_VERBOSE("initializing candidates");
0158     main.initialize(gctx, {position, direction}, BoundaryTolerance::None());
0159
0160     ACTS_VERBOSE(main.candidates().size() << " candidates remaining");
0161
0162     for (const auto& candidate : main.candidates()) {
0163       ACTS_VERBOSE(" -> " << candidate.surface().geometryId());
0164       ACTS_VERBOSE("    " << candidate.surface().toStream(gctx));
0165     }
0166
0167     if (main.currentCandidate().surface().isOnSurface(gctx, position,
0168                                                       direction)) {
0169       ACTS_VERBOSE("Already on surface at initialization, skipping candidate");
0170
0171       auto id = main.currentCandidate().surface().geometryId();
0172       csv << run << "," << position[0] << "," << position[1] << ","
0173           << position[2];
0174       csv << "," << id.volume();
0175       csv << "," << id.boundary();
0176       csv << "," << id.sensitive();
0177       csv << std::endl;
0178       if (!main.switchToNextCandidate()) {
0179         ACTS_WARNING("candidates exhausted unexpectedly");
0180         break;
0181       }
0182     }
0183
0184     auto writeIntersection = [&](const Vector3& pos, const Surface& surface) {
0185       csv << run << "," << pos[0] << "," << pos[1] << "," << pos[2];
0186       csv << "," << surface.geometryId().volume();
0187       csv << "," << surface.geometryId().boundary();
0188       csv << "," << surface.geometryId().sensitive();
0189       csv << std::endl;
0190     };
0191
0192     bool terminated = false;
0193     while (main.remainingCandidates() > 0) {
0194       const auto& candidate = main.currentCandidate();
0195
0196       ACTS_VERBOSE(candidate.position().transpose());
0197
0198       ACTS_VERBOSE("moving to position: " << position.transpose() << " (r="
0199                                           << VectorHelpers::perp(position)
0200                                           << ")");
0201
0202       Vector3 delta = candidate.position() - position;
0203
0204       std::size_t substeps =
0205           std::max(1l, std::lround(delta.norm() / 10_cm * substepsPerCm));
0206
0207       for (std::size_t j = 0; j < substeps; j++) {
0208         // position += delta / (substeps + 1);
0209         Vector3 subpos = position + dist(rng) * delta;
0210         csv << run << "," << subpos[0] << "," << subpos[1] << "," << subpos[2];
0211         csv << "," << currentVolume->geometryId().volume();
0212         csv << ",0,0";  // zero boundary and sensitive ids
0213         csv << std::endl;
0214       }
0215
0216       position = candidate.position();
0217       ACTS_VERBOSE("                 -> "
0218                    << position.transpose()
0219                    << " (r=" << VectorHelpers::perp(position) << ")");
0220
0221       writeIntersection(position, candidate.surface());
0222
0223       if (candidate.isPortalTarget()) {
0224         ACTS_VERBOSE("On portal: " << candidate.surface().toStream(gctx));
0225         currentVolume =
0226             candidate.portal().resolveVolume(gctx, position, direction).value();
0227
0228         if (currentVolume == nullptr) {
0229           ACTS_VERBOSE("switched to nullptr -> we're done");
0230           terminated = true;
0231         }
0232         break;
0233
0234       } else {
0235         ACTS_VERBOSE("Not on portal");
0236       }
0237
0238       main.switchToNextCandidate();
0239     }
0240
0241     if (terminated) {
0242       ACTS_VERBOSE("Terminate pseudo navigation");
0243       break;
0244     }
0245
0246     ACTS_VERBOSE("switched to " << currentVolume->volumeName());
0247
0248     ACTS_VERBOSE("-----");
0249   }
0250 }
0251
0252 BOOST_AUTO_TEST_CASE(NodeApiTestContainers) {
0253   // Transform3 base{AngleAxis3{30_degree, Vector3{1, 0, 0}}};
0254   Transform3 base{Transform3::Identity()};
0255
0256   std::vector<std::unique_ptr<DetectorElementBase>> detectorElements;
0257   auto makeFan = [&](const Transform3& layerBase, auto&&..., double r,
0258                      std::size_t nSensors, double thickness) {
0259     return makeFanLayer(layerBase, detectorElements, r, nSensors, thickness);
0260   };
0261
0262   Blueprint::Config cfg;
0263   cfg.envelope[AxisDirection::AxisZ] = {20_mm, 20_mm};
0264   cfg.envelope[AxisDirection::AxisR] = {0_mm, 20_mm};
0265   auto root = std::make_unique<Blueprint>(cfg);
0266
0267   root->addMaterial("GlobalMaterial", [&](MaterialDesignatorBlueprintNode&
0268                                               mat) {
0269     using enum AxisDirection;
0270     using enum AxisBoundaryType;
0271     using enum CylinderVolumeBounds::Face;
0272
0273     // Configure cylinder faces with proper binning
0274     mat.configureFace(OuterCylinder, {AxisRPhi, Bound, 20}, {AxisZ, Bound, 20});
0275     mat.configureFace(NegativeDisc, {AxisR, Bound, 15}, {AxisPhi, Bound, 25});
0276     mat.configureFace(PositiveDisc, {AxisR, Bound, 15}, {AxisPhi, Bound, 25});
0277
0278     mat.addCylinderContainer("Detector", AxisDirection::AxisR, [&](auto& det) {
0279       det.addCylinderContainer("Pixel", AxisDirection::AxisZ, [&](auto& cyl) {
0280         cyl.setAttachmentStrategy(VolumeAttachmentStrategy::Gap)
0281             .setResizeStrategy(VolumeResizeStrategy::Gap);
0282
0283         cyl.addCylinderContainer(
0284             "PixelNegativeEndcap", AxisDirection::AxisZ, [&](auto& ec) {
0285               ec.setAttachmentStrategy(VolumeAttachmentStrategy::Gap);
0286
0287               auto makeLayer = [&](const Transform3& trf, auto& layer) {
0288                 std::vector<std::shared_ptr<Surface>> surfaces;
0289                 auto layerSurfaces = makeFan(trf, 300_mm, 10, 2_mm);
0290                 std::copy(layerSurfaces.begin(), layerSurfaces.end(),
0291                           std::back_inserter(surfaces));
0292                 layerSurfaces = makeFan(trf, 500_mm, 16, 2_mm);
0293                 std::copy(layerSurfaces.begin(), layerSurfaces.end(),
0294                           std::back_inserter(surfaces));
0295
0296                 layer.setSurfaces(surfaces)
0297                     .setLayerType(LayerBlueprintNode::LayerType::Disc)
0298                     .setEnvelope(ExtentEnvelope{{
0299                         .z = {5_mm, 5_mm},
0300                         .r = {10_mm, 20_mm},
0301                     }})
0302                     .setTransform(base);
0303               };
0304
0305               ec.addLayer("PixelNeg1", [&](auto& layer) {
0306                 makeLayer(base * Translation3{Vector3{0, 0, -700_mm}}, layer);
0307               });
0308
0309               ec.addLayer("PixelNeg2", [&](auto& layer) {
0310                 makeLayer(base * Translation3{Vector3{0, 0, -500_mm}}, layer);
0311               });
0312             });
0313
0314         cyl.addCylinderContainer(
0315             "PixelBarrel", AxisDirection::AxisR, [&](auto& brl) {
0316               brl.setAttachmentStrategy(VolumeAttachmentStrategy::Gap)
0317                   .setResizeStrategy(VolumeResizeStrategy::Gap);
0318
0319               auto makeLayer = [&](const std::string& name, double r,
0320                                    std::size_t nStaves, int nSensorsPerStave) {
0321                 brl.addLayer(name, [&](auto& layer) {
0322                   std::vector<std::shared_ptr<Surface>> surfaces =
0323                       makeBarrelLayer(base, detectorElements, r, nStaves,
0324                                       nSensorsPerStave, 2.5_mm, 10_mm, 20_mm);
0325
0326                   layer.setSurfaces(surfaces)
0327                       .setLayerType(LayerBlueprintNode::LayerType::Cylinder)
0328                       .setEnvelope(ExtentEnvelope{{
0329                           .z = {5_mm, 5_mm},
0330                           .r = {1_mm, 1_mm},
0331                       }})
0332                       .setTransform(base);
0333                 });
0334               };
0335
0336               makeLayer("PixelLayer0", 30_mm, 18, 5);
0337               makeLayer("PixelLayer1", 90_mm, 30, 6);
0338
0339               brl.addStaticVolume(base,
0340                                   std::make_shared<CylinderVolumeBounds>(
0341                                       100_mm, 110_mm, 250_mm),
0342                                   "PixelSupport");
0343
0344               makeLayer("PixelLayer2", 150_mm, 40, 7);
0345               makeLayer("PixelLayer3", 250_mm, 70, 8);
0346             });
0347
0348         auto& ec =
0349             cyl.addCylinderContainer("PixelPosWrapper", AxisDirection::AxisR);
0350         ec.setResizeStrategy(VolumeResizeStrategy::Gap);
0351         ec.addStaticVolume(std::make_unique<TrackingVolume>(
0352             base * Translation3{Vector3{0, 0, 600_mm}},
0353             std::make_shared<CylinderVolumeBounds>(150_mm, 390_mm, 200_mm),
0354             "PixelPositiveEndcap"));
0355       });
0356
0357       det.addStaticVolume(
0358           base, std::make_shared<CylinderVolumeBounds>(0_mm, 23_mm, 1000_mm),
0359           "BeamPipe");
0360     });
0361   });
0362
0363   std::ofstream dot{"api_test_container.dot"};
0364   root->graphviz(dot);
0365
0366   auto trackingGeometry = root->construct({}, gctx, *logger);
0367
0368   BOOST_REQUIRE(trackingGeometry);
0369   BOOST_CHECK(trackingGeometry->geometryVersion() ==
0370               TrackingGeometry::GeometryVersion::Gen3);
0371
0372   trackingGeometry->visitVolumes([&](const TrackingVolume* volume) {
0373     std::cout << volume->volumeName() << std::endl;
0374     std::cout << " -> id: " << volume->geometryId() << std::endl;
0375     std::cout << " -> " << volume->portals().size() << " portals" << std::endl;
0376   });
0377
0378   ObjVisualization3D vis;
0379
0380   trackingGeometry->visualize(vis, gctx, {}, {});
0381
0382   vis.write("api_test_container.obj");
0383
0384   Vector3 position = Vector3::Zero();
0385   std::ofstream csv{"api_test_container.csv"};
0386   csv << "x,y,z,volume,boundary,sensitive" << std::endl;
0387
0388   std::mt19937 rnd{42};
0389
0390   std::uniform_real_distribution<> dist{-1, 1};
0391
0392   double etaWidth = 3;
0393   double thetaMin = 2 * std::atan(std::exp(-etaWidth));
0394   double thetaMax = 2 * std::atan(std::exp(etaWidth));
0395   std::uniform_real_distribution<> thetaDist{thetaMin, thetaMax};
0396
0397   using namespace UnitLiterals;
0398
0399   for (std::size_t i = 0; i < 5000; i++) {
0400     double theta = thetaDist(rnd);
0401     double phi = 2 * std::numbers::pi * dist(rnd);
0402
0403     Vector3 direction;
0404     direction[0] = std::sin(theta) * std::cos(phi);
0405     direction[1] = std::sin(theta) * std::sin(phi);
0406     direction[2] = std::cos(theta);
0407
0408     pseudoNavigation(*trackingGeometry, position, direction, csv, i, 2,
0409                      *logger->clone(std::nullopt, Logging::DEBUG));
0410   }
0411 }
0412
0413 BOOST_AUTO_TEST_CASE(NodeApiTestCuboid) {
0414   Transform3 base{Transform3::Identity()};
0415
0416   Blueprint::Config cfg;
0417   cfg.envelope[AxisDirection::AxisZ] = {20_mm, 20_mm};
0418   cfg.envelope[AxisDirection::AxisR] = {0_mm, 20_mm};
0419   auto root = std::make_unique<Blueprint>(cfg);
0420
0421   root->addMaterial("GlobalMaterial", [&](MaterialDesignatorBlueprintNode&
0422                                               mat) {
0423     using enum AxisDirection;
0424     using enum AxisBoundaryType;
0425     using enum CuboidVolumeBounds::Face;
0426
0427     // Configure valid axis combinations for each face type
0428     mat.configureFace(NegativeXFace, {AxisX, Bound, 20}, {AxisY, Bound, 20});
0429     mat.configureFace(PositiveXFace, {AxisX, Bound, 20}, {AxisY, Bound, 20});
0430     mat.configureFace(NegativeYFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0431     mat.configureFace(PositiveYFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0432     mat.configureFace(NegativeZFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0433     mat.configureFace(PositiveZFace, {AxisX, Bound, 15}, {AxisY, Bound, 25});
0434
0435     mat.addStaticVolume(
0436         base, std::make_shared<CuboidVolumeBounds>(100_mm, 100_mm, 100_mm),
0437         "TestVolume");
0438   });
0439
0440   auto trackingGeometry = root->construct({}, gctx, *logger);
0441   BOOST_REQUIRE(trackingGeometry);
0442   BOOST_CHECK(trackingGeometry->geometryVersion() ==
0443               TrackingGeometry::GeometryVersion::Gen3);
0444 }
0445
0446 BOOST_AUTO_TEST_SUITE_END();