src/Material/MaterialMapUtils.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019-2020 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/Material/MaterialMapUtils.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Material/Material.hpp"
0013 #include "Acts/Utilities/Grid.hpp"
0014 #include "Acts/Utilities/detail/Axis.hpp"
0015
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <cstddef>
0019 #include <initializer_list>
0020 #include <limits>
0021 #include <set>
0022 #include <tuple>
0023 #include <utility>
0024
0025 using Acts::VectorHelpers::perp;
0026 using Acts::VectorHelpers::phi;
0027
0028 auto Acts::materialMapperRZ(
0029     const std::function<std::size_t(std::array<std::size_t, 2> binsRZ,
0030                                     std::array<std::size_t, 2> nBinsRZ)>&
0031         materialVectorToGridMapper,
0032     std::vector<double> rPos, std::vector<double> zPos,
0033     const std::vector<Acts::Material>& material, double lengthUnit)
0034     -> MaterialMapper<Grid<Material::ParametersVector, detail::EquidistantAxis,
0035                            detail::EquidistantAxis>> {
0036   // [1] Decompose material
0037   std::vector<Material::ParametersVector> materialVector;
0038   materialVector.reserve(material.size());
0039
0040   for (const Material& mat : material) {
0041     materialVector.push_back(mat.parameters());
0042   }
0043
0044   // [2] Create Grid
0045   // sort the values
0046   std::sort(rPos.begin(), rPos.end());
0047   std::sort(zPos.begin(), zPos.end());
0048   // Get unique values
0049   rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end());
0050   zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
0051   rPos.shrink_to_fit();
0052   zPos.shrink_to_fit();
0053   // get the number of bins
0054   std::size_t nBinsR = rPos.size();
0055   std::size_t nBinsZ = zPos.size();
0056
0057   // get the minimum and maximum
0058   auto minMaxR = std::minmax_element(rPos.begin(), rPos.end());
0059   auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
0060   double rMin = *minMaxR.first;
0061   double zMin = *minMaxZ.first;
0062   double rMax = *minMaxR.second;
0063   double zMax = *minMaxZ.second;
0064   // calculate maxima (add one last bin, because bin value always corresponds to
0065   // left boundary)
0066   double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
0067   double stepR = std::fabs(rMax - rMin) / (nBinsR - 1);
0068   rMax += stepR;
0069   zMax += stepZ;
0070
0071   // Create the axis for the grid
0072   detail::EquidistantAxis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR);
0073   detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0074
0075   // Create the grid
0076   using Grid_t = Grid<Material::ParametersVector, detail::EquidistantAxis,
0077                       detail::EquidistantAxis>;
0078   Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
0079
0080   // [3] Set the material values
0081   for (std::size_t i = 1; i <= nBinsR; ++i) {
0082     for (std::size_t j = 1; j <= nBinsZ; ++j) {
0083       std::array<std::size_t, 2> nIndices = {{rPos.size(), zPos.size()}};
0084       Grid_t::index_t indices = {{i, j}};
0085       // std::vectors begin with 0 and we do not want the user needing to
0086       // take underflow or overflow bins in account this is why we need to
0087       // subtract by one
0088       grid.atLocalBins(indices) = materialVector.at(
0089           materialVectorToGridMapper({{i - 1, j - 1}}, nIndices));
0090     }
0091   }
0092   Material::ParametersVector vec;
0093   vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
0094       0., 0., 0.;
0095   grid.setExteriorBins(vec);
0096
0097   // [4] Create the transformation for the position
0098   // map (x,y,z) -> (r,z)
0099   auto transformPos = [](const Vector3& pos) {
0100     return Vector2(perp(pos), pos.z());
0101   };
0102
0103   // [5] Create the mapper & BField Service
0104   // create material mapping
0105   return MaterialMapper<Grid<Material::ParametersVector,
0106                              detail::EquidistantAxis, detail::EquidistantAxis>>(
0107       transformPos, std::move(grid));
0108 }
0109
0110 auto Acts::materialMapperXYZ(
0111     const std::function<std::size_t(std::array<std::size_t, 3> binsXYZ,
0112                                     std::array<std::size_t, 3> nBinsXYZ)>&
0113         materialVectorToGridMapper,
0114     std::vector<double> xPos, std::vector<double> yPos,
0115     std::vector<double> zPos, const std::vector<Material>& material,
0116     double lengthUnit)
0117     -> MaterialMapper<Grid<Material::ParametersVector, detail::EquidistantAxis,
0118                            detail::EquidistantAxis, detail::EquidistantAxis>> {
0119   // [1] Decompose material
0120   std::vector<Material::ParametersVector> materialVector;
0121   materialVector.reserve(material.size());
0122
0123   for (const Material& mat : material) {
0124     materialVector.push_back(mat.parameters());
0125   }
0126
0127   // [2] Create Grid
0128   // Sort the values
0129   std::sort(xPos.begin(), xPos.end());
0130   std::sort(yPos.begin(), yPos.end());
0131   std::sort(zPos.begin(), zPos.end());
0132   // Get unique values
0133   xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end());
0134   yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end());
0135   zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
0136   xPos.shrink_to_fit();
0137   yPos.shrink_to_fit();
0138   zPos.shrink_to_fit();
0139   // get the number of bins
0140   std::size_t nBinsX = xPos.size();
0141   std::size_t nBinsY = yPos.size();
0142   std::size_t nBinsZ = zPos.size();
0143
0144   // get the minimum and maximum
0145   auto minMaxX = std::minmax_element(xPos.begin(), xPos.end());
0146   auto minMaxY = std::minmax_element(yPos.begin(), yPos.end());
0147   auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
0148   // Create the axis for the grid
0149   // get minima
0150   double xMin = *minMaxX.first;
0151   double yMin = *minMaxY.first;
0152   double zMin = *minMaxZ.first;
0153   // get maxima
0154   double xMax = *minMaxX.second;
0155   double yMax = *minMaxY.second;
0156   double zMax = *minMaxZ.second;
0157   // calculate maxima (add one last bin, because bin value always corresponds to
0158   // left boundary)
0159   double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
0160   double stepY = std::fabs(yMax - yMin) / (nBinsY - 1);
0161   double stepX = std::fabs(xMax - xMin) / (nBinsX - 1);
0162   xMax += stepX;
0163   yMax += stepY;
0164   zMax += stepZ;
0165
0166   detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX);
0167   detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY);
0168   detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0169   // Create the grid
0170   using Grid_t = Grid<Material::ParametersVector, detail::EquidistantAxis,
0171                       detail::EquidistantAxis, detail::EquidistantAxis>;
0172   Grid_t grid(
0173       std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis)));
0174
0175   // [3] Set the bField values
0176   for (std::size_t i = 1; i <= nBinsX; ++i) {
0177     for (std::size_t j = 1; j <= nBinsY; ++j) {
0178       for (std::size_t k = 1; k <= nBinsZ; ++k) {
0179         Grid_t::index_t indices = {{i, j, k}};
0180         std::array<std::size_t, 3> nIndices = {
0181             {xPos.size(), yPos.size(), zPos.size()}};
0182         // std::vectors begin with 0 and we do not want the user needing to
0183         // take underflow or overflow bins in account this is why we need to
0184         // subtract by one
0185         grid.atLocalBins(indices) = materialVector.at(
0186             materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices));
0187       }
0188     }
0189   }
0190   Material::ParametersVector vec;
0191   vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
0192       0., 0., 0.;
0193   grid.setExteriorBins(vec);
0194
0195   // [4] Create the transformation for the position
0196   // map (x,y,z) -> (r,z)
0197   auto transformPos = [](const Vector3& pos) { return pos; };
0198
0199   // [5] Create the mapper & BField Service
0200   // create material mapping
0201   return MaterialMapper<
0202       Grid<Material::ParametersVector, detail::EquidistantAxis,
0203            detail::EquidistantAxis, detail::EquidistantAxis>>(transformPos,
0204                                                               std::move(grid));
0205 }