Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-15 08:15:04

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 "Acts/Definitions/Algebra.hpp"
0010 #include "Acts/Surfaces/AnnulusBounds.hpp"
0011 #include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp"
0012 #include "Acts/Utilities/VectorHelpers.hpp"
0013 
0014 #include <algorithm>
0015 #include <fstream>
0016 #include <iostream>
0017 #include <random>
0018 #include <vector>
0019 
0020 using namespace Acts;
0021 
0022 int main(int /*argc*/, char** /*argv[]*/) {
0023   std::mt19937 rng(42);
0024   std::uniform_real_distribution<double> xDist(-2, 20);
0025   std::uniform_real_distribution<double> yDist(-2, 20);
0026 
0027   std::ofstream os{"annulus.csv"};
0028   os << "x,y,inside_abs,inside_tol0,inside_tol1,inside_tol01,inside_cov"
0029      << std::endl;
0030 
0031   // === PROBLEM DATA ===
0032 
0033   // bounds definition
0034   double minRadius = 7.2;
0035   double maxRadius = 12.0;
0036   double minPhi = 0.74195;
0037   double maxPhi = 1.33970;
0038   Vector2 offset(-2., 2.);
0039   AnnulusBounds aBounds(minRadius, maxRadius, minPhi, maxPhi, offset);
0040 
0041   // helper to convert to expected local frame
0042   auto toStripFrame = [&](const Vector2& xy) -> Vector2 {
0043     auto shifted = xy + offset;
0044     double r = VectorHelpers::perp(shifted);
0045     double phi = VectorHelpers::phi(shifted);
0046     return Vector2(r, phi);
0047   };
0048 
0049   auto random_point = [&]() -> Vector2 { return {xDist(rng), yDist(rng)}; };
0050 
0051   // for covariance based check, set up one;
0052   SquareMatrix2 cov;
0053   cov << 1.0, 0, 0, 0.05;
0054 
0055   BoundaryTolerance btNone = BoundaryTolerance::None();
0056   BoundaryTolerance btInf = BoundaryTolerance::Infinite();
0057   BoundaryTolerance btAbs = BoundaryTolerance::AbsoluteEuclidean(1.0);
0058   BoundaryTolerance btCov = BoundaryTolerance::Chi2Bound(cov, 1.0);
0059 
0060   // visualization to make sense of things
0061   for (std::size_t i = 0; i < 10000; i++) {
0062     const Vector2 loc{xDist(rng), yDist(rng)};
0063     auto locPC = toStripFrame(loc);
0064     bool isInsideNone = aBounds.inside(locPC, btNone);
0065     bool isInsideInf = aBounds.inside(locPC, btInf);
0066     bool isInsideAbs = aBounds.inside(locPC, btAbs);
0067     bool isInsideCov = aBounds.inside(locPC, btCov);
0068 
0069     os << loc.x() << "," << loc.y() << "," << isInsideNone << "," << isInsideInf
0070        << "," << isInsideAbs << "," << isInsideCov << std::endl;
0071   }
0072 
0073   std::vector<std::tuple<Vector2, bool, std::string>> testPoints{{
0074       {{7.0, 6.0}, true, "center"},
0075       {{3.0, 1.0}, false, "radial far out low"},
0076       {{5.0, 4.0}, false, "radial close out low"},
0077       {{6.0, 5.0}, true, "radial close in low"},
0078       {{8.5, 7.5}, true, "radial close in high"},
0079       {{10.0, 9.0}, false, "radial close out high"},
0080       {{15.0, 15.0}, false, "radial far out high"},
0081       {{0.0, 11.0}, false, "angular far out left"},
0082       {{4.0, 9.0}, false, "angular close out left"},
0083       {{5.0, 8.5}, true, "angular close in left"},
0084       {{9.0, 5.0}, true, "angular close in right"},
0085       {{9.5, 4.5}, false, "angular close out right"},
0086       {{11.0, 0.0}, false, "angular far out right"},
0087       {{2.0, 4.0}, false, "quad low left"},
0088       {{5.0, 14.0}, false, "quad high left"},
0089       {{13.0, 6.0}, false, "quad high right"},
0090       {{7.0, 1.0}, false, "quad low right"},
0091   }};
0092 
0093   // === BENCHMARKS ===
0094 
0095   // Number of benchmark runs
0096   constexpr int NTESTS = 5'000;
0097 
0098   // Some checks are much slower, so we tune down benchmark iterations
0099   constexpr int NTESTS_SLOW = NTESTS / 10;
0100 
0101   // We use this to switch between iteration counts
0102   enum class Mode { FastOutside, SlowOutside };
0103 
0104   // Benchmark output display
0105   auto print_bench_header = [](const std::string& check_name) {
0106     std::cout << check_name << ":" << std::endl;
0107   };
0108   auto print_bench_result = [](const std::string& bench_name,
0109                                const Acts::Test::MicroBenchmarkResult& res) {
0110     std::cout << "- " << bench_name << ": " << res << std::endl;
0111   };
0112 
0113   // Benchmark runner
0114   auto run_bench = [&](auto&& iteration, int num_iters,
0115                        const std::string& bench_name) {
0116     auto bench_result = Acts::Test::microBenchmark(iteration, num_iters);
0117     print_bench_result(bench_name, bench_result);
0118   };
0119 
0120   auto run_bench_with_inputs = [&](auto&& iterationWithArg, auto&& inputs,
0121                                    const std::string& bench_name) {
0122     auto bench_result = Acts::Test::microBenchmark(iterationWithArg, inputs);
0123     print_bench_result(bench_name, bench_result);
0124   };
0125   auto run_all_benches = [&](const BoundaryTolerance& check,
0126                              const std::string& check_name, const Mode mode) {
0127     // Announce a set of benchmarks
0128     print_bench_header(check_name);
0129 
0130     // Pre-determined "interesting" test points
0131     int num_inside_points = 0;
0132     int num_outside_points = 0;
0133     switch (mode) {
0134       case Mode::FastOutside:
0135         num_inside_points = NTESTS;
0136         num_outside_points = NTESTS;
0137         break;
0138       case Mode::SlowOutside:
0139         num_inside_points = NTESTS;
0140         num_outside_points = NTESTS_SLOW;
0141       default:  // do nothing
0142         break;
0143     };
0144 
0145     for (const auto& [loc, inside, label] : testPoints) {
0146       const auto locPC = toStripFrame(loc);
0147       run_bench([&] { return aBounds.inside(locPC, check); },
0148                 inside ? num_inside_points : num_outside_points, label);
0149     }
0150 
0151     // Pre-rolled random points
0152     std::vector<Vector2> points(num_outside_points);
0153     std::generate(points.begin(), points.end(),
0154                   [&] { return toStripFrame(random_point()); });
0155     run_bench_with_inputs(
0156         [&](const auto& point) { return aBounds.inside(point, check); }, points,
0157         "Random");
0158   };
0159 
0160   // Benchmark scenarios
0161   run_all_benches(btNone, "None", Mode::FastOutside);
0162   run_all_benches(btInf, "Infinite", Mode::FastOutside);
0163   run_all_benches(btAbs, "Absolute", Mode::FastOutside);
0164   run_all_benches(btCov, "Covariance", Mode::SlowOutside);
0165 
0166   return 0;
0167 }