File indexing completed on 2025-01-18 09:13:02
0001
0002
0003
0004
0005
0006
0007
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
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
0033
0034
0035
0036
0037 double getRndDouble(std::mt19937& gen, double from, double to) {
0038 return gen() / 4294967296. * (to - from) + from;
0039 }
0040
0041
0042
0043
0044
0045
0046 int getRndInt(std::mt19937& gen, int from, int to) {
0047 return static_cast<int>(gen() / 4294967296. * (to - from) + from);
0048 }
0049
0050
0051
0052
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
0066 return {abg[0], abg[1], abg[2], delta};
0067 }
0068
0069
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.},
0103 {2.1, 10.6, 15.2}, {2.7, 19.36666, 19.5},
0104 {4.5, 34., 25.4},
0105 {-6.25, -7.9, 10.5}, {-12.8, -15.08, 11.},
0106 {-20., -22., 11.5},
0107 {-6., 8., 2.4}, {-12., 15., -3.4},
0108 {-17., 21., -8.4},
0109 {7.8, 7.8, 16.0}, {14.8, 15., 17.},
0110 {22.8, 23.3, 18.},
0111 {-5.1, 8.1, -10.}, {-1.3, 9., -10.4},
0112 {3.1, 10.1, -11.1}};
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
0144 BOOST_CHECK(vtxFound);
0145 }
0146
0147
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.},
0181 {2., 9., 15.}, {3., 19., 20.},
0182 {4.5, 34., 27.5},
0183 {-6., -8., 10.5}, {-13., -15., 11.},
0184 {-20., -22., 11.5},
0185 {7., 7., 16.0}, {14., 14., 17.},
0186 {21., 21., 18.},
0187 {-5., 8., -10.}, {-1.5, 9., -10.5},
0188 {3., 10., -11.}};
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
0220 BOOST_CHECK(vtxFound);
0221 }
0222
0223
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
0262 int nTracks = getRndInt(gen, 200, 400);
0263 for (int track = 0; track < nTracks; ++track) {
0264
0265 double posX = vtxX;
0266 double posY = vtxY;
0267 double posZ = vtxZ;
0268
0269
0270 double dirX = getRndDouble(gen, -1., 1.);
0271 double dirY = getRndDouble(gen, -1., 1.);
0272 double dirZ = getRndDouble(gen, -1., 1.);
0273
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
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
0288 double dirR2 = dirX * dirX + dirY * dirY;
0289 double D = posX * (posY + dirY) - posY * (posX + dirX);
0290
0291
0292 double r = i * 10. + getRndDouble(gen, -1., 1.);
0293
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
0300 double zDist = std::abs((x1 - posX) / dirX);
0301
0302
0303 posX = x1;
0304 posY = y1;
0305
0306 posZ += zDist * dirZ;
0307
0308 if (i == 3) {
0309
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
0352 BOOST_CHECK_EQUAL(vtxFound, nEvents);
0353 }
0354
0355
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
0394 int nTracks = getRndInt(gen, 200, 400);
0395 for (int track = 0; track < nTracks; ++track) {
0396
0397 double dirX = getRndDouble(gen, -1., 1.);
0398 double dirY = getRndDouble(gen, -1., 1.);
0399 double dirZ = getRndDouble(gen, -1., 1.);
0400
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
0408
0409 double r = rIndx * 10 + getRndDouble(gen, -1., 1.);
0410
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
0418 double zDist = std::abs((x1 - vtxX) / dirX);
0419
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
0456 BOOST_CHECK_EQUAL(vtxFound, nEvents);
0457 }