23 std::array<size_t, 2> nBinsRZ)>&
25 std::vector<double> rPos, std::vector<double> zPos,
26 std::vector<Acts::Vector2D>
bField,
double lengthUnit,
27 double BFieldUnit,
bool firstQuadrant) {
30 std::sort(rPos.begin(), rPos.end());
31 std::sort(zPos.begin(), zPos.end());
33 rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end());
34 zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
38 size_t nBinsR = rPos.size();
39 size_t nBinsZ = zPos.size();
42 auto minMaxR = std::minmax_element(rPos.begin(), rPos.end());
43 auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
44 double rMin = *minMaxR.first;
45 double zMin = *minMaxZ.first;
46 double rMax = *minMaxR.second;
47 double zMax = *minMaxZ.second;
50 double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
51 double stepR = std::fabs(rMax - rMin) / (nBinsR - 1);
55 zMin = -(*minMaxZ.second);
56 nBinsZ = 2. * nBinsZ - 1;
68 Acts::detail::EquidistantAxis>;
69 Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
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}};
80 size_t n =
std::abs(
int(j) -
int(zPos.size()));
81 Grid_t::index_t indicesFirstQuadrant = {{i - 1, n}};
83 grid.atLocalBins(indices) =
84 bField.at(localToGlobalBin(indicesFirstQuadrant, nIndices)) *
90 grid.atLocalBins(indices) =
91 bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices)) *
96 grid.setExteriorBins(Acts::Vector2D::Zero());
106 auto transformBField = [](
const Acts::Vector2D& field,
108 double r_sin_theta_2 =
pos.x() *
pos.x() +
pos.y() *
pos.y();
109 double cos_phi, sin_phi;
111 double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
112 cos_phi =
pos.x() * inv_r_sin_theta;
113 sin_phi =
pos.y() * inv_r_sin_theta;
118 return Acts::Vector3D(field.x() * cos_phi, field.x() * sin_phi, field.y());
132 std::array<size_t, 3> nBinsXYZ)>&
134 std::vector<double> xPos, std::vector<double> yPos,
135 std::vector<double> zPos, std::vector<Acts::Vector3D>
bField,
136 double lengthUnit,
double BFieldUnit,
bool firstOctant) {
139 std::sort(xPos.begin(), xPos.end());
140 std::sort(yPos.begin(), yPos.end());
141 std::sort(zPos.begin(), zPos.end());
143 xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end());
144 yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end());
145 zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
146 xPos.shrink_to_fit();
147 yPos.shrink_to_fit();
148 zPos.shrink_to_fit();
150 size_t nBinsX = xPos.size();
151 size_t nBinsY = yPos.size();
152 size_t nBinsZ = zPos.size();
155 auto minMaxX = std::minmax_element(xPos.begin(), xPos.end());
156 auto minMaxY = std::minmax_element(yPos.begin(), yPos.end());
157 auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end());
160 double xMin = *minMaxX.first;
161 double yMin = *minMaxY.first;
162 double zMin = *minMaxZ.first;
164 double xMax = *minMaxX.second;
165 double yMax = *minMaxY.second;
166 double zMax = *minMaxZ.second;
169 double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
170 double stepY = std::fabs(yMax - yMin) / (nBinsY - 1);
171 double stepX = std::fabs(xMax - xMin) / (nBinsX - 1);
178 xMin = -*minMaxX.second;
179 yMin = -*minMaxY.second;
180 zMin = -*minMaxZ.second;
181 nBinsX = 2 * nBinsX - 1;
182 nBinsY = 2 * nBinsY - 1;
183 nBinsZ = 2 * nBinsZ - 1;
185 Acts::detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit,
187 Acts::detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit,
189 Acts::detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit,
195 Acts::detail::EquidistantAxis>;
197 std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis)));
200 for (
size_t i = 1; i <= nBinsX; ++i) {
201 for (
size_t j = 1; j <= nBinsY; ++j) {
203 Grid_t::index_t indices = {{i, j,
k}};
204 std::array<size_t, 3> nIndices = {
205 {xPos.size(), yPos.size(), zPos.size()}};
210 size_t m =
std::abs(
int(i) - (
int(xPos.size())));
211 size_t n =
std::abs(
int(j) - (
int(yPos.size())));
212 size_t l =
std::abs(
int(
k) - (
int(zPos.size())));
213 Grid_t::index_t indicesFirstOctant = {{
m,
n, l}};
215 grid.atLocalBins(indices) =
216 bField.at(localToGlobalBin(indicesFirstOctant, nIndices)) *
223 grid.atLocalBins(indices) =
224 bField.at(localToGlobalBin({{i - 1, j - 1,
k - 1}}, nIndices)) *
230 grid.setExteriorBins(Acts::Vector3D::Zero());
234 auto transformPos = [](
const Acts::Vector3D&
pos) {
return pos; };
238 auto transformBField = [](
const Acts::Vector3D& field,
239 const Acts::Vector3D& ) {
return field; };
249 Acts::detail::EquidistantAxis>>
251 std::pair<double, double> zlim,
252 std::pair<size_t, size_t> nbins,
254 double rMin, rMax, zMin, zMax;
255 std::tie(rMin, rMax) = rlim;
256 std::tie(zMin, zMax) = zlim;
259 std::tie(nBinsR, nBinsZ) = nbins;
261 double stepZ =
std::abs(zMax - zMin) / (nBinsZ - 1);
262 double stepR =
std::abs(rMax - rMin) / (nBinsR - 1);
268 Acts::detail::EquidistantAxis rAxis(rMin, rMax, nBinsR);
269 Acts::detail::EquidistantAxis zAxis(zMin, zMax, nBinsZ);
274 Acts::detail::EquidistantAxis>;
275 Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
279 auto transformPos = [](
const Acts::Vector3D&
pos) {
285 auto transformBField = [](
const Acts::Vector2D& bfield,
286 const Acts::Vector3D&
pos) {
287 double r_sin_theta_2 =
pos.x() *
pos.x() +
pos.y() *
pos.y();
288 double cos_phi, sin_phi;
290 double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
291 cos_phi =
pos.x() * inv_r_sin_theta;
292 sin_phi =
pos.y() * inv_r_sin_theta;
303 for (
size_t i = 0; i <= nBinsR + 1; i++) {
304 for (
size_t j = 0; j <= nBinsZ + 1; j++) {
305 Grid_t::index_t index({i, j});
306 if (i == 0 || j == 0 || i == nBinsR + 1 || j == nBinsZ + 1) {
308 grid.atLocalBins(index) = Grid_t::value_type(0, 0);
311 Grid_t::point_t lowerLeft = grid.lowerLeftBinEdge(index);
314 grid.atLocalBins(index) =
B;