Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-16 08:03:55

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