From 748f8b6d93d3c1eed2396ca468768c48e4449fba Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Fri, 21 Aug 2020 01:04:20 +0100
Subject: [PATCH 1/2] Magnetic Field Elements, small test fixes, add reference
 for the test, move more methods to .icc or .cxx

---
 MagneticField/MagFieldElements/CMakeLists.txt |   7 +-
 .../MagFieldElements/BFieldCache.h            | 125 ++----
 .../MagFieldElements/BFieldCache.icc          |  77 +++-
 .../share/BFieldExample_test.ref              |  22 +
 .../MagFieldElements/src/BFieldCache.cxx      |  23 +
 .../test/BFieldExample_test.cxx               | 402 ++++++++++--------
 6 files changed, 368 insertions(+), 288 deletions(-)
 create mode 100644 MagneticField/MagFieldElements/share/BFieldExample_test.ref
 create mode 100644 MagneticField/MagFieldElements/src/BFieldCache.cxx

diff --git a/MagneticField/MagFieldElements/CMakeLists.txt b/MagneticField/MagFieldElements/CMakeLists.txt
index f7929d3a6cda..81af6c3cfef9 100644
--- a/MagneticField/MagFieldElements/CMakeLists.txt
+++ b/MagneticField/MagFieldElements/CMakeLists.txt
@@ -15,8 +15,7 @@ atlas_add_library( MagFieldElements
                    PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} PathResolver )
 
 		 
-atlas_add_test( BFieldExample
+atlas_add_test( BFieldExample_test
                 SOURCES  test/BFieldExample_test.cxx
-                INCLUDE_DIRS ${GTEST_INCLUDE_DIRS} #${GMOCK_INCLUDE_DIRS}
-                LINK_LIBRARIES TestTools  CxxUtils EventPrimitives ${GTEST_LIBRARIES} #${GMOCK_LIBRARIES}
-                ENVIRONMENT "JOBOPTSEARCHPATH=${CMAKE_CURRENT_SOURCE_DIR}/share" )		
+                INCLUDE_DIRS 
+                LINK_LIBRARIES MagFieldElements)		
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
index 8282eb91153e..d4895554493a 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
@@ -2,15 +2,18 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-//
-// BFieldCache.h
-//
-// Cashe of one bin of the magnetic field map.
-// Defined by ranges in z, r, phi, and the B vectors at the 8 corners of the
-// "bin".
-//
-// Masahiro Morii, Harvard University
-//
+/**
+ * BFieldCache.h
+ *
+ * Cache of one bin of the magnetic field map.
+ * Defined by ranges in z, r, phi, and the B vectors at the 8 corners of the
+ * "bin".
+ *
+ * Masahiro Morii, Harvard University
+ *
+ * RD Schaffer , Christos Anastopoulos
+ */
+
 #ifndef BFIELDCACHE_H
 #define BFIELDCACHE_H
 
@@ -22,101 +25,26 @@ class BFieldCache
 {
 public:
   // default constructor sets unphysical boundaries, so that inside() will fail
-  BFieldCache()
-    : m_zmin(0.0)
-    , m_zmax(0.0)
-    , m_rmin(0.0)
-    , m_rmax(0.0)
-    , m_phimin(0.0)
-    , m_phimax(-1.0)
-  {
-  }
+  BFieldCache() = default;
   // make this cache invalid, so that inside() will fail
-  void invalidate()
-  {
-    m_phimin = 0.0;
-    m_phimax = -1.0;
-  }
+  void invalidate();
   // set the z, r, phi range that defines the bin
   void setRange(double zmin,
                 double zmax,
                 double rmin,
                 double rmax,
                 double phimin,
-                double phimax)
-  {
-    m_zmin = zmin;
-    m_zmax = zmax;
-    m_rmin = rmin;
-    m_rmax = rmax;
-    m_phimin = phimin;
-    m_phimax = phimax;
-    m_invz = 1.0 / (zmax - zmin);
-    m_invr = 1.0 / (rmax - rmin);
-    m_invphi = 1.0 / (phimax - phimin);
-  }
-
-    // **
-    // ** The following two methods are commented out, and replaced by the third setField as optimization test
-    // **
-  // // set the field values at each corner (rescale for current scale factor)
-  // void setField(const std::array<BFieldVector<double>, 8>& field,
-  //               double scaleFactor = 1.0)
-  // {
-  //   // We pass array of 8 elements with 3 entries
-  //   // Which go to 3x8 matrix
-  //   for (int i = 0; i < 8; ++i) {
-  //     for (int j = 0; j < 3; ++j) {
-  //       m_field[j][i] = scaleFactor * field[i][j];
-  //     }
-  //   }
-  // }
-
-  // void setField(const std::array<BFieldVector<short>, 8>& field,
-  //               double scaleFactor = 1.0)
-  // {
-  //   // We pass array of 8 elements with 3 entries
-  //   // Which go to 3x8 matrix
-  //   for (int i = 0; i < 8; ++i) {
-  //     for (int j = 0; j < 3; ++j) {
-  //       m_field[j][i] = scaleFactor * field[i][j];
-  //     }
-  //   }
-  // }
+                double phimax);
 
   // set field array, filled externally
-  void setField(double field[][8]) {
-      std::copy( &field[0][0], &field[0][0] + 3*8, &m_field[0][0]);
-  }
-    
-
-  // set the field values at each corner (rescale for current scale factor)
-  void printField() 
-  {
-    // print out field values
-    std::cout << "Field at corner i, for each component j (Bz, Br, Bphi)" << std::endl;
-
-    for (int i = 0; i < 8; ++i) {
-      for (int j = 0; j < 3; ++j) {
-          std::cout << i << "," << j << ": " << m_field[j][i] << ", ";
-      }
-      std::cout << std::endl;
-    }
-  }
+  void setField(double field[][8]);
 
   // set the multiplicative factor for the field vectors
-  void setBscale(double bscale) { m_scale = bscale; }
-  float bscale() { return m_scale; }
+  void setBscale(double bscale);
+  float bscale() const;
 
   // test if (z, r, phi) is inside this bin
-  bool inside(double z, double r, double phi) const
-  {
-    if (phi < m_phimin) {
-      phi += 2.0 * M_PI;
-    }
-    return (phi >= m_phimin && phi <= m_phimax && z >= m_zmin && z <= m_zmax &&
-            r >= m_rmin && r <= m_rmax);
-  }
+  bool inside(double z, double r, double phi) const;
   // interpolate the field and return B[3].
   // also compute field derivatives if deriv[9] is given.
   void getB(const double* ATH_RESTRICT xyz,
@@ -125,10 +53,19 @@ public:
             double* ATH_RESTRICT B,
             double* ATH_RESTRICT deriv = nullptr) const;
 
+  // set the field values at each corner (rescale for current scale factor)
+  void printField() const;
+
 private:
-  double m_zmin, m_zmax;           // bin range in z
-  double m_rmin, m_rmax;           // bin range in r
-  double m_phimin, m_phimax;       // bin range in phi
+  // bin range in z
+  double m_zmin = 0.0;
+  double m_zmax = 0.0;
+  // bin range in r
+  double m_rmin = 0.0;
+  double m_rmax = 0.0;
+  // bin range in phi
+  double m_phimin = 0.0;
+  double m_phimax = -1.0;          // bin range in phi
   double m_invz, m_invr, m_invphi; // 1/(bin size) in z, r, phi
   double m_field[3][8];            // (Bz,Br,Bphi) at 8 corners of the bin
   double m_scale;                  // unit of m_field in kT
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
index 9a4e1df61592..0b81313acc86 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
@@ -2,6 +2,60 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
+inline void
+BFieldCache::invalidate()
+{
+  m_phimin = 0.0;
+  m_phimax = -1.0;
+}
+inline void
+BFieldCache::setRange(double zmin,
+                      double zmax,
+                      double rmin,
+                      double rmax,
+                      double phimin,
+                      double phimax)
+{
+  m_zmin = zmin;
+  m_zmax = zmax;
+  m_rmin = rmin;
+  m_rmax = rmax;
+  m_phimin = phimin;
+  m_phimax = phimax;
+  m_invz = 1.0 / (zmax - zmin);
+  m_invr = 1.0 / (rmax - rmin);
+  m_invphi = 1.0 / (phimax - phimin);
+}
+
+// set field array, filled externally
+inline void
+BFieldCache::setField(double field[][8])
+{
+  std::copy(&field[0][0], &field[0][0] + 3 * 8, &m_field[0][0]);
+}
+
+inline void
+BFieldCache::setBscale(double bscale)
+{
+  m_scale = bscale;
+}
+
+inline float
+BFieldCache::bscale() const
+{
+  return m_scale;
+}
+
+inline bool
+BFieldCache::inside(double z, double r, double phi) const
+{
+  if (phi < m_phimin) {
+    phi += 2.0 * M_PI;
+  }
+  return (phi >= m_phimin && phi <= m_phimax && z >= m_zmin && z <= m_zmax &&
+          r >= m_rmin && r <= m_rmax);
+}
+
 inline void
 BFieldCache::getB(const double* ATH_RESTRICT xyz,
                   double r,
@@ -54,8 +108,8 @@ BFieldCache::getB(const double* ATH_RESTRICT xyz,
 
   // compute field derivatives if requested
   if (deriv) {
-    const double sz   = m_scale * m_invz;
-    const double sr   = m_scale * m_invr;
+    const double sz = m_scale * m_invz;
+    const double sr = m_scale * m_invr;
     const double sphi = m_scale * m_invphi;
 
     std::array<double, 3> dBdz;
@@ -66,21 +120,15 @@ BFieldCache::getB(const double* ATH_RESTRICT xyz,
       const double* field = m_field[j];
       dBdz[j] =
         sz *
-        (gr * (gphi * (field[4] - field[0]) +
-               fphi * (field[5] - field[1])) +
-         fr * (gphi * (field[6] - field[2]) +
-               fphi * (field[7] - field[3])));
+        (gr * (gphi * (field[4] - field[0]) + fphi * (field[5] - field[1])) +
+         fr * (gphi * (field[6] - field[2]) + fphi * (field[7] - field[3])));
       dBdr[j] =
         sr *
-        (gz * (gphi * (field[2] - field[0]) +
-               fphi * (field[3] - field[1])) +
-         fz * (gphi * (field[6] - field[4]) +
-               fphi * (field[7] - field[5])));
+        (gz * (gphi * (field[2] - field[0]) + fphi * (field[3] - field[1])) +
+         fz * (gphi * (field[6] - field[4]) + fphi * (field[7] - field[5])));
       dBdphi[j] =
-        sphi * (gz * (gr * (field[1] - field[0]) +
-                      fr * (field[3] - field[2])) +
-                fz * (gr * (field[5] - field[4]) +
-                      fr * (field[7] - field[6])));
+        sphi * (gz * (gr * (field[1] - field[0]) + fr * (field[3] - field[2])) +
+                fz * (gr * (field[5] - field[4]) + fr * (field[7] - field[6])));
     }
     // convert to cartesian coordinates
     const double cc = c * c;
@@ -107,4 +155,3 @@ BFieldCache::getB(const double* ATH_RESTRICT xyz,
   }
 }
 
-
diff --git a/MagneticField/MagFieldElements/share/BFieldExample_test.ref b/MagneticField/MagFieldElements/share/BFieldExample_test.ref
new file mode 100644
index 000000000000..bd0c98f185b3
--- /dev/null
+++ b/MagneticField/MagFieldElements/share/BFieldExample_test.ref
@@ -0,0 +1,22 @@
+start BFieldExample test
+j, meshz 0, -1400
+j, meshz 1, -466.93
+j, meshz 2, 466.14
+j, meshz 3, 1400
+did  meshz 
+did  meshr 
+did  meshphi 
+did  field 
+did  buildLUT 
+get field std: i, bxyz 0 -2.83727e-07, 9.47007e-08, 0.00308551 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 1 -2.81403e-07, 7.49033e-08, 0.00255923 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 2 -2.79079e-07, 5.51058e-08, 0.00203296 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 3 -2.76755e-07, 3.53084e-08, 0.00150669 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 4 -2.74431e-07, 1.5511e-08, 0.000980422 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 5 -2.72107e-07, -4.28645e-09, 0.000454151 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 6 -2.69782e-07, -2.40839e-08, -7.21201e-05 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 7 -2.67458e-07, -4.38813e-08, -0.000598391 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 8 -2.65134e-07, -6.36787e-08, -0.00112466 fractional diff gt 10^-5: 0, 0, 0
+get field std: i, bxyz 9 -2.6281e-07, -8.34762e-08, -0.00165093 fractional diff gt 10^-5: 0, 0, 0
+runTest: status 0
+Test passed OK
diff --git a/MagneticField/MagFieldElements/src/BFieldCache.cxx b/MagneticField/MagFieldElements/src/BFieldCache.cxx
new file mode 100644
index 000000000000..979b01e53f7c
--- /dev/null
+++ b/MagneticField/MagFieldElements/src/BFieldCache.cxx
@@ -0,0 +1,23 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MagFieldElement/BFieldCache.h"
+#include <iostream>
+
+// set the field values at each corner (rescale for current scale factor)
+void
+BFieldCache::printField()
+{
+  // print out field values
+  std::cout << "Field at corner i, for each component j (Bz, Br, Bphi)"
+            << '\n';
+
+  for (int i = 0; i < 8; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      std::cout << i << "," << j << ": " << m_field[j][i] << ", ";
+    }
+    std::cout << '\n';
+  }
+}
+
diff --git a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
index 34ca83c9e083..d2b812ac28c2 100644
--- a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
+++ b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
@@ -2,204 +2,256 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-#include <iostream>
 #include "MagFieldElements/BFieldCache.h"
 #include "MagFieldElements/BFieldZone.h"
+#include <iostream>
 #include <unistd.h>
 
 /**
  * @brief Specimen of service caching expensive calculations.
  * It uses internally std::map which is not thread-safely expandable.
  **/
-class BFieldTest {
+class BFieldTest
+{
 public:
-    int runTest(bool doDebug = false  ) {
+  int runTest(bool doDebug = false)
+  {
+
+    // Create a field zone
+
+    double fieldz[] = {
+      19487, 19487, 19488, 19488, 19487, 19487, 19531, 19531, 19532, 19532,
+      19531, 19531, 6399,  6400,  6400,  6400,  6399,  -1561, -1561, -1560,
+      -1560, -1560, -1561, -1516, -1516, -1515, -1515, -1516, -1516, 20310,
+      20310, 20311, 20311, 20311, 20310, 20329, 20329, 20329, 20330, 20329,
+      20329, 7172,  7173,  7173,  7172,  7172,  -814,  -814,  -813,  -812,
+      -813,  -814,  -795,  -795,  -794,  -793,  -794,  -795,  20310, 20310,
+      20311, 20312, 20311, 20310, 20329, 20329, 20330, 20330, 20330, 20329,
+      7172,  7172,  7173,  7173,  7173,  7172,  -813,  -813,  -812,  -812,
+      -813,  -813,  -794,  -794,  -793,  -793,  -793,  -794,  19487, 19487,
+      19488, 19488, 19488, 19487, 19531, 19531, 19532, 19532, 19532, 19531,
+      6400,  6399,  6400,  6401,  6400,  6400,  -1561, -1561, -1560, -1559,
+      -1560, -1561, -1516, -1516, -1515, -1515, -1515, -1516, -1516, -1516
+    };
+
+    double fieldr[] = {
+      -1357, -1356, -1353, -1354, -1354, -1357, -1366, -1366, -1362, -1363,
+      -1363, -1366, -1378, -1374, -1375, -1375, -1378, -1388, -1388, -1385,
+      -1386, -1386, -1388, -1394, -1394, -1390, -1391, -1391, -1394, -318,
+      -318,  -314,  -315,  -316,  -318,  -321,  -321,  -317,  -318,  -319,
+      -321,  -325,  -321,  -322,  -322,  -325,  -328,  -328,  -324,  -325,
+      -326,  -328,  -330,  -331,  -326,  -327,  -328,  -330,  312,   312,
+      316,   315,   315,   312,   315,   315,   319,   318,   318,   315,
+      319,   318,   323,   322,   321,   319,   322,   322,   326,   325,
+      325,   322,   324,   324,   328,   327,   327,   324,   1351,  1351,
+      1356,  1354,  1354,  1351,  1360,  1360,  1365,  1363,  1363,  1360,
+      1372,  1372,  1377,  1375,  1375,  1372,  1383,  1383,  1387,  1386,
+      1386,  1383,  1388,  1388,  1393,  1391,  1391,  1388,  1388,  1388
+    };
+
+    double fieldphi[] = {
+      -2, 7,  3,  1,  6, -2, -2, 7, 3,  1,  6, -2, -2, 3, 1,  6,  -2, -2, 7, 3,
+      1,  6,  -2, -2, 7, 3,  1,  6, -2, -1, 7, 3,  1,  6, -1, -1, 7,  3,  1, 6,
+      -1, -1, 3,  1,  6, -1, -1, 7, 3,  1,  6, -1, -1, 8, 3,  1,  6,  -1, 1, 7,
+      3,  2,  6,  1,  1, 7,  3,  2, 6,  1,  1, 7,  3,  2, 6,  1,  1,  7,  3, 2,
+      6,  1,  0,  8,  3, 2,  6,  0, 2,  7,  3, 2,  6,  2, 2,  7,  3,  2,  6, 2,
+      2,  7,  3,  2,  6, 2,  2,  7, 3,  2,  6, 2,  2,  8, 3,  2,  6,  2,  2, 2
+    };
+
+    int id{ 5 };
+    double zmin{ -1400 }, zmax{ 1400 }, rmin{ 1200 }, rmax{ 1300 }, phimin{ 0 },
+      phimax{ 6.28319 }, bscale{ 1e-07 };
+    BFieldZone zone(id, zmin, zmax, rmin, rmax, phimin, phimax, bscale);
+    // add in field
+    int nmeshz{ 4 }, nmeshr{ 5 }, nmeshphi{ 6 };
+    int nfield = nmeshz * nmeshr * nmeshphi;
+
+    zone.reserve(nmeshz, nmeshr, nmeshphi);
+
+    double meshz[] = { -1400, -466.93, 466.14, 1400 };
+
+    for (int j = 0; j < nmeshz; j++) {
+      zone.appendMesh(0, meshz[j]);
+
+      if (doDebug){
+        std::cout << "j, meshz " << j << ", " << meshz[j] <<'\n';
+      }
+    }
 
+    if (doDebug){
+      std::cout << "did  meshz " << '\n';
+    }
 
-        // Create a field zone
+    double meshr[] = { 1200, 1225, 1250, 1275, 1300 };
 
+    for (int j = 0; j < nmeshr; j++) {
+      zone.appendMesh(1, meshr[j]);
+    }
 
-        double fieldz[] = { 19487, 19487, 19488, 19488, 19487, 19487, 19531, 19531, 19532, 19532, 19531, 19531,  6399,  6400,  6400,  6400,  6399, -1561, -1561, -1560, -1560, -1560, -1561, -1516, -1516, -1515, -1515, -1516, -1516, 20310, 20310, 20311, 20311, 20311, 20310, 20329, 20329, 20329, 20330, 20329, 20329,  7172,  7173,  7173,  7172,  7172,  -814,  -814,  -813,  -812,  -813,  -814,  -795,  -795,  -794,  -793,  -794,  -795, 20310, 20310, 20311, 20312, 20311, 20310, 20329, 20329, 20330, 20330, 20330, 20329,  7172,  7172,  7173,  7173,  7173,  7172,  -813,  -813,  -812,  -812,  -813,  -813,  -794,  -794,  -793,  -793,  -793,  -794, 19487, 19487, 19488, 19488, 19488, 19487, 19531, 19531, 19532, 19532, 19532, 19531,  6400,  6399,  6400,  6401,  6400,  6400, -1561, -1561, -1560, -1559, -1560, -1561, -1516, -1516, -1515, -1515, -1515, -1516, -1516, -1516 };
+    if (doDebug){
+      std::cout << "did  meshr " << '\n';
+    }
 
-        double  fieldr[]  = { -1357, -1356, -1353, -1354, -1354, -1357, -1366, -1366, -1362, -1363, -1363, -1366, -1378, -1374, -1375, -1375, -1378, -1388, -1388, -1385, -1386, -1386, -1388, -1394, -1394, -1390, -1391, -1391, -1394,  -318,  -318,  -314,  -315,  -316,  -318,  -321,  -321,  -317,  -318,  -319,  -321,  -325,  -321,  -322,  -322,  -325,  -328,  -328,  -324,  -325,  -326,  -328,  -330,  -331,  -326,  -327,  -328,  -330,   312,   312,   316,   315,   315,   312,   315,   315,   319,   318,   318,   315,   319,   318,   323,   322,   321,   319,   322,   322,   326,   325,   325,   322,   324,   324,   328,   327,   327,   324,  1351,  1351,  1356,  1354,  1354,  1351,  1360,  1360,  1365,  1363,  1363,  1360,  1372,  1372,  1377,  1375,  1375,  1372,  1383,  1383,  1387,  1386,  1386,  1383,  1388,  1388,  1393,  1391,  1391,  1388,  1388,  1388 };
+    double meshphi[] = { 0, 1.25664, 2.51327, 3.76991, 5.02655, 6.28318 };
 
-        double fieldphi[] = { -2,  7,  3,  1,  6,  -2,  -2,  7,  3,  1,  6,  -2,  -2        ,  3,  1,  6,  -2,  -2,  7,  3,  1,  6,  -2,  -2,  7,  3,  1,  6,  -2,  -1,  7,  3,  1,  6,  -1,  -1,  7,  3,  1,  6,  -1,  -1,  3,  1,  6,  -1,  -1,  7,  3,  1,  6,  -1,  -1,  8,  3,  1,  6,  -1,  1,  7,  3,  2,  6,  1,  1,  7,  3,  2,  6,  1,  1,  7,  3,  2,  6,  1,  1,  7,  3,  2,  6,  1,  0,  8,  3,  2,  6,  0,  2,  7,  3,  2,  6,  2,  2,  7,  3,  2,  6,  2,  2,  7,  3,  2,  6,  2,  2,  7,  3,  2,  6,  2,  2,  8,  3,  2,  6,  2,  2,  2 };
-            
-
-        int id{5};
-        double zmin{-1400}, zmax{1400}, rmin{1200}, rmax{1300}, phimin{0}, phimax{6.28319}, bscale{1e-07};
-        BFieldZone zone(id, zmin, zmax, rmin, rmax, phimin, phimax, bscale);
-        // add in field
-        int nmeshz{4}, nmeshr{5}, nmeshphi{6};
-        int nfield = nmeshz * nmeshr * nmeshphi;
-
-        zone.reserve(nmeshz, nmeshr, nmeshphi);
-
-        double meshz[] = { -1400, -466.93, 466.14, 1400 };
-        
-        for (int j = 0; j < nmeshz; j++) {
-            zone.appendMesh(0, meshz[j]);
-
-            if (doDebug) std::cout << "j, meshz " << j << ", " << meshz[j] << std::endl;
-            
-        }
-
-        if (doDebug) std::cout << "did  meshz " << std::endl;
-        
-        double meshr[] = { 1200, 1225, 1250, 1275, 1300 };
-        
-        for (int j = 0; j < nmeshr; j++) {
-            zone.appendMesh(1, meshr[j]);
-        }
-
-        if (doDebug) std::cout << "did  meshr " << std::endl;
-
-        double meshphi[] = { 0, 1.25664, 2.51327, 3.76991, 5.02655, 6.28318 };
-
-        for (int j = 0; j < nmeshphi; j++) {
-            zone.appendMesh(2, meshphi[j]);
-        }
-
-        if (doDebug) std::cout << "did  meshphi " << std::endl;
-
-        for (int j = 0; j < nfield; j++) {
-            BFieldVector<short> field(fieldz[j], fieldr[j], fieldphi[j]);
-            zone.appendField(field);
-        }
-        
-        if (doDebug) std::cout << "did  field " << std::endl;
-
-        // build (trivial) look up table for zone
-        zone.buildLUT();
-        
-        if (doDebug) std::cout << "did  buildLUT " << std::endl;
-
-        
-        // fill the cache, pass in current scale factor
-        BFieldCache cache3d;
-        double z{0}, r{1250}, phi{1.6};
-
-        // if (doDebug) std::cout << "do getCache old " << std::endl;
-
-        // zone.m_doNew = false;
-        
-        // zone.getCache(z, r, phi, cache3d, 1);
-
-        // cache3d.printField();
-
-        // if (doDebug) std::cout << "do getCache new " << std::endl;
-        
-        // zone.m_doNew = true;
-        
-        // zone.getCache(z, r, phi, cache3d, 1);
-
-        // cache3d.printField();
-
-
-        // get field at steps of 10 mm from 1200 to 1300
-
-        int status{0};
-
-        double z0      = z;
-        double r0      = 1200;
-        double phi0    = phi;
-        double xyz[3]  = {0, 0, 0};
-        double bxyz[3] = {0, 0, 0};
-        double bxyz_std[3][10] = { 
-            { -2.83727e-07, -2.81403e-07, -2.79079e-07, -2.76755e-07, -2.74431e-07,
-              -2.72107e-07, -2.69782e-07, -2.67458e-07, -2.65134e-07, -2.6281e-07 },
-            {  9.47007e-08,  7.49033e-08,  5.51058e-08,  3.53084e-08,  1.5511e-08,
-              -4.28645e-09, -2.40839e-08, -4.38813e-08, -6.36787e-08, -8.34762e-08 },
-            {  0.00308551,   0.00255923,   0.00203296,   0.00150669,   0.000980422,
-               0.000454151, -7.21201e-05, -0.000598391, -0.00112466,  -0.00165093 }
-        };
-        
-// get field std: i, bxyz_std 0 -2.83727e-07, 9.47007e-08, 0.00308551
-// get field new: i, bxyz_new 0 -2.83727e-07, 9.47007e-08, 0.00308551
-// get field std: i, bxyz_std 1 -2.81403e-07, 7.49033e-08, 0.00255923
-// get field new: i, bxyz_new 1 -2.81403e-07, 7.49033e-08, 0.00255923
-// get field std: i, bxyz_std 2 -2.79079e-07, 5.51058e-08, 0.00203296
-// get field new: i, bxyz_new 2 -2.79079e-07, 5.51058e-08, 0.00203296
-// get field std: i, bxyz_std 3 -2.76755e-07, 3.53084e-08, 0.00150669
-// get field new: i, bxyz_new 3 -2.76755e-07, 3.53084e-08, 0.00150669
-// get field std: i, bxyz_std 4 -2.74431e-07, 1.5511e-08, 0.000980422
-// get field new: i, bxyz_new 4 -2.74431e-07, 1.5511e-08, 0.000980422
-// get field std: i, bxyz_std 5 -2.72107e-07, -4.28645e-09, 0.000454151
-// get field new: i, bxyz_new 5 -2.72107e-07, -4.28645e-09, 0.000454151
-// get field std: i, bxyz_std 6 -2.69782e-07, -2.40839e-08, -7.21201e-05
-// get field new: i, bxyz_new 6 -2.69782e-07, -2.40839e-08, -7.21201e-05
-// get field std: i, bxyz_std 7 -2.67458e-07, -4.38813e-08, -0.000598391
-// get field new: i, bxyz_new 7 -2.67458e-07, -4.38813e-08, -0.000598391
-// get field std: i, bxyz_std 8 -2.65134e-07, -6.36787e-08, -0.00112466
-// get field new: i, bxyz_new 8 -2.65134e-07, -6.36787e-08, -0.00112466
-// get field std: i, bxyz_std 9 -2.6281e-07, -8.34762e-08, -0.00165093
-// get field new: i, bxyz_new 9 -2.6281e-07, -8.34762e-08, -0.00165093
-
-        
-        for (unsigned int i = 0; i < 10; ++i) {
-            double r1 = r0 + 5 + i*10.;
-
-            xyz[0] = r1 * cos(phi0);
-            xyz[1] = r1 * sin(phi0);
-            xyz[2] = z0;
-            
-            // do interpolation (cache3d has correct scale factor)
-            zone.getCache(z, r, phi, cache3d, 1);
-            cache3d.getB(xyz, r1, phi, bxyz, 0);
-
-            std::cout << "get field std: i, bxyz " << i << " "
-                      << bxyz[0] << ", " << bxyz[1] << ", "
-                      << bxyz[2]
-                      << " fractional diff gt 10^-5: "
-                      << int(fabs(bxyz[0] - bxyz_std[0][i])/bxyz[0] > 0.00001) << ", "
-                      << int(fabs(bxyz[1] - bxyz_std[1][i])/bxyz[1] > 0.00001) << ", "
-                      << int(fabs(bxyz[2] - bxyz_std[2][i])/bxyz[2] > 0.00001) 
-                      << std::endl;
-
-            if (fabs(bxyz[0] - bxyz_std[0][i]) > 0.00001) {
-                std::cout << "failed bz comparison - bz, bz std " << bxyz[0] << ", " << bxyz_std[0][i] << std::endl;
-                status = 1;
-            }
-            if (fabs(bxyz[1] - bxyz_std[1][i]) > 0.00001) {
-                std::cout << "failed br comparison" << std::endl;
-                std::cout << "failed br comparison - br, br std " << bxyz[1] << ", " << bxyz_std[1][i] << std::endl;
-                status = 1;
-            }
-            if (fabs(bxyz[2] - bxyz_std[2][i]) > 0.00001) {
-                std::cout << "failed bphi comparison" << std::endl;
-                std::cout << "failed bphi comparison - bphi, bphi std " << bxyz[2] << ", " << bxyz_std[2][i] << std::endl;
-                status = 1;
-            }
-            
-        }
-        
-        std::cout << "runTest: status " << status << std::endl;
-
-        return status;
-        
+    for (int j = 0; j < nmeshphi; j++) {
+      zone.appendMesh(2, meshphi[j]);
     }
 
-private:
-};
+    if (doDebug){
+      std::cout << "did  meshphi " << '\n';
+    }
 
+    for (int j = 0; j < nfield; j++) {
+      BFieldVector<short> field(fieldz[j], fieldr[j], fieldphi[j]);
+      zone.appendField(field);
+    }
 
+    if (doDebug){
+      std::cout << "did  field " << '\n';
+    }
 
+    // build (trivial) look up table for zone
+    zone.buildLUT();
 
-int main() {
+    if (doDebug){
+      std::cout << "did  buildLUT " << '\n';
+    }
 
-    std::cout << "start BFieldExample test" << std::endl;
+    // fill the cache, pass in current scale factor
+    BFieldCache cache3d;
+    double z{ 0 }, r{ 1250 }, phi{ 1.6 };
+
+    // if (doDebug) std::cout << "do getCache old " << '\n';
+    // zone.m_doNew = false;
+    // zone.getCache(z, r, phi, cache3d, 1);
+    // cache3d.printField();
+    // if (doDebug) std::cout << "do getCache new " << '\n';
+    // zone.m_doNew = true;
+    // zone.getCache(z, r, phi, cache3d, 1);
+    // cache3d.printField();
+
+    // get field at steps of 10 mm from 1200 to 1300
+    int status{ 0 };
+
+    double z0 = z;
+    double r0 = 1200;
+    double phi0 = phi;
+    double xyz[3] = { 0, 0, 0 };
+    double bxyz[3] = { 0, 0, 0 };
+    double bxyz_std[3][10] = { { -2.83727e-07,
+                                 -2.81403e-07,
+                                 -2.79079e-07,
+                                 -2.76755e-07,
+                                 -2.74431e-07,
+                                 -2.72107e-07,
+                                 -2.69782e-07,
+                                 -2.67458e-07,
+                                 -2.65134e-07,
+                                 -2.6281e-07 },
+                               { 9.47007e-08,
+                                 7.49033e-08,
+                                 5.51058e-08,
+                                 3.53084e-08,
+                                 1.5511e-08,
+                                 -4.28645e-09,
+                                 -2.40839e-08,
+                                 -4.38813e-08,
+                                 -6.36787e-08,
+                                 -8.34762e-08 },
+                               { 0.00308551,
+                                 0.00255923,
+                                 0.00203296,
+                                 0.00150669,
+                                 0.000980422,
+                                 0.000454151,
+                                 -7.21201e-05,
+                                 -0.000598391,
+                                 -0.00112466,
+                                 -0.00165093 } };
+
+    // get field std: i, bxyz_std 0 -2.83727e-07, 9.47007e-08, 0.00308551
+    // get field new: i, bxyz_new 0 -2.83727e-07, 9.47007e-08, 0.00308551
+    // get field std: i, bxyz_std 1 -2.81403e-07, 7.49033e-08, 0.00255923
+    // get field new: i, bxyz_new 1 -2.81403e-07, 7.49033e-08, 0.00255923
+    // get field std: i, bxyz_std 2 -2.79079e-07, 5.51058e-08, 0.00203296
+    // get field new: i, bxyz_new 2 -2.79079e-07, 5.51058e-08, 0.00203296
+    // get field std: i, bxyz_std 3 -2.76755e-07, 3.53084e-08, 0.00150669
+    // get field new: i, bxyz_new 3 -2.76755e-07, 3.53084e-08, 0.00150669
+    // get field std: i, bxyz_std 4 -2.74431e-07, 1.5511e-08, 0.000980422
+    // get field new: i, bxyz_new 4 -2.74431e-07, 1.5511e-08, 0.000980422
+    // get field std: i, bxyz_std 5 -2.72107e-07, -4.28645e-09, 0.000454151
+    // get field new: i, bxyz_new 5 -2.72107e-07, -4.28645e-09, 0.000454151
+    // get field std: i, bxyz_std 6 -2.69782e-07, -2.40839e-08, -7.21201e-05
+    // get field new: i, bxyz_new 6 -2.69782e-07, -2.40839e-08, -7.21201e-05
+    // get field std: i, bxyz_std 7 -2.67458e-07, -4.38813e-08, -0.000598391
+    // get field new: i, bxyz_new 7 -2.67458e-07, -4.38813e-08, -0.000598391
+    // get field std: i, bxyz_std 8 -2.65134e-07, -6.36787e-08, -0.00112466
+    // get field new: i, bxyz_new 8 -2.65134e-07, -6.36787e-08, -0.00112466
+    // get field std: i, bxyz_std 9 -2.6281e-07, -8.34762e-08, -0.00165093
+    // get field new: i, bxyz_new 9 -2.6281e-07, -8.34762e-08, -0.00165093
+
+    for (unsigned int i = 0; i < 10; ++i) {
+      double r1 = r0 + 5 + i * 10.;
+
+      xyz[0] = r1 * cos(phi0);
+      xyz[1] = r1 * sin(phi0);
+      xyz[2] = z0;
+
+      // do interpolation (cache3d has correct scale factor)
+      zone.getCache(z, r, phi, cache3d, 1);
+      cache3d.getB(xyz, r1, phi, bxyz, 0);
+
+      std::cout << "get field std: i, bxyz " << i << " " << bxyz[0] << ", "
+                << bxyz[1] << ", " << bxyz[2] << " fractional diff gt 10^-5: "
+                << int(fabs(bxyz[0] - bxyz_std[0][i]) / bxyz[0] > 0.00001)
+                << ", "
+                << int(fabs(bxyz[1] - bxyz_std[1][i]) / bxyz[1] > 0.00001)
+                << ", "
+                << int(fabs(bxyz[2] - bxyz_std[2][i]) / bxyz[2] > 0.00001)
+                << '\n';
+
+      if (fabs(bxyz[0] - bxyz_std[0][i]) > 0.00001) {
+        std::cout << "failed bz comparison - bz, bz std " << bxyz[0] << ", "
+                  << bxyz_std[0][i] << '\n';
+        status = 1;
+      }
+      if (fabs(bxyz[1] - bxyz_std[1][i]) > 0.00001) {
+        std::cout << "failed br comparison" << '\n';
+        std::cout << "failed br comparison - br, br std " << bxyz[1] << ", "
+                  << bxyz_std[1][i] << '\n';
+        status = 1;
+      }
+      if (fabs(bxyz[2] - bxyz_std[2][i]) > 0.00001) {
+        std::cout << "failed bphi comparison" << '\n';
+        std::cout << "failed bphi comparison - bphi, bphi std " << bxyz[2]
+                  << ", " << bxyz_std[2][i] << '\n';
+        status = 1;
+      }
+    }
 
+    std::cout << "runTest: status " << status << '\n';
 
-    BFieldTest bftest;
-    int status = bftest.runTest(true);
+    return status;
+  }
 
-    if (status == 0) {
-        std::cout << "Test passed OK" << std::endl;
-        return 0;
-    }
-    else {
-        std::cout << "Test failed" << std::endl;
-        return 1;
-    }
-    
+private:
+};
+
+int
+main()
+{
+
+  std::cout << "start BFieldExample test" << '\n';
+
+  BFieldTest bftest;
+  int status = bftest.runTest(true);
+
+  if (status == 0) {
+    std::cout << "Test passed OK" << '\n';
+    return 0;
+  } else {
+    std::cout << "Test failed" << '\n';
+    return 1;
+  }
 }
-- 
GitLab


From 330461b3dd200622f9dfb4f3883f1782187c1548 Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Fri, 21 Aug 2020 01:59:33 +0100
Subject: [PATCH 2/2] Fix compilation of printField

---
 MagneticField/MagFieldElements/src/BFieldCache.cxx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/MagneticField/MagFieldElements/src/BFieldCache.cxx b/MagneticField/MagFieldElements/src/BFieldCache.cxx
index 979b01e53f7c..6e1c2871078e 100644
--- a/MagneticField/MagFieldElements/src/BFieldCache.cxx
+++ b/MagneticField/MagFieldElements/src/BFieldCache.cxx
@@ -2,12 +2,12 @@
   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
-#include "MagFieldElement/BFieldCache.h"
+#include "MagFieldElements/BFieldCache.h"
 #include <iostream>
 
 // set the field values at each corner (rescale for current scale factor)
 void
-BFieldCache::printField()
+BFieldCache::printField() const
 {
   // print out field values
   std::cout << "Field at corner i, for each component j (Bz, Br, Bphi)"
-- 
GitLab