Core/Surfaces/SurfaceArrayTests.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/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Surfaces/CylinderSurface.hpp"
0014 #include "Acts/Surfaces/PlanarBounds.hpp"
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/Surfaces/RectangleBounds.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Surfaces/SurfaceArray.hpp"
0019 #include "Acts/Utilities/Axis.hpp"
0020 #include "Acts/Utilities/AxisDefinitions.hpp"
0021 #include "Acts/Utilities/Helpers.hpp"
0022
0023 #include <algorithm>
0024 #include <cmath>
0025 #include <cstddef>
0026 #include <fstream>
0027 #include <iomanip>
0028 #include <iostream>
0029 #include <memory>
0030 #include <numbers>
0031 #include <sstream>
0032 #include <string>
0033 #include <tuple>
0034 #include <utility>
0035 #include <vector>
0036
0037 #include <boost/format.hpp>
0038
0039 using Acts::VectorHelpers::phi;
0040
0041 using namespace Acts;
0042
0043 namespace ActsTests {
0044
0045 // Create a test context
0046 GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0047
0048 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
0049 struct SurfaceArrayFixture {
0050   std::vector<std::shared_ptr<const Surface>> m_surfaces;
0051
0052   SurfaceArrayFixture() { BOOST_TEST_MESSAGE("setup fixture"); }
0053   ~SurfaceArrayFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
0054
0055   SrfVec fullPhiTestSurfacesEC(std::size_t n = 10, double shift = 0,
0056                                double zbase = 0, double r = 10) {
0057     SrfVec res;
0058
0059     double phiStep = 2 * std::numbers::pi / n;
0060     for (std::size_t i = 0; i < n; ++i) {
0061       double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
0062
0063       Transform3 trans;
0064       trans.setIdentity();
0065       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0066       trans.translate(Vector3(r, 0, z));
0067
0068       auto bounds = std::make_shared<const RectangleBounds>(2, 1);
0069       std::shared_ptr<const Surface> srf =
0070           Surface::makeShared<PlaneSurface>(trans, bounds);
0071
0072       res.push_back(srf);
0073       m_surfaces.push_back(
0074           std::move(srf));  // keep shared, will get destroyed at the end
0075     }
0076
0077     return res;
0078   }
0079
0080   SrfVec fullPhiTestSurfacesBRL(int n = 10, double shift = 0, double zbase = 0,
0081                                 double incl = std::numbers::pi / 9.,
0082                                 double w = 2, double h = 1.5) {
0083     SrfVec res;
0084
0085     double phiStep = 2 * std::numbers::pi / n;
0086     for (int i = 0; i < n; ++i) {
0087       double z = zbase;
0088
0089       Transform3 trans;
0090       trans.setIdentity();
0091       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0092       trans.translate(Vector3(10, 0, z));
0093       trans.rotate(Eigen::AngleAxisd(incl, Vector3(0, 0, 1)));
0094       trans.rotate(Eigen::AngleAxisd(std::numbers::pi / 2., Vector3(0, 1, 0)));
0095
0096       auto bounds = std::make_shared<const RectangleBounds>(w, h);
0097       std::shared_ptr<const Surface> srf =
0098           Surface::makeShared<PlaneSurface>(trans, bounds);
0099
0100       res.push_back(srf);
0101       m_surfaces.push_back(
0102           std::move(srf));  // keep shared, will get destroyed at the end
0103     }
0104
0105     return res;
0106   }
0107
0108   SrfVec straightLineSurfaces(
0109       std::size_t n = 10., double step = 3, const Vector3& origin = {0, 0, 1.5},
0110       const Transform3& pretrans = Transform3::Identity(),
0111       const Vector3& dir = {0, 0, 1}) {
0112     SrfVec res;
0113     for (std::size_t i = 0; i < n; ++i) {
0114       Transform3 trans;
0115       trans.setIdentity();
0116       trans.translate(origin + dir * step * i);
0117       // trans.rotate(AngleAxis3(std::numbers::pi/9., Vector3(0, 0, 1)));
0118       trans.rotate(AngleAxis3(std::numbers::pi / 2., Vector3(1, 0, 0)));
0119       trans = trans * pretrans;
0120
0121       auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
0122
0123       std::shared_ptr<const Surface> srf =
0124           Surface::makeShared<PlaneSurface>(trans, bounds);
0125
0126       res.push_back(srf);
0127       m_surfaces.push_back(
0128           std::move(srf));  // keep shared, will get destroyed at the end
0129     }
0130
0131     return res;
0132   }
0133
0134   SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
0135     double z0 = -(nZ - 1) * w;
0136     SrfVec res;
0137
0138     for (int i = 0; i < nZ; i++) {
0139       double z = i * w * 2 + z0;
0140       SrfVec ring =
0141           fullPhiTestSurfacesBRL(nPhi, 0, z, std::numbers::pi / 9., w, h);
0142       res.insert(res.end(), ring.begin(), ring.end());
0143     }
0144
0145     return res;
0146   }
0147
0148   void draw_surfaces(const SrfVec& surfaces, const std::string& fname) {
0149     std::ofstream os;
0150     os.open(fname);
0151
0152     os << std::fixed << std::setprecision(4);
0153
0154     std::size_t nVtx = 0;
0155     for (const auto& srfx : surfaces) {
0156       std::shared_ptr<const PlaneSurface> srf =
0157           std::dynamic_pointer_cast<const PlaneSurface>(srfx);
0158       const PlanarBounds* bounds =
0159           dynamic_cast<const PlanarBounds*>(&srf->bounds());
0160
0161       for (const auto& vtxloc : bounds->vertices()) {
0162         Vector3 vtx = srf->localToGlobalTransform(tgContext) *
0163                       Vector3(vtxloc.x(), vtxloc.y(), 0);
0164         os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0165       }
0166
0167       // connect them
0168       os << "f";
0169       for (std::size_t i = 1; i <= bounds->vertices().size(); ++i) {
0170         os << " " << nVtx + i;
0171       }
0172       os << "\n";
0173
0174       nVtx += bounds->vertices().size();
0175     }
0176
0177     os.close();
0178   }
0179 };
0180
0181 BOOST_AUTO_TEST_SUITE(SurfacesSuite)
0182
0183 BOOST_FIXTURE_TEST_CASE(SurfaceArray_create, SurfaceArrayFixture) {
0184   GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0185
0186   SrfVec brl = makeBarrel(30, 7, 2, 1);
0187   std::vector<const Surface*> brlRaw = unpackSmartPointers(brl);
0188   draw_surfaces(brl, "SurfaceArray_create_BRL_1.obj");
0189
0190   Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0191       -std::numbers::pi, std::numbers::pi, 30u);
0192   Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-14, 14, 7u);
0193
0194   double R = 10;
0195   auto itransform = [R](const Vector2& loc) {
0196     return Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0197   };
0198
0199   auto cylinder =
0200       Surface::makeShared<CylinderSurface>(Transform3::Identity(), R, 10);
0201   SurfaceArray sa(tgContext, brl, cylinder, 1., std::tuple{phiAxis, zAxis});
0202
0203   // let's see if we can access all surfaces
0204   sa.toStream(tgContext, std::cout);
0205
0206   // the bin a surface's own centre falls into has to hold it
0207   for (const auto& srf : brl) {
0208     const Vector3 ctr = srf->referencePosition(tgContext, AxisDirection::AxisR);
0209     const Vector3 normal = srf->normal(tgContext, ctr, Vector3::UnitZ());
0210     const auto binContent = sa.at(tgContext, ctr, normal);
0211
0212     BOOST_CHECK(std::ranges::find(binContent, srf.get()) != binContent.end());
0213   }
0214
0215   const auto neighbors = sa.neighbors(tgContext, itransform(Vector2(0, 0)),
0216                                       itransform(Vector2(0, 0)).normalized());
0217   BOOST_CHECK_EQUAL(neighbors.size(), 16u);
0218 }
0219
0220 BOOST_AUTO_TEST_CASE(SurfaceArray_singleElement) {
0221   const double w = 3;
0222   const double h = 4;
0223   const auto bounds = std::make_shared<const RectangleBounds>(w, h);
0224   auto srf = Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0225
0226   SurfaceArray sa(srf);
0227
0228   const auto binContent =
0229       sa.at(tgContext, Vector3(42, 42, 42), Vector3::UnitX());
0230   BOOST_CHECK_EQUAL(binContent.size(), 1u);
0231   BOOST_CHECK_EQUAL(binContent[0], srf.get());
0232   BOOST_CHECK_EQUAL(sa.surfaces().size(), 1u);
0233   BOOST_CHECK_EQUAL(sa.surfaces().at(0), srf.get());
0234 }
0235
0236 BOOST_AUTO_TEST_CASE(SurfaceArrayToStreamPreservesStreamState) {
0237   const auto bounds = std::make_shared<const RectangleBounds>(3., 4.);
0238   auto surface =
0239       Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0240   SurfaceArray surfaceArray(surface);
0241
0242   std::ostringstream stream;
0243   stream << std::scientific << std::showpos << std::setfill('#')
0244          << std::setprecision(3);
0245   stream.width(17);
0246
0247   const auto flags = stream.flags();
0248   const auto precision = stream.precision();
0249   const auto width = stream.width();
0250   const auto fill = stream.fill();
0251
0252   surfaceArray.toStream(tgContext, stream);
0253
0254   BOOST_CHECK(stream.flags() == flags);
0255   BOOST_CHECK_EQUAL(stream.precision(), precision);
0256   BOOST_CHECK_EQUAL(stream.width(), width);
0257   BOOST_CHECK_EQUAL(stream.fill(), fill);
0258 }
0259
0260 BOOST_AUTO_TEST_SUITE_END()
0261
0262 }  // namespace ActsTests