Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:10:22

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018 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 http://mozilla.org/MPL/2.0/.
0008 
0009 #include "Acts/Definitions/TrackParametrization.hpp"
0010 #include "Acts/Definitions/Units.hpp"
0011 #include "Acts/MagneticField/BFieldMapUtils.hpp"
0012 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0013 #include "Acts/MagneticField/SolenoidBField.hpp"
0014 #include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp"
0015 #include "Acts/Utilities/VectorHelpers.hpp"
0016 
0017 #include <chrono>
0018 #include <fstream>
0019 #include <iostream>
0020 #include <random>
0021 #include <string>
0022 
0023 using namespace Acts::UnitLiterals;
0024 
0025 int main(int argc, char* argv[]) {
0026   std::size_t iters_map = 5e2;
0027   std::size_t iters_solenoid = 3;
0028   std::size_t runs_solenoid = 1000;
0029   if (argc >= 2) {
0030     iters_map = std::stoi(argv[1]);
0031   }
0032   if (argc >= 3) {
0033     iters_solenoid = std::stoi(argv[2]);
0034   }
0035   if (argc >= 4) {
0036     runs_solenoid = std::stoi(argv[3]);
0037   }
0038 
0039   const double L = 5.8_m;
0040   const double R = (2.56 + 2.46) * 0.5 * 0.5_m;
0041   const std::size_t nCoils = 1154;
0042   const double bMagCenter = 2_T;
0043   const std::size_t nBinsR = 150;
0044   const std::size_t nBinsZ = 200;
0045 
0046   double rMin = -0.1;
0047   double rMax = R * 2.;
0048   double zMin = 2 * (-L / 2.);
0049   double zMax = 2 * (L / 2.);
0050 
0051   Acts::SolenoidBField bSolenoidField({R, L, nCoils, bMagCenter});
0052   std::cout << "Building interpolated field map" << std::endl;
0053   auto bFieldMap = Acts::solenoidFieldMap({rMin, rMax}, {zMin, zMax},
0054                                           {nBinsR, nBinsZ}, bSolenoidField);
0055   Acts::MagneticFieldContext mctx{};
0056 
0057   std::minstd_rand rng;
0058   std::uniform_real_distribution<double> zDist(1.5 * (-L / 2.), 1.5 * L / 2.);
0059   std::uniform_real_distribution<double> rDist(0, R * 1.5);
0060   std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
0061   auto genPos = [&]() -> Acts::Vector3 {
0062     const double z = zDist(rng), r = rDist(rng), phi = phiDist(rng);
0063     return {r * std::cos(phi), r * std::sin(phi), z};
0064   };
0065 
0066   std::ofstream os{"bfield_bench.csv"};
0067 
0068   auto csv = [&](const std::string& name, auto res) {
0069     os << name << "," << res.run_timings.size() << "," << res.iters_per_run
0070        << "," << res.totalTime().count() << "," << res.runTimeMedian().count()
0071        << "," << 1.96 * res.runTimeError().count() << ","
0072        << res.iterTimeAverage().count() << ","
0073        << 1.96 * res.iterTimeError().count();
0074 
0075     os << std::endl;
0076   };
0077 
0078   os << "name,runs,iters,total_time,run_time_median,run_time_error,iter_"
0079         "time_average,iter_time_error"
0080      << std::endl;
0081 
0082   // SolenoidBField lookup is so slow that the cost of generating a random field
0083   // lookup position is negligible in comparison...
0084   std::cout << "Benchmarking random SolenoidBField lookup: " << std::flush;
0085   const auto solenoid_result = Acts::Test::microBenchmark(
0086       [&] { return bSolenoidField.getField(genPos()); }, iters_solenoid,
0087       runs_solenoid);
0088   std::cout << solenoid_result << std::endl;
0089   csv("solenoid", solenoid_result);
0090 
0091   // ...but for interpolated B-field map, the overhead of a field lookup is
0092   // comparable to that of generating a random position, so we must be more
0093   // careful. Hence we do two microbenchmarks which represent a kind of
0094   // lower and upper bound on field lookup performance.
0095   //
0096   // - The first benchmark operates at constant position, so it measures only
0097   //   field lookup overhead but has unrealistically good cache locality. In
0098   //   that sense, it provides a lower bound of field lookup performance.
0099   std::cout << "Benchmarking interpolated field lookup: " << std::flush;
0100   const auto fixedPos = genPos();
0101   const auto map_fixed_nocache_result = Acts::Test::microBenchmark(
0102       [&] { return bFieldMap.getField(fixedPos); }, iters_map);
0103   std::cout << map_fixed_nocache_result << std::endl;
0104   csv("interp_nocache_fixed", map_fixed_nocache_result);
0105 
0106   std::cout << "Benchmarking random interpolated field lookup: " << std::flush;
0107 
0108   // - The second benchmark generates random positions, so it is biased by the
0109   //   cost of random position generation and has unrealistically bad cache
0110   //   locality, but provides an upper bound of field lookup performance.
0111   const auto map_rand_result = Acts::Test::microBenchmark(
0112       [&] { return bFieldMap.getField(genPos()); }, iters_map);
0113   std::cout << map_rand_result << std::endl;
0114   csv("interp_nocache_random", map_rand_result);
0115 
0116   // - This variation of the first benchmark uses a fixed position again, but
0117   //   uses the cache infrastructure to evaluate how much of an impact it has on
0118   //   performance in this scenario. We expect this to improve performance as
0119   //   the cache will always be valid for the fixed point.
0120   {
0121     std::cout << "Benchmarking cached interpolated field lookup: "
0122               << std::flush;
0123     auto cache = bFieldMap.makeCache(mctx);
0124     const auto map_cached_result_cache = Acts::Test::microBenchmark(
0125         [&] { return bFieldMap.getField(fixedPos, cache).value(); }, iters_map);
0126     std::cout << map_cached_result_cache << std::endl;
0127     csv("interp_cache_fixed", map_cached_result_cache);
0128   }
0129 
0130   // - This variation of the second benchmark again generates random positions
0131   //   and uses the cache infrastructure to evaluate the impact on performance.
0132   //   We expect this to deteriorate performance, as the cache will most likely
0133   //   be invalid and need to be recreated, on top of the underlying lookup.
0134   {
0135     std::cout << "Benchmarking cached random interpolated field lookup: "
0136               << std::flush;
0137     auto cache2 = bFieldMap.makeCache(mctx);
0138     const auto map_rand_result_cache = Acts::Test::microBenchmark(
0139         [&] { return bFieldMap.getField(genPos(), cache2).value(); },
0140         iters_map);
0141     std::cout << map_rand_result_cache << std::endl;
0142     csv("interp_cache_random", map_rand_result_cache);
0143   }
0144 
0145   // - The fourth benchmark tests a more 'realistic' access pattern than fixed
0146   //   or random positions: it advances along a straight line (which is close to
0147   //   a slightly curved line which happens in particle propagation). This
0148   //   instance does not use the cache infrastructure, so is effectively close
0149   //   to the random points benchmark, although positions are not really random.
0150   {
0151     std::cout << "Benchmarking advancing interpolated field lookup: "
0152               << std::flush;
0153     Acts::Vector3 pos{0, 0, 0};
0154     Acts::Vector3 dir{};
0155     dir.setRandom();
0156     double h = 1e-3;
0157     std::vector<Acts::Vector3> steps;
0158     steps.reserve(iters_map);
0159     for (std::size_t i = 0; i < iters_map; i++) {
0160       pos += dir * h;
0161       double z = pos[Acts::eFreePos2];
0162       if (Acts::VectorHelpers::perp(pos) > rMax || z >= zMax || z < zMin) {
0163         break;
0164       }
0165       steps.push_back(pos);
0166     }
0167     const auto map_adv_result = Acts::Test::microBenchmark(
0168         [&](const auto& s) { return bFieldMap.getField(s); }, steps);
0169     std::cout << map_adv_result << std::endl;
0170     csv("interp_nocache_adv", map_adv_result);
0171 
0172     // - This variation of the fourth benchmark advances in a straight line, but
0173     //   also uses the cache infrastructure. As subsequent positions are close
0174     //   to one another, the cache will be valid for a certain number of points,
0175     //   before becoming invalid. This means we expect performance to improve
0176     //   over the uncached straight line advance.
0177 
0178     std::cout << "Benchmarking cached advancing interpolated field lookup: "
0179               << std::flush;
0180     auto cache = bFieldMap.makeCache(mctx);
0181     const auto map_adv_result_cache = Acts::Test::microBenchmark(
0182         [&](const auto& s) { return bFieldMap.getField(s, cache).value(); },
0183         steps);
0184     std::cout << map_adv_result_cache << std::endl;
0185     csv("interp_cache_adv", map_adv_result_cache);
0186   }
0187 }