Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:27

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