athena/src/FieldMapBrBz.cpp

0001 #include <DD4hep/DetFactoryHelper.h>
0002 #include <DD4hep/FieldTypes.h>
0003 #include <DD4hep/Printout.h>
0004 #include <XML/Utilities.h>
0005
0006 #include <cstdlib>
0007 #include <filesystem>
0008 #include <fstream>
0009 #include <iostream>
0010 #include <sstream>
0011 #include <stdexcept>
0012 #include <string>
0013 #include <tuple>
0014 namespace fs = std::filesystem;
0015
0016 #include "FileLoaderHelper.h"
0017
0018 using namespace dd4hep;
0019
0020
0021 // implementation of the field map
0022 class FieldMapBrBz : public dd4hep::CartesianField::Object
0023 {
0024 public:
0025     FieldMapBrBz(const std::string &field_type = "magnetic");
0026     void Configure(double rmin, double rmax, double rstep, double zmin, double zmax, double zstep);
0027     void LoadMap(const std::string &map_file, double scale);
0028     void GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz);
0029     void SetTransform(const Transform3D &tr) { trans = tr; trans_inv = tr.Inverse(); }
0030
0031     virtual void fieldComponents(const double *pos, double *field);
0032
0033
0034 private:
0035     Transform3D trans, trans_inv;
0036     double rmin, rmax, rstep, zmin, zmax, zstep;
0037     std::vector<std::vector<std::vector<double>>> Bvals;
0038 };
0039
0040 // constructor
0041 FieldMapBrBz::FieldMapBrBz(const std::string &field_type)
0042 {
0043     std::string ftype = field_type;
0044     for (auto &c : ftype) { c = tolower(c); }
0045
0046     // set type
0047     if (ftype == "magnetic") {
0048         type = CartesianField::MAGNETIC;
0049     } else if (ftype == "electric") {
0050         type = CartesianField::ELECTRIC;
0051     } else {
0052         type = CartesianField::UNKNOWN;
0053         std::cout << "FieldMapBrBz Warning: Unknown field type " << field_type << "!" << std::endl;
0054     }
0055 }
0056
0057 void FieldMapBrBz::Configure(double r1, double r2, double rs, double z1, double z2, double zs)
0058 {
0059     rmin = r1;
0060     rmax = r2;
0061     rstep = rs;
0062     zmin = z1;
0063     zmax = z2;
0064     zstep = zs;
0065
0066     int nr = int((r2 - r1)/rs) + 2;
0067     int nz = int((z2 - z1)/zs) + 2;
0068
0069     Bvals.resize(nr);
0070     for (auto &B2 : Bvals) {
0071         B2.resize(nz);
0072         for (auto &B : B2) {
0073             B.resize(2, 0.);
0074         }
0075     }
0076 }
0077
0078 void FieldMapBrBz::GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz)
0079 {
0080     // boundary check
0081     if (r > rmax || r < rmin || z > zmax || z < zmin) {
0082         ir = -1;
0083         iz = -1;
0084         return;
0085     }
0086     // get indices
0087     double idr, idz;
0088     dr = std::modf((r - rmin)/rstep, &idr);
0089     dz = std::modf((z - zmin)/zstep, &idz);
0090
0091     ir = static_cast<int>(idr);
0092     iz = static_cast<int>(idz);
0093 }
0094
0095 // load data
0096 void FieldMapBrBz::LoadMap(const std::string &map_file, double scale)
0097 {
0098     std::string line;
0099     std::ifstream input(map_file);
0100     if (!input) {
0101         std::cout << "FieldMapBrBz Error: file \"" << map_file << "\" cannot be read." << std::endl;
0102     }
0103
0104     double r, z, br, bz;
0105     int ir, iz;
0106     double dr, dz;
0107     while (std::getline(input, line).good()) {
0108         std::istringstream iss(line);
0109         iss >> r >> z >> br >> bz;
0110         GetIndices(r, z, ir, iz, dr, dz);
0111         if (ir < 0 || iz < 0) {
0112             std::cout << "FieldMapBrBz Warning: coordinates out of range ("
0113                       << r << ", " << z << "), skipped it." << std::endl;
0114         } else {
0115             Bvals[ir][iz] = {br*scale, bz*scale};
0116             // ROOT::Math::XYZPoint p(r, 0, z);
0117             // std::cout << p << " -> " << trans*p << std::endl;
0118             // std::cout << ir << ", " << iz << ", " << br << ", " << bz << std::endl;
0119         }
0120     }
0121 }
0122
0123 // get field components
0124 void FieldMapBrBz::fieldComponents(const double *pos, double *field)
0125 {
0126     // coordinate conversion
0127     auto p = trans_inv*ROOT::Math::XYZPoint(pos[0], pos[1], pos[2]);
0128
0129     // coordinates conversion
0130     const double r = sqrt(p.x()*p.x() + p.y()*p.y());
0131     const double z = p.z();
0132     const double phi = atan2(p.y(), p.x());
0133
0134     int ir, iz;
0135     double dr, dz;
0136     GetIndices(r, z, ir, iz, dr, dz);
0137
0138     // out of the range
0139     if (ir < 0 || iz < 0) { return; }
0140
0141     // p1    p3
0142     //    p
0143     // p0    p2
0144     auto &p0 = Bvals[ir][iz];
0145     auto &p1 = Bvals[ir][iz + 1];
0146     auto &p2 = Bvals[ir + 1][iz];
0147     auto &p3 = Bvals[ir + 1][iz + 1];
0148
0149     // linear interpolation
0150     double Br = p0[0] * (1-dr) * (1-dz)
0151               + p1[0] * (1-dr) *    dz
0152               + p2[0] *    dr  * (1-dz)
0153               + p3[0] *    dr  *    dz;
0154
0155     double Bz = p0[1] * (1-dr) * (1-dz)
0156               + p1[1] * (1-dr) *    dz
0157               + p2[1] *    dr  * (1-dz)
0158               + p3[1] *    dr  *    dz;
0159
0160     // convert Br Bz to Bx By Bz
0161     auto B = trans*ROOT::Math::XYZPoint(Br*sin(phi), Br*cos(phi), Bz);
0162     field[0] += B.x()*tesla;
0163     field[1] += B.y()*tesla;
0164     field[2] += B.z()*tesla;
0165     return;
0166 }
0167
0168
0169 // assign the field map to CartesianField
0170 static Ref_t create_field_map_brbz(Detector & /*lcdd*/, xml::Handle_t handle)
0171 {
0172     xml_comp_t x_par(handle);
0173
0174     if (!x_par.hasAttr(_Unicode(field_map))) {
0175         throw std::runtime_error("FieldMapBrBz Error: must have an xml attribute \"field_map\" for the field map.");
0176     }
0177
0178     CartesianField field;
0179     std::string field_type = x_par.attr<std::string>(_Unicode(field_type));
0180
0181     // dimensions
0182     xml_comp_t x_dim = x_par.dimensions();
0183
0184     // min, max, step
0185     xml_comp_t r_dim = x_dim.child(_Unicode(transverse));
0186     xml_comp_t z_dim = x_dim.child(_Unicode(longitudinal));
0187
0188     std::string field_map_file = x_par.attr<std::string>(_Unicode(field_map));
0189     std::string field_map_url = x_par.attr<std::string>(_Unicode(url));
0190
0191     EnsureFileFromURLExists(field_map_url,field_map_file);
0192
0193     double field_map_scale = x_par.attr<double>(_Unicode(scale));
0194
0195     if( !fs::exists(fs::path(field_map_file))  ) {
0196         printout(ERROR, "FieldMapBrBz", "file " + field_map_file + " does not exist");
0197         printout(ERROR, "FieldMapBrBz", "use a FileLoader plugin before the field element");
0198         std::_Exit(EXIT_FAILURE);
0199     }
0200
0201     auto map = new FieldMapBrBz(field_type);
0202     map->Configure(r_dim.rmin(), r_dim.rmax(), r_dim.step(), z_dim.zmin(), z_dim.zmax(), z_dim.step());
0203
0204     // translation, rotation
0205     static double deg2r = ROOT::Math::Pi()/180.;
0206     RotationZYX rot(0., 0., 0.);
0207     if (x_dim.hasChild(_Unicode(rotation))) {
0208         xml_comp_t rot_dim = x_dim.child(_Unicode(rotation));
0209         rot = RotationZYX(rot_dim.z()*deg2r, rot_dim.y()*deg2r, rot_dim.x()*deg2r);
0210     }
0211
0212     Translation3D trans(0., 0., 0.);
0213     if (x_dim.hasChild(_Unicode(translation))) {
0214         xml_comp_t trans_dim = x_dim.child(_Unicode(translation));
0215         trans = Translation3D(trans_dim.x(), trans_dim.y(), trans_dim.z());
0216     }
0217     map->SetTransform(trans*rot);
0218
0219     map->LoadMap(field_map_file, field_map_scale);
0220     field.assign(map, x_par.nameStr(), "FieldMapBrBz");
0221
0222     return field;
0223 }
0224
0225 DECLARE_XMLELEMENT(FieldMapBrBz, create_field_map_brbz)
0226