src/Material/MaterialGridHelper.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 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/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::fabs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1);
0036   double stepAxis2 = std::fabs(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::fabs(maxAxis1 - minAxis1) / std::max(nBinsAxis1 - 1, std::size_t(1));
0068   double stepAxis2 =
0069       std::fabs(maxAxis2 - minAxis2) / std::max(nBinsAxis2 - 1, std::size_t(1));
0070   double stepAxis3 =
0071       std::fabs(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::BinningValue& type) {
0088   std::function<double(Acts::Vector3)> transfoGlobalToLocal;
0089
0090   switch (type) {
0091     case Acts::binX:
0092       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0093         return (pos.x());
0094       };
0095       break;
0096
0097     case Acts::binY:
0098       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0099         return (pos.y());
0100       };
0101       break;
0102
0103     case Acts::binR:
0104       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0105         return (Acts::VectorHelpers::perp(pos));
0106       };
0107       break;
0108
0109     case Acts::binPhi:
0110       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0111         return (Acts::VectorHelpers::phi(pos));
0112       };
0113       break;
0114
0115     case Acts::binZ:
0116       transfoGlobalToLocal = [](const Acts::Vector3& pos) -> double {
0117         return (pos.z());
0118       };
0119       break;
0120
0121       // case Acts::binRPhi:
0122       // case Acts::binEta:
0123       // case Acts::binH:
0124       // case Acts::binMag:
0125       // case Acts::binValues:
0126     default:
0127       throw std::invalid_argument("Incorrect bin, should be x,y,z,r,phi");
0128   }
0129   return transfoGlobalToLocal;
0130 }
0131
0132 Acts::Grid2D Acts::createGrid2D(
0133     const Acts::BinUtility& bins,
0134     std::function<Acts::Vector2(Acts::Vector3)>& transfoGlobalToLocal) {
0135   auto bu = bins.binningData();
0136
0137   bool isCartesian = false;
0138   bool isCylindrical = false;
0139
0140   for (std::size_t b = 0; b < bu.size(); b++) {
0141     if (bu[b].binvalue == Acts::binX || bu[b].binvalue == Acts::binY) {
0142       isCartesian = true;
0143     }
0144     if (bu[b].binvalue == Acts::binR || bu[b].binvalue == Acts::binPhi) {
0145       isCylindrical = true;
0146     }
0147   }
0148   if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
0149     throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
0150   }
0151
0152   // First we need to create the 2 axis
0153   MaterialGridAxisData gridAxis1{bu[0].min, bu[0].max, bu[0].bins()};
0154   MaterialGridAxisData gridAxis2{bu[1].min, bu[1].max, bu[1].bins()};
0155
0156   std::function<double(Acts::Vector3)> coord1 =
0157       globalToLocalFromBin(bu[0].binvalue);
0158   std::function<double(Acts::Vector3)> coord2 =
0159       globalToLocalFromBin(bu[1].binvalue);
0160   Transform3 transfo = bins.transform().inverse();
0161   transfoGlobalToLocal = [coord1, coord2,
0162                           transfo](Acts::Vector3 pos) -> Acts::Vector2 {
0163     pos = transfo * pos;
0164     return {coord1(pos), coord2(pos)};
0165   };
0166   return Acts::createGrid(gridAxis1, gridAxis2);
0167 }
0168
0169 Acts::Grid3D Acts::createGrid3D(
0170     const Acts::BinUtility& bins,
0171     std::function<Acts::Vector3(Acts::Vector3)>& transfoGlobalToLocal) {
0172   auto bu = bins.binningData();
0173   // First we need to create the 3 axis
0174
0175   bool isCartesian = false;
0176   bool isCylindrical = false;
0177
0178   for (std::size_t b = 0; b < bu.size(); b++) {
0179     if (bu[b].binvalue == Acts::binX || bu[b].binvalue == Acts::binY) {
0180       isCartesian = true;
0181     }
0182     if (bu[b].binvalue == Acts::binR || bu[b].binvalue == Acts::binPhi) {
0183       isCylindrical = true;
0184     }
0185   }
0186   if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
0187     throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
0188   }
0189
0190   MaterialGridAxisData gridAxis1{bu[0].min, bu[0].max, bu[0].bins()};
0191
0192   MaterialGridAxisData gridAxis2{bu[1].min, bu[1].max, bu[1].bins()};
0193
0194   MaterialGridAxisData gridAxis3{bu[2].min, bu[2].max, bu[2].bins()};
0195
0196   std::function<double(Acts::Vector3)> coord1 =
0197       globalToLocalFromBin(bu[0].binvalue);
0198   std::function<double(Acts::Vector3)> coord2 =
0199       globalToLocalFromBin(bu[1].binvalue);
0200   std::function<double(Acts::Vector3)> coord3 =
0201       globalToLocalFromBin(bu[2].binvalue);
0202   Transform3 transfo = bins.transform().inverse();
0203
0204   transfoGlobalToLocal = [coord1, coord2, coord3,
0205                           transfo](Acts::Vector3 pos) -> Acts::Vector3 {
0206     pos = transfo * pos;
0207     return {coord1(pos), coord2(pos), coord3(pos)};
0208   };
0209   return Acts::createGrid(gridAxis1, gridAxis2, gridAxis3);
0210 }
0211
0212 Acts::MaterialGrid2D Acts::mapMaterialPoints(Acts::Grid2D& grid) {
0213   // Build material grid
0214   // Re-build the axes
0215   Acts::Grid2D::point_t min = grid.minPosition();
0216   Acts::Grid2D::point_t max = grid.maxPosition();
0217   Acts::Grid2D::index_t nBins = grid.numLocalBins();
0218
0219   Acts::EAxis axis1(min[0], max[0], nBins[0]);
0220   Acts::EAxis axis2(min[1], max[1], nBins[1]);
0221
0222   // Fill the material Grid by averaging the material in the 2D grid
0223   Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
0224   for (std::size_t index = 0; index < grid.size(); index++) {
0225     mGrid.at(index) = grid.at(index).average().parameters();
0226   }
0227
0228   return mGrid;
0229 }
0230
0231 Acts::MaterialGrid3D Acts::mapMaterialPoints(Acts::Grid3D& grid) {
0232   // Build material grid
0233   // Re-build the axes
0234   Acts::Grid3D::point_t min = grid.minPosition();
0235   Acts::Grid3D::point_t max = grid.maxPosition();
0236   Acts::Grid3D::index_t nBins = grid.numLocalBins();
0237
0238   Acts::EAxis axis1(min[0], max[0], nBins[0]);
0239   Acts::EAxis axis2(min[1], max[1], nBins[1]);
0240   Acts::EAxis axis3(min[2], max[2], nBins[2]);
0241
0242   // Fill the material Grid by averaging the material in the 3D grid
0243   Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
0244   for (std::size_t index = 0; index < grid.size(); index++) {
0245     mGrid.at(index) = grid.at(index).average().parameters();
0246   }
0247   return mGrid;
0248 }