Fatras/Digitization/SegmentizerTests.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/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Surfaces/DiscBounds.hpp"
0015 #include "Acts/Surfaces/DiscSurface.hpp"
0016 #include "Acts/Surfaces/PlanarBounds.hpp"
0017 #include "Acts/Surfaces/PlaneSurface.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/RectangleBounds.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Utilities/IAxis.hpp"
0022 #include "Acts/Utilities/IMultiAxis.hpp"
0023 #include "ActsFatras/Digitization/Segmentizer.hpp"
0024
0025 #include <cmath>
0026 #include <fstream>
0027 #include <memory>
0028 #include <string>
0029 #include <tuple>
0030 #include <utility>
0031 #include <vector>
0032
0033 #include "DigitizationCsvOutput.hpp"
0034 #include "PlanarSurfaceTestBeds.hpp"
0035
0036 namespace bdata = boost::unit_test::data;
0037
0038 using namespace Acts;
0039 using namespace ActsFatras;
0040
0041 namespace ActsTests {
0042
0043 BOOST_AUTO_TEST_SUITE(DigitizationSuite)
0044
0045 BOOST_AUTO_TEST_CASE(SegmentizerCartesian) {
0046   auto geoCtx = GeometryContext::dangerouslyDefaultConstruct();
0047
0048   auto rectangleBounds = std::make_shared<RectangleBounds>(1., 1.);
0049   auto planeSurface = Surface::makeShared<PlaneSurface>(Transform3::Identity(),
0050                                                         rectangleBounds);
0051
0052   // The segmentation
0053   const auto segmentation = IMultiAxis::create(
0054       *IAxis::createEquidistant(AxisBoundaryType::Bound, -1., 1., 20,
0055                                 AxisDirection::AxisX),
0056       *IAxis::createEquidistant(AxisBoundaryType::Bound, -1., 1., 20,
0057                                 AxisDirection::AxisY));
0058
0059   Segmentizer cl;
0060
0061   // Test: Normal hit into the surface
0062   Vector2 nPosition(0.37, 0.76);
0063   auto nSegments =
0064       cl.segments(geoCtx, *planeSurface, *segmentation, {nPosition, nPosition});
0065   BOOST_CHECK_EQUAL(nSegments.size(), 1);
0066   BOOST_CHECK_EQUAL(nSegments[0].bin[0], 13);
0067   BOOST_CHECK_EQUAL(nSegments[0].bin[1], 17);
0068
0069   // Test: Inclined hit into the surface - negative x direction
0070   Vector2 ixPositionS(0.37, 0.76);
0071   Vector2 ixPositionE(0.02, 0.73);
0072   auto ixSegments = cl.segments(geoCtx, *planeSurface, *segmentation,
0073                                 {ixPositionS, ixPositionE});
0074   BOOST_CHECK_EQUAL(ixSegments.size(), 4);
0075
0076   // Test: Inclined hit into the surface - positive y direction
0077   Vector2 iyPositionS(0.37, 0.76);
0078   Vector2 iyPositionE(0.39, 0.91);
0079   auto iySegments = cl.segments(geoCtx, *planeSurface, *segmentation,
0080                                 {iyPositionS, iyPositionE});
0081   BOOST_CHECK_EQUAL(iySegments.size(), 3);
0082
0083   // Test: Inclined hit into the surface - x/y direction
0084   Vector2 ixyPositionS(-0.27, 0.76);
0085   Vector2 ixyPositionE(-0.02, -0.73);
0086   auto ixySegments = cl.segments(geoCtx, *planeSurface, *segmentation,
0087                                  {ixyPositionS, ixyPositionE});
0088   BOOST_CHECK_EQUAL(ixySegments.size(), 18);
0089 }
0090
0091 BOOST_AUTO_TEST_CASE(SegmentizerPolarRadial) {
0092   auto geoCtx = GeometryContext::dangerouslyDefaultConstruct();
0093
0094   auto radialBounds = std::make_shared<const RadialBounds>(5., 10., 0.25, 0.);
0095   auto radialDisc =
0096       Surface::makeShared<DiscSurface>(Transform3::Identity(), radialBounds);
0097
0098   // The segmentation
0099   const auto segmentation = IMultiAxis::create(
0100       *IAxis::createEquidistant(AxisBoundaryType::Bound, 5., 10., 2,
0101                                 AxisDirection::AxisR),
0102       *IAxis::createEquidistant(AxisBoundaryType::Bound, -0.25, 0.25, 250,
0103                                 AxisDirection::AxisPhi));
0104
0105   Segmentizer cl;
0106
0107   // Test: Normal hit into the surface
0108   Vector2 nPosition(6.76, 0.5);
0109   auto nSegments =
0110       cl.segments(geoCtx, *radialDisc, *segmentation, {nPosition, nPosition});
0111   BOOST_CHECK_EQUAL(nSegments.size(), 1);
0112   BOOST_CHECK_EQUAL(nSegments[0].bin[0], 0);
0113   BOOST_CHECK_EQUAL(nSegments[0].bin[1], 161);
0114
0115   // Test: now opver more phi strips
0116   Vector2 sPositionS(6.76, 0.5);
0117   Vector2 sPositionE(7.03, -0.3);
0118   auto sSegment =
0119       cl.segments(geoCtx, *radialDisc, *segmentation, {sPositionS, sPositionE});
0120   BOOST_CHECK_EQUAL(sSegment.size(), 59);
0121
0122   // Test: jump over R boundary, but stay in phi bin
0123   sPositionS = Vector2(6.76, 0.);
0124   sPositionE = Vector2(7.83, 0.);
0125   sSegment =
0126       cl.segments(geoCtx, *radialDisc, *segmentation, {sPositionS, sPositionE});
0127   BOOST_CHECK_EQUAL(sSegment.size(), 2);
0128 }
0129
0130 /// Unit test for testing the Segmentizer
0131 BOOST_DATA_TEST_CASE(
0132     RandomSegmentizerTest,
0133     bdata::random((
0134         bdata::engine = std::mt19937(), bdata::seed = 1,
0135         bdata::distribution = std::uniform_real_distribution<double>(0., 1.))) ^
0136         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 2,
0137                        bdata::distribution =
0138                            std::uniform_real_distribution<double>(0., 1.))) ^
0139         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0140                        bdata::distribution =
0141                            std::uniform_real_distribution<double>(0., 1.))) ^
0142         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0143                        bdata::distribution =
0144                            std::uniform_real_distribution<double>(0., 1.))) ^
0145         bdata::xrange(25),
0146     startR0, startR1, endR0, endR1, index) {
0147   auto geoCtx = GeometryContext::dangerouslyDefaultConstruct();
0148   Segmentizer cl;
0149
0150   // Test beds with random numbers generated inside
0151   PlanarSurfaceTestBeds pstd;
0152   auto testBeds = pstd(1.);
0153
0154   DigitizationCsvOutput csvHelper;
0155
0156   for (const auto& tb : testBeds) {
0157     const auto& name = std::get<0>(tb);
0158     const auto* surface = (std::get<1>(tb)).get();
0159     const auto& segmentation = std::get<2>(tb);
0160     const auto& randomizer = std::get<3>(tb);
0161
0162     if (index == 0) {
0163       std::ofstream shape;
0164       std::ofstream grid;
0165       const auto centerXY = surface->center(geoCtx).segment<2>(0);
0166       // 0 - write the shape
0167       shape.open("Segmentizer" + name + "Borders.csv");
0168       if (surface->type() == Surface::Plane) {
0169         const auto* pBounds =
0170             static_cast<const PlanarBounds*>(&(surface->bounds()));
0171         csvHelper.writePolygon(shape, pBounds->vertices(1), -centerXY);
0172       } else if (surface->type() == Surface::Disc) {
0173         const auto* dBounds =
0174             static_cast<const DiscBounds*>(&(surface->bounds()));
0175         csvHelper.writePolygon(shape, dBounds->vertices(72), -centerXY);
0176       }
0177       // 1 - write the grid
0178       grid.open("Segmentizer" + name + "Grid.csv");
0179       const IAxis& axis0 = segmentation->getAxis(0);
0180       const IAxis& axis1 = segmentation->getAxis(1);
0181       if (axis0.getDirection() == AxisDirection::AxisX &&
0182           axis1.getDirection() == AxisDirection::AxisY) {
0183         double bxmin = axis0.getMin();
0184         double bxmax = axis0.getMax();
0185         double bymin = axis1.getMin();
0186         double bymax = axis1.getMax();
0187         const std::vector<double> xboundaries = axis0.getBinEdges();
0188         const std::vector<double> yboundaries = axis1.getBinEdges();
0189         for (const double xval : xboundaries) {
0190           csvHelper.writeLine(grid, {xval, bymin}, {xval, bymax});
0191         }
0192         for (const double yval : yboundaries) {
0193           csvHelper.writeLine(grid, {bxmin, yval}, {bxmax, yval});
0194         }
0195       } else if (axis0.getDirection() == AxisDirection::AxisR &&
0196                  axis1.getDirection() == AxisDirection::AxisPhi) {
0197         double brmin = axis0.getMin();
0198         double brmax = axis0.getMax();
0199         double bphimin = axis1.getMin();
0200         double bphimax = axis1.getMax();
0201         const std::vector<double> rboundaries = axis0.getBinEdges();
0202         const std::vector<double> phiboundaries = axis1.getBinEdges();
0203         for (const double r : rboundaries) {
0204           csvHelper.writeArc(grid, r, bphimin, bphimax);
0205         }
0206         for (const double phi : phiboundaries) {
0207           double cphi = std::cos(phi);
0208           double sphi = std::sin(phi);
0209           csvHelper.writeLine(grid, {brmin * cphi, brmin * sphi},
0210                               {brmax * cphi, brmax * sphi});
0211         }
0212       }
0213     }
0214
0215     auto start = randomizer(startR0, startR1);
0216     auto end = randomizer(endR0, endR1);
0217
0218     std::ofstream segments;
0219     segments.open("Segmentizer" + name + "Segments_n" + std::to_string(index) +
0220                   ".csv");
0221
0222     std::ofstream cluster;
0223     cluster.open("Segmentizer" + name + "Cluster_n" + std::to_string(index) +
0224                  ".csv");
0225
0226     /// Run the Segmentizer
0227     auto cSegments = cl.segments(geoCtx, *surface, *segmentation, {start, end});
0228
0229     for (const auto& cs : cSegments) {
0230       csvHelper.writeLine(segments, cs.path2D[0], cs.path2D[1]);
0231     }
0232
0233     segments.close();
0234     cluster.close();
0235   }
0236 }
0237
0238 BOOST_AUTO_TEST_SUITE_END()
0239
0240 }  // namespace ActsTests