MagneticFieldMap3D.cxx

#include "SimCore/MagneticFieldMap3D.h"
 
#include "Framework/Exception/Exception.h"
 
// STL
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
 
// Geant4
#include "G4SystemOfUnits.hh"
 
using namespace std;
 
namespace simcore {
MagneticFieldMap3D::MagneticFieldMap3D(const char* filename, double xOffset,
                                       double yOffset, double zOffset)
    : nx_(0),
      ny_(0),
      nz_(0),
      x_offset_(xOffset),
      y_offset_(yOffset),
      z_offset_(zOffset),
      invert_x_(false),
      invert_y_(false),
      invert_z_(false) {
  ifstream file(filename);  // Open the file for reading.
 
  // Throw an error if file does not exist.
  if (!file.good()) {
    EXCEPTION_RAISE("FileDNE", "The field map file '" + std::string(filename) +
                                   "' does not exist!");
  }
 
  ldmx_log(trace)
      << "-----------------------------------------------------------";
  ldmx_log(trace) << "    Magnetic Field Map 3D";
  ldmx_log(trace)
      << "-----------------------------------------------------------";
 
  ldmx_log(trace) << "Reading the field grid from " << filename << " ... ";
  ldmx_log(trace) << "  Offsets: " << xOffset << " " << yOffset << " "
                  << zOffset;
 
  // Ignore first blank line
  char buffer[256];
  file.getline(buffer, 256);
 
  // Read table dimensions
  file >> nx_ >> ny_ >> nz_;  // Note dodgy order
 
  ldmx_log(trace) << "  Number of values: " << nx_ << " " << ny_ << " " << nz_;
 
  // Set up storage space for table
  x_field_.resize(nx_);
  y_field_.resize(nx_);
  z_field_.resize(nx_);
  int ix, iy, iz;
  for (ix = 0; ix < nx_; ix++) {
    x_field_[ix].resize(ny_);
    y_field_[ix].resize(ny_);
    z_field_[ix].resize(ny_);
    for (iy = 0; iy < ny_; iy++) {
      x_field_[ix][iy].resize(nz_);
      y_field_[ix][iy].resize(nz_);
      z_field_[ix][iy].resize(nz_);
    }
  }
 
  // Ignore other header information
  // The first line whose second character is '0' is considered to
  // be the last line of the header.
  do {
    file.getline(buffer, 256);
  } while (buffer[1] != '0');
 
  // Read in the data
  double xval, yval, zval, bx, by, bz;
  for (ix = 0; ix < nx_; ix++) {
    for (iy = 0; iy < ny_; iy++) {
      for (iz = 0; iz < nz_; iz++) {
        file >> xval >> yval >> zval >> bx >> by >> bz;
        if (ix == 0 && iy == 0 && iz == 0) {
          minx_ = xval;
          miny_ = yval;
          minz_ = zval;
        }
        x_field_[ix][iy][iz] = bx;
        y_field_[ix][iy][iz] = by;
        z_field_[ix][iy][iz] = bz;
      }
    }
  }
  file.close();
 
  maxx_ = xval;
  maxy_ = yval;
  maxz_ = zval;
 
  ldmx_log(trace) << "  ... done reading ";
  ldmx_log(trace) << "Read values of field from file " << filename;
  ldmx_log(trace) << "  Assumed the order: x_, y_, z_, Bx, By, Bz";
  ldmx_log(trace) << "  Min values: " << minx_ << " " << miny_ << " " << minz_
                  << " mm ";
  ldmx_log(trace) << "  Max values: " << maxx_ << " " << maxy_ << " " << maxz_
                  << " mm ";
  ldmx_log(trace) << "  Field offsets: " << x_offset_ << " " << y_offset_ << " "
                  << z_offset_ << " mm ";
 
  // Should really check that the limits are not the wrong way around.
  if (maxx_ < minx_) {
    swap(maxx_, minx_);
    invert_x_ = true;
  }
  if (maxy_ < miny_) {
    swap(maxy_, miny_);
    invert_y_ = true;
  }
  if (maxz_ < minz_) {
    swap(maxz_, minz_);
    invert_z_ = true;
  }
 
  ldmx_log(trace) << "After reordering if necessary";
  ldmx_log(trace) << "  Min values: " << minx_ << " " << miny_ << " " << minz_
                  << " mm ";
  ldmx_log(trace) << "  Max values: " << maxx_ << " " << maxy_ << " " << maxz_
                  << " mm ";
  ;
 
  dx_ = maxx_ - minx_;
  dy_ = maxy_ - miny_;
  dz_ = maxz_ - minz_;
 
  ldmx_log(trace) << "  Range of values: " << dx_ << " " << dy_ << " " << dz_
                  << " mm";
  ldmx_log(trace) << "Done loading field map";
  ldmx_log(trace)
      << "-----------------------------------------------------------";
}
 
void MagneticFieldMap3D::GetFieldValue(const double point[4],
                                       double* bfield) const {
  double x = point[0] - x_offset_;
  double y = point[1] - y_offset_;
  double z = point[2] - z_offset_;
  double eps = 1E-6;
 
  // Check that the point is within the defined region
  if ((x >= minx_ && (x < maxx_ - eps)) && (y >= miny_ && (y < maxy_ - eps)) &&
      (z >= minz_ && (z < maxz_ - eps))) {
    // Position of given point within region, normalized to the range
    // [0,1]
    double xfraction = (x - minx_) / dx_;
    double yfraction = (y - miny_) / dy_;
    double zfraction = (z - minz_) / dz_;
 
    if (invert_x_) {
      xfraction = 1 - xfraction;
    }
    if (invert_y_) {
      yfraction = 1 - yfraction;
    }
    if (invert_z_) {
      zfraction = 1 - zfraction;
    }
 
    // Need addresses of these to pass to modf below.
    // modf uses its second argument as an OUTPUT argument.
    double xdindex, ydindex, zdindex;
 
    // Position of the point within the cuboid defined by the
    // nearest surrounding tabulated points
    double xlocal = (std::modf(xfraction * (nx_ - 1), &xdindex));
    double ylocal = (std::modf(yfraction * (ny_ - 1), &ydindex));
    double zlocal = (std::modf(zfraction * (nz_ - 1), &zdindex));
 
    // The indices of the nearest tabulated point whose coordinates
    // are all less than those of the given point
    int xindex = static_cast<int>(xdindex);
    int yindex = static_cast<int>(ydindex);
    int zindex = static_cast<int>(zdindex);
 
#ifdef DEBUG_INTERPOLATING_FIELD
    ldmx_log(trace) << "Local x_,y_,z_: " << xlocal << " " << ylocal << " "
                    << zlocal;
    ldmx_log(trace) << "Index x_,y_,z_: " << xindex << " " << yindex << " "
                    << zindex;
    double valx0z0, mulx0z0, valx1z0, mulx1z0;
    double valx0z1, mulx0z1, valx1z1, mulx1z1;
    valx0z0 = table[xindex][0][zindex];
    mulx0z0 = (1 - xlocal) * (1 - zlocal);
    valx1z0 = table[xindex + 1][0][zindex];
    mulx1z0 = xlocal * (1 - zlocal);
    valx0z1 = table[xindex][0][zindex + 1];
    mulx0z1 = (1 - xlocal) * zlocal;
    valx1z1 = table[xindex + 1][0][zindex + 1];
    mulx1z1 = xlocal * zlocal;
#endif
 
    // Full 3-dimensional version
    bfield[0] =
        x_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
            (1 - zlocal) +
        x_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
            zlocal +
        x_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
            (1 - zlocal) +
        x_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
            zlocal +
        x_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
            (1 - zlocal) +
        x_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
            zlocal +
        x_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
            (1 - zlocal) +
        x_field_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
    bfield[1] =
        y_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
            (1 - zlocal) +
        y_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
            zlocal +
        y_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
            (1 - zlocal) +
        y_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
            zlocal +
        y_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
            (1 - zlocal) +
        y_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
            zlocal +
        y_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
            (1 - zlocal) +
        y_field_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
    bfield[2] =
        z_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
            (1 - zlocal) +
        z_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
            zlocal +
        z_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
            (1 - zlocal) +
        z_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
            zlocal +
        z_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
            (1 - zlocal) +
        z_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
            zlocal +
        z_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
            (1 - zlocal) +
        z_field_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
 
  } else {
    bfield[0] = 0.0;
    bfield[1] = 0.0;
    bfield[2] = 0.0;
  }
}
 
}  // namespace simcore