Skip to content

Commit 33d63aa

Browse files
wdconincwdconinc
committed
feat: convert FieldMapB to using covfie
1 parent 27bf021 commit 33d63aa

File tree

2 files changed

+39
-160
lines changed

2 files changed

+39
-160
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@ include(GNUInstallDirs)
4646
# Dependencies
4747
find_package(DD4hep 1.27 REQUIRED COMPONENTS DDCore DDRec)
4848
find_package(fmt REQUIRED)
49+
find_package(covfie REQUIRED)
4950

5051
#-----------------------------------------------------------------------------------
5152
set(a_lib_name ${PROJECT_NAME})

‎src/FieldMapB.cpp‎

Lines changed: 38 additions & 160 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,14 @@
55
#include <DD4hep/FieldTypes.h>
66
#include <DD4hep/Printout.h>
77
#include <XML/Utilities.h>
8+
#include <covfie/core/algebra/affine.hpp>
9+
#include <covfie/core/backend/primitive/array.hpp>
10+
#include <covfie/core/backend/transformer/affine.hpp>
11+
#include <covfie/core/backend/transformer/linear.hpp>
12+
#include <covfie/core/backend/transformer/strided.hpp>
13+
#include <covfie/core/field.hpp>
14+
#include <covfie/core/field_view.hpp>
15+
#include <covfie/core/parameter_pack.hpp>
816

917
#include <cstdlib>
1018
#include <filesystem>
@@ -19,6 +27,14 @@ namespace fs = std::filesystem;
1927

2028
#include "FileLoaderHelper.h"
2129

30+
using fieldBrBz_t =
31+
covfie::field<covfie::backend::affine<covfie::backend::linear<covfie::backend::strided<
32+
covfie::vector::size2, covfie::backend::array<covfie::vector::float2>>>>>;
33+
34+
using fieldBxByBz_t =
35+
covfie::field<covfie::backend::affine<covfie::backend::linear<covfie::backend::strided<
36+
covfie::vector::size3, covfie::backend::array<covfie::vector::float3>>>>>;
37+
2238
using namespace dd4hep;
2339

2440
// Two coordinate options:
@@ -40,11 +56,9 @@ class FieldMapB : public dd4hep::CartesianField::Object {
4056
public:
4157
FieldMapB(const std::string& field_type_str = "magnetic",
4258
const std::string& coord_type_str = "BrBz");
43-
void Configure(std::vector<xml_comp_t> dimensions);
59+
4460
void LoadMap(const std::string& map_file, float scale);
45-
bool GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR, float* deltaZ);
46-
bool GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ, float* deltaX,
47-
float* deltaY, float* deltaZ);
61+
4862
void SetCoordTranslation(const Transform3D& tr) {
4963
coordTranslate = tr;
5064
coordTranslate_inv = tr.Inverse();
@@ -62,12 +76,10 @@ class FieldMapB : public dd4hep::CartesianField::Object {
6276
std::optional<Transform3D> coordTranslate_inv{std::nullopt};
6377
std::optional<Transform3D> fieldRot{std::nullopt}; // field rotation
6478
std::optional<Transform3D> fieldRot_inv{std::nullopt};
65-
std::vector<float> steps, mins, maxs; // B map cell info
66-
int ir, ix, iy, iz; // lookup indices
67-
float dr, dx, dy, dz; // deltas for interpolation
68-
std::vector<std::vector<std::array<float, 2>>> Bvals_RZ; // B map values: {R}, {Z}, {Br,Bz}
69-
std::vector<std::vector<std::vector<std::array<float, 3>>>>
70-
Bvals_XYZ; // B map values: {X}, {Y}, {Z}, {Bx,By,Bz}
79+
std::vector<float> steps, mins, maxs; // B map cell info
80+
81+
std::shared_ptr<fieldBrBz_t> fieldBrBz;
82+
std::shared_ptr<fieldBxByBz_t> fieldBxByBz;
7183
};
7284

7385
// constructor
@@ -94,120 +106,21 @@ FieldMapB::FieldMapB(const std::string& field_type_str, const std::string& coord
94106
}
95107
}
96108

97-
// fill field vector
98-
void FieldMapB::Configure(std::vector<xml_comp_t> dimensions) {
99-
// Fill vectors with step size, min, max for each dimension
100-
for (auto el : dimensions) {
101-
steps.push_back(el.step());
102-
mins.push_back(getAttrOrDefault<float>(el, _Unicode(min), 0));
103-
maxs.push_back(getAttrOrDefault<float>(el, _Unicode(max), 0));
104-
}
105-
106-
if (fieldCoord == FieldCoord::BrBz) {
107-
// N bins increased by 1 beyond grid size to account for edge cases at upper limits
108-
int nr = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
109-
int nz = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;
110-
111-
Bvals_RZ.resize(nr);
112-
for (auto& B2 : Bvals_RZ) {
113-
B2.resize(nz);
114-
}
115-
} else {
116-
// N bins increased by 1 beyond grid size to account for edge cases at upper limits
117-
int nx = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
118-
int ny = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;
119-
int nz = std::roundf((maxs[2] - mins[2]) / steps[2]) + 2;
120-
121-
Bvals_XYZ.resize(nx);
122-
for (auto& B3 : Bvals_XYZ) {
123-
B3.resize(ny);
124-
for (auto& B2 : B3) {
125-
B2.resize(nz);
126-
}
127-
}
128-
}
129-
}
130-
131-
// get RZ cell indices corresponding to point of interest
132-
bool FieldMapB::GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR, float* deltaZ) {
133-
// boundary check
134-
if (R > maxs[0] || R < mins[0] || Z > maxs[1] || Z < mins[1]) {
135-
return false;
136-
}
137-
138-
// get indices
139-
float idx_1_f, idx_2_f;
140-
*deltaR = std::modf((R - mins[0]) / steps[0], &idx_1_f);
141-
*deltaZ = std::modf((Z - mins[1]) / steps[1], &idx_2_f);
142-
*idxR = static_cast<int>(idx_1_f);
143-
*idxZ = static_cast<int>(idx_2_f);
144-
145-
return true;
146-
}
147-
148-
// get XYZ cell indices corresponding to point of interest
149-
bool FieldMapB::GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ,
150-
float* deltaX, float* deltaY, float* deltaZ) {
151-
// boundary check
152-
if (X > maxs[0] || X < mins[0] || Y > maxs[1] || Y < mins[1] || Z > maxs[2] || Z < mins[2]) {
153-
return false;
154-
}
155-
156-
// get indices
157-
float idx_1_f, idx_2_f, idx_3_f;
158-
*deltaX = std::modf((X - mins[0]) / steps[0], &idx_1_f);
159-
*deltaY = std::modf((Y - mins[1]) / steps[1], &idx_2_f);
160-
*deltaZ = std::modf((Z - mins[2]) / steps[2], &idx_3_f);
161-
*idxX = static_cast<int>(idx_1_f);
162-
*idxY = static_cast<int>(idx_2_f);
163-
*idxZ = static_cast<int>(idx_3_f);
164-
165-
return true;
166-
}
167-
168109
// load data
169-
void FieldMapB::LoadMap(const std::string& map_file, float scale) {
110+
void FieldMapB::LoadMap(const std::string& map_file, float) {
170111
std::string line;
171112
std::ifstream input(map_file);
172113
if (!input) {
173114
printout(ERROR, "FieldMapB", "FieldMapB Error: file " + map_file + " cannot be read.");
174115
}
175116

176-
std::vector<float> coord = {};
177-
std::vector<float> Bcomp = {};
178-
179-
while (std::getline(input, line).good()) {
180-
std::istringstream iss(line);
181-
182-
coord.clear();
183-
Bcomp.clear();
184-
185-
if (fieldCoord == FieldCoord::BrBz) {
186-
coord.resize(2);
187-
Bcomp.resize(2);
188-
iss >> coord[0] >> coord[1] >> Bcomp[0] >> Bcomp[1];
189-
190-
if (!GetIndices(coord[0], coord[1], &ir, &iz, &dr, &dz)) {
191-
printout(WARNING, "FieldMapB", "coordinates out of range, skipped it.");
192-
} else { // scale field
193-
Bvals_RZ[ir][iz] = {Bcomp[0] * scale * float(tesla), Bcomp[1] * scale * float(tesla)};
194-
}
195-
} else {
196-
coord.resize(3);
197-
Bcomp.resize(3);
198-
iss >> coord[0] >> coord[1] >> coord[2] >> Bcomp[0] >> Bcomp[1] >> Bcomp[2];
199-
200-
if (!GetIndices(coord[0], coord[1], coord[2], &ix, &iy, &iz, &dx, &dy, &dz)) {
201-
printout(WARNING, "FieldMap", "coordinates out of range, skipped it.");
202-
} else { // scale and rotate B field vector
203-
auto B = ROOT::Math::XYZPoint(Bcomp[0], Bcomp[1], Bcomp[2]);
204-
B *= scale * float(tesla);
205-
if (fieldRot.has_value()) {
206-
B = fieldRot.value() * B;
207-
}
208-
Bvals_XYZ[ix][iy][iz] = {float(B.x()), float(B.y()), float(B.z())};
209-
}
210-
}
117+
switch (fieldCoord) {
118+
case FieldCoord::BrBz:
119+
fieldBrBz = std::make_shared<fieldBrBz_t>(input);
120+
break;
121+
case FieldCoord::BxByBz:
122+
fieldBxByBz = std::make_shared<fieldBxByBz_t>(input);
123+
break;
211124
}
212125
}
213126

@@ -223,65 +136,31 @@ void FieldMapB::fieldComponents(const double* pos, double* field) {
223136
const float r = sqrt(p.x() * p.x() + p.y() * p.y());
224137
const float z = p.z();
225138

226-
if (!GetIndices(r, z, &ir, &iz, &dr, &dz)) {
227-
// out of range
228-
return;
229-
}
230-
231-
// p1 p3
232-
// p
233-
// p0 p2
234-
auto& p0 = Bvals_RZ[ir][iz];
235-
auto& p1 = Bvals_RZ[ir][iz + 1];
236-
auto& p2 = Bvals_RZ[ir + 1][iz];
237-
auto& p3 = Bvals_RZ[ir + 1][iz + 1];
238-
239-
// Bilinear interpolation
240-
float Br = p0[0] * (1 - dr) * (1 - dz) + p1[0] * (1 - dr) * dz + p2[0] * dr * (1 - dz) +
241-
p3[0] * dr * dz;
242-
243-
float Bz = p0[1] * (1 - dr) * (1 - dz) + p1[1] * (1 - dr) * dz + p2[1] * dr * (1 - dz) +
244-
p3[1] * dr * dz;
139+
fieldBrBz_t::view_t view(*fieldBrBz);
140+
fieldBrBz_t::output_t b = view.at(r, z);
141+
float Br = b[0];
142+
float Bz = b[1];
245143

246144
// convert Br Bz to Bx By Bz and rotate field
247145
const float phi = atan2(p.y(), p.x());
248146
auto B = fieldRot.has_value()
249147
? fieldRot.value() * ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz)
250148
: ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz);
149+
251150
field[0] += B.x();
252151
field[1] += B.y();
253152
field[2] += B.z();
153+
254154
} else { // BxByBz
255155

256-
if (!GetIndices(p.x(), p.y(), p.z(), &ix, &iy, &iz, &dx, &dy, &dz)) {
257-
return; // out of range
258-
}
259-
260-
float b[3] = {0};
261-
for (int comp = 0; comp < 3; comp++) { // field component loop
262-
// Trilinear interpolation
263-
// First along X, along 4 lines
264-
float b00 = Bvals_XYZ[ix][iy][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz][comp] * dx;
265-
float b01 =
266-
Bvals_XYZ[ix][iy][iz + 1][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz + 1][comp] * dx;
267-
float b10 =
268-
Bvals_XYZ[ix][iy + 1][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy + 1][iz][comp] * dx;
269-
float b11 = Bvals_XYZ[ix][iy + 1][iz + 1][comp] * (1 - dx) +
270-
Bvals_XYZ[ix + 1][iy + 1][iz + 1][comp] * dx;
271-
// Next along Y, along 2 lines
272-
float b0 = b00 * (1 - dy) + b10 * dy;
273-
float b1 = b01 * (1 - dy) + b11 * dy;
274-
// Finally along Z
275-
b[comp] = b0 * (1 - dz) + b1 * dz;
276-
}
156+
fieldBxByBz_t::view_t view(*fieldBxByBz);
157+
fieldBxByBz_t::output_t b = view.at(p.x(), p.y(), p.z());
277158

278159
// field rotation done in LoadMap()
279160
field[0] += b[0];
280161
field[1] += b[1];
281162
field[2] += b[2];
282163
}
283-
284-
return;
285164
}
286165

287166
// assign the field map to CartesianField
@@ -331,7 +210,6 @@ static Ref_t create_field_map_b(Detector& /*lcdd*/, xml::Handle_t handle) {
331210
}
332211

333212
auto map = new FieldMapB(field_type, coord_type);
334-
map->Configure(dimensions);
335213

336214
// translation
337215
if (x_dim.hasChild(_Unicode(rotationField))) {

0 commit comments

Comments
 (0)