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