From bf997eb640537c9d91f0ee7f54e95216a680a98a Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Mon, 2 Dec 2024 13:36:44 +0100 Subject: [PATCH] refactor: modernise map utils (#3615) ## Summary by CodeRabbit - **New Features** - Introduced a new function to compute minimum and maximum values along with unique bin counts for input vectors. - **Refactor** - Enhanced the handling of magnetic field data by improving type usage and logic in multiple functions. - Streamlined calculations for minimum and maximum values using the new utility function. - Simplified the material mapping functions to improve maintainability and reduce complexity. - **Bug Fixes** - Adjusted error handling and control flow for grid value settings based on magnetic field data. --- .../Acts/MagneticField/BFieldMapUtils.hpp | 15 +- Core/include/Acts/Utilities/Helpers.hpp | 45 +++ Core/src/MagneticField/BFieldMapUtils.cpp | 269 +++++++----------- Core/src/Material/MaterialMapUtils.cpp | 98 ++----- 4 files changed, 175 insertions(+), 252 deletions(-) diff --git a/Core/include/Acts/MagneticField/BFieldMapUtils.hpp b/Core/include/Acts/MagneticField/BFieldMapUtils.hpp index 09ac0035d9b..4aa1934699b 100644 --- a/Core/include/Acts/MagneticField/BFieldMapUtils.hpp +++ b/Core/include/Acts/MagneticField/BFieldMapUtils.hpp @@ -75,7 +75,7 @@ fieldMapRZ(const std::function binsRZ, std::array nBinsRZ)>& localToGlobalBin, std::vector rPos, std::vector zPos, - std::vector bField, + const std::vector& bField, double lengthUnit = UnitConstants::mm, double BFieldUnit = UnitConstants::T, bool firstQuadrant = false); @@ -137,7 +137,7 @@ fieldMapXYZ( std::array nBinsXYZ)>& localToGlobalBin, std::vector xPos, std::vector yPos, - std::vector zPos, std::vector bField, + std::vector zPos, const std::vector& bField, double lengthUnit = UnitConstants::mm, double BFieldUnit = UnitConstants::T, bool firstOctant = false); @@ -145,17 +145,18 @@ fieldMapXYZ( /// creates a field mapper by sampling grid points from the analytical /// solenoid field. /// -/// @param rlim pair of r bounds -/// @param zlim pair of z bounds -/// @param nbins pair of bin counts +/// @param rLim pair of r bounds +/// @param zLim pair of z bounds +/// @param nBins pair of bin counts /// @param field the solenoid field instance /// /// @return A field map instance for use in interpolation. Acts::InterpolatedBFieldMap< Acts::Grid, Acts::Axis>> -solenoidFieldMap(std::pair rlim, std::pair zlim, - std::pair nbins, +solenoidFieldMap(const std::pair& rLim, + const std::pair& zLim, + const std::pair& nBins, const SolenoidBField& field); } // namespace Acts diff --git a/Core/include/Acts/Utilities/Helpers.hpp b/Core/include/Acts/Utilities/Helpers.hpp index 44bf23d9157..56d9b36481f 100644 --- a/Core/include/Acts/Utilities/Helpers.hpp +++ b/Core/include/Acts/Utilities/Helpers.hpp @@ -221,4 +221,49 @@ struct overloaded : Ts... { template overloaded(Ts...) -> overloaded; +namespace detail { + +/// Computes the minimum, maximum, and bin count for a given vector of values. +/// +/// This function processes a vector of doubles to compute: +/// - The minimum value (@c xMin) +/// - The maximum value (@c xMax), adjusted to include an additional bin +/// - The bin count (@c xBinCount) based on the number of unique values +/// +/// The computation is performed as follows: +/// 1. Sorts the input vector using @c std::ranges::sort to prepare for uniqueness. +/// 2. Determines the number of unique values using @c std::unique and calculates the bin count. +/// 3. Calculates the minimum and maximum using @c std::ranges::minmax. +/// 4. Adjusts the maximum to include an additional bin by adding the bin step +/// size. +/// +/// @param xPos A reference to a vector of doubles. +/// @return A tuple containing: +/// - The minimum value (double) +/// - The adjusted maximum value (double) +/// - The bin count (std::size_t) +/// +/// @note The vector xPos will be modified during the call. +inline auto getMinMaxAndBinCount(std::vector& xPos) { + // sort the values for unique() + std::ranges::sort(xPos); + + // get the number of bins over unique values + auto it = std::unique(xPos.begin(), xPos.end()); + const std::size_t xBinCount = std::distance(xPos.begin(), it); + + // get the minimum and maximum + auto [xMin, xMax] = std::ranges::minmax(xPos); + + // calculate maxima (add one last bin, because bin value always corresponds to + // left boundary) + const double stepX = (xMax - xMin) / static_cast(xBinCount - 1); + xMax += stepX; + + // Return all values as a tuple + return std::make_tuple(xMin, xMax, xBinCount); +} + +} // namespace detail + } // namespace Acts diff --git a/Core/src/MagneticField/BFieldMapUtils.cpp b/Core/src/MagneticField/BFieldMapUtils.cpp index b5b92191589..f5b1b410ffe 100644 --- a/Core/src/MagneticField/BFieldMapUtils.cpp +++ b/Core/src/MagneticField/BFieldMapUtils.cpp @@ -12,12 +12,14 @@ #include "Acts/MagneticField/SolenoidBField.hpp" #include "Acts/Utilities/Axis.hpp" #include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Utilities/VectorHelpers.hpp" #include "Acts/Utilities/detail/grid_helper.hpp" #include #include +#include #include #include #include @@ -35,100 +37,73 @@ Acts::fieldMapRZ( std::array nBinsRZ)>& localToGlobalBin, std::vector rPos, std::vector zPos, - std::vector bField, double lengthUnit, double BFieldUnit, + const std::vector& bField, double lengthUnit, double BFieldUnit, bool firstQuadrant) { // [1] Create Grid - // sort the values - std::ranges::sort(rPos); - std::ranges::sort(zPos); - // Get unique values - rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - rPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsR = rPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum. We just sorted the vectors, so these are just - // the first and last elements. - double rMin = rPos[0]; - double zMin = zPos[0]; - double rMax = rPos[nBinsR - 1]; - double zMax = zPos[nBinsZ - 1]; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; - zMax += stepZ; + const auto [rMin, rMax, rBinCount] = detail::getMinMaxAndBinCount(rPos); + auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos); + + const std::size_t nBinsR = rBinCount; + std::size_t nBinsZ = zBinCount; + if (firstQuadrant) { zMin = -zPos[nBinsZ - 1]; - nBinsZ = static_cast(2. * nBinsZ - 1); + nBinsZ = 2 * nBinsZ - 1; } // Create the axis for the grid - Acts::Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); - Acts::Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); + Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); + Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); // Create the grid - Grid grid(Type, std::move(rAxis), std::move(zAxis)); + Grid grid(Type, std::move(rAxis), std::move(zAxis)); using Grid_t = decltype(grid); // [2] Set the bField values + const std::array nIndices = {{rBinCount, zBinCount}}; for (std::size_t i = 1; i <= nBinsR; ++i) { for (std::size_t j = 1; j <= nBinsZ; ++j) { - std::array nIndices = {{rPos.size(), zPos.size()}}; Grid_t::index_t indices = {{i, j}}; + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one if (firstQuadrant) { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one - std::size_t n = - std::abs(static_cast(j) - static_cast(zPos.size())); - Grid_t::index_t indicesFirstQuadrant = {{i - 1, n}}; + std::size_t n = std::abs(static_cast(j) - + static_cast(zBinCount)); grid.atLocalBins(indices) = - bField.at(localToGlobalBin(indicesFirstQuadrant, nIndices)) * - BFieldUnit; + bField.at(localToGlobalBin({{i - 1, n}}, nIndices)) * BFieldUnit; } else { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one grid.atLocalBins(indices) = bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices)) * BFieldUnit; } } } - grid.setExteriorBins(Acts::Vector2::Zero()); + grid.setExteriorBins(Vector2::Zero()); - // [3] Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { - return Acts::Vector2(perp(pos), pos.z()); + // [3] Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { + return Vector2(perp(pos), pos.z()); }; - // [4] Create the transformation for the bfield - // map (Br,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector2& field, - const Acts::Vector3& pos) { - double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y(); - double cos_phi = 0, sin_phi = 0; - if (r_sin_theta_2 > std::numeric_limits::min()) { - double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2); - cos_phi = pos.x() * inv_r_sin_theta; - sin_phi = pos.y() * inv_r_sin_theta; - } else { - cos_phi = 1.; - sin_phi = 0.; + // [4] Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector2& field, const Vector3& pos) { + const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y(); + double cosPhi = 1.; + double sinPhi = 0.; + + if (rSinTheta2 > std::numeric_limits::min()) { + const double invRsinTheta = 1. / std::sqrt(rSinTheta2); + cosPhi = pos.x() * invRsinTheta; + sinPhi = pos.y() * invRsinTheta; } - return Acts::Vector3(field.x() * cos_phi, field.x() * sin_phi, field.y()); + + return Vector3(field.x() * cosPhi, field.x() * sinPhi, field.y()); }; - // [5] Create the mapper & BField Service - // create field mapping - return Acts::InterpolatedBFieldMap( + // [5] Create the mapper & BField Service create field mapping + return InterpolatedBFieldMap( {transformPos, transformBField, std::move(grid)}); } @@ -141,88 +116,58 @@ Acts::fieldMapXYZ( std::array nBinsXYZ)>& localToGlobalBin, std::vector xPos, std::vector yPos, - std::vector zPos, std::vector bField, + std::vector zPos, const std::vector& bField, double lengthUnit, double BFieldUnit, bool firstOctant) { // [1] Create Grid - // Sort the values - std::ranges::sort(xPos); - std::ranges::sort(yPos); - std::ranges::sort(zPos); - // Get unique values - xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end()); - yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - xPos.shrink_to_fit(); - yPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsX = xPos.size(); - std::size_t nBinsY = yPos.size(); - std::size_t nBinsZ = zPos.size(); + auto [xMin, xMax, xBinCount] = detail::getMinMaxAndBinCount(xPos); + auto [yMin, yMax, yBinCount] = detail::getMinMaxAndBinCount(yPos); + auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos); - // Create the axis for the grid - // get minima and maximia. We just sorted the vectors, so these are just the - // first and last elements. - double xMin = xPos[0]; - double yMin = yPos[0]; - double zMin = zPos[0]; - // get maxima - double xMax = xPos[nBinsX - 1]; - double yMax = yPos[nBinsY - 1]; - double zMax = zPos[nBinsZ - 1]; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::abs(yMax - yMin) / (nBinsY - 1); - double stepX = std::abs(xMax - xMin) / (nBinsX - 1); - xMax += stepX; - yMax += stepY; - zMax += stepZ; + std::size_t nBinsX = xBinCount; + std::size_t nBinsY = yBinCount; + std::size_t nBinsZ = zBinCount; - // If only the first octant is given if (firstOctant) { xMin = -xPos[nBinsX - 1]; - yMin = -yPos[nBinsY - 1]; - zMin = -zPos[nBinsZ - 1]; nBinsX = 2 * nBinsX - 1; + yMin = -yPos[nBinsY - 1]; nBinsY = 2 * nBinsY - 1; + zMin = -zPos[nBinsZ - 1]; nBinsZ = 2 * nBinsZ - 1; } - Acts::Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); - Acts::Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); - Acts::Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); + + Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); + Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); + Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); // Create the grid Grid grid(Type, std::move(xAxis), std::move(yAxis), std::move(zAxis)); using Grid_t = decltype(grid); // [2] Set the bField values + const std::array nIndices = { + {xBinCount, yBinCount, zBinCount}}; + + auto calcAbsDiff = [](std::size_t val, std::size_t binCount) { + return std::abs(static_cast(val) - + static_cast(binCount)); + }; + for (std::size_t i = 1; i <= nBinsX; ++i) { for (std::size_t j = 1; j <= nBinsY; ++j) { for (std::size_t k = 1; k <= nBinsZ; ++k) { Grid_t::index_t indices = {{i, j, k}}; - std::array nIndices = { - {xPos.size(), yPos.size(), zPos.size()}}; + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one if (firstOctant) { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one - std::size_t m = - std::abs(static_cast(i) - (static_cast(xPos.size()))); - std::size_t n = - std::abs(static_cast(j) - (static_cast(yPos.size()))); - std::size_t l = - std::abs(static_cast(k) - (static_cast(zPos.size()))); - Grid_t::index_t indicesFirstOctant = {{m, n, l}}; + const std::size_t l = calcAbsDiff(i, xBinCount); + const std::size_t m = calcAbsDiff(j, yBinCount); + const std::size_t n = calcAbsDiff(k, zBinCount); grid.atLocalBins(indices) = - bField.at(localToGlobalBin(indicesFirstOctant, nIndices)) * - BFieldUnit; - + bField.at(localToGlobalBin({{l, m, n}}, nIndices)) * BFieldUnit; } else { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one grid.atLocalBins(indices) = bField.at(localToGlobalBin({{i - 1, j - 1, k - 1}}, nIndices)) * BFieldUnit; @@ -230,74 +175,67 @@ Acts::fieldMapXYZ( } } } - grid.setExteriorBins(Acts::Vector3::Zero()); + grid.setExteriorBins(Vector3::Zero()); - // [3] Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { return pos; }; + // [3] Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { return pos; }; - // [4] Create the transformation for the bfield - // map (Bx,By,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector3& field, - const Acts::Vector3& /*pos*/) { return field; }; + // [4] Create the transformation for the BField map (Bx,By,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector3& field, const Vector3& /*pos*/) { + return field; + }; - // [5] Create the mapper & BField Service - // create field mapping - return Acts::InterpolatedBFieldMap( + // [5] Create the mapper & BField Service create field mapping + return InterpolatedBFieldMap( {transformPos, transformBField, std::move(grid)}); } Acts::InterpolatedBFieldMap< Acts::Grid, Acts::Axis>> -Acts::solenoidFieldMap(std::pair rlim, - std::pair zlim, - std::pair nbins, +Acts::solenoidFieldMap(const std::pair& rLim, + const std::pair& zLim, + const std::pair& nBins, const SolenoidBField& field) { - auto [rMin, rMax] = rlim; - auto [zMin, zMax] = zlim; - const auto [nBinsR, nBinsZ] = nbins; + auto [rMin, rMax] = rLim; + auto [zMin, zMax] = zLim; + const auto [nBinsR, nBinsZ] = nBins; double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; zMax += stepZ; // Create the axis for the grid - Acts::Axis rAxis(rMin, rMax, nBinsR); - Acts::Axis zAxis(zMin, zMax, nBinsZ); + Axis rAxis(rMin, rMax, nBinsR); + Axis zAxis(zMin, zMax, nBinsZ); // Create the grid - Grid grid(Type, std::move(rAxis), std::move(zAxis)); + Grid grid(Type, std::move(rAxis), std::move(zAxis)); using Grid_t = decltype(grid); - // Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { - return Acts::Vector2(perp(pos), pos.z()); + // Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { + return Vector2(perp(pos), pos.z()); }; - // Create the transformation for the bfield - // map (Br,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector2& bfield, - const Acts::Vector3& pos) { - double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y(); - double cos_phi = 0, sin_phi = 0; - if (r_sin_theta_2 > std::numeric_limits::min()) { - double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2); - cos_phi = pos.x() * inv_r_sin_theta; - sin_phi = pos.y() * inv_r_sin_theta; - } else { - cos_phi = 1.; - sin_phi = 0.; + // Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector2& bField, const Vector3& pos) { + const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y(); + double cosPhi = 1.; + double sinPhi = 0.; + + if (rSinTheta2 > std::numeric_limits::min()) { + const double invRsinTheta = 1. / std::sqrt(rSinTheta2); + cosPhi = pos.x() * invRsinTheta; + sinPhi = pos.y() * invRsinTheta; } - return Acts::Vector3(bfield.x() * cos_phi, bfield.x() * sin_phi, - bfield.y()); + + return Vector3(bField.x() * cosPhi, bField.x() * sinPhi, bField.y()); }; - // iterate over all bins, set their value to the solenoid value - // at their lower left position + // iterate over all bins, set their value to the solenoid value at their lower + // left position for (std::size_t i = 0; i <= nBinsR + 1; i++) { for (std::size_t j = 0; j <= nBinsZ + 1; j++) { Grid_t::index_t index({i, j}); @@ -314,9 +252,8 @@ Acts::solenoidFieldMap(std::pair rlim, } } - // Create the mapper & BField Service - // create field mapping - Acts::InterpolatedBFieldMap map( + // Create the mapper & BField Service create field mapping + InterpolatedBFieldMap map( {transformPos, transformBField, std::move(grid)}); return map; } diff --git a/Core/src/Material/MaterialMapUtils.cpp b/Core/src/Material/MaterialMapUtils.cpp index b51b1a3e223..8f4d23cd71e 100644 --- a/Core/src/Material/MaterialMapUtils.cpp +++ b/Core/src/Material/MaterialMapUtils.cpp @@ -12,6 +12,7 @@ #include "Acts/Material/Material.hpp" #include "Acts/Utilities/Axis.hpp" #include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Helpers.hpp" #include #include @@ -43,31 +44,8 @@ auto Acts::materialMapperRZ( } // [2] Create Grid - // sort the values - std::ranges::sort(rPos); - std::ranges::sort(zPos); - // Get unique values - rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - rPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsR = rPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum - auto minMaxR = std::minmax_element(rPos.begin(), rPos.end()); - auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end()); - double rMin = *minMaxR.first; - double zMin = *minMaxZ.first; - double rMax = *minMaxR.second; - double zMax = *minMaxZ.second; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; - zMax += stepZ; + const auto [rMin, rMax, nBinsR] = detail::getMinMaxAndBinCount(rPos); + const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos); // Create the axis for the grid Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); @@ -79,13 +57,13 @@ auto Acts::materialMapperRZ( using Grid_t = decltype(grid); // [3] Set the material values + const std::array nIndices = {{nBinsR, nBinsZ}}; for (std::size_t i = 1; i <= nBinsR; ++i) { for (std::size_t j = 1; j <= nBinsZ; ++j) { - std::array nIndices = {{rPos.size(), zPos.size()}}; Grid_t::index_t indices = {{i, j}}; - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one grid.atLocalBins(indices) = materialVector.at( materialVectorToGridMapper({{i - 1, j - 1}}, nIndices)); } @@ -95,14 +73,12 @@ auto Acts::materialMapperRZ( 0., 0., 0.; grid.setExteriorBins(vec); - // [4] Create the transformation for the position - // map (x,y,z) -> (r,z) + // [4] Create the transformation for the position map (x,y,z) -> (r,z) auto transformPos = [](const Vector3& pos) { return Vector2(perp(pos), pos.z()); }; - // [5] Create the mapper & BField Service - // create material mapping + // [5] Create the mapper & BField Service create material mapping return MaterialMapper(transformPos, std::move(grid)); } @@ -125,44 +101,11 @@ auto Acts::materialMapperXYZ( } // [2] Create Grid - // Sort the values - std::ranges::sort(xPos); - std::ranges::sort(yPos); - std::ranges::sort(zPos); - // Get unique values - xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end()); - yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - xPos.shrink_to_fit(); - yPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsX = xPos.size(); - std::size_t nBinsY = yPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum - auto minMaxX = std::minmax_element(xPos.begin(), xPos.end()); - auto minMaxY = std::minmax_element(yPos.begin(), yPos.end()); - auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end()); - // Create the axis for the grid - // get minima - double xMin = *minMaxX.first; - double yMin = *minMaxY.first; - double zMin = *minMaxZ.first; - // get maxima - double xMax = *minMaxX.second; - double yMax = *minMaxY.second; - double zMax = *minMaxZ.second; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::abs(yMax - yMin) / (nBinsY - 1); - double stepX = std::abs(xMax - xMin) / (nBinsX - 1); - xMax += stepX; - yMax += stepY; - zMax += stepZ; + const auto [xMin, xMax, nBinsX] = detail::getMinMaxAndBinCount(xPos); + const auto [yMin, yMax, nBinsY] = detail::getMinMaxAndBinCount(yPos); + const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos); + // Create the axis for the grid Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); @@ -172,15 +115,14 @@ auto Acts::materialMapperXYZ( using Grid_t = decltype(grid); // [3] Set the bField values + const std::array nIndices = {{nBinsX, nBinsY, nBinsZ}}; for (std::size_t i = 1; i <= nBinsX; ++i) { for (std::size_t j = 1; j <= nBinsY; ++j) { for (std::size_t k = 1; k <= nBinsZ; ++k) { Grid_t::index_t indices = {{i, j, k}}; - std::array nIndices = { - {xPos.size(), yPos.size(), zPos.size()}}; - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one grid.atLocalBins(indices) = materialVector.at( materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices)); } @@ -191,11 +133,9 @@ auto Acts::materialMapperXYZ( 0., 0., 0.; grid.setExteriorBins(vec); - // [4] Create the transformation for the position - // map (x,y,z) -> (r,z) + // [4] Create the transformation for the position map (x,y,z) -> (r,z) auto transformPos = [](const Vector3& pos) { return pos; }; - // [5] Create the mapper & BField Service - // create material mapping + // [5] Create the mapper & BField Service create material mapping return MaterialMapper(transformPos, std::move(grid)); }