ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MaterialGridHelper.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MaterialGridHelper.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 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 
11 Acts::Grid2D Acts::createGrid(std::array<double, 3> gridAxis1,
12  std::array<double, 3> gridAxis2) {
13  // get the number of bins
14  size_t nBinsAxis1 = gridAxis1[2];
15  size_t nBinsAxis2 = gridAxis2[2];
16 
17  // get the minimum and maximum
18  double minAxis1 = gridAxis1[0];
19  double minAxis2 = gridAxis2[0];
20  double maxAxis1 = gridAxis1[1];
21  double maxAxis2 = gridAxis2[1];
22  // calculate maxima (add one last bin, because bin value always corresponds
23  // to
24  // left boundary)
25  double stepAxis1 = std::fabs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1);
26  double stepAxis2 = std::fabs(maxAxis2 - minAxis2) / (nBinsAxis2 - 1);
27  maxAxis1 += stepAxis1;
28  maxAxis2 += stepAxis2;
29 
30  // Create the axis for the grid
31  Acts::EAxis axis1(minAxis1, maxAxis1, nBinsAxis1);
32  Acts::EAxis axis2(minAxis2, maxAxis2, nBinsAxis2);
33 
34  // The material mapping grid
35  return Acts::Grid2D(std::make_tuple(std::move(axis1), std::move(axis2)));
36 }
37 
38 Acts::Grid3D Acts::createGrid(std::array<double, 3> gridAxis1,
39  std::array<double, 3> gridAxis2,
40  std::array<double, 3> gridAxis3) {
41  // get the number of bins
42  size_t nBinsAxis1 = gridAxis1[2];
43  size_t nBinsAxis2 = gridAxis2[2];
44  size_t nBinsAxis3 = gridAxis3[2];
45 
46  // get the minimum and maximum
47  double minAxis1 = gridAxis1[0];
48  double minAxis2 = gridAxis2[0];
49  double minAxis3 = gridAxis3[0];
50  double maxAxis1 = gridAxis1[1];
51  double maxAxis2 = gridAxis2[1];
52  double maxAxis3 = gridAxis3[1];
53  // calculate maxima (add one last bin, because bin value always corresponds
54  // to
55  // left boundary)
56  double stepAxis1 = std::fabs(maxAxis1 - minAxis1) / (nBinsAxis1 - 1);
57  double stepAxis2 = std::fabs(maxAxis2 - minAxis2) / (nBinsAxis2 - 1);
58  double stepAxis3 = std::fabs(maxAxis3 - minAxis3) / (nBinsAxis3 - 1);
59  maxAxis1 += stepAxis1;
60  maxAxis2 += stepAxis2;
61  maxAxis3 += stepAxis3;
62 
63  // Create the axis for the grid
64  Acts::EAxis axis1(minAxis1, maxAxis1, nBinsAxis1);
65  Acts::EAxis axis2(minAxis2, maxAxis2, nBinsAxis2);
66  Acts::EAxis axis3(minAxis3, maxAxis3, nBinsAxis3);
67 
68  // The material mapping grid
69  return Acts::Grid3D(
70  std::make_tuple(std::move(axis1), std::move(axis2), std::move(axis3)));
71 }
72 
73 std::function<double(Acts::Vector3D)> Acts::globalToLocalFromBin(
74  Acts::BinningValue& type) {
75  std::function<double(Acts::Vector3D)> transfoGlobalToLocal;
76 
77  switch (type) {
78  case Acts::binX:
79  transfoGlobalToLocal = [](Acts::Vector3D pos) -> double {
80  return (pos.x());
81  };
82  break;
83 
84  case Acts::binY:
85  transfoGlobalToLocal = [](Acts::Vector3D pos) -> double {
86  return (pos.y());
87  };
88  break;
89 
90  case Acts::binR:
91  transfoGlobalToLocal = [](Acts::Vector3D pos) -> double {
93  };
94  break;
95 
96  case Acts::binPhi:
97  transfoGlobalToLocal = [](Acts::Vector3D pos) -> double {
98  return (Acts::VectorHelpers::phi(pos));
99  };
100  break;
101 
102  case Acts::binZ:
103  transfoGlobalToLocal = [](Acts::Vector3D pos) -> double {
104  return (pos.z());
105  };
106  break;
107 
108  case Acts::binRPhi:
109  case Acts::binEta:
110  case Acts::binH:
111  case Acts::binMag:
112  case Acts::binValues:
113  throw std::invalid_argument("Incorrect bin, should be x,y,z,r,phi,z");
114  break;
115  }
116  return transfoGlobalToLocal;
117 }
118 
120  const Acts::BinUtility& bins,
121  std::function<Acts::Vector2D(Acts::Vector3D)>& transfoGlobalToLocal) {
122  auto bu = bins.binningData();
123  // First we nee to create the 2 axis
124  std::array<double, 3> gridAxis1;
125  std::array<double, 3> gridAxis2;
126 
127  bool isCartesian = false;
128  bool isCylindrical = false;
129 
130  for (size_t b = 0; b < bu.size(); b++) {
131  if (bu[b].binvalue == Acts::binX || bu[b].binvalue == Acts::binY) {
132  isCartesian = true;
133  }
134  if (bu[b].binvalue == Acts::binR || bu[b].binvalue == Acts::binPhi) {
135  isCylindrical = true;
136  }
137  }
138  if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
139  throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
140  }
141 
142  gridAxis1[0] = bu[0].min;
143  gridAxis1[1] = bu[0].max;
144  gridAxis1[2] = bu[0].bins();
145 
146  gridAxis2[0] = bu[1].min;
147  gridAxis2[1] = bu[1].max;
148  gridAxis2[2] = bu[1].bins();
149 
150  std::function<double(Acts::Vector3D)> coord1 =
151  globalToLocalFromBin(bu[0].binvalue);
152  std::function<double(Acts::Vector3D)> coord2 =
153  globalToLocalFromBin(bu[1].binvalue);
154 
155  transfoGlobalToLocal = [coord1,
156  coord2](Acts::Vector3D pos) -> Acts::Vector2D {
157  return {coord1(pos), coord2(pos)};
158  };
159  return (Acts::createGrid(std::move(gridAxis1), std::move(gridAxis2)));
160 }
161 
163  const Acts::BinUtility& bins,
164  std::function<Acts::Vector3D(Acts::Vector3D)>& transfoGlobalToLocal) {
165  auto bu = bins.binningData();
166  // First we nee to create the 3 axis
167  std::array<double, 3> gridAxis1;
168  std::array<double, 3> gridAxis2;
169  std::array<double, 3> gridAxis3;
170 
171  bool isCartesian = false;
172  bool isCylindrical = false;
173 
174  for (size_t b = 0; b < bu.size(); b++) {
175  if (bu[b].binvalue == Acts::binX || bu[b].binvalue == Acts::binY) {
176  isCartesian = true;
177  }
178  if (bu[b].binvalue == Acts::binR || bu[b].binvalue == Acts::binPhi) {
179  isCylindrical = true;
180  }
181  }
182  if (!(isCartesian || isCylindrical) || (isCylindrical && isCartesian)) {
183  throw std::invalid_argument("Incorrect bin, should be x,y,z or r,phi,z");
184  }
185 
186  gridAxis1[0] = bu[0].min;
187  gridAxis1[1] = bu[0].max;
188  gridAxis1[2] = bu[0].bins();
189 
190  gridAxis2[0] = bu[1].min;
191  gridAxis2[1] = bu[1].max;
192  gridAxis2[2] = bu[1].bins();
193 
194  gridAxis3[0] = bu[2].min;
195  gridAxis3[1] = bu[2].max;
196  gridAxis3[2] = bu[2].bins();
197 
198  std::function<double(Acts::Vector3D)> coord1 =
199  globalToLocalFromBin(bu[0].binvalue);
200  std::function<double(Acts::Vector3D)> coord2 =
201  globalToLocalFromBin(bu[1].binvalue);
202  std::function<double(Acts::Vector3D)> coord3 =
203  globalToLocalFromBin(bu[2].binvalue);
204 
205  transfoGlobalToLocal = [coord1, coord2,
206  coord3](Acts::Vector3D pos) -> Acts::Vector3D {
207  return {coord1(pos), coord2(pos), coord3(pos)};
208  };
209  return (Acts::createGrid(std::move(gridAxis1), std::move(gridAxis2),
210  std::move(gridAxis3)));
211 }
212 
214  Acts::Grid2D& grid, const Acts::RecordedMaterialPoint& mPoints,
215  std::function<Acts::Vector2D(Acts::Vector3D)>& transfoGlobalToLocal) {
216  // Walk over each point
217  for (const auto& rm : mPoints) {
218  // Search for fitting grid point and accumulate
219  Acts::Grid2D::index_t index =
220  grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(rm.second));
221  grid.atLocalBins(index).accumulate(rm.first);
222  }
223 
224  // Build material grid
225  // Re-build the axes
228  Acts::Grid2D::index_t nBins = grid.numLocalBins();
229 
230  Acts::EAxis axis1(min[0], max[0], nBins[0]);
231  Acts::EAxis axis2(min[1], max[1], nBins[1]);
232 
233  // Build the grid and fill it with data
234  Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
235  for (size_t index = 0; index < grid.size(); index++) {
236  mGrid.at(index) = grid.at(index).average().classificationNumbers();
237  }
238 
239  return mGrid;
240 }
241 
243  Acts::Grid3D& grid, const Acts::RecordedMaterialPoint& mPoints,
244  std::function<Acts::Vector3D(Acts::Vector3D)>& transfoGlobalToLocal) {
245  // Walk over each point
246  for (const auto& rm : mPoints) {
247  // Search for fitting grid point and accumulate
248  Acts::Grid3D::index_t index =
249  grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(rm.second));
250  grid.atLocalBins(index).accumulate(rm.first);
251  }
252 
253  // Build material grid
254  // Re-build the axes
257  Acts::Grid3D::index_t nBins = grid.numLocalBins();
258 
259  Acts::EAxis axis1(min[0], max[0], nBins[0]);
260  Acts::EAxis axis2(min[1], max[1], nBins[1]);
261  Acts::EAxis axis3(min[2], max[2], nBins[2]);
262 
263  // Build the grid and fill it with data
264  Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
265  for (size_t index = 0; index < grid.size(); index++) {
266  mGrid.at(index) = grid.at(index).average().classificationNumbers();
267  }
268  return mGrid;
269 }