Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:13:02

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/Vertexing/SingleSeedVertexFinder.hpp"
0014 
0015 #include <cmath>
0016 #include <iostream>
0017 #include <random>
0018 #include <vector>
0019 
0020 /// @brief SpacePoint definition to be used for the unit tests. Implements all the relevant methods.
0021 struct SpacePoint4SSVFT {
0022   SpacePoint4SSVFT(double x, double y, double z) : m_x(x), m_y(y), m_z(z) {}
0023   double m_x;
0024   double m_y;
0025   double m_z;
0026   double x() const { return m_x; }
0027   double y() const { return m_y; }
0028   double z() const { return m_z; }
0029   double r() const { return std::sqrt(m_x * m_x + m_y * m_y); }
0030 };
0031 
0032 /// @brief Provides random double number between $from and $to
0033 /// @param gen random number generator
0034 /// @param from lower threshold
0035 /// @param to upper threshold
0036 /// @return random number in [from,to)
0037 double getRndDouble(std::mt19937& gen, double from, double to) {
0038   return gen() / 4294967296. * (to - from) + from;
0039 }
0040 
0041 /// @brief Provides random integer number between $from and $to
0042 /// @param gen random number generator
0043 /// @param from lower threshold
0044 /// @param to upper threshold
0045 /// @return random number in [from,to)
0046 int getRndInt(std::mt19937& gen, int from, int to) {
0047   return static_cast<int>(gen() / 4294967296. * (to - from) + from);
0048 }
0049 
0050 /// @brief Calculates equation of the plane (alpha*x + beta*y + gamma*z + delta = 0), given the three points
0051 /// @param a,b,c The three points
0052 /// @return Parameters of the plane {alpha,beta,gamma,delta}
0053 std::vector<double> makePlaneFromTriplet(SpacePoint4SSVFT aa,
0054                                          SpacePoint4SSVFT bb,
0055                                          SpacePoint4SSVFT cc) {
0056   Acts::Vector3 a{aa.x(), aa.y(), aa.z()};
0057   Acts::Vector3 b{bb.x(), bb.y(), bb.z()};
0058   Acts::Vector3 c{cc.x(), cc.y(), cc.z()};
0059 
0060   Acts::Vector3 ba = b - a, ca = c - a;
0061 
0062   Acts::Vector3 abg = ba.cross(ca).normalized();
0063   double delta = -1. * abg.dot(a);
0064 
0065   // plane (alpha*x + beta*y + gamma*z + delta = 0)
0066   return {abg[0], abg[1], abg[2], delta};
0067 }
0068 
0069 /// @brief Unit test for SingleSeedVertexFinder. Fits a set of the spacepoints with planes and compare the result to the easy-to-calculate expected result
0070 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_small_planes_test) {
0071   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
0072   singleSeedVtxCfg.maxPhideviation = 0.2;
0073   singleSeedVtxCfg.maxXYdeviation = 0.1;
0074   singleSeedVtxCfg.maxXYZdeviation = 0.1;
0075   singleSeedVtxCfg.minTheta = 0.5;
0076   singleSeedVtxCfg.rMinNear = 6.f * Acts::UnitConstants::mm;
0077   singleSeedVtxCfg.rMaxNear = 14.9f * Acts::UnitConstants::mm;
0078   singleSeedVtxCfg.rMinMiddle = 15.f * Acts::UnitConstants::mm;
0079   singleSeedVtxCfg.rMaxMiddle = 24.9f * Acts::UnitConstants::mm;
0080   singleSeedVtxCfg.rMinFar = 25.f * Acts::UnitConstants::mm;
0081   singleSeedVtxCfg.rMaxFar = 44.9f * Acts::UnitConstants::mm;
0082   singleSeedVtxCfg.numPhiSlices = 10;
0083   singleSeedVtxCfg.useFracPhiSlices = 1.0;
0084   singleSeedVtxCfg.numZSlices = 10;
0085   singleSeedVtxCfg.useFracZSlices = 1.0;
0086   singleSeedVtxCfg.maxAbsZ = 50. * Acts::UnitConstants::mm;
0087   singleSeedVtxCfg.maxZPosition = 20.f * Acts::UnitConstants::mm;
0088   singleSeedVtxCfg.maxRPosition = 5.f * Acts::UnitConstants::mm;
0089   singleSeedVtxCfg.minimalizeWRT = "planes";
0090   singleSeedVtxCfg.maxIterations = 3;
0091   singleSeedVtxCfg.removeFraction = 0.3;
0092   singleSeedVtxCfg.minVtxShift = 0.3f * Acts::UnitConstants::mm;
0093   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
0094       singleSeedVtxCfg);
0095 
0096   double vtxX = 1.;
0097   double vtxY = -1.;
0098   double vtxZ = 10;
0099 
0100   std::vector<std::vector<double>> positions = {
0101       {10.8, 0., 7.},      {20.9, 0.7, 4.},
0102       {30.5, 1.9, 1.},  // plane ((-0.25,-0.25,-0.9), 9.0)
0103       {2.1, 10.6, 15.2},   {2.7, 19.36666, 19.5},
0104       {4.5, 34., 25.4},  // plane ((0.8, -0.3, 0.5), -6.1)
0105       {-6.25, -7.9, 10.5}, {-12.8, -15.08, 11.},
0106       {-20., -22., 11.5},  // plane ((-0.02,-0.05,-0.98), 9.77)
0107       {-6., 8., 2.4},      {-12., 15., -3.4},
0108       {-17., 21., -8.4},  // plane ((0.1,0.5,0.5), -4.6)
0109       {7.8, 7.8, 16.0},    {14.8, 15., 17.},
0110       {22.8, 23.3, 18.},  // this will be removed after iteration
0111       {-5.1, 8.1, -10.},   {-1.3, 9., -10.4},
0112       {3.1, 10.1, -11.1}};  // this will not form a triplet
0113 
0114   std::vector<SpacePoint4SSVFT> inputSpacepoints;
0115   for (auto pos : positions) {
0116     inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
0117   }
0118 
0119   auto t1 = std::chrono::high_resolution_clock::now();
0120   auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
0121   auto t2 = std::chrono::high_resolution_clock::now();
0122 
0123   bool vtxFound = false;
0124   if (vtx.ok()) {
0125     std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
0126               << " ms at x = " << vtx.value()[0] << "mm, y = " << vtx.value()[1]
0127               << "mm, z = " << vtx.value()[2] << "mm" << std::endl;
0128     std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0129               << "mm, z = " << vtxZ << "mm" << std::endl;
0130 
0131     if (std::abs(vtxX - vtx.value()[0]) < 1e-3 &&
0132         std::abs(vtxY - vtx.value()[1]) < 1e-3 &&
0133         std::abs(vtxZ - vtx.value()[2]) < 1e-3) {
0134       vtxFound = true;
0135     }
0136   } else {
0137     std::cout << "Not found a vertex in the event after "
0138               << (t2 - t1).count() / 1e6 << " ms" << std::endl;
0139     std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0140               << "mm, z = " << vtxZ << "mm" << std::endl;
0141   }
0142 
0143   // check if found vertex has compatible position
0144   BOOST_CHECK(vtxFound);
0145 }
0146 
0147 /// @brief Unit test for SingleSeedVertexFinder. Fits a set of the spacepoints with straigh lines and compare the result to the easy-to-calculate expected result
0148 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_small_rays_test) {
0149   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
0150   singleSeedVtxCfg.maxPhideviation = 0.2;
0151   singleSeedVtxCfg.maxXYdeviation = 0.1;
0152   singleSeedVtxCfg.maxXYZdeviation = 0.1;
0153   singleSeedVtxCfg.minTheta = 0.5;
0154   singleSeedVtxCfg.rMinNear = 6.f * Acts::UnitConstants::mm;
0155   singleSeedVtxCfg.rMaxNear = 14.9f * Acts::UnitConstants::mm;
0156   singleSeedVtxCfg.rMinMiddle = 15.f * Acts::UnitConstants::mm;
0157   singleSeedVtxCfg.rMaxMiddle = 24.9f * Acts::UnitConstants::mm;
0158   singleSeedVtxCfg.rMinFar = 25.f * Acts::UnitConstants::mm;
0159   singleSeedVtxCfg.rMaxFar = 44.9f * Acts::UnitConstants::mm;
0160   singleSeedVtxCfg.numPhiSlices = 10;
0161   singleSeedVtxCfg.useFracPhiSlices = 1.0;
0162   singleSeedVtxCfg.numZSlices = 10;
0163   singleSeedVtxCfg.useFracZSlices = 1.0;
0164   singleSeedVtxCfg.maxAbsZ = 50. * Acts::UnitConstants::mm;
0165   singleSeedVtxCfg.maxZPosition = 20.f * Acts::UnitConstants::mm;
0166   singleSeedVtxCfg.maxRPosition = 5.f * Acts::UnitConstants::mm;
0167   singleSeedVtxCfg.minimalizeWRT = "rays";
0168   singleSeedVtxCfg.maxIterations = 1;
0169   singleSeedVtxCfg.removeFraction = 0.3;
0170   singleSeedVtxCfg.minVtxShift = 0.3f * Acts::UnitConstants::mm;
0171   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
0172       singleSeedVtxCfg);
0173 
0174   double vtxX = 1.;
0175   double vtxY = -1.;
0176   double vtxZ = 10;
0177 
0178   std::vector<std::vector<double>> positions = {
0179       {11., 0., 7.},      {21., 1., 4.},
0180       {31., 2., 1.},  // this ray is Ok
0181       {2., 9., 15.},      {3., 19., 20.},
0182       {4.5, 34., 27.5},  // this ray is Ok
0183       {-6., -8., 10.5},   {-13., -15., 11.},
0184       {-20., -22., 11.5},  // this ray is Ok
0185       {7., 7., 16.0},     {14., 14., 17.},
0186       {21., 21., 18.},  // this will be removed after iteration
0187       {-5., 8., -10.},    {-1.5, 9., -10.5},
0188       {3., 10., -11.}};  // this will not form a triplet
0189 
0190   std::vector<SpacePoint4SSVFT> inputSpacepoints;
0191   for (auto pos : positions) {
0192     inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
0193   }
0194 
0195   auto t1 = std::chrono::high_resolution_clock::now();
0196   auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
0197   auto t2 = std::chrono::high_resolution_clock::now();
0198 
0199   bool vtxFound = false;
0200   if (vtx.ok()) {
0201     std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
0202               << " ms at x = " << vtx.value()[0] << "mm, y = " << vtx.value()[1]
0203               << "mm, z = " << vtx.value()[2] << "mm" << std::endl;
0204     std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0205               << "mm, z = " << vtxZ << "mm" << std::endl;
0206 
0207     if (std::abs(vtxX - vtx.value()[0]) < 1e-3 &&
0208         std::abs(vtxY - vtx.value()[1]) < 1e-3 &&
0209         std::abs(vtxZ - vtx.value()[2]) < 1e-3) {
0210       vtxFound = true;
0211     }
0212   } else {
0213     std::cout << "Not found a vertex in the event after "
0214               << (t2 - t1).count() / 1e6 << " ms" << std::endl;
0215     std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0216               << "mm, z = " << vtxZ << "mm" << std::endl;
0217   }
0218 
0219   // check if found vertex has compatible position
0220   BOOST_CHECK(vtxFound);
0221 }
0222 
0223 /// @brief Unit test for SingleSeedVertexFinder. Generates real-looking sets of the spacepoints, then fits them with planes, finds a vertex, and then verifies the reconstructed vertex is actually near the original one
0224 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_planes_test) {
0225   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
0226   singleSeedVtxCfg.maxPhideviation = 0.1;
0227   singleSeedVtxCfg.maxXYdeviation = 0.1;
0228   singleSeedVtxCfg.maxXYZdeviation = 0.1;
0229   singleSeedVtxCfg.minTheta = 0.5;
0230   singleSeedVtxCfg.rMinNear = 8.f * Acts::UnitConstants::mm;
0231   singleSeedVtxCfg.rMaxNear = 12.f * Acts::UnitConstants::mm;
0232   singleSeedVtxCfg.rMinMiddle = 18.f * Acts::UnitConstants::mm;
0233   singleSeedVtxCfg.rMaxMiddle = 22.f * Acts::UnitConstants::mm;
0234   singleSeedVtxCfg.rMinFar = 28.f * Acts::UnitConstants::mm;
0235   singleSeedVtxCfg.rMaxFar = 32.f * Acts::UnitConstants::mm;
0236   singleSeedVtxCfg.numPhiSlices = 60;
0237   singleSeedVtxCfg.useFracPhiSlices = 0.8;
0238   singleSeedVtxCfg.numZSlices = 150;
0239   singleSeedVtxCfg.useFracZSlices = 0.8;
0240   singleSeedVtxCfg.maxAbsZ = 75. * Acts::UnitConstants::mm;
0241   singleSeedVtxCfg.maxZPosition = 15.f * Acts::UnitConstants::mm;
0242   singleSeedVtxCfg.maxRPosition = 2.5f * Acts::UnitConstants::mm;
0243   singleSeedVtxCfg.minimalizeWRT = "planes";
0244   singleSeedVtxCfg.maxIterations = 20;
0245   singleSeedVtxCfg.removeFraction = 0.1;
0246   singleSeedVtxCfg.minVtxShift = 0.05f * Acts::UnitConstants::mm;
0247   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
0248       singleSeedVtxCfg);
0249 
0250   std::mt19937 gen(299792458);
0251 
0252   int vtxFound = 0;
0253   int nEvents = 5;
0254   for (int event = 0; event < nEvents; event++) {
0255     double vtxX = getRndDouble(gen, -1., 1.);
0256     double vtxY = getRndDouble(gen, -1., 1.);
0257     double vtxZ = getRndDouble(gen, -10., 10.);
0258 
0259     std::vector<SpacePoint4SSVFT> inputSpacepoints;
0260 
0261     // make straight lines originating from the given vertex
0262     int nTracks = getRndInt(gen, 200, 400);
0263     for (int track = 0; track < nTracks; ++track) {
0264       // initial position of the track
0265       double posX = vtxX;
0266       double posY = vtxY;
0267       double posZ = vtxZ;
0268 
0269       // initial direction of the track
0270       double dirX = getRndDouble(gen, -1., 1.);
0271       double dirY = getRndDouble(gen, -1., 1.);
0272       double dirZ = getRndDouble(gen, -1., 1.);
0273       // rotation of the track
0274       double theta = getRndDouble(gen, 0.03, 0.09) *
0275                      (getRndDouble(gen, -1., 1.) > 0 ? 1 : -1);
0276       double sgn = std::copysign(1., dirY);
0277 
0278       for (int i = 1; i <= 3; ++i) {
0279         if (i != 1) {
0280           // rotate direction, not for the first time
0281           double dirXtmp = cos(theta) * dirX - sin(theta) * dirY;
0282           double dirYtmp = sin(theta) * dirX + cos(theta) * dirY;
0283           dirX = dirXtmp;
0284           dirY = dirYtmp;
0285         }
0286 
0287         // double sgn=std::copysign(1.,dirY);
0288         double dirR2 = dirX * dirX + dirY * dirY;
0289         double D = posX * (posY + dirY) - posY * (posX + dirX);
0290         // add some smearing to the layers
0291         // layers are (9-11), (19-21), and (29-31)
0292         double r = i * 10. + getRndDouble(gen, -1., 1.);
0293         // intersection of the layer and the straigh line
0294         double x1 =
0295             (D * dirY + sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) / dirR2;
0296         double y1 =
0297             (-D * dirX + std::abs(dirY) * std::sqrt(r * r * dirR2 - D * D)) /
0298             dirR2;
0299         // how many units from the vertex to the intersection
0300         double zDist = std::abs((x1 - posX) / dirX);
0301 
0302         // position of the new spacepoint
0303         posX = x1;
0304         posY = y1;
0305         // use the same amount of units for distance in Z
0306         posZ += zDist * dirZ;
0307 
0308         if (i == 3) {
0309           // move z position, so the vertex will be part of the plane
0310           auto abgd = makePlaneFromTriplet({vtxX, vtxY, vtxZ},
0311                                            inputSpacepoints.rbegin()[1],
0312                                            inputSpacepoints.rbegin()[0]);
0313           posZ = -1. * (abgd[0] * posX + abgd[1] * posY + abgd[3]) / abgd[2];
0314         }
0315 
0316         inputSpacepoints.emplace_back(posX, posY, posZ);
0317       }
0318     }
0319 
0320     auto t1 = std::chrono::high_resolution_clock::now();
0321     auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
0322     auto t2 = std::chrono::high_resolution_clock::now();
0323 
0324     if (vtx.ok()) {
0325       std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
0326                 << " ms at x = " << vtx.value()[0]
0327                 << "mm, y = " << vtx.value()[1] << "mm, z = " << vtx.value()[2]
0328                 << "mm" << std::endl;
0329       std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0330                 << "mm, z = " << vtxZ << "mm" << std::endl;
0331       std::cout << "Difference is in x = " << std::abs(vtxX - vtx.value()[0])
0332                 << "mm, y = " << std::abs(vtxY - vtx.value()[1])
0333                 << "mm, z = " << std::abs(vtxZ - vtx.value()[2]) << "mm"
0334                 << std::endl;
0335       if (std::abs(vtxX - vtx.value()[0]) < 2.0 &&
0336           std::abs(vtxY - vtx.value()[1]) < 2.0 &&
0337           std::abs(vtxZ - vtx.value()[2]) < 0.3) {
0338         ++vtxFound;
0339       }
0340     } else {
0341       std::cout << "Not found a vertex in the event after "
0342                 << (t2 - t1).count() / 1e6 << " ms" << std::endl;
0343       std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0344                 << "mm, z = " << vtxZ << "mm" << std::endl;
0345     }
0346   }
0347 
0348   std::cout << "Found " << vtxFound << " out of " << nEvents << " vertices"
0349             << std::endl;
0350 
0351   // check if all vertices have compatible positions
0352   BOOST_CHECK_EQUAL(vtxFound, nEvents);
0353 }
0354 
0355 /// @brief Unit test for SingleSeedVertexFinder. Generates real-looking sets of the spacepoints, then fits them with rays, finds a vertex, and then verifies the reconstructed vertex is actually near the original one
0356 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_rays_test) {
0357   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
0358   singleSeedVtxCfg.maxPhideviation = 0.1;
0359   singleSeedVtxCfg.maxXYdeviation = 0.1;
0360   singleSeedVtxCfg.maxXYZdeviation = 0.1;
0361   singleSeedVtxCfg.minTheta = 0.5;
0362   singleSeedVtxCfg.rMinNear = 8.f * Acts::UnitConstants::mm;
0363   singleSeedVtxCfg.rMaxNear = 12.f * Acts::UnitConstants::mm;
0364   singleSeedVtxCfg.rMinMiddle = 18.f * Acts::UnitConstants::mm;
0365   singleSeedVtxCfg.rMaxMiddle = 22.f * Acts::UnitConstants::mm;
0366   singleSeedVtxCfg.rMinFar = 28.f * Acts::UnitConstants::mm;
0367   singleSeedVtxCfg.rMaxFar = 32.f * Acts::UnitConstants::mm;
0368   singleSeedVtxCfg.numPhiSlices = 60;
0369   singleSeedVtxCfg.useFracPhiSlices = 0.8;
0370   singleSeedVtxCfg.numZSlices = 150;
0371   singleSeedVtxCfg.useFracZSlices = 0.8;
0372   singleSeedVtxCfg.maxAbsZ = 75. * Acts::UnitConstants::mm;
0373   singleSeedVtxCfg.maxZPosition = 15.f * Acts::UnitConstants::mm;
0374   singleSeedVtxCfg.maxRPosition = 2.5f * Acts::UnitConstants::mm;
0375   singleSeedVtxCfg.minimalizeWRT = "rays";
0376   singleSeedVtxCfg.maxIterations = 10;
0377   singleSeedVtxCfg.removeFraction = 0.1;
0378   singleSeedVtxCfg.minVtxShift = 0.05f * Acts::UnitConstants::mm;
0379   Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
0380       singleSeedVtxCfg);
0381 
0382   std::mt19937 gen(299792458);
0383 
0384   int vtxFound = 0;
0385   int nEvents = 5;
0386   for (int event = 0; event < nEvents; event++) {
0387     double vtxX = getRndDouble(gen, -1., 1.);
0388     double vtxY = getRndDouble(gen, -1., 1.);
0389     double vtxZ = getRndDouble(gen, -10., 10.);
0390 
0391     std::vector<SpacePoint4SSVFT> inputSpacepoints;
0392 
0393     // make straight lines originating from the given vertex
0394     int nTracks = getRndInt(gen, 200, 400);
0395     for (int track = 0; track < nTracks; ++track) {
0396       // direction of the track
0397       double dirX = getRndDouble(gen, -1., 1.);
0398       double dirY = getRndDouble(gen, -1., 1.);
0399       double dirZ = getRndDouble(gen, -1., 1.);
0400       // use upper or lower intersection?
0401       int part = (getRndDouble(gen, -1., 1.) > 0 ? 1 : -1);
0402 
0403       for (int rIndx = 1; rIndx <= 3; rIndx += 1) {
0404         double sgn = std::copysign(1., dirY);
0405         double dirR2 = dirX * dirX + dirY * dirY;
0406         double D = vtxX * (vtxY + dirY) - vtxY * (vtxX + dirX);
0407         // add some smearing to the layers
0408         // layers are (9-11), (19-21), and (29-31)
0409         double r = rIndx * 10 + getRndDouble(gen, -1., 1.);
0410         // intersection of the layer and the straigh line
0411         double x1 =
0412             (D * dirY + part * sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) /
0413             dirR2;
0414         double y1 = (-D * dirX +
0415                      part * std::abs(dirY) * std::sqrt(r * r * dirR2 - D * D)) /
0416                     dirR2;
0417         // how many units from the vertex to the intersection
0418         double zDist = std::abs((x1 - vtxX) / dirX);
0419         // use the same amount of units for distance in Z
0420         inputSpacepoints.emplace_back(x1, y1, zDist * dirZ + vtxZ);
0421       }
0422     }
0423 
0424     auto t1 = std::chrono::high_resolution_clock::now();
0425     auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
0426     auto t2 = std::chrono::high_resolution_clock::now();
0427 
0428     if (vtx.ok()) {
0429       std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
0430                 << " ms at x = " << vtx.value()[0]
0431                 << "mm, y = " << vtx.value()[1] << "mm, z = " << vtx.value()[2]
0432                 << "mm" << std::endl;
0433       std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0434                 << "mm, z = " << vtxZ << "mm" << std::endl;
0435       std::cout << "Difference is in x = " << std::abs(vtxX - vtx.value()[0])
0436                 << "mm, y = " << std::abs(vtxY - vtx.value()[1])
0437                 << "mm, z = " << std::abs(vtxZ - vtx.value()[2]) << "mm"
0438                 << std::endl;
0439       if (std::abs(vtxX - vtx.value()[0]) < 0.3 &&
0440           std::abs(vtxY - vtx.value()[1]) < 0.3 &&
0441           std::abs(vtxZ - vtx.value()[2]) < 0.3) {
0442         ++vtxFound;
0443       }
0444     } else {
0445       std::cout << "Not found a vertex in the event after "
0446                 << (t2 - t1).count() / 1e6 << " ms" << std::endl;
0447       std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
0448                 << "mm, z = " << vtxZ << "mm" << std::endl;
0449     }
0450   }
0451 
0452   std::cout << "Found " << vtxFound << " out of " << nEvents << " vertices"
0453             << std::endl;
0454 
0455   // check if all vertices have compatible positions
0456   BOOST_CHECK_EQUAL(vtxFound, nEvents);
0457 }