Skip to content

feat: convert FieldMapB to using covfie #934

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
+39 −160
Open

feat: convert FieldMapB to using covfie #934

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ include(GNUInstallDirs)
# Dependencies
find_package(DD4hep 1.27 REQUIRED COMPONENTS DDCore DDRec)
find_package(fmt REQUIRED)
find_package(covfie REQUIRED)

#-----------------------------------------------------------------------------------
set(a_lib_name ${PROJECT_NAME})
Expand Down
198 changes: 38 additions & 160 deletions src/FieldMapB.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,14 @@
#include <DD4hep/FieldTypes.h>
#include <DD4hep/Printout.h>
#include <XML/Utilities.h>
#include <covfie/core/algebra/affine.hpp>
#include <covfie/core/backend/primitive/array.hpp>
#include <covfie/core/backend/transformer/affine.hpp>
#include <covfie/core/backend/transformer/linear.hpp>
#include <covfie/core/backend/transformer/strided.hpp>
#include <covfie/core/field.hpp>
#include <covfie/core/field_view.hpp>
#include <covfie/core/parameter_pack.hpp>

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

#include "FileLoaderHelper.h"

using fieldBrBz_t =
covfie::field<covfie::backend::affine<covfie::backend::linear<covfie::backend::strided<
covfie::vector::size2, covfie::backend::array<covfie::vector::float2>>>>>;

using fieldBxByBz_t =
covfie::field<covfie::backend::affine<covfie::backend::linear<covfie::backend::strided<
covfie::vector::size3, covfie::backend::array<covfie::vector::float3>>>>>;

using namespace dd4hep;

// Two coordinate options:
Expand All @@ -40,11 +56,9 @@ class FieldMapB : public dd4hep::CartesianField::Object {
public:
FieldMapB(const std::string& field_type_str = "magnetic",
const std::string& coord_type_str = "BrBz");
void Configure(std::vector<xml_comp_t> dimensions);

void LoadMap(const std::string& map_file, float scale);
bool GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR, float* deltaZ);
bool GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ, float* deltaX,
float* deltaY, float* deltaZ);

void SetCoordTranslation(const Transform3D& tr) {
coordTranslate = tr;
coordTranslate_inv = tr.Inverse();
Expand All @@ -62,12 +76,10 @@ class FieldMapB : public dd4hep::CartesianField::Object {
std::optional<Transform3D> coordTranslate_inv{std::nullopt};
std::optional<Transform3D> fieldRot{std::nullopt}; // field rotation
std::optional<Transform3D> fieldRot_inv{std::nullopt};
std::vector<float> steps, mins, maxs; // B map cell info
int ir, ix, iy, iz; // lookup indices
float dr, dx, dy, dz; // deltas for interpolation
std::vector<std::vector<std::array<float, 2>>> Bvals_RZ; // B map values: {R}, {Z}, {Br,Bz}
std::vector<std::vector<std::vector<std::array<float, 3>>>>
Bvals_XYZ; // B map values: {X}, {Y}, {Z}, {Bx,By,Bz}
std::vector<float> steps, mins, maxs; // B map cell info

std::shared_ptr<fieldBrBz_t> fieldBrBz;
std::shared_ptr<fieldBxByBz_t> fieldBxByBz;
};

// constructor
Expand All @@ -94,120 +106,21 @@ FieldMapB::FieldMapB(const std::string& field_type_str, const std::string& coord
}
}

// fill field vector
void FieldMapB::Configure(std::vector<xml_comp_t> dimensions) {
// Fill vectors with step size, min, max for each dimension
for (auto el : dimensions) {
steps.push_back(el.step());
mins.push_back(getAttrOrDefault<float>(el, _Unicode(min), 0));
maxs.push_back(getAttrOrDefault<float>(el, _Unicode(max), 0));
}

if (fieldCoord == FieldCoord::BrBz) {
// N bins increased by 1 beyond grid size to account for edge cases at upper limits
int nr = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
int nz = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;

Bvals_RZ.resize(nr);
for (auto& B2 : Bvals_RZ) {
B2.resize(nz);
}
} else {
// N bins increased by 1 beyond grid size to account for edge cases at upper limits
int nx = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
int ny = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;
int nz = std::roundf((maxs[2] - mins[2]) / steps[2]) + 2;

Bvals_XYZ.resize(nx);
for (auto& B3 : Bvals_XYZ) {
B3.resize(ny);
for (auto& B2 : B3) {
B2.resize(nz);
}
}
}
}

// get RZ cell indices corresponding to point of interest
bool FieldMapB::GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR, float* deltaZ) {
// boundary check
if (R > maxs[0] || R < mins[0] || Z > maxs[1] || Z < mins[1]) {
return false;
}

// get indices
float idx_1_f, idx_2_f;
*deltaR = std::modf((R - mins[0]) / steps[0], &idx_1_f);
*deltaZ = std::modf((Z - mins[1]) / steps[1], &idx_2_f);
*idxR = static_cast<int>(idx_1_f);
*idxZ = static_cast<int>(idx_2_f);

return true;
}

// get XYZ cell indices corresponding to point of interest
bool FieldMapB::GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ,
float* deltaX, float* deltaY, float* deltaZ) {
// boundary check
if (X > maxs[0] || X < mins[0] || Y > maxs[1] || Y < mins[1] || Z > maxs[2] || Z < mins[2]) {
return false;
}

// get indices
float idx_1_f, idx_2_f, idx_3_f;
*deltaX = std::modf((X - mins[0]) / steps[0], &idx_1_f);
*deltaY = std::modf((Y - mins[1]) / steps[1], &idx_2_f);
*deltaZ = std::modf((Z - mins[2]) / steps[2], &idx_3_f);
*idxX = static_cast<int>(idx_1_f);
*idxY = static_cast<int>(idx_2_f);
*idxZ = static_cast<int>(idx_3_f);

return true;
}

// load data
void FieldMapB::LoadMap(const std::string& map_file, float scale) {
void FieldMapB::LoadMap(const std::string& map_file, float) {
std::string line;
std::ifstream input(map_file);
if (!input) {
printout(ERROR, "FieldMapB", "FieldMapB Error: file " + map_file + " cannot be read.");
}

std::vector<float> coord = {};
std::vector<float> Bcomp = {};

while (std::getline(input, line).good()) {
std::istringstream iss(line);

coord.clear();
Bcomp.clear();

if (fieldCoord == FieldCoord::BrBz) {
coord.resize(2);
Bcomp.resize(2);
iss >> coord[0] >> coord[1] >> Bcomp[0] >> Bcomp[1];

if (!GetIndices(coord[0], coord[1], &ir, &iz, &dr, &dz)) {
printout(WARNING, "FieldMapB", "coordinates out of range, skipped it.");
} else { // scale field
Bvals_RZ[ir][iz] = {Bcomp[0] * scale * float(tesla), Bcomp[1] * scale * float(tesla)};
}
} else {
coord.resize(3);
Bcomp.resize(3);
iss >> coord[0] >> coord[1] >> coord[2] >> Bcomp[0] >> Bcomp[1] >> Bcomp[2];

if (!GetIndices(coord[0], coord[1], coord[2], &ix, &iy, &iz, &dx, &dy, &dz)) {
printout(WARNING, "FieldMap", "coordinates out of range, skipped it.");
} else { // scale and rotate B field vector
auto B = ROOT::Math::XYZPoint(Bcomp[0], Bcomp[1], Bcomp[2]);
B *= scale * float(tesla);
if (fieldRot.has_value()) {
B = fieldRot.value() * B;
}
Bvals_XYZ[ix][iy][iz] = {float(B.x()), float(B.y()), float(B.z())};
}
}
switch (fieldCoord) {
case FieldCoord::BrBz:
fieldBrBz = std::make_shared<fieldBrBz_t>(input);
break;
case FieldCoord::BxByBz:
fieldBxByBz = std::make_shared<fieldBxByBz_t>(input);
break;
}
}

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

if (!GetIndices(r, z, &ir, &iz, &dr, &dz)) {
// out of range
return;
}

// p1 p3
// p
// p0 p2
auto& p0 = Bvals_RZ[ir][iz];
auto& p1 = Bvals_RZ[ir][iz + 1];
auto& p2 = Bvals_RZ[ir + 1][iz];
auto& p3 = Bvals_RZ[ir + 1][iz + 1];

// Bilinear interpolation
float Br = p0[0] * (1 - dr) * (1 - dz) + p1[0] * (1 - dr) * dz + p2[0] * dr * (1 - dz) +
p3[0] * dr * dz;

float Bz = p0[1] * (1 - dr) * (1 - dz) + p1[1] * (1 - dr) * dz + p2[1] * dr * (1 - dz) +
p3[1] * dr * dz;
fieldBrBz_t::view_t view(*fieldBrBz);
fieldBrBz_t::output_t b = view.at(r, z);
float Br = b[0];
float Bz = b[1];

// convert Br Bz to Bx By Bz and rotate field
const float phi = atan2(p.y(), p.x());
auto B = fieldRot.has_value()
? fieldRot.value() * ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz)
: ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz);

field[0] += B.x();
field[1] += B.y();
field[2] += B.z();

} else { // BxByBz

if (!GetIndices(p.x(), p.y(), p.z(), &ix, &iy, &iz, &dx, &dy, &dz)) {
return; // out of range
}

float b[3] = {0};
for (int comp = 0; comp < 3; comp++) { // field component loop
// Trilinear interpolation
// First along X, along 4 lines
float b00 = Bvals_XYZ[ix][iy][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz][comp] * dx;
float b01 =
Bvals_XYZ[ix][iy][iz + 1][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz + 1][comp] * dx;
float b10 =
Bvals_XYZ[ix][iy + 1][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy + 1][iz][comp] * dx;
float b11 = Bvals_XYZ[ix][iy + 1][iz + 1][comp] * (1 - dx) +
Bvals_XYZ[ix + 1][iy + 1][iz + 1][comp] * dx;
// Next along Y, along 2 lines
float b0 = b00 * (1 - dy) + b10 * dy;
float b1 = b01 * (1 - dy) + b11 * dy;
// Finally along Z
b[comp] = b0 * (1 - dz) + b1 * dz;
}
fieldBxByBz_t::view_t view(*fieldBxByBz);
fieldBxByBz_t::output_t b = view.at(p.x(), p.y(), p.z());

// field rotation done in LoadMap()
field[0] += b[0];
field[1] += b[1];
field[2] += b[2];
}

return;
}

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

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

// translation
if (x_dim.hasChild(_Unicode(rotationField))) {
Expand Down
Loading