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