Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:11:25 

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/MagneticField/BFieldMapUtils.hpp"
0010 
0011 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0012 #include "Acts/MagneticField/SolenoidBField.hpp"
0013 #include "Acts/Utilities/Axis.hpp"
0014 #include "Acts/Utilities/Grid.hpp"
0015 #include "Acts/Utilities/Helpers.hpp"
0016 #include "Acts/Utilities/Result.hpp"
0017 #include "Acts/Utilities/VectorHelpers.hpp"
0018 #include "Acts/Utilities/detail/grid_helper.hpp"
0019 
0020 #include <algorithm>
0021 #include <cmath>
0022 #include <cstddef>
0023 #include <cstdlib>
0024 #include <initializer_list>
0025 #include <limits>
0026 #include <set>
0027 #include <tuple>
0028 
0029 using Acts::VectorHelpers::perp;
0030 using Acts::VectorHelpers::phi;
0031 
0032 Acts::InterpolatedBFieldMap<
0033     Acts::Grid<Acts::Vector2, Acts::Axis<Acts::AxisType::Equidistant>,
0034                Acts::Axis<Acts::AxisType::Equidistant>>>
0035 Acts::fieldMapRZ(
0036     const std::function<std::size_t(std::array<std::size_t, 2> binsRZ,
0037                                     std::array<std::size_t, 2> nBinsRZ)>&
0038         localToGlobalBin,
0039     std::vector<double> rPos, std::vector<double> zPos,
0040     const std::vector<Vector2>& bField, double lengthUnit, double BFieldUnit,
0041     bool firstQuadrant) {
0042   // [1] Create Grid
0043   const auto [rMin, rMax, rBinCount] = detail::getMinMaxAndBinCount(rPos);
0044   auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos);
0045 
0046   const std::size_t nBinsR = rBinCount;
0047   std::size_t nBinsZ = zBinCount;
0048 
0049   if (firstQuadrant) {
0050     zMin = -zPos[nBinsZ - 1];
0051     nBinsZ = 2 * nBinsZ - 1;
0052   }
0053 
0054   // Create the axis for the grid
0055   Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR);
0056   Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0057 
0058   // Create the grid
0059   Grid grid(Type<Vector2>, std::move(rAxis), std::move(zAxis));
0060   using Grid_t = decltype(grid);
0061 
0062   // [2] Set the bField values
0063   const std::array<std::size_t, 2> nIndices = {{rBinCount, zBinCount}};
0064   for (std::size_t i = 1; i <= nBinsR; ++i) {
0065     for (std::size_t j = 1; j <= nBinsZ; ++j) {
0066       Grid_t::index_t indices = {{i, j}};
0067       // std::vectors begin with 0 and we do not want the user needing to take
0068       // underflow or overflow bins in account this is why we need to subtract
0069       // by one
0070       if (firstQuadrant) {
0071         std::size_t n = std::abs(static_cast<std::ptrdiff_t>(j) -
0072                                  static_cast<std::ptrdiff_t>(zBinCount));
0073 
0074         grid.atLocalBins(indices) =
0075             bField.at(localToGlobalBin({{i - 1, n}}, nIndices)) * BFieldUnit;
0076       } else {
0077         grid.atLocalBins(indices) =
0078             bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices)) *
0079             BFieldUnit;
0080       }
0081     }
0082   }
0083   grid.setExteriorBins(Vector2::Zero());
0084 
0085   // [3] Create the transformation for the position map (x,y,z) -> (r,z)
0086   auto transformPos = [](const Vector3& pos) {
0087     return Vector2(perp(pos), pos.z());
0088   };
0089 
0090   // [4] Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz)
0091   auto transformBField = [](const Vector2& field, const Vector3& pos) {
0092     const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y();
0093     double cosPhi = 1.;
0094     double sinPhi = 0.;
0095 
0096     if (rSinTheta2 > std::numeric_limits<double>::min()) {
0097       const double invRsinTheta = 1. / std::sqrt(rSinTheta2);
0098       cosPhi = pos.x() * invRsinTheta;
0099       sinPhi = pos.y() * invRsinTheta;
0100     }
0101 
0102     return Vector3(field.x() * cosPhi, field.x() * sinPhi, field.y());
0103   };
0104 
0105   // [5] Create the mapper & BField Service create field mapping
0106   return InterpolatedBFieldMap<Grid_t>(
0107       {transformPos, transformBField, std::move(grid)});
0108 }
0109 
0110 Acts::InterpolatedBFieldMap<
0111     Acts::Grid<Acts::Vector3, Acts::Axis<Acts::AxisType::Equidistant>,
0112                Acts::Axis<Acts::AxisType::Equidistant>,
0113                Acts::Axis<Acts::AxisType::Equidistant>>>
0114 Acts::fieldMapXYZ(
0115     const std::function<std::size_t(std::array<std::size_t, 3> binsXYZ,
0116                                     std::array<std::size_t, 3> nBinsXYZ)>&
0117         localToGlobalBin,
0118     std::vector<double> xPos, std::vector<double> yPos,
0119     std::vector<double> zPos, const std::vector<Vector3>& bField,
0120     double lengthUnit, double BFieldUnit, bool firstOctant) {
0121   // [1] Create Grid
0122   auto [xMin, xMax, xBinCount] = detail::getMinMaxAndBinCount(xPos);
0123   auto [yMin, yMax, yBinCount] = detail::getMinMaxAndBinCount(yPos);
0124   auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos);
0125 
0126   std::size_t nBinsX = xBinCount;
0127   std::size_t nBinsY = yBinCount;
0128   std::size_t nBinsZ = zBinCount;
0129 
0130   if (firstOctant) {
0131     xMin = -xPos[nBinsX - 1];
0132     nBinsX = 2 * nBinsX - 1;
0133     yMin = -yPos[nBinsY - 1];
0134     nBinsY = 2 * nBinsY - 1;
0135     zMin = -zPos[nBinsZ - 1];
0136     nBinsZ = 2 * nBinsZ - 1;
0137   }
0138 
0139   Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX);
0140   Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY);
0141   Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0142   // Create the grid
0143   Grid grid(Type<Vector3>, std::move(xAxis), std::move(yAxis),
0144             std::move(zAxis));
0145   using Grid_t = decltype(grid);
0146 
0147   // [2] Set the bField values
0148   const std::array<std::size_t, 3> nIndices = {
0149       {xBinCount, yBinCount, zBinCount}};
0150 
0151   auto calcAbsDiff = [](std::size_t val, std::size_t binCount) {
0152     return std::abs(static_cast<std::ptrdiff_t>(val) -
0153                     static_cast<std::ptrdiff_t>(binCount));
0154   };
0155 
0156   for (std::size_t i = 1; i <= nBinsX; ++i) {
0157     for (std::size_t j = 1; j <= nBinsY; ++j) {
0158       for (std::size_t k = 1; k <= nBinsZ; ++k) {
0159         Grid_t::index_t indices = {{i, j, k}};
0160         // std::vectors begin with 0 and we do not want the user needing to take
0161         // underflow or overflow bins in account this is why we need to subtract
0162         // by one
0163         if (firstOctant) {
0164           const std::size_t l = calcAbsDiff(i, xBinCount);
0165           const std::size_t m = calcAbsDiff(j, yBinCount);
0166           const std::size_t n = calcAbsDiff(k, zBinCount);
0167 
0168           grid.atLocalBins(indices) =
0169               bField.at(localToGlobalBin({{l, m, n}}, nIndices)) * BFieldUnit;
0170         } else {
0171           grid.atLocalBins(indices) =
0172               bField.at(localToGlobalBin({{i - 1, j - 1, k - 1}}, nIndices)) *
0173               BFieldUnit;
0174         }
0175       }
0176     }
0177   }
0178   grid.setExteriorBins(Vector3::Zero());
0179 
0180   // [3] Create the transformation for the position map (x,y,z) -> (r,z)
0181   auto transformPos = [](const Vector3& pos) { return pos; };
0182 
0183   // [4] Create the transformation for the BField map (Bx,By,Bz) -> (Bx,By,Bz)
0184   auto transformBField = [](const Vector3& field, const Vector3& /*pos*/) {
0185     return field;
0186   };
0187 
0188   // [5] Create the mapper & BField Service create field mapping
0189   return InterpolatedBFieldMap<Grid_t>(
0190       {transformPos, transformBField, std::move(grid)});
0191 }
0192 
0193 Acts::InterpolatedBFieldMap<
0194     Acts::Grid<Acts::Vector2, Acts::Axis<Acts::AxisType::Equidistant>,
0195                Acts::Axis<Acts::AxisType::Equidistant>>>
0196 Acts::solenoidFieldMap(const std::pair<double, double>& rLim,
0197                        const std::pair<double, double>& zLim,
0198                        const std::pair<std::size_t, std::size_t>& nBins,
0199                        const SolenoidBField& field) {
0200   auto [rMin, rMax] = rLim;
0201   auto [zMin, zMax] = zLim;
0202   const auto [nBinsR, nBinsZ] = nBins;
0203 
0204   double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1);
0205   double stepR = std::abs(rMax - rMin) / (nBinsR - 1);
0206   rMax += stepR;
0207   zMax += stepZ;
0208 
0209   // Create the axis for the grid
0210   Axis rAxis(rMin, rMax, nBinsR);
0211   Axis zAxis(zMin, zMax, nBinsZ);
0212 
0213   // Create the grid
0214   Grid grid(Type<Vector2>, std::move(rAxis), std::move(zAxis));
0215   using Grid_t = decltype(grid);
0216 
0217   // Create the transformation for the position map (x,y,z) -> (r,z)
0218   auto transformPos = [](const Vector3& pos) {
0219     return Vector2(perp(pos), pos.z());
0220   };
0221 
0222   // Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz)
0223   auto transformBField = [](const Vector2& bField, const Vector3& pos) {
0224     const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y();
0225     double cosPhi = 1.;
0226     double sinPhi = 0.;
0227 
0228     if (rSinTheta2 > std::numeric_limits<double>::min()) {
0229       const double invRsinTheta = 1. / std::sqrt(rSinTheta2);
0230       cosPhi = pos.x() * invRsinTheta;
0231       sinPhi = pos.y() * invRsinTheta;
0232     }
0233 
0234     return Vector3(bField.x() * cosPhi, bField.x() * sinPhi, bField.y());
0235   };
0236 
0237   // iterate over all bins, set their value to the solenoid value at their lower
0238   // left position
0239   for (std::size_t i = 0; i <= nBinsR + 1; i++) {
0240     for (std::size_t j = 0; j <= nBinsZ + 1; j++) {
0241       Grid_t::index_t index({i, j});
0242       if (i == 0 || j == 0 || i == nBinsR + 1 || j == nBinsZ + 1) {
0243         // under or overflow bin, set zero
0244         grid.atLocalBins(index) = Grid_t::value_type(0, 0);
0245       } else {
0246         // regular bin, get lower left boundary
0247         Grid_t::point_t lowerLeft = grid.lowerLeftBinEdge(index);
0248         // do lookup
0249         Vector2 B = field.getField(Vector2(lowerLeft[0], lowerLeft[1]));
0250         grid.atLocalBins(index) = B;
0251       }
0252     }
0253   }
0254 
0255   // Create the mapper & BField Service create field mapping
0256   InterpolatedBFieldMap<Grid_t> map(
0257       {transformPos, transformBField, std::move(grid)});
0258   return map;
0259 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!