Core/Vertexing/HoughVertexFinderTests.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/Definitions/Units.hpp"
0013 #include "Acts/Utilities/Result.hpp"
0014 #include "Acts/Vertexing/HoughVertexFinder.hpp"
0015 #include "Acts/Vertexing/Vertex.hpp"
0016
0017 #include <cmath>
0018 #include <random>
0019 #include <stdexcept>
0020 #include <vector>
0021
0022 using namespace Acts;
0023
0024 namespace ActsTests {
0025
0026 /// @brief SpacePoint definition to be used for the unit tests. Implements all the relevant methods.
0027 struct SpacePoint4HVFT {
0028   SpacePoint4HVFT(double x, double y, double z) : m_x(x), m_y(y), m_z(z) {}
0029   double m_x;
0030   double m_y;
0031   double m_z;
0032   double x() const { return m_x; }
0033   double y() const { return m_y; }
0034   double z() const { return m_z; }
0035   double r() const { return std::sqrt(m_x * m_x + m_y * m_y); }
0036 };
0037
0038 /// @brief Provides random double number between $from and $to
0039 /// @param gen random number generator
0040 /// @param from lower threshold
0041 /// @param to upper threshold
0042 /// @return random number in [from,to)
0043 double getRndDouble(std::mt19937& gen, double from, double to) {
0044   std::uniform_real_distribution<double> dist(from, to);
0045   return dist(gen);
0046 }
0047
0048 /// @brief Provides random integer number between $from and $to
0049 /// @param gen random number generator
0050 /// @param from lower threshold
0051 /// @param to upper threshold
0052 /// @return random number in [from,to]
0053 int getRndInt(std::mt19937& gen, int from, int to) {
0054   std::uniform_int_distribution<int> dist(from, to);
0055   return dist(gen);
0056 }
0057
0058 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0059
0060 /// @brief Unit test for HoughVertexFinder. Compare the result to the easy-to-calculate expected result
0061 BOOST_AUTO_TEST_CASE(hough_vertex_finder_small_test) {
0062   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0063   houghVtxCfg.targetSPs = 1000;
0064   houghVtxCfg.minAbsEta = 0.3;
0065   houghVtxCfg.maxAbsEta = 3.0;
0066   houghVtxCfg.minHits = 3;
0067   houghVtxCfg.fillNeighbours = 0;
0068   houghVtxCfg.absEtaRanges = std::vector<double>({3.0});
0069   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0070   houghVtxCfg.rangeIterZ = std::vector<double>({100.05 * UnitConstants::mm});
0071   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001});
0072   houghVtxCfg.nBinsCotThetaIterZ = std::vector<unsigned int>({1000});
0073   houghVtxCfg.binsCotThetaDecrease = 1.0;
0074   houghVtxCfg.peakWidth = 1;
0075   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0076   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0077   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0078
0079   HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(std::move(houghVtxCfg));
0080
0081   double vtxX = 0., vtxY = 0., vtxZ = 20.;
0082
0083   std::vector<std::vector<double>> positions = {
0084       {10., 0., 25.},   {20., 0., 30.},   {30., 0., 35.},      // track 1
0085       {0., 5., 19.},    {0., 10., 18.},   {0, 15., 17.},       // track 2
0086       {-6., -4., 22.5}, {-12., -8., 25.}, {-18., -12., 27.5},  // track 3
0087       {-8., 2., 23.5},  {-16., 4., 27.},  {-24., 6., 30.5}};   // track 4
0088
0089   std::vector<SpacePoint4HVFT> inputSpacepoints;
0090   for (auto pos : positions) {
0091     inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
0092   }
0093
0094   auto vtx = houghVertexFinder.find(inputSpacepoints);
0095
0096   bool vtxFound = false;
0097   if (vtx.ok()) {
0098     // check if the found vertex has a compatible position
0099     if (std::abs(vtxX - vtx.value()[0]) < 1e-3 &&
0100         std::abs(vtxY - vtx.value()[1]) < 1e-3 &&
0101         std::abs(vtxZ - vtx.value()[2]) < 1e-3) {
0102       vtxFound = true;
0103     }
0104   }
0105
0106   BOOST_CHECK(vtxFound);
0107 }
0108
0109 /// @brief Unit test for HoughVertexFinder. Generates real-looking sets of the spacepoints, then finds a vertex, and then verifies the reconstructed vertex is actually near the original one
0110 BOOST_AUTO_TEST_CASE(hough_vertex_finder_full_test) {
0111   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0112   houghVtxCfg.targetSPs = 1000;
0113   houghVtxCfg.minAbsEta = 0.3;
0114   houghVtxCfg.maxAbsEta = 3.0;
0115   houghVtxCfg.minHits = 3;
0116   houghVtxCfg.fillNeighbours = 0;
0117   houghVtxCfg.absEtaRanges = std::vector<double>({3.0});
0118   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0119   houghVtxCfg.rangeIterZ = std::vector<double>({100.05 * UnitConstants::mm});
0120   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001});
0121   houghVtxCfg.nBinsCotThetaIterZ = std::vector<unsigned int>({1000});
0122   houghVtxCfg.binsCotThetaDecrease = 1.35;
0123   houghVtxCfg.peakWidth = 3;
0124   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0125   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0126   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0127
0128   HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(houghVtxCfg);
0129
0130   std::mt19937 gen(299792458);
0131
0132   int vtxFound = 0;
0133   int nEvents = 20;
0134   for (int event = 0; event < nEvents; event++) {
0135     double vtxX = getRndDouble(gen, -0.1, 0.1);
0136     double vtxY = getRndDouble(gen, -0.1, 0.1);
0137     double vtxZ = getRndDouble(gen, -50., 50.);
0138
0139     std::vector<SpacePoint4HVFT> inputSpacepoints;
0140
0141     // make straight lines originating from the given vertex
0142     int nTracks = getRndInt(gen, 200, 1000);
0143     for (int track = 0; track < nTracks; ++track) {
0144       // direction of the track
0145       double dirX = getRndDouble(gen, -1., 1.);
0146       double dirY = getRndDouble(gen, -1., 1.);
0147       double dirZ = getRndDouble(gen, -1., 1.);
0148       // use upper or lower intersection?
0149       int part = getRndInt(gen, 0, 1) * 2 - 1;
0150
0151       for (int rIndx = 1; rIndx <= 3; rIndx += 1) {
0152         double sgn = std::copysign(1., dirY);
0153         double dirR2 = dirX * dirX + dirY * dirY;
0154         double D = vtxX * (vtxY + dirY) - vtxY * (vtxX + dirX);
0155         // add some smearing to the layers
0156         // layers are (9-11), (19-21), and (29-31)
0157         double r = rIndx * 10 + getRndDouble(gen, -1., 1.);
0158         // intersection of the layer and the straigh line
0159         double x1 =
0160             (D * dirY + part * sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) /
0161             dirR2;
0162         double y1 = (-D * dirX +
0163                      part * std::abs(dirY) * std::sqrt(r * r * dirR2 - D * D)) /
0164                     dirR2;
0165         // how many units from the vertex to the intersection
0166         double zDist = std::abs((x1 - vtxX) / dirX);
0167         // use the same amount of units for distance in Z
0168         inputSpacepoints.emplace_back(x1, y1, zDist * dirZ + vtxZ);
0169       }
0170     }
0171
0172     auto vtx = houghVertexFinder.find(inputSpacepoints);
0173
0174     if (vtx.ok()) {
0175       // check if the found vertex has a compatible position
0176       if (std::abs(vtxX - vtx.value()[0]) < 0.2 &&
0177           std::abs(vtxY - vtx.value()[1]) < 0.2 &&
0178           std::abs(vtxZ - vtx.value()[2]) < 0.2) {
0179         ++vtxFound;
0180       }
0181     }
0182   }
0183
0184   BOOST_CHECK_EQUAL(vtxFound, nEvents);
0185 }
0186
0187 /// @brief Unit test for HoughVertexFinder. Provides no input spacepoints
0188 BOOST_AUTO_TEST_CASE(hough_vertex_finder_empty_test) {
0189   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0190   HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(std::move(houghVtxCfg));
0191
0192   // no input spacepoints
0193   std::vector<SpacePoint4HVFT> inputSpacepoints;
0194
0195   auto vtx = houghVertexFinder.find(inputSpacepoints);
0196
0197   bool vtxFound = false;
0198   if (vtx.ok()) {
0199     vtxFound = true;
0200   }
0201
0202   BOOST_CHECK(!vtxFound);
0203 }
0204
0205 /// @brief Unit test for HoughVertexFinder. Does not provides enough spacepoints
0206 BOOST_AUTO_TEST_CASE(hough_vertex_finder_insufficient_test) {
0207   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0208   houghVtxCfg.targetSPs = 1000;
0209   houghVtxCfg.minAbsEta = 0.3;
0210   houghVtxCfg.maxAbsEta = 3.0;
0211   houghVtxCfg.minHits = 3;  // requires 3 spacepoints per track
0212   houghVtxCfg.fillNeighbours = 0;
0213   houghVtxCfg.absEtaRanges = std::vector<double>({3.0});
0214   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0215   houghVtxCfg.rangeIterZ = std::vector<double>({100.05 * UnitConstants::mm});
0216   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001});
0217   houghVtxCfg.nBinsCotThetaIterZ = std::vector<unsigned int>({1000});
0218   houghVtxCfg.binsCotThetaDecrease = 1.0;
0219   houghVtxCfg.peakWidth = 1;
0220   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0221   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0222   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0223
0224   HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(std::move(houghVtxCfg));
0225
0226   // only 2 spacepoints per track provided
0227   std::vector<std::vector<double>> positions = {
0228       {10., 0., 25.},   {20., 0., 30.},     // track 1
0229       {0., 5., 19.},    {0., 10., 18.},     // track 2
0230       {-6., -4., 22.5}, {-12., -8., 25.}};  // track 3
0231
0232   std::vector<SpacePoint4HVFT> inputSpacepoints;
0233   for (auto pos : positions) {
0234     inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
0235   }
0236
0237   auto vtx = houghVertexFinder.find(inputSpacepoints);
0238
0239   bool vtxFound = false;
0240   if (vtx.ok()) {
0241     vtxFound = true;
0242   }
0243
0244   BOOST_CHECK(!vtxFound);
0245 }
0246
0247 /// @brief Unit test for HoughVertexFinder. Misconfigured #1
0248 BOOST_AUTO_TEST_CASE(hough_vertex_finder_misconfig1_test) {
0249   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0250   houghVtxCfg.targetSPs = 1000;
0251   houghVtxCfg.minAbsEta = 0.3;
0252   houghVtxCfg.maxAbsEta = 3.0;
0253   houghVtxCfg.minHits = 3;
0254   houghVtxCfg.fillNeighbours = 0;
0255   houghVtxCfg.absEtaRanges = std::vector<double>({3.0, 4.0});  // misconfigured
0256   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0257   houghVtxCfg.rangeIterZ = std::vector<double>({100.05 * UnitConstants::mm});
0258   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001});
0259   houghVtxCfg.nBinsCotThetaIterZ = std::vector<unsigned int>({1000});
0260   houghVtxCfg.binsCotThetaDecrease = 1.0;
0261   houghVtxCfg.peakWidth = 1;
0262   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0263   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0264   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0265
0266   BOOST_CHECK_THROW(HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(
0267                         std::move(houghVtxCfg)),
0268                     std::invalid_argument);
0269 }
0270
0271 /// @brief Unit test for HoughVertexFinder. Misconfigured #2
0272 BOOST_AUTO_TEST_CASE(hough_vertex_finder_misconfig2_test) {
0273   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0274   houghVtxCfg.targetSPs = 1000;
0275   houghVtxCfg.minAbsEta = 0.3;
0276   houghVtxCfg.maxAbsEta = 3.0;
0277   houghVtxCfg.minHits = 3;
0278   houghVtxCfg.fillNeighbours = 0;
0279   houghVtxCfg.absEtaRanges = std::vector<double>({3.0});
0280   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0281   houghVtxCfg.rangeIterZ = std::vector<double>(
0282       {100.05 * UnitConstants::mm, 100.05 * UnitConstants::mm});
0283   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001});  // misconfigured
0284   houghVtxCfg.nBinsCotThetaIterZ = std::vector<unsigned int>({1000, 1000});
0285   houghVtxCfg.binsCotThetaDecrease = 1.0;
0286   houghVtxCfg.peakWidth = 1;
0287   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0288   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0289   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0290
0291   BOOST_CHECK_THROW(HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(
0292                         std::move(houghVtxCfg)),
0293                     std::invalid_argument);
0294 }
0295
0296 /// @brief Unit test for HoughVertexFinder. Misconfigured #3
0297 BOOST_AUTO_TEST_CASE(hough_vertex_finder_misconfig3_test) {
0298   HoughVertexFinder<SpacePoint4HVFT>::Config houghVtxCfg;
0299   houghVtxCfg.targetSPs = 1000;
0300   houghVtxCfg.minAbsEta = 0.3;
0301   houghVtxCfg.maxAbsEta = 3.0;
0302   houghVtxCfg.minHits = 3;
0303   houghVtxCfg.fillNeighbours = 0;
0304   houghVtxCfg.absEtaRanges = std::vector<double>({3.0});
0305   houghVtxCfg.absEtaFractions = std::vector<double>({1.0});
0306   houghVtxCfg.rangeIterZ = std::vector<double>(
0307       {100.05 * UnitConstants::mm, 100.05 * UnitConstants::mm});
0308   houghVtxCfg.nBinsZIterZ = std::vector<unsigned int>({2001, 2001});
0309   houghVtxCfg.nBinsCotThetaIterZ =
0310       std::vector<unsigned int>({1000});  // misconfigured
0311   houghVtxCfg.binsCotThetaDecrease = 1.0;
0312   houghVtxCfg.peakWidth = 1;
0313   houghVtxCfg.defVtxPosition[0] = 0. * UnitConstants::mm;
0314   houghVtxCfg.defVtxPosition[1] = 0. * UnitConstants::mm;
0315   houghVtxCfg.defVtxPosition[2] = 0. * UnitConstants::mm;
0316
0317   BOOST_CHECK_THROW(HoughVertexFinder<SpacePoint4HVFT> houghVertexFinder(
0318                         std::move(houghVtxCfg)),
0319                     std::invalid_argument);
0320 }
0321
0322 BOOST_AUTO_TEST_SUITE_END()
0323
0324 }  // namespace ActsTests