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