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