From 6b53d42384ba2994b2d282fd6dc3cd22b6642b7e Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Wed, 9 Sep 2020 02:21:30 +0100
Subject: [PATCH 1/2] start moving methods to .icc

---
 .../MagFieldElements/BFieldMesh.h             |  65 ++--
 .../MagFieldElements/BFieldMesh.icc           | 336 +++++++++---------
 .../MagFieldElements/BFieldMeshZR.h           |  68 ++--
 .../MagFieldElements/BFieldMeshZR.icc         |  86 +++++
 .../MagFieldElements/BFieldZone.h             |  15 +-
 .../MagFieldElements/BFieldZone.icc           |  37 ++
 6 files changed, 356 insertions(+), 251 deletions(-)

diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
index b6dc1f34ffce..b3b166f30fb4 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 1abe3474d8d6..c88219a89c5d 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 b79841a62b3f..b82c6bd275b0 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 781dbf7dc885..1b76cf8b1dee 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 57fbb3ea6233..888e3f8d059a 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 071400e3729d..ddb5022730ae 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
 {
-- 
GitLab


From b25c36af47430679346bb5b19f7e63125b7c51a2 Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Wed, 9 Sep 2020 02:45:37 +0100
Subject: [PATCH 2/2] start moving methods to .icc

---
 .../MagFieldElements/AtlasFieldCache.h        |   1 -
 .../MagFieldElements/BFieldMesh.h             |  65 ++--
 .../MagFieldElements/BFieldMesh.icc           | 345 +++++++++---------
 .../MagFieldElements/BFieldMeshZR.h           |   4 +-
 .../MagFieldElements/src/AtlasFieldCache.cxx  |   5 +-
 .../MagFieldElements/src/AtlasFieldMap.cxx    |  16 -
 6 files changed, 213 insertions(+), 223 deletions(-)

diff --git a/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h
index 7161f00aa43e..5878c8b06909 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h
@@ -23,7 +23,6 @@
 #include "MagFieldElements/BFieldMeshZR.h"
 #include "MagFieldElements/BFieldZone.h"
 
-#include <iostream>
 #include <memory>
 namespace MagField {
 
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
index b3b166f30fb4..b6dc1f34ffce 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
@@ -39,24 +39,37 @@ public:
              double rmax,
              double phimin,
              double phimax,
-             double bscale);
+             double bscale)
+    : m_scale(bscale)
+    , m_nomScale(bscale)
+  {
+    m_min = { zmin, rmin, phimin };
+    m_max = { zmax, rmax, phimax };
+  }
   // set ranges
   void setRange(double zmin,
                 double zmax,
                 double rmin,
                 double rmax,
                 double phimin,
-                double phimax);
+                double phimax)
+  {
+    m_min = { zmin, rmin, phimin };
+    m_max = { zmax, rmax, phimax };
+  }
   // set bscale
-  void setBscale(double bscale);
+  void setBscale(double bscale) { m_scale = m_nomScale = bscale; }
   // scale bscale by a factor
-  void scaleBscale(double factor);
+  void scaleBscale(double factor) { m_scale = factor * m_nomScale; }
   // allocate space to vectors
   void reserve(int nz, int nr, int nphi, int nfield);
-  void reserve(int nz, int nr, int nphi);
+  void reserve(int nz, int nr, int nphi)
+  {
+    reserve(nz, nr, nphi, nz * nr * nphi);
+  }
   // add elements to vectors
-  void appendMesh(int i, double mesh);
-  void appendField(const BFieldVector<T>& field);
+  void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); }
+  void appendField(const BFieldVector<T>& field) { m_field.push_back(field); }
   // build Look Up Table
   void buildLUT();
   // test if a point is inside this zone
@@ -67,35 +80,41 @@ 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;
-  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;
+  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; }
   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 7b72c0939438..b43cfc406b05 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
@@ -3,21 +3,22 @@
 */
 
 #include "CxxUtils/vec.h"
+//
+// Reserve space in the vectors to avoid unnecessary memory re-allocations.
+//
 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)
+void
+BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield)
 {
-  m_min = { zmin, rmin, phimin };
-  m_max = { zmax, rmax, phimax };
+  m_mesh[0].reserve(nz);
+  m_mesh[1].reserve(nr);
+  m_mesh[2].reserve(nphi);
+  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,37 +33,9 @@ BFieldMesh<T>::inside(double z, double r, double phi) const
           r >= rmin() && r <= rmax());
 }
 
-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
-}
-
+//
+// Find and return the cache of the bin containing (z,r,phi)
+//
 template<class T>
 void
 BFieldMesh<T>::getCache(double z,
@@ -100,7 +73,6 @@ BFieldMesh<T>::getCache(double z,
   // store the bin edges
   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
   const double sf = scaleFactor;
@@ -145,153 +117,166 @@ BFieldMesh<T>::getCache(double z,
   cache.setBscale(m_scale);
 }
 
+//
+// Compute the magnetic field (and the derivatives) without caching
+//
 template<class T>
 void
-BFieldMesh<T>::setRange(double zmin,
-                        double zmax,
-                        double rmin,
-                        double rmax,
-                        double phimin,
-                        double phimax)
-{
-  m_min = { zmin, rmin, phimin };
-  m_max = { zmax, rmax, phimax };
-}
-
-template<class T>
-void
-BFieldMesh<T>::setBscale(double bscale)
-{
-  m_scale = m_nomScale = bscale;
-}
-
-template<class T>
-void
-BFieldMesh<T>::scaleBscale(double factor)
+BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz,
+                    double* ATH_RESTRICT B,
+                    double* ATH_RESTRICT deriv) const
 {
-  m_scale = factor * m_nomScale;
-}
+  // 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];
 
-template<class T>
-void
-BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield)
-{
-  m_mesh[0].reserve(nz);
-  m_mesh[1].reserve(nr);
-  m_mesh[2].reserve(nphi);
-  m_field.reserve(nfield);
-}
+  // 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>::reserve(int nz, int nr, int nphi)
-{
-  reserve(nz, nr, nphi, nz * nr * nphi);
-}
-
-template<class T>
-void
-BFieldMesh<T>::appendMesh(int i, double mesh)
-{
-  m_mesh[i].push_back(mesh);
+    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];
+  }
 }
 
+//
+// Construct the look-up table to accelerate bin-finding.
+//
 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
+BFieldMesh<T>::buildLUT()
 {
-  return m_scale;
+  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>
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
index b82c6bd275b0..1d78885d5312 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h
@@ -40,8 +40,8 @@ public:
                        BFieldCacheZR& cache,
                        double scaleFactor = 1.0) const;
   // accessors
-  double min(size_t i);
-  double max(size_t i);
+  double min(size_t i) const;
+  double max(size_t i) const;
   double zmin() const;
   double zmax() const;
   double rmin() const;
diff --git a/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx b/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx
index a3927c4553d7..d9e66122f27f 100644
--- a/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx
+++ b/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx
@@ -11,6 +11,9 @@
 //
 #include "MagFieldElements/AtlasFieldCache.h"
 
+#include <cmath>
+#include <iostream>
+
 void
 MagField::AtlasFieldCache::getField(const double* ATH_RESTRICT xyz,
                                     double* ATH_RESTRICT bxyz,
@@ -89,7 +92,7 @@ MagField::AtlasFieldCache::getFieldZR(const double* ATH_RESTRICT xyz,
   const double x = xyz[0];
   const double y = xyz[1];
   const double z = xyz[2];
-  const double r = sqrt(x * x + y * y);
+  const double r = std::sqrt(x * x + y * y);
 
   // test if the cache was initialized and the ZR cache is valid for current
   // position
diff --git a/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx b/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx
index 1956e5f35d9b..d1f928528cd0 100644
--- a/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx
+++ b/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx
@@ -89,46 +89,30 @@ MagField::AtlasFieldMap::initializeMap(TFile* rootfile,
   // }
   int id;
   double zmin;
-
   double zmax;
-
   double rmin;
-
   double rmax;
-
   double phimin;
-
   double phimax;
   double bscale;
   int ncond;
   bool* finite;
   double* p1x;
-
   double* p1y;
-
   double* p1z;
-
   double* p2x;
-
   double* p2y;
-
   double* p2z;
   double* curr;
   int nmeshz;
-
   int nmeshr;
-
   int nmeshphi;
   double* meshz;
-
   double* meshr;
-
   double* meshphi;
   int nfield;
   short* fieldz;
-
   short* fieldr;
-
   short* fieldphi;
   // define the fixed-sized branches first
   tree->SetBranchAddress("id", &id);
-- 
GitLab