diff --git a/CMakeLists.txt b/CMakeLists.txt index 49e5778d3f..8395556a57 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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}) diff --git a/src/FieldMapB.cpp b/src/FieldMapB.cpp index b58fc1eafe..ec6ee70d19 100644 --- a/src/FieldMapB.cpp +++ b/src/FieldMapB.cpp @@ -5,6 +5,14 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include #include #include @@ -19,6 +27,14 @@ namespace fs = std::filesystem; #include "FileLoaderHelper.h" +using fieldBrBz_t = + covfie::field>>>>; + +using fieldBxByBz_t = + covfie::field>>>>; + using namespace dd4hep; // Two coordinate options: @@ -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 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(); @@ -62,12 +76,10 @@ class FieldMapB : public dd4hep::CartesianField::Object { std::optional coordTranslate_inv{std::nullopt}; std::optional fieldRot{std::nullopt}; // field rotation std::optional fieldRot_inv{std::nullopt}; - std::vector steps, mins, maxs; // B map cell info - int ir, ix, iy, iz; // lookup indices - float dr, dx, dy, dz; // deltas for interpolation - std::vector>> Bvals_RZ; // B map values: {R}, {Z}, {Br,Bz} - std::vector>>> - Bvals_XYZ; // B map values: {X}, {Y}, {Z}, {Bx,By,Bz} + std::vector steps, mins, maxs; // B map cell info + + std::shared_ptr fieldBrBz; + std::shared_ptr fieldBxByBz; }; // constructor @@ -94,120 +106,21 @@ FieldMapB::FieldMapB(const std::string& field_type_str, const std::string& coord } } -// fill field vector -void FieldMapB::Configure(std::vector dimensions) { - // Fill vectors with step size, min, max for each dimension - for (auto el : dimensions) { - steps.push_back(el.step()); - mins.push_back(getAttrOrDefault(el, _Unicode(min), 0)); - maxs.push_back(getAttrOrDefault(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(idx_1_f); - *idxZ = static_cast(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(idx_1_f); - *idxY = static_cast(idx_2_f); - *idxZ = static_cast(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 coord = {}; - std::vector 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(input); + break; + case FieldCoord::BxByBz: + fieldBxByBz = std::make_shared(input); + break; } } @@ -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 @@ -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))) {