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