diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
index b6dc1f34ffce2e01150893a173971af2c1929dbc..b3b166f30fb499a94e43b3803716fc76ea0922b7 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
@@ -39,37 +39,24 @@ public:
              double rmax,
              double phimin,
              double phimax,
-             double bscale)
-    : m_scale(bscale)
-    , m_nomScale(bscale)
-  {
-    m_min = { zmin, rmin, phimin };
-    m_max = { zmax, rmax, phimax };
-  }
+             double bscale);
   // set ranges
   void setRange(double zmin,
                 double zmax,
                 double rmin,
                 double rmax,
                 double phimin,
-                double phimax)
-  {
-    m_min = { zmin, rmin, phimin };
-    m_max = { zmax, rmax, phimax };
-  }
+                double phimax);
   // set bscale
-  void setBscale(double bscale) { m_scale = m_nomScale = bscale; }
+  void setBscale(double bscale);
   // scale bscale by a factor
-  void scaleBscale(double factor) { m_scale = factor * m_nomScale; }
+  void scaleBscale(double factor);
   // allocate space to vectors
   void reserve(int nz, int nr, int nphi, int nfield);
-  void reserve(int nz, int nr, int nphi)
-  {
-    reserve(nz, nr, nphi, nz * nr * nphi);
-  }
+  void reserve(int nz, int nr, int nphi);
   // add elements to vectors
-  void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); }
-  void appendField(const BFieldVector<T>& field) { m_field.push_back(field); }
+  void appendMesh(int i, double mesh);
+  void appendField(const BFieldVector<T>& field);
   // build Look Up Table
   void buildLUT();
   // test if a point is inside this zone
@@ -80,41 +67,35 @@ public:
                 double phi,
                 BFieldCache& cache,
                 double scaleFactor = 1.0) const;
-  // get the B field
-  void getB(const double* ATH_RESTRICT xyz,
-            double* ATH_RESTRICT B,
-            double* ATH_RESTRICT deriv = nullptr) const;
   // accessors
-  double min(size_t i) const { return m_min[i]; }
-  double max(size_t i) const { return m_max[i]; }
-  double zmin() const { return m_min[0]; }
-  double zmax() const { return m_max[0]; }
-  double rmin() const { return m_min[1]; }
-  double rmax() const { return m_max[1]; }
-  double phimin() const { return m_min[2]; }
-  double phimax() const { return m_max[2]; }
-  unsigned nmesh(size_t i) const { return m_mesh[i].size(); }
-  double mesh(size_t i, size_t j) const { return m_mesh[i][j]; }
-  unsigned nfield() const { return m_field.size(); }
-  const BFieldVector<T>& field(size_t i) const { return m_field[i]; }
-  double bscale() const { return m_scale; }
+  double min(size_t i) const;
+  double max(size_t i) const;
+  double zmin() const;
+  double zmax() const;
+  double rmin() const;
+  double rmax() const;
+  double phimin() const;
+  double phimax() const;
+  unsigned nmesh(size_t i) const;
+  double mesh(size_t i, size_t j) const;
+  unsigned nfield() const;
+  const BFieldVector<T>& field(size_t i) const;
+  double bscale() const;
   int memSize() const;
 
 protected:
   std::array<double, 3> m_min;
   std::array<double, 3> m_max;
-  std::array<std::vector<double>,3> m_mesh;
+  std::array<std::vector<double>, 3> m_mesh;
 
 private:
   std::vector<BFieldVector<T>> m_field;
   double m_scale = 1.0;
   double m_nomScale; // nominal m_scale from the map
-
   // look-up table and related variables
-  std::array<std::vector<int>,3> m_LUT;
-  std::array<double,3> m_invUnit; // inverse unit size in the LUT
+  std::array<std::vector<int>, 3> m_LUT;
+  std::array<double, 3> m_invUnit; // inverse unit size in the LUT
   int m_roff, m_zoff;
-
 };
 #include "MagFieldElements/BFieldMesh.icc"
 #endif
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
index 1abe3474d8d6176af368125686d235dd71f3c59f..c88219a89c5d457b4538453b3edb0ea9cf4b4e32 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
@@ -2,9 +2,22 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-//
-// Reserve space in the vectors to avoid unnecessary memory re-allocations.
-//
+// constructor
+template<class T>
+BFieldMesh<T>::BFieldMesh(double zmin,
+                          double zmax,
+                          double rmin,
+                          double rmax,
+                          double phimin,
+                          double phimax,
+                          double bscale)
+  : m_scale(bscale)
+  , m_nomScale(bscale)
+{
+  m_min = { zmin, rmin, phimin };
+  m_max = { zmax, rmax, phimax };
+}
+
 template<class T>
 void
 BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield)
@@ -15,9 +28,6 @@ BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield)
   m_field.reserve(nfield);
 }
 
-//
-// Test if a point (z,r,phi) is inside this mesh region.
-//
 template<class T>
 bool
 BFieldMesh<T>::inside(double z, double r, double phi) const
@@ -32,9 +42,37 @@ BFieldMesh<T>::inside(double z, double r, double phi) const
           r >= rmin() && r <= rmax());
 }
 
-//
-// Find and return the cache of the bin containing (z,r,phi)
-//
+template<class T>
+void
+BFieldMesh<T>::buildLUT()
+{
+  for (int j = 0; j < 3; ++j) { // z, r, phi
+    // align the m_mesh edges to m_min/m_max
+    m_mesh[j].front() = m_min[j];
+    m_mesh[j].back() = m_max[j];
+    // determine the unit size, q, to be used in the LUTs
+    const double width = m_mesh[j].back() - m_mesh[j].front();
+    double q(width);
+    for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) {
+      q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]);
+    }
+    // find the number of units in the LUT
+    int n = int(width / q) + 1;
+    q = width / (n + 0.5);
+    m_invUnit[j] = 1.0 / q; // new unit size
+    ++n;
+    int m = 0;                    // mesh number
+    for (int i = 0; i < n; ++i) { // LUT index
+      if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) {
+        m++;
+      }
+      m_LUT[j].push_back(m);
+    }
+  }
+  m_roff = m_mesh[2].size();          // index offset for incrementing r by 1
+  m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1
+}
+
 template<class T>
 void
 BFieldMesh<T>::getCache(double z,
@@ -70,8 +108,8 @@ BFieldMesh<T>::getCache(double z,
     ++iphi;
   }
   // store the bin edges
-  cache.setRange(mz[iz], mz[iz + 1], mr[ir], mr[ir + 1], mphi[iphi], mphi[iphi + 1]);
-
+  cache.setRange(
+    mz[iz], mz[iz + 1], mr[ir], mr[ir + 1], mphi[iphi], mphi[iphi + 1]);
 
   // store the B field at the 8 corners in cache
   const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner
@@ -105,164 +143,143 @@ BFieldMesh<T>::getCache(double z,
   cache.setBscale(m_scale);
 }
 
-//
-// Compute the magnetic field (and the derivatives) without caching
-//
 template<class T>
 void
-BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz,
-                    double* ATH_RESTRICT B,
-                    double* ATH_RESTRICT deriv) const
+BFieldMesh<T>::setRange(double zmin,
+                        double zmax,
+                        double rmin,
+                        double rmax,
+                        double phimin,
+                        double phimax)
 {
-  // cylindrical coordinates
-  const double x = xyz[0];
-  const double y = xyz[1];
-  const double z = xyz[2];
-  const double r = sqrt(x * x + y * y);
-  double phi = std::atan2(y, x);
-  if (phi < phimin()) {
-    phi += 2.0 * M_PI;
-  }
-  // is it inside this map?
-  if (!inside(z, r, phi)) { // no
-    B[0] = B[1] = B[2] = 0.0;
-    if (deriv) {
-      for (int i = 0; i < 9; ++i) {
-        deriv[i] = 0.0;
-      }
-    }
-    return;
-  }
-  // find the bin
-  // z
-  const std::vector<double>& mz(m_mesh[0]);
-  int iz = int((z - zmin()) * m_invUnit[0]); // index to LUT
-  iz = m_LUT[0][iz];                         // tentative mesh index from LUT
-  if (z > mz[iz + 1]) {
-    ++iz;
-  }
-  // r
-  const std::vector<double>& mr(m_mesh[1]);
-  int ir = int((r - rmin()) * m_invUnit[1]); // index to LUT
-  ir = m_LUT[1][ir];                         // tentative mesh index from LUT
-  if (r > mr[ir + 1]) {
-    ++ir;
-  }
-  // phi
-  const std::vector<double>& mphi(m_mesh[2]);
-  int iphi = int((phi - phimin()) * m_invUnit[2]); // index to LUT
-  iphi = m_LUT[2][iphi]; // tentative mesh index from LUT
-  if (phi > mphi[iphi + 1]) {
-    ++iphi;
-  }
-  // get the B field at the 8 corners
-  int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner
-  const std::array<BFieldVector<T>, 8> field = { m_field[im0],
-                                           m_field[im0 + 1],
-                                           m_field[im0 + m_roff],
-                                           m_field[im0 + m_roff + 1],
-                                           m_field[im0 + m_zoff],
-                                           m_field[im0 + m_zoff + 1],
-                                           m_field[im0 + m_zoff + m_roff],
-                                           m_field[im0 + m_zoff + m_roff + 1] };
-  // fractional position inside this mesh
-  const double fz = (z - mz[iz]) / (mz[iz + 1] - mz[iz]);
-  const double gz = 1.0 - fz;
-  const double fr = (r - mr[ir]) / (mr[ir + 1] - mr[ir]);
-  const double gr = 1.0 - fr;
-  const double fphi = (phi - mphi[iphi]) / (mphi[iphi + 1] - mphi[iphi]);
-  const double gphi = 1.0 - fphi;
-  const double scale = m_scale;
-  // interpolate field values in z, r, phi
-  double Bzrphi[3];
-  for (int i = 0; i < 3; ++i) { // z, r, phi
-    Bzrphi[i] = scale * (gz * (gr * (gphi * field[0][i] + fphi * field[1][i]) +
-                               fr * (gphi * field[2][i] + fphi * field[3][i])) +
-                         fz * (gr * (gphi * field[4][i] + fphi * field[5][i]) +
-                               fr * (gphi * field[6][i] + fphi * field[7][i])));
-  }
-  // convert (Bz,Br,Bphi) to (Bx,By,Bz)
-  const double c = (r > 0.0) ? x / r : cos(mphi[iphi]);
-  const double s = (r > 0.0) ? y / r : sin(mphi[iphi]);
-  B[0] = Bzrphi[1] * c - Bzrphi[2] * s;
-  B[1] = Bzrphi[1] * s + Bzrphi[2] * c;
-  B[2] = Bzrphi[0];
+  m_min = { zmin, rmin, phimin };
+  m_max = { zmax, rmax, phimax };
+}
 
-  // compute field derivatives if requested
-  if (deriv) {
-    const double sz = m_scale / (mz[iz + 1] - mz[iz]);
-    const double sr = m_scale / (mr[ir + 1] - mr[ir]);
-    const double sphi = m_scale / (mphi[iphi + 1] - mphi[iphi]);
+template<class T>
+void
+BFieldMesh<T>::setBscale(double bscale)
+{
+  m_scale = m_nomScale = bscale;
+}
 
-    std::array<double,3> dBdz;
-    std::array<double,3> dBdr;
-    std::array<double,3> dBdphi;
-    for (int j = 0; j < 3; ++j) { // Bz, Br, Bphi components
-      dBdz[j] = sz * (gr * (gphi * (field[4][j] - field[0][j]) +
-                            fphi * (field[5][j] - field[1][j])) +
-                      fr * (gphi * (field[6][j] - field[2][j]) +
-                            fphi * (field[7][j] - field[3][j])));
-      dBdr[j] = sr * (gz * (gphi * (field[2][j] - field[0][j]) +
-                            fphi * (field[3][j] - field[1][j])) +
-                      fz * (gphi * (field[6][j] - field[4][j]) +
-                            fphi * (field[7][j] - field[5][j])));
-      dBdphi[j] = sphi * (gz * (gr * (field[1][j] - field[0][j]) +
-                                fr * (field[3][j] - field[2][j])) +
-                          fz * (gr * (field[5][j] - field[4][j]) +
-                                fr * (field[7][j] - field[6][j])));
-    }
-    // convert to cartesian coordinates
-    const double cc = c * c;
-    const double cs = c * s;
-    const double ss = s * s;
-    deriv[0] = cc * dBdr[1] - cs * dBdr[2] - cs * dBdphi[1] / r +
-               ss * dBdphi[2] / r + s * B[1] / r;
-    deriv[1] = cs * dBdr[1] - ss * dBdr[2] + cc * dBdphi[1] / r -
-               cs * dBdphi[2] / r - c * B[1] / r;
-    deriv[2] = c * dBdz[1] - s * dBdz[2];
-    deriv[3] = cs * dBdr[1] + cc * dBdr[2] - ss * dBdphi[1] / r -
-               cs * dBdphi[2] / r - s * B[0] / r;
-    deriv[4] = ss * dBdr[1] + cs * dBdr[2] + cs * dBdphi[1] / r +
-               cc * dBdphi[2] / r + c * B[0] / r;
-    deriv[5] = s * dBdz[1] + c * dBdz[2];
-    deriv[6] = c * dBdr[0] - s * dBdphi[0] / r;
-    deriv[7] = s * dBdr[1] + c * dBdphi[0] / r;
-    deriv[8] = dBdz[0];
-  }
+template<class T>
+void
+BFieldMesh<T>::scaleBscale(double factor)
+{
+  m_scale = factor * m_nomScale;
 }
 
-//
-// Construct the look-up table to accelerate bin-finding.
-//
 template<class T>
 void
-BFieldMesh<T>::buildLUT()
+BFieldMesh<T>::reserve(int nz, int nr, int nphi)
 {
-  for (int j = 0; j < 3; ++j) { // z, r, phi
-    // align the m_mesh edges to m_min/m_max
-    m_mesh[j].front() = m_min[j];
-    m_mesh[j].back() = m_max[j];
-    // determine the unit size, q, to be used in the LUTs
-    const double width = m_mesh[j].back() - m_mesh[j].front();
-    double q(width);
-    for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) {
-      q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]);
-    }
-    // find the number of units in the LUT
-    int n = int(width / q) + 1;
-    q = width / (n + 0.5);
-    m_invUnit[j] = 1.0 / q; // new unit size
-    ++n;
-    int m = 0;                    // mesh number
-    for (int i = 0; i < n; ++i) { // LUT index
-      if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) {
-        m++;
-      }
-      m_LUT[j].push_back(m);
-    }
-  }
-  m_roff = m_mesh[2].size();          // index offset for incrementing r by 1
-  m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1
+  reserve(nz, nr, nphi, nz * nr * nphi);
+}
+
+template<class T>
+void
+BFieldMesh<T>::appendMesh(int i, double mesh)
+{
+  m_mesh[i].push_back(mesh);
+}
+
+template<class T>
+void
+BFieldMesh<T>::appendField(const BFieldVector<T>& field)
+{
+  m_field.push_back(field);
+}
+
+template<class T>
+double
+BFieldMesh<T>::min(size_t i) const
+{
+  return m_min[i];
+}
+
+template<class T>
+double
+BFieldMesh<T>::max(size_t i) const
+{
+  return m_max[i];
+}
+
+template<class T>
+double
+BFieldMesh<T>::zmin() const
+{
+  return m_min[0];
+}
+
+template<class T>
+double
+BFieldMesh<T>::zmax() const
+{
+  return m_max[0];
+}
+
+template<class T>
+double
+BFieldMesh<T>::rmin() const
+{
+  return m_min[1];
+}
+
+template<class T>
+double
+BFieldMesh<T>::rmax() const
+{
+  return m_max[1];
+}
+
+template<class T>
+double
+BFieldMesh<T>::phimin() const
+{
+  return m_min[2];
+}
+
+template<class T>
+double
+BFieldMesh<T>::phimax() const
+{
+  return m_max[2];
+}
+
+template<class T>
+unsigned
+BFieldMesh<T>::nmesh(size_t i) const
+{
+  return m_mesh[i].size();
+}
+
+template<class T>
+double
+BFieldMesh<T>::mesh(size_t i, size_t j) const
+{
+  return m_mesh[i][j];
+}
+
+template<class T>
+unsigned
+BFieldMesh<T>::nfield() const
+{
+  return m_field.size();
+}
+
+template<class T>
+const BFieldVector<T>&
+BFieldMesh<T>::field(size_t i) const
+{
+  return m_field[i];
+}
+
+template<class T>
+double
+BFieldMesh<T>::bscale() const
+{
+  return m_scale;
 }
 
 template<class T>
@@ -280,4 +297,3 @@ BFieldMesh<T>::memSize() const
   return size;
 }
 
-
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
index b79841a62b3f3f510298d69a0adad48ab901afb4..b82c6bd275b0a17eee7170b980f7f721553a6faf 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
@@ -2,73 +2,63 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-//
-// BFieldMeshZR.h
-//
-// A 2-dim z-r mesh inside the solenoid field map
-//
-// Masahiro Morii, Harvard University
-//
+/*
+ * BFieldMeshZR.h
+ *
+ * A 2-dim z-r mesh inside the solenoid field map
+ *
+ * Masahiro Morii, Harvard University
+ *
+ * AthenaMT : RD Schaffer , Christos Anastopoulos
+ */
 #ifndef BFIELDMESHZR_H
 #define BFIELDMESHZR_H
 
 #include "MagFieldElements/BFieldCacheZR.h"
 #include "MagFieldElements/BFieldVectorZR.h"
+#include <array>
 #include <cmath>
 #include <vector>
-#include <array>
 
 class BFieldMeshZR
 {
 public:
   // constructor
-  BFieldMeshZR(double zmin, double zmax, double rmin, double rmax)
-  {
-    m_min = { zmin, rmin };
-    m_max = { zmax, rmax };
-  }
+  BFieldMeshZR(double zmin, double zmax, double rmin, double rmax);
   // allocate space to vectors
-  void reserve(int nz, int nr)
-  {
-    m_mesh[0].reserve(nz);
-    m_mesh[1].reserve(nr);
-    m_field.reserve(nz * nr);
-  }
+  void reserve(int nz, int nr);
   // add elements to vectors
-  void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); }
-  void appendField(const BFieldVectorZR& field) { m_field.push_back(field); }
+  void appendMesh(int i, double mesh);
+  void appendField(const BFieldVectorZR& field);
   // build LUT
   void buildLUT();
   // test if a point is inside this zone
-  inline bool inside(double z, double r) const
-  {
-    return (z >= zmin() && z <= zmax() && r >= rmin() && r <= rmax());
-  }
+  bool inside(double z, double r) const;
   // find the bin
   inline void getCache(double z,
                        double r,
                        BFieldCacheZR& cache,
                        double scaleFactor = 1.0) const;
   // accessors
-  double min(size_t i) const { return m_min[i]; }
-  double max(size_t i) const { return m_max[i]; }
-  double zmin() const { return m_min[0]; }
-  double zmax() const { return m_max[0]; }
-  double rmin() const { return m_min[1]; }
-  double rmax() const { return m_max[1]; }
-  unsigned nmesh(size_t i) const { return m_mesh[i].size(); }
-  double mesh(size_t i, size_t j) const { return m_mesh[i][j]; }
-  unsigned nfield() const { return m_field.size(); }
-  const BFieldVectorZR& field(size_t i) const { return m_field[i]; }
+  double min(size_t i);
+  double max(size_t i);
+  double zmin() const;
+  double zmax() const;
+  double rmin() const;
+  double rmax() const;
+  unsigned nmesh(size_t i) const;
+  double mesh(size_t i, size_t j) const;
+  unsigned nfield() const;
+  const BFieldVectorZR& field(size_t i) const;
   int memSize() const;
 
 private:
-  std::array<double,2> m_min;
-  std::array<double,2> m_max;
-  std::array<std::vector<double>,2> m_mesh;
+  std::array<double, 2> m_min;
+  std::array<double, 2> m_max;
+  std::array<std::vector<double>, 2> m_mesh;
   std::vector<BFieldVectorZR> m_field;
   // look-up table and related variables
-  std::array<std::vector<int>,2> m_LUT;
+  std::array<std::vector<int>, 2> m_LUT;
   std::array<double, 2> m_invUnit; // inverse unit size in the LUT
   int m_zoff;
 };
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc
index 781dbf7dc885eef916b3f8c87d55c4da9b8d454e..1b76cf8b1dee79300de6ca089fd736c41efe4b60 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc
@@ -36,3 +36,89 @@ BFieldMeshZR::getCache(double z,
   cache.setField(3, m_field[im0 + m_zoff + 1], scaleFactor);
 }
 
+inline BFieldMeshZR::BFieldMeshZR(double zmin,
+                                  double zmax,
+                                  double rmin,
+                                  double rmax)
+{
+  m_min = { zmin, rmin };
+  m_max = { zmax, rmax };
+}
+
+inline void
+BFieldMeshZR::reserve(int nz, int nr)
+{
+  m_mesh[0].reserve(nz);
+  m_mesh[1].reserve(nr);
+  m_field.reserve(nz * nr);
+}
+// add elements to vectors
+inline void
+BFieldMeshZR::appendMesh(int i, double mesh)
+{
+  m_mesh[i].push_back(mesh);
+}
+inline void
+BFieldMeshZR::appendField(const BFieldVectorZR& field)
+{
+  m_field.push_back(field);
+}
+
+// test if a point is inside this zone
+inline bool
+BFieldMeshZR::inside(double z, double r) const
+{
+  return (z >= zmin() && z <= zmax() && r >= rmin() && r <= rmax());
+}
+
+inline double
+BFieldMeshZR::min(size_t i) const
+{
+  return m_min[i];
+}
+inline double
+BFieldMeshZR::max(size_t i) const
+{
+  return m_max[i];
+}
+inline double
+BFieldMeshZR::zmin() const
+{
+  return m_min[0];
+}
+inline double
+BFieldMeshZR::zmax() const
+{
+  return m_max[0];
+}
+inline double
+BFieldMeshZR::rmin() const
+{
+  return m_min[1];
+}
+inline double
+BFieldMeshZR::rmax() const
+{
+  return m_max[1];
+}
+inline unsigned
+BFieldMeshZR::nmesh(size_t i) const
+{
+  return m_mesh[i].size();
+}
+inline double
+BFieldMeshZR::mesh(size_t i, size_t j) const
+{
+  return m_mesh[i][j];
+}
+inline unsigned
+BFieldMeshZR::nfield() const
+{
+  return m_field.size();
+}
+inline const BFieldVectorZR&
+BFieldMeshZR::field(size_t i) const
+{
+  return m_field[i];
+}
+
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h
index 57fbb3ea62338aefee086f46ac501fa3e69da6ee..888e3f8d059a56fdba94c138a94b79fced504544 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h
@@ -29,12 +29,7 @@ public:
              double rmax,
              double phimin,
              double phimax,
-             double scale)
-    : BFieldMesh<short>(zmin, zmax, rmin, rmax, phimin, phimax, scale)
-    , m_id(id)
-  {
-    ;
-  }
+             double scale);
   // add elements to vectors
   void appendCond(const BFieldCond& cond);
   // compute Biot-Savart magnetic field and add to B[3]
@@ -45,10 +40,10 @@ public:
   // Scaling is done in cachec
   void scaleField(double factor);
   // accessors
-  int id() const { return m_id; }
-  unsigned ncond() const { return m_cond.size(); }
-  const BFieldCond& cond(int i) const { return m_cond[i]; }
-  const std::vector<BFieldCond>* condVector() const { return &m_cond; }
+  int id() const;
+  unsigned ncond() const;
+  const BFieldCond& cond(int i) const;
+  const std::vector<BFieldCond>* condVector() const;
   int memSize() const;
   // adjust the min/max edges to a new value
   void adjustMin(int i, double x);
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc
index 071400e3729dfef1ee7ac92c3d500d1653f4fa56..ddb5022730ae8913d9cdcbe186f528d46e80aa14 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc
@@ -1,6 +1,19 @@
 /*
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
+
+inline BFieldZone::BFieldZone(int id,
+                              double zmin,
+                              double zmax,
+                              double rmin,
+                              double rmax,
+                              double phimin,
+                              double phimax,
+                              double scale)
+  : BFieldMesh<short>(zmin, zmax, rmin, rmax, phimin, phimax, scale)
+  , m_id(id)
+{}
+
 inline void
 BFieldZone::appendCond(const BFieldCond& cond)
 {
@@ -26,6 +39,30 @@ BFieldZone::addBiotSavart(const double* ATH_RESTRICT xyz,
   }
 }
 
+inline int
+BFieldZone::id() const
+{
+  return m_id;
+}
+
+inline unsigned
+BFieldZone::ncond() const
+{
+  return m_cond.size();
+}
+
+inline const BFieldCond&
+BFieldZone::cond(int i) const
+{
+  return m_cond[i];
+}
+
+inline const std::vector<BFieldCond>*
+BFieldZone::condVector() const
+{
+  return &m_cond;
+}
+
 inline int
 BFieldZone::memSize() const
 {