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