From 64693a44b0dd7b4a685528197d4a2506fdfdd88f Mon Sep 17 00:00:00 2001 From: Christos Anastopoulos <christos.anastopoulos@cern.ch> Date: Tue, 25 Aug 2020 15:44:19 +0000 Subject: [PATCH] MagFieldElement further follow ATLAS style on inline conventions --- MagneticField/MagFieldElements/CMakeLists.txt | 1 - .../MagFieldElements/BFieldCache.h | 5 +- .../MagFieldElements/BFieldCache.icc | 100 +---------------- .../MagFieldElements/BFieldCacheZR.h | 34 +----- .../MagFieldElements/BFieldCacheZR.icc | 95 ++++++---------- .../MagFieldElements/src/BFieldCache.cxx | 104 ++++++++++++++++-- .../MagFieldElements/src/BFieldCacheZR.cxx | 70 ++++++++++++ .../test/BFieldExample_test.cxx | 9 -- 8 files changed, 204 insertions(+), 214 deletions(-) create mode 100644 MagneticField/MagFieldElements/src/BFieldCacheZR.cxx diff --git a/MagneticField/MagFieldElements/CMakeLists.txt b/MagneticField/MagFieldElements/CMakeLists.txt index 81af6c3cfef9..7ed582eeedd9 100644 --- a/MagneticField/MagFieldElements/CMakeLists.txt +++ b/MagneticField/MagFieldElements/CMakeLists.txt @@ -17,5 +17,4 @@ atlas_add_library( MagFieldElements atlas_add_test( BFieldExample_test SOURCES test/BFieldExample_test.cxx - INCLUDE_DIRS LINK_LIBRARIES MagFieldElements) diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h index d4895554493a..7ba7c5ffed86 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h @@ -19,7 +19,6 @@ #include "CxxUtils/restrict.h" #include "MagFieldElements/BFieldVector.h" -#include <cmath> class BFieldCache { @@ -28,6 +27,7 @@ public: BFieldCache() = default; // make this cache invalid, so that inside() will fail void invalidate(); + // set the z, r, phi range that defines the bin void setRange(double zmin, double zmax, @@ -53,9 +53,6 @@ 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: // bin range in z double m_zmin = 0.0; diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc index 0b81313acc86..d9528098145a 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc @@ -1,7 +1,7 @@ /* Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ - +#include <cmath> inline void BFieldCache::invalidate() { @@ -56,102 +56,4 @@ BFieldCache::inside(double z, double r, double phi) const r >= m_rmin && r <= m_rmax); } -inline void -BFieldCache::getB(const double* ATH_RESTRICT xyz, - double r, - double phi, - double* ATH_RESTRICT B, - double* ATH_RESTRICT deriv) const -{ - - const double x = xyz[0]; - const double y = xyz[1]; - const double z = xyz[2]; - - // make sure phi is inside [m_phimin,m_phimax] - if (phi < m_phimin) { - phi += 2 * M_PI; - } - // fractional position inside this bin - const double fz = (z - m_zmin) * m_invz; - const double gz = 1.0 - fz; - const double fr = (r - m_rmin) * m_invr; - const double gr = 1.0 - fr; - const double fphi = (phi - m_phimin) * m_invphi; - 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 components - const double* field = m_field[i]; - Bzrphi[i] = scale * (gz * (gr * (gphi * field[0] + fphi * field[1]) + - fr * (gphi * field[2] + fphi * field[3])) + - fz * (gr * (gphi * field[4] + fphi * field[5]) + - fr * (gphi * field[6] + fphi * field[7]))); - } - // convert (Bz,Br,Bphi) to (Bx,By,Bz) - double invr; - double c; - double s; - if (r > 0.0) { - invr = 1.0 / r; - c = x * invr; - s = y * invr; - } else { - invr = 0.0; - c = cos(m_phimin); - s = sin(m_phimin); - } - B[0] = Bzrphi[1] * c - Bzrphi[2] * s; - B[1] = Bzrphi[1] * s + Bzrphi[2] * c; - B[2] = Bzrphi[0]; - - // compute field derivatives if requested - if (deriv) { - 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; - std::array<double, 3> dBdr; - std::array<double, 3> dBdphi; - - for (int j = 0; j < 3; ++j) { // Bz, Br, Bphi components - 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]))); - 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]))); - 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]))); - } - // convert to cartesian coordinates - const double cc = c * c; - const double cs = c * s; - const double ss = s * s; - const double ccinvr = cc * invr; - const double csinvr = cs * invr; - const double ssinvr = ss * invr; - const double sinvr = s * invr; - const double cinvr = c * invr; - deriv[0] = cc * dBdr[1] - cs * dBdr[2] - csinvr * dBdphi[1] + - ssinvr * dBdphi[2] + sinvr * B[1]; - deriv[1] = cs * dBdr[1] - ss * dBdr[2] + ccinvr * dBdphi[1] - - csinvr * dBdphi[2] - cinvr * B[1]; - deriv[2] = c * dBdz[1] - s * dBdz[2]; - deriv[3] = cs * dBdr[1] + cc * dBdr[2] - ssinvr * dBdphi[1] - - csinvr * dBdphi[2] - sinvr * B[0]; - deriv[4] = ss * dBdr[1] + cs * dBdr[2] + csinvr * dBdphi[1] + - ccinvr * dBdphi[2] + cinvr * B[0]; - deriv[5] = s * dBdz[1] + c * dBdz[2]; - deriv[6] = c * dBdr[0] - sinvr * dBdphi[0]; - deriv[7] = s * dBdr[0] + cinvr * dBdphi[0]; - deriv[8] = dBdz[0]; - } -} diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.h index 56056e02e5d2..9597916e5ca0 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.h @@ -22,40 +22,16 @@ class BFieldCacheZR { public: // default constructor sets unphysical boundaries, so that inside() will fail - BFieldCacheZR() - : m_zmin(0.0) - , m_zmax(-1.0) - , m_rmin(0.0) - , m_rmax(-1.0) - {} + BFieldCacheZR(); // invalidate this cache, so that inside() will fail - void invalidate() - { - m_rmin = 0.0; - m_rmax = -1.0; - } + void invalidate(); // set the z, r range that defines the bin - void setRange(double zmin, double zmax, double rmin, double rmax) - { - m_zmin = zmin; - m_zmax = zmax; - m_rmin = rmin; - m_rmax = rmax; - m_invz = 1.0 / (zmax - zmin); - m_invr = 1.0 / (rmax - rmin); - } + void setRange(double zmin, double zmax, double rmin, double rmax); // set the field values at each corner (rescale for current scale factor) - void setField(int i, const BFieldVectorZR& field, double scaleFactor = 1.0) - { - m_field[0][i] = scaleFactor * field[0]; - m_field[1][i] = scaleFactor * field[1]; - } + void setField(int i, const BFieldVectorZR& field, double scaleFactor = 1.0); // set the multiplicative factor for the field vectors // test if (z, r) is inside this bin - bool inside(double z, double r) const - { - return (z >= m_zmin && z <= m_zmax && r >= m_rmin && r <= m_rmax); - } + bool inside(double z, double r) const; // interpolate the field and return B[3]. // also compute field derivatives if deriv[9] is given. void getB(const double* ATH_RESTRICT xyz, diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.icc index b90a35375dd5..16aded6a769c 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCacheZR.icc @@ -1,70 +1,41 @@ /* Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ - +inline BFieldCacheZR::BFieldCacheZR() + : m_zmin(0.0) + , m_zmax(-1.0) + , m_rmin(0.0) + , m_rmax(-1.0) +{} inline void -BFieldCacheZR::getB(const double* ATH_RESTRICT xyz, - double r, - double* ATH_RESTRICT B, - double* ATH_RESTRICT deriv) const +BFieldCacheZR::invalidate() { - const double x = xyz[0]; - const double y = xyz[1]; - const double z = xyz[2]; - // fractional position inside this bin - const double fz = (z - m_zmin) * m_invz; - const double gz = 1.0 - fz; - const double fr = (r - m_rmin) * m_invr; - const double gr = 1.0 - fr; - // interpolate field values in z, r - double Bzr[2]; - for (int i = 0; i < 2; ++i) { // z, r components - const double* field = m_field[i]; - Bzr[i] = gz * (gr * field[0] + fr * field[1]) + - fz * (gr * field[2] + fr * field[3]); - } - // convert (Bz,Br) to (Bx,By,Bz) - double invr; - if (r > 0.0) { - invr = 1.0 / r; - } else { - invr = 0.0; - } - const double c(x * invr); - const double s(y * invr); - B[0] = Bzr[1] * c; - B[1] = Bzr[1] * s; - B[2] = Bzr[0]; - - // compute field derivatives if requested - if (deriv) { - std::array<double, 2> dBdz; - std::array<double, 2> dBdr; - for (int j = 0; j < 2; ++j) { // Bz, Br components - const double* field = m_field[j]; - dBdz[j] = - m_invz * (gr * (field[2] - field[0]) + - fr * (field[3] - field[1])); - dBdr[j] = - m_invr * (gz * (field[1] - field[0]) + - fz * (field[3] - field[2])); - } - // convert to cartesian coordinates - const double cc = c * c; - const double cs = c * s; - const double ss = s * s; - const double sinvr = s * invr; - const double cinvr = c * invr; - deriv[0] = cc * dBdr[1] + sinvr * B[1]; - deriv[1] = cs * dBdr[1] - cinvr * B[1]; - deriv[2] = c * dBdz[1]; - deriv[3] = cs * dBdr[1] - sinvr * B[0]; - deriv[4] = ss * dBdr[1] + cinvr * B[0]; - deriv[5] = s * dBdz[1]; - deriv[6] = c * dBdr[0]; - deriv[7] = s * dBdr[0]; - deriv[8] = dBdz[0]; - } + m_rmin = 0.0; + m_rmax = -1.0; +} +inline void +BFieldCacheZR::setRange(double zmin, double zmax, double rmin, double rmax) +{ + m_zmin = zmin; + m_zmax = zmax; + m_rmin = rmin; + m_rmax = rmax; + m_invz = 1.0 / (zmax - zmin); + m_invr = 1.0 / (rmax - rmin); +} +// set the field values at each corner (rescale for current scale factor) +inline void +BFieldCacheZR::setField(int i, + const BFieldVectorZR& field, + double scaleFactor) +{ + m_field[0][i] = scaleFactor * field[0]; + m_field[1][i] = scaleFactor * field[1]; +} +inline bool +BFieldCacheZR::inside(double z, double r) const +{ + return (z >= m_zmin && z <= m_zmax && r >= m_rmin && r <= m_rmax); } diff --git a/MagneticField/MagFieldElements/src/BFieldCache.cxx b/MagneticField/MagFieldElements/src/BFieldCache.cxx index 6e1c2871078e..4a24558c24e7 100644 --- a/MagneticField/MagFieldElements/src/BFieldCache.cxx +++ b/MagneticField/MagFieldElements/src/BFieldCache.cxx @@ -3,21 +3,105 @@ */ #include "MagFieldElements/BFieldCache.h" -#include <iostream> +#include <cmath> -// set the field values at each corner (rescale for current scale factor) void -BFieldCache::printField() const +BFieldCache::getB(const double* ATH_RESTRICT xyz, + double r, + double phi, + double* ATH_RESTRICT B, + double* ATH_RESTRICT deriv) const { - // 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] << ", "; + const double x = xyz[0]; + const double y = xyz[1]; + const double z = xyz[2]; + + // make sure phi is inside [m_phimin,m_phimax] + if (phi < m_phimin) { + phi += 2 * M_PI; + } + // fractional position inside this bin + const double fz = (z - m_zmin) * m_invz; + const double gz = 1.0 - fz; + const double fr = (r - m_rmin) * m_invr; + const double gr = 1.0 - fr; + const double fphi = (phi - m_phimin) * m_invphi; + 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 components + const double* field = m_field[i]; + Bzrphi[i] = scale * (gz * (gr * (gphi * field[0] + fphi * field[1]) + + fr * (gphi * field[2] + fphi * field[3])) + + fz * (gr * (gphi * field[4] + fphi * field[5]) + + fr * (gphi * field[6] + fphi * field[7]))); + } + // convert (Bz,Br,Bphi) to (Bx,By,Bz) + double invr; + double c; + double s; + if (r > 0.0) { + invr = 1.0 / r; + c = x * invr; + s = y * invr; + } else { + invr = 0.0; + c = cos(m_phimin); + s = sin(m_phimin); + } + B[0] = Bzrphi[1] * c - Bzrphi[2] * s; + B[1] = Bzrphi[1] * s + Bzrphi[2] * c; + B[2] = Bzrphi[0]; + + // compute field derivatives if requested + if (deriv) { + 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; + std::array<double, 3> dBdr; + std::array<double, 3> dBdphi; + + for (int j = 0; j < 3; ++j) { // Bz, Br, Bphi components + 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]))); + 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]))); + 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]))); } - std::cout << '\n'; + // convert to cartesian coordinates + const double cc = c * c; + const double cs = c * s; + const double ss = s * s; + const double ccinvr = cc * invr; + const double csinvr = cs * invr; + const double ssinvr = ss * invr; + const double sinvr = s * invr; + const double cinvr = c * invr; + deriv[0] = cc * dBdr[1] - cs * dBdr[2] - csinvr * dBdphi[1] + + ssinvr * dBdphi[2] + sinvr * B[1]; + deriv[1] = cs * dBdr[1] - ss * dBdr[2] + ccinvr * dBdphi[1] - + csinvr * dBdphi[2] - cinvr * B[1]; + deriv[2] = c * dBdz[1] - s * dBdz[2]; + deriv[3] = cs * dBdr[1] + cc * dBdr[2] - ssinvr * dBdphi[1] - + csinvr * dBdphi[2] - sinvr * B[0]; + deriv[4] = ss * dBdr[1] + cs * dBdr[2] + csinvr * dBdphi[1] + + ccinvr * dBdphi[2] + cinvr * B[0]; + deriv[5] = s * dBdz[1] + c * dBdz[2]; + deriv[6] = c * dBdr[0] - sinvr * dBdphi[0]; + deriv[7] = s * dBdr[0] + cinvr * dBdphi[0]; + deriv[8] = dBdz[0]; } } + diff --git a/MagneticField/MagFieldElements/src/BFieldCacheZR.cxx b/MagneticField/MagFieldElements/src/BFieldCacheZR.cxx new file mode 100644 index 000000000000..12240d28d525 --- /dev/null +++ b/MagneticField/MagFieldElements/src/BFieldCacheZR.cxx @@ -0,0 +1,70 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#include "MagFieldElements/BFieldCacheZR.h" +#include <cmath> + +void +BFieldCacheZR::getB(const double* ATH_RESTRICT xyz, + double r, + double* ATH_RESTRICT B, + double* ATH_RESTRICT deriv) const +{ + const double x = xyz[0]; + const double y = xyz[1]; + const double z = xyz[2]; + // fractional position inside this bin + const double fz = (z - m_zmin) * m_invz; + const double gz = 1.0 - fz; + const double fr = (r - m_rmin) * m_invr; + const double gr = 1.0 - fr; + // interpolate field values in z, r + double Bzr[2]; + for (int i = 0; i < 2; ++i) { // z, r components + const double* field = m_field[i]; + Bzr[i] = gz * (gr * field[0] + fr * field[1]) + + fz * (gr * field[2] + fr * field[3]); + } + // convert (Bz,Br) to (Bx,By,Bz) + double invr; + if (r > 0.0) { + invr = 1.0 / r; + } else { + invr = 0.0; + } + const double c(x * invr); + const double s(y * invr); + B[0] = Bzr[1] * c; + B[1] = Bzr[1] * s; + B[2] = Bzr[0]; + + // compute field derivatives if requested + if (deriv) { + std::array<double, 2> dBdz; + std::array<double, 2> dBdr; + for (int j = 0; j < 2; ++j) { // Bz, Br components + const double* field = m_field[j]; + dBdz[j] = + m_invz * (gr * (field[2] - field[0]) + fr * (field[3] - field[1])); + dBdr[j] = + m_invr * (gz * (field[1] - field[0]) + fz * (field[3] - field[2])); + } + // convert to cartesian coordinates + const double cc = c * c; + const double cs = c * s; + const double ss = s * s; + const double sinvr = s * invr; + const double cinvr = c * invr; + deriv[0] = cc * dBdr[1] + sinvr * B[1]; + deriv[1] = cs * dBdr[1] - cinvr * B[1]; + deriv[2] = c * dBdz[1]; + deriv[3] = cs * dBdr[1] - sinvr * B[0]; + deriv[4] = ss * dBdr[1] + cinvr * B[0]; + deriv[5] = s * dBdz[1]; + deriv[6] = c * dBdr[0]; + deriv[7] = s * dBdr[0]; + deriv[8] = dBdz[0]; + } +} + diff --git a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx index d2b812ac28c2..7eec7f7db067 100644 --- a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx +++ b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx @@ -122,15 +122,6 @@ public: 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 }; -- GitLab