diff --git a/MagneticField/MagFieldElements/CMakeLists.txt b/MagneticField/MagFieldElements/CMakeLists.txt
index d24244fcfb5aba4d9e7d09679e0fc93b42d3acaf..f7929d3a6cda021157f6b392bb73e39e158fa708 100644
--- a/MagneticField/MagFieldElements/CMakeLists.txt
+++ b/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
+ 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" )
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
index e42b999ef3bebe16e2df4f9564f8d7de6a9fe526..8282eb91153ed3a7297593e6b7da13787a262b5d 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
@@ -55,28 +55,52 @@ public:
m_invr = 1.0 / (rmax - rmin);
m_invphi = 1.0 / (phimax - phimin);
}
- // 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];
- }
- }
+
+ // **
+ // ** 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];
+ // }
+ // }
+ // }
+
+ // set field array, filled externally
+ void setField(double field[][8]) {
+ std::copy( &field[0][0], &field[0][0] + 3*8, &m_field[0][0]);
}
+
- void setField(const std::array<BFieldVector<short>, 8>& field,
- double scaleFactor = 1.0)
+ // set the field values at each corner (rescale for current scale factor)
+ void printField()
{
- // We pass array of 8 elements with 3 entries
- // Which go to 3x8 matrix
+ // 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) {
- m_field[j][i] = scaleFactor * field[i][j];
+ std::cout << i << "," << j << ": " << m_field[j][i] << ", ";
}
+ std::cout << std::endl;
}
}
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
index 9c85785a4af2c2a59ca69d5a3409b8986c81660d..7424d87ab3de46217aa4186211d851693c20b545 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
@@ -112,6 +112,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
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
index 810c9739fdc2588d80b5cf1a200d2b29049ca1a2..114bffa4fdcebe9052d81b1141c4203aea1d765b 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
@@ -70,21 +70,60 @@ 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
- 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]
+
+
+ // ** The following field initialization is commented out, and replaced by subsequent matrix
+ // ** initialization as optimization test
+
+
+ // 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]
+ // };
+
+ // cache.setField(field, scaleFactor);
+
+ 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(field, scaleFactor);
+ cache.setField(field1);
+
// store the B scale
cache.setBscale(m_scale);
}
diff --git a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..34ca83c9e083321c520689ea9a72d2aeb8aab821
--- /dev/null
+++ b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
@@ -0,0 +1,205 @@
+/*
+ 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[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;
+
+ }
+
+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;
+ }
+
+}