Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MaterialMapUtils.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MaterialMapUtils.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
15 
16 #include <algorithm>
17 #include <cmath>
18 #include <cstddef>
19 #include <initializer_list>
20 #include <limits>
21 #include <set>
22 #include <tuple>
23 #include <utility>
24 
27 
29  const std::function<size_t(std::array<size_t, 2> binsRZ,
30  std::array<size_t, 2> nBinsRZ)>&
31  materialVectorToGridMapper,
32  std::vector<double> rPos, std::vector<double> zPos,
33  const std::vector<Acts::Material>& material, double lengthUnit)
34  -> MaterialMapper<
36  detail::EquidistantAxis>> {
37  // [1] Decompose material
38  std::vector<Material::ParametersVector> materialVector;
39  materialVector.reserve(material.size());
40 
41  for (const Material& mat : material) {
42  materialVector.push_back(mat.parameters());
43  }
44 
45  // [2] Create Grid
46  // sort the values
47  std::sort(rPos.begin(), rPos.end());
48  std::sort(zPos.begin(), zPos.end());
49  // Get unique values
50  rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end());
51  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
52  rPos.shrink_to_fit();
53  zPos.shrink_to_fit();
54  // get the number of bins
55  size_t nBinsR = rPos.size();
56  size_t nBinsZ = zPos.size();
57 
58  // get the minimum and maximum
59  auto minMaxR = std::minmax_element(rPos.begin(), rPos.end());
60  auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
61  double rMin = *minMaxR.first;
62  double zMin = *minMaxZ.first;
63  double rMax = *minMaxR.second;
64  double zMax = *minMaxZ.second;
65  // calculate maxima (add one last bin, because bin value always corresponds to
66  // left boundary)
67  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
68  double stepR = std::fabs(rMax - rMin) / (nBinsR - 1);
69  rMax += stepR;
70  zMax += stepZ;
71 
72  // Create the axis for the grid
73  detail::EquidistantAxis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR);
74  detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
75 
76  // Create the grid
77  using Grid_t = detail::Grid<Material::ParametersVector,
78  detail::EquidistantAxis, detail::EquidistantAxis>;
79  Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
80 
81  // [3] Set the material values
82  for (size_t i = 1; i <= nBinsR; ++i) {
83  for (size_t j = 1; j <= nBinsZ; ++j) {
84  std::array<size_t, 2> nIndices = {{rPos.size(), zPos.size()}};
85  Grid_t::index_t indices = {{i, j}};
86  // std::vectors begin with 0 and we do not want the user needing to
87  // take underflow or overflow bins in account this is why we need to
88  // subtract by one
89  grid.atLocalBins(indices) = materialVector.at(
90  materialVectorToGridMapper({{i - 1, j - 1}}, nIndices));
91  }
92  }
93  Material::ParametersVector vec;
94  vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
95  0., 0., 0.;
96  grid.setExteriorBins(vec);
97 
98  // [4] Create the transformation for the position
99  // map (x,y,z) -> (r,z)
100  auto transformPos = [](const Vector3& pos) {
101  return Vector2(perp(pos), pos.z());
102  };
103 
104  // [5] Create the mapper & BField Service
105  // create material mapping
106  return MaterialMapper<
107  detail::Grid<Material::ParametersVector, detail::EquidistantAxis,
108  detail::EquidistantAxis>>(transformPos, std::move(grid));
109 }
110 
112  const std::function<size_t(std::array<size_t, 3> binsXYZ,
113  std::array<size_t, 3> nBinsXYZ)>&
114  materialVectorToGridMapper,
115  std::vector<double> xPos, std::vector<double> yPos,
116  std::vector<double> zPos, const std::vector<Material>& material,
117  double lengthUnit)
118  -> MaterialMapper<
119  detail::Grid<Material::ParametersVector, detail::EquidistantAxis,
120  detail::EquidistantAxis, detail::EquidistantAxis>> {
121  // [1] Decompose material
122  std::vector<Material::ParametersVector> materialVector;
123  materialVector.reserve(material.size());
124 
125  for (const Material& mat : material) {
126  materialVector.push_back(mat.parameters());
127  }
128 
129  // [2] Create Grid
130  // Sort the values
131  std::sort(xPos.begin(), xPos.end());
132  std::sort(yPos.begin(), yPos.end());
133  std::sort(zPos.begin(), zPos.end());
134  // Get unique values
135  xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end());
136  yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end());
137  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
138  xPos.shrink_to_fit();
139  yPos.shrink_to_fit();
140  zPos.shrink_to_fit();
141  // get the number of bins
142  size_t nBinsX = xPos.size();
143  size_t nBinsY = yPos.size();
144  size_t nBinsZ = zPos.size();
145 
146  // get the minimum and maximum
147  auto minMaxX = std::minmax_element(xPos.begin(), xPos.end());
148  auto minMaxY = std::minmax_element(yPos.begin(), yPos.end());
149  auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
150  // Create the axis for the grid
151  // get minima
152  double xMin = *minMaxX.first;
153  double yMin = *minMaxY.first;
154  double zMin = *minMaxZ.first;
155  // get maxima
156  double xMax = *minMaxX.second;
157  double yMax = *minMaxY.second;
158  double zMax = *minMaxZ.second;
159  // calculate maxima (add one last bin, because bin value always corresponds to
160  // left boundary)
161  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
162  double stepY = std::fabs(yMax - yMin) / (nBinsY - 1);
163  double stepX = std::fabs(xMax - xMin) / (nBinsX - 1);
164  xMax += stepX;
165  yMax += stepY;
166  zMax += stepZ;
167 
168  detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX);
169  detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY);
170  detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
171  // Create the grid
172  using Grid_t =
173  detail::Grid<Material::ParametersVector, detail::EquidistantAxis,
174  detail::EquidistantAxis, detail::EquidistantAxis>;
175  Grid_t grid(
176  std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis)));
177 
178  // [3] Set the bField values
179  for (size_t i = 1; i <= nBinsX; ++i) {
180  for (size_t j = 1; j <= nBinsY; ++j) {
181  for (size_t k = 1; k <= nBinsZ; ++k) {
182  Grid_t::index_t indices = {{i, j, k}};
183  std::array<size_t, 3> nIndices = {
184  {xPos.size(), yPos.size(), zPos.size()}};
185  // std::vectors begin with 0 and we do not want the user needing to
186  // take underflow or overflow bins in account this is why we need to
187  // subtract by one
188  grid.atLocalBins(indices) = materialVector.at(
189  materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices));
190  }
191  }
192  }
193  Material::ParametersVector vec;
194  vec << std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
195  0., 0., 0.;
196  grid.setExteriorBins(vec);
197 
198  // [4] Create the transformation for the position
199  // map (x,y,z) -> (r,z)
200  auto transformPos = [](const Vector3& pos) { return pos; };
201 
202  // [5] Create the mapper & BField Service
203  // create material mapping
204  return MaterialMapper<
205  detail::Grid<Material::ParametersVector, detail::EquidistantAxis,
206  detail::EquidistantAxis, detail::EquidistantAxis>>(
207  transformPos, std::move(grid));
208 }