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