src/Material/MaterialMapUtils.cpp

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/Material/MaterialMapUtils.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Material/Material.hpp"
0013 #include "Acts/Utilities/Axis.hpp"
0014 #include "Acts/Utilities/Grid.hpp"
0015 #include "Acts/Utilities/Helpers.hpp"
0016
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <cstddef>
0020 #include <initializer_list>
0021 #include <limits>
0022 #include <set>
0023 #include <tuple>
0024 #include <utility>
0025
0026 using Acts::VectorHelpers::perp;
0027 using Acts::VectorHelpers::phi;
0028
0029 auto Acts::materialMapperRZ(
0030     const std::function<std::size_t(std::array<std::size_t, 2> binsRZ,
0031                                     std::array<std::size_t, 2> nBinsRZ)>&
0032         materialVectorToGridMapper,
0033     std::vector<double> rPos, std::vector<double> zPos,
0034     const std::vector<Acts::Material>& material, double lengthUnit)
0035     -> MaterialMapper<
0036         Grid<Material::ParametersVector, Axis<Acts::AxisType::Equidistant>,
0037              Axis<Acts::AxisType::Equidistant>>> {
0038   // [1] Decompose material
0039   std::vector<Material::ParametersVector> materialVector;
0040   materialVector.reserve(material.size());
0041
0042   for (const Material& mat : material) {
0043     materialVector.push_back(mat.parameters());
0044   }
0045
0046   // [2] Create Grid
0047   const auto [rMin, rMax, nBinsR] = detail::getMinMaxAndBinCount(rPos);
0048   const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos);
0049
0050   // Create the axis for the grid
0051   Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR);
0052   Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0053
0054   // Create the grid
0055   Grid grid(Type<Material::ParametersVector>, std::move(rAxis),
0056             std::move(zAxis));
0057   using Grid_t = decltype(grid);
0058
0059   // [3] Set the material values
0060   const std::array<std::size_t, 2> nIndices = {{nBinsR, nBinsZ}};
0061   for (std::size_t i = 1; i <= nBinsR; ++i) {
0062     for (std::size_t j = 1; j <= nBinsZ; ++j) {
0063       Grid_t::index_t indices = {{i, j}};
0064       // std::vectors begin with 0 and we do not want the user needing to take
0065       // underflow or overflow bins in account this is why we need to subtract
0066       // by one
0067       grid.atLocalBins(indices) = materialVector.at(
0068           materialVectorToGridMapper({{i - 1, j - 1}}, nIndices));
0069     }
0070   }
0071   Material::ParametersVector vec;
0072   vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
0073       0., 0., 0.;
0074   grid.setExteriorBins(vec);
0075
0076   // [4] Create the transformation for the position map (x,y,z) -> (r,z)
0077   auto transformPos = [](const Vector3& pos) {
0078     return Vector2(perp(pos), pos.z());
0079   };
0080
0081   // [5] Create the mapper & BField Service create material mapping
0082   return MaterialMapper(transformPos, std::move(grid));
0083 }
0084
0085 auto Acts::materialMapperXYZ(
0086     const std::function<std::size_t(std::array<std::size_t, 3> binsXYZ,
0087                                     std::array<std::size_t, 3> nBinsXYZ)>&
0088         materialVectorToGridMapper,
0089     std::vector<double> xPos, std::vector<double> yPos,
0090     std::vector<double> zPos, const std::vector<Material>& material,
0091     double lengthUnit)
0092     -> MaterialMapper<Grid<
0093         Material::ParametersVector, Axis<Acts::AxisType::Equidistant>,
0094         Axis<Acts::AxisType::Equidistant>, Axis<Acts::AxisType::Equidistant>>> {
0095   // [1] Decompose material
0096   std::vector<Material::ParametersVector> materialVector;
0097   materialVector.reserve(material.size());
0098
0099   for (const Material& mat : material) {
0100     materialVector.push_back(mat.parameters());
0101   }
0102
0103   // [2] Create Grid
0104   const auto [xMin, xMax, nBinsX] = detail::getMinMaxAndBinCount(xPos);
0105   const auto [yMin, yMax, nBinsY] = detail::getMinMaxAndBinCount(yPos);
0106   const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos);
0107
0108   // Create the axis for the grid
0109   Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX);
0110   Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY);
0111   Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
0112   // Create the grid
0113   Grid grid(Type<Material::ParametersVector>, std::move(xAxis),
0114             std::move(yAxis), std::move(zAxis));
0115   using Grid_t = decltype(grid);
0116
0117   // [3] Set the bField values
0118   const std::array<std::size_t, 3> nIndices = {{nBinsX, nBinsY, nBinsZ}};
0119   for (std::size_t i = 1; i <= nBinsX; ++i) {
0120     for (std::size_t j = 1; j <= nBinsY; ++j) {
0121       for (std::size_t k = 1; k <= nBinsZ; ++k) {
0122         Grid_t::index_t indices = {{i, j, k}};
0123         // std::vectors begin with 0 and we do not want the user needing to take
0124         // underflow or overflow bins in account this is why we need to subtract
0125         // by one
0126         grid.atLocalBins(indices) = materialVector.at(
0127             materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices));
0128       }
0129     }
0130   }
0131   Material::ParametersVector vec;
0132   vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
0133       0., 0., 0.;
0134   grid.setExteriorBins(vec);
0135
0136   // [4] Create the transformation for the position map (x,y,z) -> (r,z)
0137   auto transformPos = [](const Vector3& pos) { return pos; };
0138
0139   // [5] Create the mapper & BField Service create material mapping
0140   return MaterialMapper(transformPos, std::move(grid));
0141 }