adding in unit test for field optimizations. This version is a detailed...

adding in unit test for field optimizations. This version is a detailed framework for testing, which will be simplified in further commits.

MagneticField/MagFieldElements/CMakeLists.txt

@@ -13,3 +13,10 @@ atlas_add_library( MagFieldElements
                    INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
                    LINK_LIBRARIES CxxUtils EventPrimitives GaudiKernel
                    PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} PathResolver )
+
+		 
+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" )		

MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h

@@ -80,6 +80,26 @@ public:
     }
   }
 
+  // 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;
+    }
+  }
+
   // set the multiplicative factor for the field vectors
   void setBscale(double bscale) { m_scale = bscale; }
   float bscale() { return m_scale; }

MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h

@@ -98,6 +98,9 @@ public:
   double bscale() const { return m_scale; }
   int memSize() const;
 
+  bool m_doNew{false};
+
+    
 protected:
   std::array<double, 3> m_min;
   std::array<double, 3> m_max;
@@ -112,6 +115,7 @@ private:
   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

MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc

@@ -70,10 +70,14 @@ 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
   const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner
+
+  if (!m_doNew) {
+
   std::array<BFieldVector<T>, 8> field = {
     m_field[im0],
     m_field[im0 + 1],
@@ -84,7 +88,50 @@ BFieldMesh<T>::getCache(double z,
     m_field[im0 + m_zoff + m_roff],
     m_field[im0 + m_zoff + m_roff + 1]
   };
+
   cache.setField(field, scaleFactor);
+
+  }
+  else {
+      
+  // cache.printField();
+
+  const double sf = scaleFactor;
+
+  double field1[3][8] = {
+      { sf * m_field[im0][0],                        
+        sf * m_field[im0 + 1][0],                    
+        sf * m_field[im0 + m_roff][0],               
+        sf * m_field[im0 + m_roff + 1][0],           
+        sf * m_field[im0 + m_zoff][0],               
+        sf * m_field[im0 + m_zoff + 1][0],           
+        sf * m_field[im0 + m_zoff + m_roff][0],      
+        sf * m_field[im0 + m_zoff + m_roff + 1][0] },
+      { sf * m_field[im0][1],                      
+        sf * m_field[im0 + 1][1],                  
+        sf * m_field[im0 + m_roff][1],             
+        sf * m_field[im0 + m_roff + 1][1],         
+        sf * m_field[im0 + m_zoff][1],             
+        sf * m_field[im0 + m_zoff + 1][1],         
+        sf * m_field[im0 + m_zoff + m_roff][1],    
+        sf * m_field[im0 + m_zoff + m_roff + 1][1] },
+      { sf * m_field[im0][2],                     
+        sf * m_field[im0 + 1][2],                 
+        sf * m_field[im0 + m_roff][2],            
+        sf * m_field[im0 + m_roff + 1][2],        
+        sf * m_field[im0 + m_zoff][2],            
+        sf * m_field[im0 + m_zoff + 1][2],        
+        sf * m_field[im0 + m_zoff + m_roff][2],   
+        sf * m_field[im0 + m_zoff + m_roff + 1][2] }
+  };
+  cache.setField(field1);
+  
+  // cache.printField();
+
+  }
+  
+
+  
   // store the B scale
   cache.setBscale(m_scale);
 }

MagneticField/MagFieldElements/test/BFieldExample_test.cxx

+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <iostream>
+#include "MagFieldElements/BFieldCache.h"
+#include "MagFieldElements/BFieldZone.h"
+#include <unistd.h>
+
+/**
+ * @brief Specimen of service caching expensive calculations.
+ * It uses internally std::map which is not thread-safely expandable.
+ **/
+class BFieldTest {
+public:
+    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] << 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_std[3] = {0, 0, 0};
+        double bxyz_new[3] = {0, 0, 0};
+        
+        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.m_doNew = false;
+            zone.getCache(z, r, phi, cache3d, 1);
+            cache3d.getB(xyz, r1, phi, bxyz_std, 0);
+
+            std::cout << "get field std: i, bxyz_std " << i << " "
+                      << bxyz_std[0] << ", " << bxyz_std[1] << ", "
+                      << bxyz_std[2] << std::endl;
+
+            zone.m_doNew = true;
+            zone.getCache(z, r, phi, cache3d, 1);
+            cache3d.getB(xyz, r1, phi, bxyz_new, 0);
+
+            std::cout << "get field new: i, bxyz_new " << i << " "
+                      << bxyz_new[0] << ", " << bxyz_new[1] << ", "
+                      << bxyz_new[2] << std::endl;
+
+            if (bxyz_std[0] != bxyz_new[0]) {
+                std::cout << "failed bz comparison" << std::endl;
+                status = 1;
+            }
+            if (bxyz_std[1] != bxyz_new[1]) {
+                std::cout << "failed br comparison" << std::endl;
+                status = 1;
+            }
+            if (bxyz_std[2] != bxyz_new[2]) {
+                std::cout << "failed bphi comparison" << std::endl;
+                status = 1;
+            }
+            
+        }
+        
+        std::cout << "runTest: status " << status << std::endl;
+
+        return status;
+        
+    }
+
+private:
+};
+
+
+
+
+int main() {
+
+    std::cout << "start BFieldExample test" << std::endl;
+
+
+    BFieldTest bftest;
+    int status = bftest.runTest(true);
+
+    if (status == 0) {
+        std::cout << "Test passed OK" << std::endl;
+        return 0;
+    }
+    else {
+        std::cout << "Test failed" << std::endl;
+        return 1;
+    }
+    
+}