Fatras/Digitization/SegmentizerTests.cpp

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