ECCE @ EIC 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 #include <iostream>
11 #include <limits>
15 
18 
20  const std::function<size_t(std::array<size_t, 2> binsRZ,
21  std::array<size_t, 2> nBinsRZ)>&
22  materialVectorToGridMapper,
23  std::vector<double> rPos, std::vector<double> zPos,
24  std::vector<Acts::Material> material, double lengthUnit)
26  detail::EquidistantAxis>> {
27  // [1] Decompose material
28  std::vector<ActsVectorF<5>> materialVector;
29  materialVector.reserve(material.size());
30 
31  for (Material& mat : material) {
32  materialVector.push_back(mat.classificationNumbers());
33  }
34 
35  // [2] Create Grid
36  // sort the values
37  std::sort(rPos.begin(), rPos.end());
38  std::sort(zPos.begin(), zPos.end());
39  // Get unique values
40  rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end());
41  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
42  rPos.shrink_to_fit();
43  zPos.shrink_to_fit();
44  // get the number of bins
45  size_t nBinsR = rPos.size();
46  size_t nBinsZ = zPos.size();
47 
48  // get the minimum and maximum
49  auto minMaxR = std::minmax_element(rPos.begin(), rPos.end());
50  auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
51  double rMin = *minMaxR.first;
52  double zMin = *minMaxZ.first;
53  double rMax = *minMaxR.second;
54  double zMax = *minMaxZ.second;
55  // calculate maxima (add one last bin, because bin value always corresponds to
56  // left boundary)
57  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
58  double stepR = std::fabs(rMax - rMin) / (nBinsR - 1);
59  rMax += stepR;
60  zMax += stepZ;
61 
62  // Create the axis for the grid
63  detail::EquidistantAxis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR);
64  detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
65 
66  // Create the grid
68  detail::EquidistantAxis>;
69  Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
70 
71  // [3] Set the material values
72  for (size_t i = 1; i <= nBinsR; ++i) {
73  for (size_t j = 1; j <= nBinsZ; ++j) {
74  std::array<size_t, 2> nIndices = {{rPos.size(), zPos.size()}};
75  Grid_t::index_t indices = {{i, j}};
76  // std::vectors begin with 0 and we do not want the user needing to
77  // take underflow or overflow bins in account this is why we need to
78  // subtract by one
79  grid.atLocalBins(indices) = materialVector.at(
80  materialVectorToGridMapper({{i - 1, j - 1}}, nIndices));
81  }
82  }
83  ActsVectorF<5> vec;
85  0., 0., 0.;
86  grid.setExteriorBins(vec);
87 
88  // [4] Create the transformation for the position
89  // map (x,y,z) -> (r,z)
90  auto transformPos = [](const Vector3D& pos) {
91  return Vector2D(perp(pos), pos.z());
92  };
93 
94  // [5] Create the mapper & BField Service
95  // create material mapping
97  detail::EquidistantAxis>>(transformPos,
98  std::move(grid));
99 }
100 
102  const std::function<size_t(std::array<size_t, 3> binsXYZ,
103  std::array<size_t, 3> nBinsXYZ)>&
104  materialVectorToGridMapper,
105  std::vector<double> xPos, std::vector<double> yPos,
106  std::vector<double> zPos, std::vector<Material> material, double lengthUnit)
107  -> MaterialMapper<
109  detail::EquidistantAxis, detail::EquidistantAxis>> {
110  // [1] Decompose material
111  std::vector<ActsVectorF<5>> materialVector;
112  materialVector.reserve(material.size());
113 
114  for (Material& mat : material) {
115  materialVector.push_back(mat.classificationNumbers());
116  }
117 
118  // [2] Create Grid
119  // Sort the values
120  std::sort(xPos.begin(), xPos.end());
121  std::sort(yPos.begin(), yPos.end());
122  std::sort(zPos.begin(), zPos.end());
123  // Get unique values
124  xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end());
125  yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end());
126  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
127  xPos.shrink_to_fit();
128  yPos.shrink_to_fit();
129  zPos.shrink_to_fit();
130  // get the number of bins
131  size_t nBinsX = xPos.size();
132  size_t nBinsY = yPos.size();
133  size_t nBinsZ = zPos.size();
134 
135  // get the minimum and maximum
136  auto minMaxX = std::minmax_element(xPos.begin(), xPos.end());
137  auto minMaxY = std::minmax_element(yPos.begin(), yPos.end());
138  auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
139  // Create the axis for the grid
140  // get minima
141  double xMin = *minMaxX.first;
142  double yMin = *minMaxY.first;
143  double zMin = *minMaxZ.first;
144  // get maxima
145  double xMax = *minMaxX.second;
146  double yMax = *minMaxY.second;
147  double zMax = *minMaxZ.second;
148  // calculate maxima (add one last bin, because bin value always corresponds to
149  // left boundary)
150  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
151  double stepY = std::fabs(yMax - yMin) / (nBinsY - 1);
152  double stepX = std::fabs(xMax - xMin) / (nBinsX - 1);
153  xMax += stepX;
154  yMax += stepY;
155  zMax += stepZ;
156 
157  detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX);
158  detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY);
159  detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ);
160  // Create the grid
161  using Grid_t = detail::Grid<ActsVectorF<5>, detail::EquidistantAxis,
162  detail::EquidistantAxis, detail::EquidistantAxis>;
163  Grid_t grid(
164  std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis)));
165 
166  // [3] Set the bField values
167  for (size_t i = 1; i <= nBinsX; ++i) {
168  for (size_t j = 1; j <= nBinsY; ++j) {
169  for (size_t k = 1; k <= nBinsZ; ++k) {
170  Grid_t::index_t indices = {{i, j, k}};
171  std::array<size_t, 3> nIndices = {
172  {xPos.size(), yPos.size(), zPos.size()}};
173  // std::vectors begin with 0 and we do not want the user needing to
174  // take underflow or overflow bins in account this is why we need to
175  // subtract by one
176  grid.atLocalBins(indices) = materialVector.at(
177  materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices));
178  }
179  }
180  }
181  ActsVectorF<5> vec;
183  0., 0., 0.;
184  grid.setExteriorBins(vec);
185 
186  // [4] Create the transformation for the position
187  // map (x,y,z) -> (r,z)
188  auto transformPos = [](const Vector3D& pos) { return pos; };
189 
190  // [5] Create the mapper & BField Service
191  // create material mapping
192  return MaterialMapper<
193  detail::Grid<ActsVectorF<5>, detail::EquidistantAxis,
194  detail::EquidistantAxis, detail::EquidistantAxis>>(
195  transformPos, std::move(grid));
196 }