src/Material/MaterialGridHelper.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/MaterialGridHelper.hpp"
0010
0011 #include "Acts/Utilities/BinningData.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <memory>
0016 #include <stdexcept>
0017 #include <tuple>
0018 #include <utility>
0019 #include <vector>
0020
0021 Acts::Grid2D Acts::createGrid(Acts::MaterialGridAxisData gridAxis1,
0022                               Acts::MaterialGridAxisData gridAxis2) {
0023   // get the number of bins
0024   std::size_t nBinsAxis1 = std::get<2>(gridAxis1);
0025   std::size_t nBinsAxis2 = std::get<2>(gridAxis2);
0026
0027   // get the minimum and maximum
0028   double minAxis1 = std::get<0>(gridAxis1);
0029   double minAxis2 = std::get<0>(gridAxis2);
0030   double maxAxis1 = std::get<1>(gridAxis1);
0031   double maxAxis2 = std::get<1>(gridAxis2);
0032   // calculate maxima (add one last bin, because bin value always corresponds
0033   // to
0034   // left boundary)
0035   double stepAxis1 = std::abs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1);
0036   double stepAxis2 = std::abs(maxAxis2 - minAxis2) / (nBinsAxis2 - 1);
0037   maxAxis1 += stepAxis1;
0038   maxAxis2 += stepAxis2;
0039
0040   // Create the axis for the grid
0041   Acts::EAxis axis1(minAxis1, maxAxis1, nBinsAxis1);
0042   Acts::EAxis axis2(minAxis2, maxAxis2, nBinsAxis2);
0043
0044   // The material mapping grid
0045   return Acts::Grid2D(std::make_tuple(std::move(axis1), std::move(axis2)));
0046 }
0047
0048 Acts::Grid3D Acts::createGrid(Acts::MaterialGridAxisData gridAxis1,
0049                               Acts::MaterialGridAxisData gridAxis2,
0050                               Acts::MaterialGridAxisData gridAxis3) {
0051   // get the number of bins
0052   std::size_t nBinsAxis1 = std::get<2>(gridAxis1);
0053   std::size_t nBinsAxis2 = std::get<2>(gridAxis2);
0054   std::size_t nBinsAxis3 = std::get<2>(gridAxis3);
0055
0056   // get the minimum and maximum
0057   double minAxis1 = std::get<0>(gridAxis1);
0058   double minAxis2 = std::get<0>(gridAxis2);
0059   double minAxis3 = std::get<0>(gridAxis3);
0060   double maxAxis1 = std::get<1>(gridAxis1);
0061   double maxAxis2 = std::get<1>(gridAxis2);
0062   double maxAxis3 = std::get<1>(gridAxis3);
0063   // calculate maxima (add one last bin, because bin value always corresponds
0064   // to
0065   // left boundary)
0066   double stepAxis1 =
0067       std::abs(maxAxis1 - minAxis1) / std::max(nBinsAxis1 - 1, std::size_t{1});
0068   double stepAxis2 =
0069       std::abs(maxAxis2 - minAxis2) / std::max(nBinsAxis2 - 1, std::size_t{1});
0070   double stepAxis3 =
0071       std::abs(maxAxis3 - minAxis3) / std::max(nBinsAxis3 - 1, std::size_t{1});
0072   maxAxis1 += stepAxis1;
0073   maxAxis2 += stepAxis2;
0074   maxAxis3 += stepAxis3;
0075
0076   // Create the axis for the grid
0077   Acts::EAxis axis1(minAxis1, maxAxis1, nBinsAxis1);
0078   Acts::EAxis axis2(minAxis2, maxAxis2, nBinsAxis2);
0079   Acts::EAxis axis3(minAxis3, maxAxis3, nBinsAxis3);
0080
0081   // The material mapping grid
0082   return Acts::Grid3D(
0083       std::make_tuple(std::move(axis1), std::move(axis2), std::move(axis3)));
0084 }
0085
0086 std::function<double(Acts::Vector3)> Acts::globalToLocalFromBin(
0087     Acts::AxisDirection& type) {
0088   std::function<double(Acts::Vector3)> transfoGlobalToLocal;
0089
0090   switch (type) {
0091     case Acts::AxisDirection::AxisX:
0092       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0093         return (pos.x());
0094       };
0095       break;
0096
0097     case Acts::AxisDirection::AxisY:
0098       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0099         return (pos.y());
0100       };
0101       break;
0102
0103     case Acts::AxisDirection::AxisR:
0104       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0105         return (Acts::VectorHelpers::perp(pos));
0106       };
0107       break;
0108
0109     case Acts::AxisDirection::AxisPhi:
0110       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0111         return (Acts::VectorHelpers::phi(pos));
0112       };
0113       break;
0114
0115     case Acts::AxisDirection::AxisZ:
0116       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0117         return (pos.z());
0118       };
0119       break;
0120
0121       // case Acts::AxisDirection::AxisRPhi:
0122       // case Acts::AxisDirection::AxisEta:
0123       // case Acts::AxisDirection::AxisTheta:
0124       // case Acts::AxisDirection::AxisMag:
0125     default:
0126       throw std::invalid_argument("Incorrect bin, should be x,y,z,r,phi");
0127   }
0128   return transfoGlobalToLocal;
0129 }
0130
0131 Acts::Grid2D Acts::createGrid2D(
0132     const Acts::BinUtility& bins,
0133     std::function<Acts::Vector2(Acts::Vector3)>& transfoGlobalToLocal) {
0134   auto bu = bins.binningData();
0135
0136   bool isCartesian = false;
0137   bool isCylindrical = false;
0138
0139   for (std::size_t b = 0; b < bu.size(); b++) {
0140     if (bu[b].binvalue == Acts::AxisDirection::AxisX ||
0141         bu[b].binvalue == Acts::AxisDirection::AxisY) {
0142       isCartesian = true;
0143     }
0144     if (bu[b].binvalue == Acts::AxisDirection::AxisR ||
0145         bu[b].binvalue == Acts::AxisDirection::AxisPhi) {
0146       isCylindrical = true;
0147     }
0148   }
0149   if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
0150     throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
0151   }
0152
0153   // First we need to create the 2 axis
0154   MaterialGridAxisData gridAxis1{bu[0].min, bu[0].max, bu[0].bins()};
0155   MaterialGridAxisData gridAxis2{bu[1].min, bu[1].max, bu[1].bins()};
0156
0157   std::function<double(Acts::Vector3)> coord1 =
0158       globalToLocalFromBin(bu[0].binvalue);
0159   std::function<double(Acts::Vector3)> coord2 =
0160       globalToLocalFromBin(bu[1].binvalue);
0161   Transform3 transfo = bins.transform().inverse();
0162   transfoGlobalToLocal = [coord1, coord2,
0163                           transfo](Acts::Vector3 pos) -> Acts::Vector2 {
0164     pos = transfo * pos;
0165     return {coord1(pos), coord2(pos)};
0166   };
0167   return Acts::createGrid(gridAxis1, gridAxis2);
0168 }
0169
0170 Acts::Grid3D Acts::createGrid3D(
0171     const Acts::BinUtility& bins,
0172     std::function<Acts::Vector3(Acts::Vector3)>& transfoGlobalToLocal) {
0173   auto bu = bins.binningData();
0174   // First we need to create the 3 axis
0175
0176   bool isCartesian = false;
0177   bool isCylindrical = false;
0178
0179   for (std::size_t b = 0; b < bu.size(); b++) {
0180     if (bu[b].binvalue == Acts::AxisDirection::AxisX ||
0181         bu[b].binvalue == Acts::AxisDirection::AxisY) {
0182       isCartesian = true;
0183     }
0184     if (bu[b].binvalue == Acts::AxisDirection::AxisR ||
0185         bu[b].binvalue == Acts::AxisDirection::AxisPhi) {
0186       isCylindrical = true;
0187     }
0188   }
0189   if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
0190     throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
0191   }
0192
0193   MaterialGridAxisData gridAxis1{bu[0].min, bu[0].max, bu[0].bins()};
0194
0195   MaterialGridAxisData gridAxis2{bu[1].min, bu[1].max, bu[1].bins()};
0196
0197   MaterialGridAxisData gridAxis3{bu[2].min, bu[2].max, bu[2].bins()};
0198
0199   std::function<double(Acts::Vector3)> coord1 =
0200       globalToLocalFromBin(bu[0].binvalue);
0201   std::function<double(Acts::Vector3)> coord2 =
0202       globalToLocalFromBin(bu[1].binvalue);
0203   std::function<double(Acts::Vector3)> coord3 =
0204       globalToLocalFromBin(bu[2].binvalue);
0205   Transform3 transfo = bins.transform().inverse();
0206
0207   transfoGlobalToLocal = [coord1, coord2, coord3,
0208                           transfo](Acts::Vector3 pos) -> Acts::Vector3 {
0209     pos = transfo * pos;
0210     return {coord1(pos), coord2(pos), coord3(pos)};
0211   };
0212   return Acts::createGrid(gridAxis1, gridAxis2, gridAxis3);
0213 }
0214
0215 Acts::MaterialGrid2D Acts::mapMaterialPoints(Acts::Grid2D& grid) {
0216   // Build material grid
0217   // Re-build the axes
0218   Acts::Grid2D::point_t min = grid.minPosition();
0219   Acts::Grid2D::point_t max = grid.maxPosition();
0220   Acts::Grid2D::index_t nBins = grid.numLocalBins();
0221
0222   Acts::EAxis axis1(min[0], max[0], nBins[0]);
0223   Acts::EAxis axis2(min[1], max[1], nBins[1]);
0224
0225   // Fill the material Grid by averaging the material in the 2D grid
0226   Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
0227   for (std::size_t index = 0; index < grid.size(); index++) {
0228     mGrid.at(index) = grid.at(index).average().parameters();
0229   }
0230
0231   return mGrid;
0232 }
0233
0234 Acts::MaterialGrid3D Acts::mapMaterialPoints(Acts::Grid3D& grid) {
0235   // Build material grid
0236   // Re-build the axes
0237   Acts::Grid3D::point_t min = grid.minPosition();
0238   Acts::Grid3D::point_t max = grid.maxPosition();
0239   Acts::Grid3D::index_t nBins = grid.numLocalBins();
0240
0241   Acts::EAxis axis1(min[0], max[0], nBins[0]);
0242   Acts::EAxis axis2(min[1], max[1], nBins[1]);
0243   Acts::EAxis axis3(min[2], max[2], nBins[2]);
0244
0245   // Fill the material Grid by averaging the material in the 3D grid
0246   Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
0247   for (std::size_t index = 0; index < grid.size(); index++) {
0248     mGrid.at(index) = grid.at(index).average().parameters();
0249   }
0250   return mGrid;
0251 }