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