epic/src/FieldMapB.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 - 2026 Wouter Deconinck, Dhevan Gangadharan, Aranya Giri
0003
0004 #include <DD4hep/DetFactoryHelper.h>
0005 #include <DD4hep/FieldTypes.h>
0006 #include <DD4hep/Printout.h>
0007 #include <XML/Utilities.h>
0008
0009 #include <cstdlib>
0010 #include <cmath>
0011 #include <filesystem>
0012 #include <fstream>
0013 #include <iostream>
0014 #include <optional>
0015 #include <sstream>
0016 #include <stdexcept>
0017 #include <string>
0018 #include <tuple>
0019 namespace fs = std::filesystem;
0020
0021 #include "FileLoaderHelper.h"
0022
0023 using namespace dd4hep;
0024
0025 //  Two coordinate options:
0026 //  coord_type="BrBz"
0027 //    expected fields contained within <dimensions>:
0028 //          <R step="" min="" max=""/>
0029 //          <Z step="" min="" max=""/>
0030 //  coord_type="BxByBz"
0031 //    expected fields contained within <dimensions>:
0032 //          <X step="" min="" max=""/>
0033 //          <Y step="" min="" max=""/>
0034 //          <Z step="" min="" max=""/>
0035
0036 // implementation of the field map
0037 class FieldMapB : public dd4hep::CartesianField::Object {
0038
0039   enum FieldCoord { BrBz, BxByBz };
0040
0041 public:
0042   FieldMapB(const std::string& field_type_str = "magnetic",
0043             const std::string& coord_type_str = "BrBz");
0044   void Configure(std::vector<xml_comp_t> dimensions);
0045   void LoadMap(const std::string& map_file, float scale);
0046   bool GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR, float* deltaZ) const;
0047   bool GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ, float* deltaX,
0048                   float* deltaY, float* deltaZ) const;
0049   void SetCoordTranslation(const Transform3D& tr) {
0050     coordTranslate     = tr;
0051     coordTranslate_inv = tr.Inverse();
0052   }
0053   void SetFieldRotation(const Transform3D& tr) {
0054     fieldRot     = tr;
0055     fieldRot_inv = tr.Inverse();
0056   }
0057
0058   virtual void fieldComponents(const double* pos, double* field);
0059
0060 private:
0061   FieldCoord fieldCoord;                                   // field coordinate type
0062   std::optional<Transform3D> coordTranslate{std::nullopt}; // coord translation
0063   std::optional<Transform3D> coordTranslate_inv{std::nullopt};
0064   std::optional<Transform3D> fieldRot{std::nullopt}; // field rotation
0065   std::optional<Transform3D> fieldRot_inv{std::nullopt};
0066   std::vector<float> steps, mins, maxs;                    // B map cell info
0067   std::vector<std::vector<std::array<float, 2>>> Bvals_RZ; // B map values:  {R}, {Z}, {Br,Bz}
0068   std::vector<std::vector<std::vector<std::array<float, 3>>>>
0069       Bvals_XYZ; // B map values: {X}, {Y}, {Z}, {Bx,By,Bz}
0070 };
0071
0072 // constructor
0073 FieldMapB::FieldMapB(const std::string& field_type_str, const std::string& coord_type_str) {
0074   std::string ftype = field_type_str;
0075   for (auto& c : ftype) {
0076     c = tolower(c);
0077   }
0078
0079   // set type
0080   if (ftype == "magnetic") {
0081     field_type = CartesianField::MAGNETIC;
0082   } else if (ftype == "electric") {
0083     field_type = CartesianField::ELECTRIC;
0084   } else {
0085     field_type = CartesianField::UNKNOWN;
0086     printout(ERROR, "FieldMap", "Unknown field type " + ftype);
0087   }
0088
0089   if (coord_type_str.compare("BrBz") == 0) { // BrBz
0090     fieldCoord = FieldCoord::BrBz;
0091   } else { // BxByBz
0092     fieldCoord = FieldCoord::BxByBz;
0093   }
0094 }
0095
0096 // fill field vector
0097 void FieldMapB::Configure(std::vector<xml_comp_t> dimensions) {
0098   // Fill vectors with step size, min, max for each dimension
0099   for (auto el : dimensions) {
0100     steps.push_back(el.step());
0101     mins.push_back(getAttrOrDefault<float>(el, _Unicode(min), 0));
0102     maxs.push_back(getAttrOrDefault<float>(el, _Unicode(max), 0));
0103   }
0104
0105   if (fieldCoord == FieldCoord::BrBz) {
0106     // N bins increased by 1 beyond grid size to account for edge cases at upper limits
0107     int nr = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
0108     int nz = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;
0109
0110     Bvals_RZ.resize(nr);
0111     for (auto& B2 : Bvals_RZ) {
0112       B2.resize(nz);
0113     }
0114   } else {
0115     // N bins increased by 1 beyond grid size to account for edge cases at upper limits
0116     int nx = std::roundf((maxs[0] - mins[0]) / steps[0]) + 2;
0117     int ny = std::roundf((maxs[1] - mins[1]) / steps[1]) + 2;
0118     int nz = std::roundf((maxs[2] - mins[2]) / steps[2]) + 2;
0119
0120     Bvals_XYZ.resize(nx);
0121     for (auto& B3 : Bvals_XYZ) {
0122       B3.resize(ny);
0123       for (auto& B2 : B3) {
0124         B2.resize(nz);
0125       }
0126     }
0127   }
0128 }
0129
0130 // get RZ cell indices corresponding to point of interest
0131 bool FieldMapB::GetIndices(float R, float Z, int* idxR, int* idxZ, float* deltaR,
0132                            float* deltaZ) const {
0133   // boundary check
0134   if (R > maxs[0] || R < mins[0] || Z > maxs[1] || Z < mins[1]) {
0135     return false;
0136   }
0137
0138   // get indices
0139   float idx_1_f, idx_2_f;
0140   *deltaR = std::modf((R - mins[0]) / steps[0], &idx_1_f);
0141   *deltaZ = std::modf((Z - mins[1]) / steps[1], &idx_2_f);
0142   *idxR   = static_cast<int>(idx_1_f);
0143   *idxZ   = static_cast<int>(idx_2_f);
0144
0145   return true;
0146 }
0147
0148 // get XYZ cell indices corresponding to point of interest
0149 bool FieldMapB::GetIndices(float X, float Y, float Z, int* idxX, int* idxY, int* idxZ,
0150                            float* deltaX, float* deltaY, float* deltaZ) const {
0151   // boundary check
0152   if (X > maxs[0] || X < mins[0] || Y > maxs[1] || Y < mins[1] || Z > maxs[2] || Z < mins[2]) {
0153     return false;
0154   }
0155
0156   // get indices
0157   float idx_1_f, idx_2_f, idx_3_f;
0158   *deltaX = std::modf((X - mins[0]) / steps[0], &idx_1_f);
0159   *deltaY = std::modf((Y - mins[1]) / steps[1], &idx_2_f);
0160   *deltaZ = std::modf((Z - mins[2]) / steps[2], &idx_3_f);
0161   *idxX   = static_cast<int>(idx_1_f);
0162   *idxY   = static_cast<int>(idx_2_f);
0163   *idxZ   = static_cast<int>(idx_3_f);
0164
0165   return true;
0166 }
0167
0168 // load data
0169 void FieldMapB::LoadMap(const std::string& map_file, float scale) {
0170   std::string line;
0171   std::ifstream input(map_file);
0172   if (!input) {
0173     printout(ERROR, "FieldMapB", "FieldMapB Error: file " + map_file + " cannot be read.");
0174   }
0175
0176   std::vector<float> coord = {};
0177   std::vector<float> Bcomp = {};
0178
0179   while (std::getline(input, line).good()) {
0180     std::istringstream iss(line);
0181
0182     coord.clear();
0183     Bcomp.clear();
0184
0185     if (fieldCoord == FieldCoord::BrBz) {
0186       coord.resize(2);
0187       Bcomp.resize(2);
0188       iss >> coord[0] >> coord[1] >> Bcomp[0] >> Bcomp[1];
0189
0190       int ir, iz;
0191       float dr, dz;
0192       if (!GetIndices(coord[0], coord[1], &ir, &iz, &dr, &dz)) {
0193         printout(WARNING, "FieldMapB", "coordinates out of range, skipped it.");
0194       } else { // scale field
0195         Bvals_RZ[ir][iz] = {Bcomp[0] * scale * float(tesla), Bcomp[1] * scale * float(tesla)};
0196       }
0197     } else {
0198       coord.resize(3);
0199       Bcomp.resize(3);
0200       iss >> coord[0] >> coord[1] >> coord[2] >> Bcomp[0] >> Bcomp[1] >> Bcomp[2];
0201
0202       int ix, iy, iz;
0203       float dx, dy, dz;
0204       if (!GetIndices(coord[0], coord[1], coord[2], &ix, &iy, &iz, &dx, &dy, &dz)) {
0205         printout(WARNING, "FieldMap", "coordinates out of range, skipped it.");
0206       } else { // scale and rotate B field vector
0207         auto B = ROOT::Math::XYZPoint(Bcomp[0], Bcomp[1], Bcomp[2]);
0208         B *= scale * float(tesla);
0209         if (fieldRot.has_value()) {
0210           B = fieldRot.value() * B;
0211         }
0212         Bvals_XYZ[ix][iy][iz] = {float(B.x()), float(B.y()), float(B.z())};
0213       }
0214     }
0215   }
0216 }
0217
0218 // get field components
0219 void FieldMapB::fieldComponents(const double* pos, double* field) {
0220   if (std::isnan(pos[0]) || std::isnan(pos[1]) || std::isnan(pos[2])) {
0221     printout(ERROR, "FieldMapB", "fieldComponents called for pos = (%f, %f, %f)", pos[0], pos[1],
0222              pos[2]);
0223     std::_Exit(EXIT_FAILURE);
0224   }
0225
0226   // coordinate conversion
0227   auto p = coordTranslate_inv.has_value()
0228                ? coordTranslate_inv.value() * ROOT::Math::XYZPoint(pos[0], pos[1], pos[2])
0229                : ROOT::Math::XYZPoint(pos[0], pos[1], pos[2]);
0230
0231   if (fieldCoord == FieldCoord::BrBz) {
0232     // coordinates conversion
0233     const float r = sqrt(p.x() * p.x() + p.y() * p.y());
0234     const float z = p.z();
0235
0236     int ir, iz;
0237     float dr, dz;
0238     if (!GetIndices(r, z, &ir, &iz, &dr, &dz)) {
0239       // out of range
0240       return;
0241     }
0242
0243     // p1    p3
0244     //    p
0245     // p0    p2
0246     auto& p0 = Bvals_RZ[ir][iz];
0247     auto& p1 = Bvals_RZ[ir][iz + 1];
0248     auto& p2 = Bvals_RZ[ir + 1][iz];
0249     auto& p3 = Bvals_RZ[ir + 1][iz + 1];
0250
0251     // Bilinear interpolation
0252     float Br = p0[0] * (1 - dr) * (1 - dz) + p1[0] * (1 - dr) * dz + p2[0] * dr * (1 - dz) +
0253                p3[0] * dr * dz;
0254
0255     float Bz = p0[1] * (1 - dr) * (1 - dz) + p1[1] * (1 - dr) * dz + p2[1] * dr * (1 - dz) +
0256                p3[1] * dr * dz;
0257
0258     // convert Br Bz to Bx By Bz and rotate field
0259     const float phi = atan2(p.y(), p.x());
0260     auto B = fieldRot.has_value()
0261                  ? fieldRot.value() * ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz)
0262                  : ROOT::Math::XYZPoint(Br * cos(phi), Br * sin(phi), Bz);
0263     field[0] += B.x();
0264     field[1] += B.y();
0265     field[2] += B.z();
0266   } else { // BxByBz
0267
0268     int ix, iy, iz;
0269     float dx, dy, dz;
0270     if (!GetIndices(p.x(), p.y(), p.z(), &ix, &iy, &iz, &dx, &dy, &dz)) {
0271       return; // out of range
0272     }
0273
0274     float b[3] = {0};
0275     for (int comp = 0; comp < 3; comp++) { // field component loop
0276       // Trilinear interpolation
0277       // First along X, along 4 lines
0278       float b00 = Bvals_XYZ[ix][iy][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz][comp] * dx;
0279       float b01 =
0280           Bvals_XYZ[ix][iy][iz + 1][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy][iz + 1][comp] * dx;
0281       float b10 =
0282           Bvals_XYZ[ix][iy + 1][iz][comp] * (1 - dx) + Bvals_XYZ[ix + 1][iy + 1][iz][comp] * dx;
0283       float b11 = Bvals_XYZ[ix][iy + 1][iz + 1][comp] * (1 - dx) +
0284                   Bvals_XYZ[ix + 1][iy + 1][iz + 1][comp] * dx;
0285       // Next along Y, along 2 lines
0286       float b0 = b00 * (1 - dy) + b10 * dy;
0287       float b1 = b01 * (1 - dy) + b11 * dy;
0288       // Finally along Z
0289       b[comp] = b0 * (1 - dz) + b1 * dz;
0290     }
0291
0292     // field rotation done in LoadMap()
0293     field[0] += b[0];
0294     field[1] += b[1];
0295     field[2] += b[2];
0296   }
0297
0298   return;
0299 }
0300
0301 // assign the field map to CartesianField
0302 static Ref_t create_field_map_b(Detector& /*lcdd*/, xml::Handle_t handle) {
0303   xml_comp_t x_par(handle);
0304
0305   if (!x_par.hasAttr(_Unicode(field_map))) {
0306     throw std::runtime_error(
0307         "FieldMapB Error: must have an xml attribute \"field_map\" for the field map.");
0308   }
0309
0310   CartesianField field;
0311   std::string field_type = x_par.attr<std::string>(_Unicode(field_type));
0312
0313   std::string coord_type = x_par.attr<std::string>(_Unicode(coord_type));
0314
0315   // dimensions
0316   xml_comp_t x_dim = x_par.dimensions();
0317
0318   // vector of dimension parameters: step, min, max
0319   std::vector<xml_comp_t> dimensions;
0320
0321   if (coord_type.compare("BrBz") == 0) {
0322     dimensions.push_back(x_dim.child(_Unicode(R)));
0323     dimensions.push_back(x_dim.child(_Unicode(Z)));
0324   } else if (coord_type.compare("BxByBz") == 0) {
0325     dimensions.push_back(x_dim.child(_Unicode(X)));
0326     dimensions.push_back(x_dim.child(_Unicode(Y)));
0327     dimensions.push_back(x_dim.child(_Unicode(Z)));
0328   } else {
0329     printout(ERROR, "FieldMapB", "Coordinate type: " + coord_type + ", is not BrBz nor BxByBz");
0330     std::_Exit(EXIT_FAILURE);
0331   }
0332
0333   std::string field_map_file  = x_par.attr<std::string>(_Unicode(field_map));
0334   std::string field_map_url   = x_par.attr<std::string>(_Unicode(url));
0335   std::string field_map_cache = getAttrOrDefault<std::string>(x_par, _Unicode(cache), "");
0336
0337   EnsureFileFromURLExists(field_map_url, field_map_file, field_map_cache);
0338
0339   float field_map_scale = x_par.attr<float>(_Unicode(scale));
0340
0341   if (!fs::exists(fs::path(field_map_file))) {
0342     printout(ERROR, "FieldMapB", "file " + field_map_file + " does not exist");
0343     printout(ERROR, "FieldMapB", "use a FileLoader plugin before the field element");
0344     std::_Exit(EXIT_FAILURE);
0345   }
0346
0347   auto map = new FieldMapB(field_type, coord_type);
0348   map->Configure(dimensions);
0349
0350   // translation
0351   if (x_dim.hasChild(_Unicode(rotationField))) {
0352     static float deg2r = ROOT::Math::Pi() / 180.;
0353     xml_comp_t rot_dim = x_dim.child(_Unicode(rotationField));
0354     RotationZYX rot(rot_dim.z() * deg2r, rot_dim.y() * deg2r, rot_dim.x() * deg2r);
0355     map->SetFieldRotation(Transform3D(rot));
0356   }
0357   // rotation
0358   if (x_dim.hasChild(_Unicode(translationCoord))) {
0359     xml_comp_t trans_dim = x_dim.child(_Unicode(translationCoord));
0360     Translation3D trans(trans_dim.x(), trans_dim.y(), trans_dim.z());
0361     map->SetCoordTranslation(Transform3D(trans));
0362   }
0363
0364   map->LoadMap(field_map_file, field_map_scale);
0365   field.assign(map, x_par.nameStr(), "FieldMapB");
0366
0367   return field;
0368 }
0369
0370 DECLARE_XMLELEMENT(epic_FieldMapB, create_field_map_b)