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