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