From 6b53d42384ba2994b2d282fd6dc3cd22b6642b7e Mon Sep 17 00:00:00 2001 From: christos <christos@cern.ch> Date: Wed, 9 Sep 2020 02:21:30 +0100 Subject: [PATCH 1/2] start moving methods to .icc --- .../MagFieldElements/BFieldMesh.h | 65 ++-- .../MagFieldElements/BFieldMesh.icc | 336 +++++++++--------- .../MagFieldElements/BFieldMeshZR.h | 68 ++-- .../MagFieldElements/BFieldMeshZR.icc | 86 +++++ .../MagFieldElements/BFieldZone.h | 15 +- .../MagFieldElements/BFieldZone.icc | 37 ++ 6 files changed, 356 insertions(+), 251 deletions(-) diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h index b6dc1f34ffce..b3b166f30fb4 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h @@ -39,37 +39,24 @@ public: double rmax, double phimin, double phimax, - double bscale) - : m_scale(bscale) - , m_nomScale(bscale) - { - m_min = { zmin, rmin, phimin }; - m_max = { zmax, rmax, phimax }; - } + double bscale); // set ranges void setRange(double zmin, double zmax, double rmin, double rmax, double phimin, - double phimax) - { - m_min = { zmin, rmin, phimin }; - m_max = { zmax, rmax, phimax }; - } + double phimax); // set bscale - void setBscale(double bscale) { m_scale = m_nomScale = bscale; } + void setBscale(double bscale); // scale bscale by a factor - void scaleBscale(double factor) { m_scale = factor * m_nomScale; } + void scaleBscale(double factor); // allocate space to vectors void reserve(int nz, int nr, int nphi, int nfield); - void reserve(int nz, int nr, int nphi) - { - reserve(nz, nr, nphi, nz * nr * nphi); - } + void reserve(int nz, int nr, int nphi); // add elements to vectors - void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); } - void appendField(const BFieldVector<T>& field) { m_field.push_back(field); } + void appendMesh(int i, double mesh); + void appendField(const BFieldVector<T>& field); // build Look Up Table void buildLUT(); // test if a point is inside this zone @@ -80,41 +67,35 @@ public: double phi, BFieldCache& cache, double scaleFactor = 1.0) const; - // get the B field - void getB(const double* ATH_RESTRICT xyz, - double* ATH_RESTRICT B, - double* ATH_RESTRICT deriv = nullptr) const; // accessors - double min(size_t i) const { return m_min[i]; } - double max(size_t i) const { return m_max[i]; } - double zmin() const { return m_min[0]; } - double zmax() const { return m_max[0]; } - double rmin() const { return m_min[1]; } - double rmax() const { return m_max[1]; } - double phimin() const { return m_min[2]; } - double phimax() const { return m_max[2]; } - unsigned nmesh(size_t i) const { return m_mesh[i].size(); } - double mesh(size_t i, size_t j) const { return m_mesh[i][j]; } - unsigned nfield() const { return m_field.size(); } - const BFieldVector<T>& field(size_t i) const { return m_field[i]; } - double bscale() const { return m_scale; } + double min(size_t i) const; + double max(size_t i) const; + double zmin() const; + double zmax() const; + double rmin() const; + double rmax() const; + double phimin() const; + double phimax() const; + unsigned nmesh(size_t i) const; + double mesh(size_t i, size_t j) const; + unsigned nfield() const; + const BFieldVector<T>& field(size_t i) const; + double bscale() const; int memSize() const; protected: std::array<double, 3> m_min; std::array<double, 3> m_max; - std::array<std::vector<double>,3> m_mesh; + std::array<std::vector<double>, 3> m_mesh; private: std::vector<BFieldVector<T>> m_field; double m_scale = 1.0; double m_nomScale; // nominal m_scale from the map - // look-up table and related variables - std::array<std::vector<int>,3> m_LUT; - std::array<double,3> m_invUnit; // inverse unit size in the LUT + 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 1abe3474d8d6..c88219a89c5d 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc @@ -2,9 +2,22 @@ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ -// -// Reserve space in the vectors to avoid unnecessary memory re-allocations. -// +// constructor +template<class T> +BFieldMesh<T>::BFieldMesh(double zmin, + double zmax, + double rmin, + double rmax, + double phimin, + double phimax, + double bscale) + : m_scale(bscale) + , m_nomScale(bscale) +{ + m_min = { zmin, rmin, phimin }; + m_max = { zmax, rmax, phimax }; +} + template<class T> void BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield) @@ -15,9 +28,6 @@ BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield) m_field.reserve(nfield); } -// -// Test if a point (z,r,phi) is inside this mesh region. -// template<class T> bool BFieldMesh<T>::inside(double z, double r, double phi) const @@ -32,9 +42,37 @@ BFieldMesh<T>::inside(double z, double r, double phi) const r >= rmin() && r <= rmax()); } -// -// Find and return the cache of the bin containing (z,r,phi) -// +template<class T> +void +BFieldMesh<T>::buildLUT() +{ + for (int j = 0; j < 3; ++j) { // z, r, phi + // align the m_mesh edges to m_min/m_max + m_mesh[j].front() = m_min[j]; + m_mesh[j].back() = m_max[j]; + // determine the unit size, q, to be used in the LUTs + const double width = m_mesh[j].back() - m_mesh[j].front(); + double q(width); + for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) { + q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]); + } + // find the number of units in the LUT + int n = int(width / q) + 1; + q = width / (n + 0.5); + m_invUnit[j] = 1.0 / q; // new unit size + ++n; + int m = 0; // mesh number + for (int i = 0; i < n; ++i) { // LUT index + if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) { + m++; + } + m_LUT[j].push_back(m); + } + } + m_roff = m_mesh[2].size(); // index offset for incrementing r by 1 + m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1 +} + template<class T> void BFieldMesh<T>::getCache(double z, @@ -70,8 +108,8 @@ 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 in cache const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner @@ -105,164 +143,143 @@ BFieldMesh<T>::getCache(double z, cache.setBscale(m_scale); } -// -// Compute the magnetic field (and the derivatives) without caching -// template<class T> void -BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz, - double* ATH_RESTRICT B, - double* ATH_RESTRICT deriv) const +BFieldMesh<T>::setRange(double zmin, + double zmax, + double rmin, + double rmax, + double phimin, + double phimax) { - // cylindrical coordinates - const double x = xyz[0]; - const double y = xyz[1]; - const double z = xyz[2]; - const double r = sqrt(x * x + y * y); - double phi = std::atan2(y, x); - if (phi < phimin()) { - phi += 2.0 * M_PI; - } - // is it inside this map? - if (!inside(z, r, phi)) { // no - B[0] = B[1] = B[2] = 0.0; - if (deriv) { - for (int i = 0; i < 9; ++i) { - deriv[i] = 0.0; - } - } - return; - } - // find the bin - // z - const std::vector<double>& mz(m_mesh[0]); - int iz = int((z - zmin()) * m_invUnit[0]); // index to LUT - iz = m_LUT[0][iz]; // tentative mesh index from LUT - if (z > mz[iz + 1]) { - ++iz; - } - // r - const std::vector<double>& mr(m_mesh[1]); - int ir = int((r - rmin()) * m_invUnit[1]); // index to LUT - ir = m_LUT[1][ir]; // tentative mesh index from LUT - if (r > mr[ir + 1]) { - ++ir; - } - // phi - const std::vector<double>& mphi(m_mesh[2]); - int iphi = int((phi - phimin()) * m_invUnit[2]); // index to LUT - iphi = m_LUT[2][iphi]; // tentative mesh index from LUT - if (phi > mphi[iphi + 1]) { - ++iphi; - } - // get the B field at the 8 corners - int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner - const 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] }; - // fractional position inside this mesh - const double fz = (z - mz[iz]) / (mz[iz + 1] - mz[iz]); - const double gz = 1.0 - fz; - const double fr = (r - mr[ir]) / (mr[ir + 1] - mr[ir]); - const double gr = 1.0 - fr; - const double fphi = (phi - mphi[iphi]) / (mphi[iphi + 1] - mphi[iphi]); - 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 - Bzrphi[i] = scale * (gz * (gr * (gphi * field[0][i] + fphi * field[1][i]) + - fr * (gphi * field[2][i] + fphi * field[3][i])) + - fz * (gr * (gphi * field[4][i] + fphi * field[5][i]) + - fr * (gphi * field[6][i] + fphi * field[7][i]))); - } - // convert (Bz,Br,Bphi) to (Bx,By,Bz) - const double c = (r > 0.0) ? x / r : cos(mphi[iphi]); - const double s = (r > 0.0) ? y / r : sin(mphi[iphi]); - B[0] = Bzrphi[1] * c - Bzrphi[2] * s; - B[1] = Bzrphi[1] * s + Bzrphi[2] * c; - B[2] = Bzrphi[0]; + m_min = { zmin, rmin, phimin }; + m_max = { zmax, rmax, phimax }; +} - // compute field derivatives if requested - if (deriv) { - const double sz = m_scale / (mz[iz + 1] - mz[iz]); - const double sr = m_scale / (mr[ir + 1] - mr[ir]); - const double sphi = m_scale / (mphi[iphi + 1] - mphi[iphi]); +template<class T> +void +BFieldMesh<T>::setBscale(double bscale) +{ + m_scale = m_nomScale = bscale; +} - 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 - dBdz[j] = sz * (gr * (gphi * (field[4][j] - field[0][j]) + - fphi * (field[5][j] - field[1][j])) + - fr * (gphi * (field[6][j] - field[2][j]) + - fphi * (field[7][j] - field[3][j]))); - dBdr[j] = sr * (gz * (gphi * (field[2][j] - field[0][j]) + - fphi * (field[3][j] - field[1][j])) + - fz * (gphi * (field[6][j] - field[4][j]) + - fphi * (field[7][j] - field[5][j]))); - dBdphi[j] = sphi * (gz * (gr * (field[1][j] - field[0][j]) + - fr * (field[3][j] - field[2][j])) + - fz * (gr * (field[5][j] - field[4][j]) + - fr * (field[7][j] - field[6][j]))); - } - // convert to cartesian coordinates - const double cc = c * c; - const double cs = c * s; - const double ss = s * s; - deriv[0] = cc * dBdr[1] - cs * dBdr[2] - cs * dBdphi[1] / r + - ss * dBdphi[2] / r + s * B[1] / r; - deriv[1] = cs * dBdr[1] - ss * dBdr[2] + cc * dBdphi[1] / r - - cs * dBdphi[2] / r - c * B[1] / r; - deriv[2] = c * dBdz[1] - s * dBdz[2]; - deriv[3] = cs * dBdr[1] + cc * dBdr[2] - ss * dBdphi[1] / r - - cs * dBdphi[2] / r - s * B[0] / r; - deriv[4] = ss * dBdr[1] + cs * dBdr[2] + cs * dBdphi[1] / r + - cc * dBdphi[2] / r + c * B[0] / r; - deriv[5] = s * dBdz[1] + c * dBdz[2]; - deriv[6] = c * dBdr[0] - s * dBdphi[0] / r; - deriv[7] = s * dBdr[1] + c * dBdphi[0] / r; - deriv[8] = dBdz[0]; - } +template<class T> +void +BFieldMesh<T>::scaleBscale(double factor) +{ + m_scale = factor * m_nomScale; } -// -// Construct the look-up table to accelerate bin-finding. -// template<class T> void -BFieldMesh<T>::buildLUT() +BFieldMesh<T>::reserve(int nz, int nr, int nphi) { - for (int j = 0; j < 3; ++j) { // z, r, phi - // align the m_mesh edges to m_min/m_max - m_mesh[j].front() = m_min[j]; - m_mesh[j].back() = m_max[j]; - // determine the unit size, q, to be used in the LUTs - const double width = m_mesh[j].back() - m_mesh[j].front(); - double q(width); - for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) { - q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]); - } - // find the number of units in the LUT - int n = int(width / q) + 1; - q = width / (n + 0.5); - m_invUnit[j] = 1.0 / q; // new unit size - ++n; - int m = 0; // mesh number - for (int i = 0; i < n; ++i) { // LUT index - if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) { - m++; - } - m_LUT[j].push_back(m); - } - } - m_roff = m_mesh[2].size(); // index offset for incrementing r by 1 - m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1 + reserve(nz, nr, nphi, nz * nr * nphi); +} + +template<class T> +void +BFieldMesh<T>::appendMesh(int i, double mesh) +{ + m_mesh[i].push_back(mesh); +} + +template<class T> +void +BFieldMesh<T>::appendField(const BFieldVector<T>& field) +{ + m_field.push_back(field); +} + +template<class T> +double +BFieldMesh<T>::min(size_t i) const +{ + return m_min[i]; +} + +template<class T> +double +BFieldMesh<T>::max(size_t i) const +{ + return m_max[i]; +} + +template<class T> +double +BFieldMesh<T>::zmin() const +{ + return m_min[0]; +} + +template<class T> +double +BFieldMesh<T>::zmax() const +{ + return m_max[0]; +} + +template<class T> +double +BFieldMesh<T>::rmin() const +{ + return m_min[1]; +} + +template<class T> +double +BFieldMesh<T>::rmax() const +{ + return m_max[1]; +} + +template<class T> +double +BFieldMesh<T>::phimin() const +{ + return m_min[2]; +} + +template<class T> +double +BFieldMesh<T>::phimax() const +{ + return m_max[2]; +} + +template<class T> +unsigned +BFieldMesh<T>::nmesh(size_t i) const +{ + return m_mesh[i].size(); +} + +template<class T> +double +BFieldMesh<T>::mesh(size_t i, size_t j) const +{ + return m_mesh[i][j]; +} + +template<class T> +unsigned +BFieldMesh<T>::nfield() const +{ + return m_field.size(); +} + +template<class T> +const BFieldVector<T>& +BFieldMesh<T>::field(size_t i) const +{ + return m_field[i]; +} + +template<class T> +double +BFieldMesh<T>::bscale() const +{ + return m_scale; } template<class T> @@ -280,4 +297,3 @@ BFieldMesh<T>::memSize() const return size; } - diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h index b79841a62b3f..b82c6bd275b0 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h @@ -2,73 +2,63 @@ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ -// -// BFieldMeshZR.h -// -// A 2-dim z-r mesh inside the solenoid field map -// -// Masahiro Morii, Harvard University -// +/* + * BFieldMeshZR.h + * + * A 2-dim z-r mesh inside the solenoid field map + * + * Masahiro Morii, Harvard University + * + * AthenaMT : RD Schaffer , Christos Anastopoulos + */ #ifndef BFIELDMESHZR_H #define BFIELDMESHZR_H #include "MagFieldElements/BFieldCacheZR.h" #include "MagFieldElements/BFieldVectorZR.h" +#include <array> #include <cmath> #include <vector> -#include <array> class BFieldMeshZR { public: // constructor - BFieldMeshZR(double zmin, double zmax, double rmin, double rmax) - { - m_min = { zmin, rmin }; - m_max = { zmax, rmax }; - } + BFieldMeshZR(double zmin, double zmax, double rmin, double rmax); // allocate space to vectors - void reserve(int nz, int nr) - { - m_mesh[0].reserve(nz); - m_mesh[1].reserve(nr); - m_field.reserve(nz * nr); - } + void reserve(int nz, int nr); // add elements to vectors - void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); } - void appendField(const BFieldVectorZR& field) { m_field.push_back(field); } + void appendMesh(int i, double mesh); + void appendField(const BFieldVectorZR& field); // build LUT void buildLUT(); // test if a point is inside this zone - inline bool inside(double z, double r) const - { - return (z >= zmin() && z <= zmax() && r >= rmin() && r <= rmax()); - } + bool inside(double z, double r) const; // find the bin inline void getCache(double z, double r, BFieldCacheZR& cache, double scaleFactor = 1.0) const; // accessors - double min(size_t i) const { return m_min[i]; } - double max(size_t i) const { return m_max[i]; } - double zmin() const { return m_min[0]; } - double zmax() const { return m_max[0]; } - double rmin() const { return m_min[1]; } - double rmax() const { return m_max[1]; } - unsigned nmesh(size_t i) const { return m_mesh[i].size(); } - double mesh(size_t i, size_t j) const { return m_mesh[i][j]; } - unsigned nfield() const { return m_field.size(); } - const BFieldVectorZR& field(size_t i) const { return m_field[i]; } + double min(size_t i); + double max(size_t i); + double zmin() const; + double zmax() const; + double rmin() const; + double rmax() const; + unsigned nmesh(size_t i) const; + double mesh(size_t i, size_t j) const; + unsigned nfield() const; + const BFieldVectorZR& field(size_t i) const; int memSize() const; private: - std::array<double,2> m_min; - std::array<double,2> m_max; - std::array<std::vector<double>,2> m_mesh; + std::array<double, 2> m_min; + std::array<double, 2> m_max; + std::array<std::vector<double>, 2> m_mesh; std::vector<BFieldVectorZR> m_field; // look-up table and related variables - std::array<std::vector<int>,2> m_LUT; + std::array<std::vector<int>, 2> m_LUT; std::array<double, 2> m_invUnit; // inverse unit size in the LUT int m_zoff; }; diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc index 781dbf7dc885..1b76cf8b1dee 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.icc @@ -36,3 +36,89 @@ BFieldMeshZR::getCache(double z, cache.setField(3, m_field[im0 + m_zoff + 1], scaleFactor); } +inline BFieldMeshZR::BFieldMeshZR(double zmin, + double zmax, + double rmin, + double rmax) +{ + m_min = { zmin, rmin }; + m_max = { zmax, rmax }; +} + +inline void +BFieldMeshZR::reserve(int nz, int nr) +{ + m_mesh[0].reserve(nz); + m_mesh[1].reserve(nr); + m_field.reserve(nz * nr); +} +// add elements to vectors +inline void +BFieldMeshZR::appendMesh(int i, double mesh) +{ + m_mesh[i].push_back(mesh); +} +inline void +BFieldMeshZR::appendField(const BFieldVectorZR& field) +{ + m_field.push_back(field); +} + +// test if a point is inside this zone +inline bool +BFieldMeshZR::inside(double z, double r) const +{ + return (z >= zmin() && z <= zmax() && r >= rmin() && r <= rmax()); +} + +inline double +BFieldMeshZR::min(size_t i) const +{ + return m_min[i]; +} +inline double +BFieldMeshZR::max(size_t i) const +{ + return m_max[i]; +} +inline double +BFieldMeshZR::zmin() const +{ + return m_min[0]; +} +inline double +BFieldMeshZR::zmax() const +{ + return m_max[0]; +} +inline double +BFieldMeshZR::rmin() const +{ + return m_min[1]; +} +inline double +BFieldMeshZR::rmax() const +{ + return m_max[1]; +} +inline unsigned +BFieldMeshZR::nmesh(size_t i) const +{ + return m_mesh[i].size(); +} +inline double +BFieldMeshZR::mesh(size_t i, size_t j) const +{ + return m_mesh[i][j]; +} +inline unsigned +BFieldMeshZR::nfield() const +{ + return m_field.size(); +} +inline const BFieldVectorZR& +BFieldMeshZR::field(size_t i) const +{ + return m_field[i]; +} + diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h index 57fbb3ea6233..888e3f8d059a 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.h @@ -29,12 +29,7 @@ public: double rmax, double phimin, double phimax, - double scale) - : BFieldMesh<short>(zmin, zmax, rmin, rmax, phimin, phimax, scale) - , m_id(id) - { - ; - } + double scale); // add elements to vectors void appendCond(const BFieldCond& cond); // compute Biot-Savart magnetic field and add to B[3] @@ -45,10 +40,10 @@ public: // Scaling is done in cachec void scaleField(double factor); // accessors - int id() const { return m_id; } - unsigned ncond() const { return m_cond.size(); } - const BFieldCond& cond(int i) const { return m_cond[i]; } - const std::vector<BFieldCond>* condVector() const { return &m_cond; } + int id() const; + unsigned ncond() const; + const BFieldCond& cond(int i) const; + const std::vector<BFieldCond>* condVector() const; int memSize() const; // adjust the min/max edges to a new value void adjustMin(int i, double x); diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc index 071400e3729d..ddb5022730ae 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldZone.icc @@ -1,6 +1,19 @@ /* Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ + +inline BFieldZone::BFieldZone(int id, + double zmin, + double zmax, + double rmin, + double rmax, + double phimin, + double phimax, + double scale) + : BFieldMesh<short>(zmin, zmax, rmin, rmax, phimin, phimax, scale) + , m_id(id) +{} + inline void BFieldZone::appendCond(const BFieldCond& cond) { @@ -26,6 +39,30 @@ BFieldZone::addBiotSavart(const double* ATH_RESTRICT xyz, } } +inline int +BFieldZone::id() const +{ + return m_id; +} + +inline unsigned +BFieldZone::ncond() const +{ + return m_cond.size(); +} + +inline const BFieldCond& +BFieldZone::cond(int i) const +{ + return m_cond[i]; +} + +inline const std::vector<BFieldCond>* +BFieldZone::condVector() const +{ + return &m_cond; +} + inline int BFieldZone::memSize() const { -- GitLab From b25c36af47430679346bb5b19f7e63125b7c51a2 Mon Sep 17 00:00:00 2001 From: christos <christos@cern.ch> Date: Wed, 9 Sep 2020 02:45:37 +0100 Subject: [PATCH 2/2] start moving methods to .icc --- .../MagFieldElements/AtlasFieldCache.h | 1 - .../MagFieldElements/BFieldMesh.h | 65 ++-- .../MagFieldElements/BFieldMesh.icc | 345 +++++++++--------- .../MagFieldElements/BFieldMeshZR.h | 4 +- .../MagFieldElements/src/AtlasFieldCache.cxx | 5 +- .../MagFieldElements/src/AtlasFieldMap.cxx | 16 - 6 files changed, 213 insertions(+), 223 deletions(-) diff --git a/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h index 7161f00aa43e..5878c8b06909 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h +++ b/MagneticField/MagFieldElements/MagFieldElements/AtlasFieldCache.h @@ -23,7 +23,6 @@ #include "MagFieldElements/BFieldMeshZR.h" #include "MagFieldElements/BFieldZone.h" -#include <iostream> #include <memory> namespace MagField { diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h index b3b166f30fb4..b6dc1f34ffce 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h @@ -39,24 +39,37 @@ public: double rmax, double phimin, double phimax, - double bscale); + double bscale) + : m_scale(bscale) + , m_nomScale(bscale) + { + m_min = { zmin, rmin, phimin }; + m_max = { zmax, rmax, phimax }; + } // set ranges void setRange(double zmin, double zmax, double rmin, double rmax, double phimin, - double phimax); + double phimax) + { + m_min = { zmin, rmin, phimin }; + m_max = { zmax, rmax, phimax }; + } // set bscale - void setBscale(double bscale); + void setBscale(double bscale) { m_scale = m_nomScale = bscale; } // scale bscale by a factor - void scaleBscale(double factor); + void scaleBscale(double factor) { m_scale = factor * m_nomScale; } // allocate space to vectors void reserve(int nz, int nr, int nphi, int nfield); - void reserve(int nz, int nr, int nphi); + void reserve(int nz, int nr, int nphi) + { + reserve(nz, nr, nphi, nz * nr * nphi); + } // add elements to vectors - void appendMesh(int i, double mesh); - void appendField(const BFieldVector<T>& field); + void appendMesh(int i, double mesh) { m_mesh[i].push_back(mesh); } + void appendField(const BFieldVector<T>& field) { m_field.push_back(field); } // build Look Up Table void buildLUT(); // test if a point is inside this zone @@ -67,35 +80,41 @@ public: double phi, BFieldCache& cache, double scaleFactor = 1.0) const; + // get the B field + void getB(const double* ATH_RESTRICT xyz, + double* ATH_RESTRICT B, + double* ATH_RESTRICT deriv = nullptr) const; // accessors - double min(size_t i) const; - double max(size_t i) const; - double zmin() const; - double zmax() const; - double rmin() const; - double rmax() const; - double phimin() const; - double phimax() const; - unsigned nmesh(size_t i) const; - double mesh(size_t i, size_t j) const; - unsigned nfield() const; - const BFieldVector<T>& field(size_t i) const; - double bscale() const; + double min(size_t i) const { return m_min[i]; } + double max(size_t i) const { return m_max[i]; } + double zmin() const { return m_min[0]; } + double zmax() const { return m_max[0]; } + double rmin() const { return m_min[1]; } + double rmax() const { return m_max[1]; } + double phimin() const { return m_min[2]; } + double phimax() const { return m_max[2]; } + unsigned nmesh(size_t i) const { return m_mesh[i].size(); } + double mesh(size_t i, size_t j) const { return m_mesh[i][j]; } + unsigned nfield() const { return m_field.size(); } + const BFieldVector<T>& field(size_t i) const { return m_field[i]; } + double bscale() const { return m_scale; } int memSize() const; protected: std::array<double, 3> m_min; std::array<double, 3> m_max; - std::array<std::vector<double>, 3> m_mesh; + std::array<std::vector<double>,3> m_mesh; private: std::vector<BFieldVector<T>> m_field; double m_scale = 1.0; double m_nomScale; // nominal m_scale from the map + // look-up table and related variables - std::array<std::vector<int>, 3> m_LUT; - std::array<double, 3> m_invUnit; // inverse unit size in the LUT + 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 7b72c0939438..b43cfc406b05 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc @@ -3,21 +3,22 @@ */ #include "CxxUtils/vec.h" +// +// Reserve space in the vectors to avoid unnecessary memory re-allocations. +// template<class T> -BFieldMesh<T>::BFieldMesh(double zmin, - double zmax, - double rmin, - double rmax, - double phimin, - double phimax, - double bscale) - : m_scale(bscale) - , m_nomScale(bscale) +void +BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield) { - m_min = { zmin, rmin, phimin }; - m_max = { zmax, rmax, phimax }; + m_mesh[0].reserve(nz); + m_mesh[1].reserve(nr); + m_mesh[2].reserve(nphi); + m_field.reserve(nfield); } +// +// Test if a point (z,r,phi) is inside this mesh region. +// template<class T> bool BFieldMesh<T>::inside(double z, double r, double phi) const @@ -32,37 +33,9 @@ BFieldMesh<T>::inside(double z, double r, double phi) const r >= rmin() && r <= rmax()); } -template<class T> -void -BFieldMesh<T>::buildLUT() -{ - for (int j = 0; j < 3; ++j) { // z, r, phi - // align the m_mesh edges to m_min/m_max - m_mesh[j].front() = m_min[j]; - m_mesh[j].back() = m_max[j]; - // determine the unit size, q, to be used in the LUTs - const double width = m_mesh[j].back() - m_mesh[j].front(); - double q(width); - for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) { - q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]); - } - // find the number of units in the LUT - int n = int(width / q) + 1; - q = width / (n + 0.5); - m_invUnit[j] = 1.0 / q; // new unit size - ++n; - int m = 0; // mesh number - for (int i = 0; i < n; ++i) { // LUT index - if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) { - m++; - } - m_LUT[j].push_back(m); - } - } - m_roff = m_mesh[2].size(); // index offset for incrementing r by 1 - m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1 -} - +// +// Find and return the cache of the bin containing (z,r,phi) +// template<class T> void BFieldMesh<T>::getCache(double z, @@ -100,7 +73,6 @@ BFieldMesh<T>::getCache(double z, // store the bin edges 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 in cache const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner const double sf = scaleFactor; @@ -145,153 +117,166 @@ BFieldMesh<T>::getCache(double z, cache.setBscale(m_scale); } +// +// Compute the magnetic field (and the derivatives) without caching +// template<class T> void -BFieldMesh<T>::setRange(double zmin, - double zmax, - double rmin, - double rmax, - double phimin, - double phimax) -{ - m_min = { zmin, rmin, phimin }; - m_max = { zmax, rmax, phimax }; -} - -template<class T> -void -BFieldMesh<T>::setBscale(double bscale) -{ - m_scale = m_nomScale = bscale; -} - -template<class T> -void -BFieldMesh<T>::scaleBscale(double factor) +BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz, + double* ATH_RESTRICT B, + double* ATH_RESTRICT deriv) const { - m_scale = factor * m_nomScale; -} + // cylindrical coordinates + const double x = xyz[0]; + const double y = xyz[1]; + const double z = xyz[2]; + const double r = sqrt(x * x + y * y); + double phi = std::atan2(y, x); + if (phi < phimin()) { + phi += 2.0 * M_PI; + } + // is it inside this map? + if (!inside(z, r, phi)) { // no + B[0] = B[1] = B[2] = 0.0; + if (deriv) { + for (int i = 0; i < 9; ++i) { + deriv[i] = 0.0; + } + } + return; + } + // find the bin + // z + const std::vector<double>& mz(m_mesh[0]); + int iz = int((z - zmin()) * m_invUnit[0]); // index to LUT + iz = m_LUT[0][iz]; // tentative mesh index from LUT + if (z > mz[iz + 1]) { + ++iz; + } + // r + const std::vector<double>& mr(m_mesh[1]); + int ir = int((r - rmin()) * m_invUnit[1]); // index to LUT + ir = m_LUT[1][ir]; // tentative mesh index from LUT + if (r > mr[ir + 1]) { + ++ir; + } + // phi + const std::vector<double>& mphi(m_mesh[2]); + int iphi = int((phi - phimin()) * m_invUnit[2]); // index to LUT + iphi = m_LUT[2][iphi]; // tentative mesh index from LUT + if (phi > mphi[iphi + 1]) { + ++iphi; + } + // get the B field at the 8 corners + int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner + const 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] + }; + // fractional position inside this mesh + const double fz = (z - mz[iz]) / (mz[iz + 1] - mz[iz]); + const double gz = 1.0 - fz; + const double fr = (r - mr[ir]) / (mr[ir + 1] - mr[ir]); + const double gr = 1.0 - fr; + const double fphi = (phi - mphi[iphi]) / (mphi[iphi + 1] - mphi[iphi]); + 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 + Bzrphi[i] = scale * (gz * (gr * (gphi * field[0][i] + fphi * field[1][i]) + + fr * (gphi * field[2][i] + fphi * field[3][i])) + + fz * (gr * (gphi * field[4][i] + fphi * field[5][i]) + + fr * (gphi * field[6][i] + fphi * field[7][i]))); + } + // convert (Bz,Br,Bphi) to (Bx,By,Bz) + const double c = (r > 0.0) ? x / r : cos(mphi[iphi]); + const double s = (r > 0.0) ? y / r : sin(mphi[iphi]); + B[0] = Bzrphi[1] * c - Bzrphi[2] * s; + B[1] = Bzrphi[1] * s + Bzrphi[2] * c; + B[2] = Bzrphi[0]; -template<class T> -void -BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield) -{ - m_mesh[0].reserve(nz); - m_mesh[1].reserve(nr); - m_mesh[2].reserve(nphi); - m_field.reserve(nfield); -} + // compute field derivatives if requested + if (deriv) { + const double sz = m_scale / (mz[iz + 1] - mz[iz]); + const double sr = m_scale / (mr[ir + 1] - mr[ir]); + const double sphi = m_scale / (mphi[iphi + 1] - mphi[iphi]); -template<class T> -void -BFieldMesh<T>::reserve(int nz, int nr, int nphi) -{ - reserve(nz, nr, nphi, nz * nr * nphi); -} - -template<class T> -void -BFieldMesh<T>::appendMesh(int i, double mesh) -{ - m_mesh[i].push_back(mesh); + 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 + dBdz[j] = sz * (gr * (gphi * (field[4][j] - field[0][j]) + + fphi * (field[5][j] - field[1][j])) + + fr * (gphi * (field[6][j] - field[2][j]) + + fphi * (field[7][j] - field[3][j]))); + dBdr[j] = sr * (gz * (gphi * (field[2][j] - field[0][j]) + + fphi * (field[3][j] - field[1][j])) + + fz * (gphi * (field[6][j] - field[4][j]) + + fphi * (field[7][j] - field[5][j]))); + dBdphi[j] = sphi * (gz * (gr * (field[1][j] - field[0][j]) + + fr * (field[3][j] - field[2][j])) + + fz * (gr * (field[5][j] - field[4][j]) + + fr * (field[7][j] - field[6][j]))); + } + // convert to cartesian coordinates + const double cc = c * c; + const double cs = c * s; + const double ss = s * s; + deriv[0] = cc * dBdr[1] - cs * dBdr[2] - cs * dBdphi[1] / r + + ss * dBdphi[2] / r + s * B[1] / r; + deriv[1] = cs * dBdr[1] - ss * dBdr[2] + cc * dBdphi[1] / r - + cs * dBdphi[2] / r - c * B[1] / r; + deriv[2] = c * dBdz[1] - s * dBdz[2]; + deriv[3] = cs * dBdr[1] + cc * dBdr[2] - ss * dBdphi[1] / r - + cs * dBdphi[2] / r - s * B[0] / r; + deriv[4] = ss * dBdr[1] + cs * dBdr[2] + cs * dBdphi[1] / r + + cc * dBdphi[2] / r + c * B[0] / r; + deriv[5] = s * dBdz[1] + c * dBdz[2]; + deriv[6] = c * dBdr[0] - s * dBdphi[0] / r; + deriv[7] = s * dBdr[1] + c * dBdphi[0] / r; + deriv[8] = dBdz[0]; + } } +// +// Construct the look-up table to accelerate bin-finding. +// template<class T> void -BFieldMesh<T>::appendField(const BFieldVector<T>& field) -{ - m_field.push_back(field); -} - -template<class T> -double -BFieldMesh<T>::min(size_t i) const -{ - return m_min[i]; -} - -template<class T> -double -BFieldMesh<T>::max(size_t i) const -{ - return m_max[i]; -} - -template<class T> -double -BFieldMesh<T>::zmin() const -{ - return m_min[0]; -} - -template<class T> -double -BFieldMesh<T>::zmax() const -{ - return m_max[0]; -} - -template<class T> -double -BFieldMesh<T>::rmin() const -{ - return m_min[1]; -} - -template<class T> -double -BFieldMesh<T>::rmax() const -{ - return m_max[1]; -} - -template<class T> -double -BFieldMesh<T>::phimin() const -{ - return m_min[2]; -} - -template<class T> -double -BFieldMesh<T>::phimax() const -{ - return m_max[2]; -} - -template<class T> -unsigned -BFieldMesh<T>::nmesh(size_t i) const -{ - return m_mesh[i].size(); -} - -template<class T> -double -BFieldMesh<T>::mesh(size_t i, size_t j) const -{ - return m_mesh[i][j]; -} - -template<class T> -unsigned -BFieldMesh<T>::nfield() const -{ - return m_field.size(); -} - -template<class T> -const BFieldVector<T>& -BFieldMesh<T>::field(size_t i) const -{ - return m_field[i]; -} - -template<class T> -double -BFieldMesh<T>::bscale() const +BFieldMesh<T>::buildLUT() { - return m_scale; + for (int j = 0; j < 3; ++j) { // z, r, phi + // align the m_mesh edges to m_min/m_max + m_mesh[j].front() = m_min[j]; + m_mesh[j].back() = m_max[j]; + // determine the unit size, q, to be used in the LUTs + const double width = m_mesh[j].back() - m_mesh[j].front(); + double q(width); + for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) { + q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]); + } + // find the number of units in the LUT + int n = int(width / q) + 1; + q = width / (n + 0.5); + m_invUnit[j] = 1.0 / q; // new unit size + ++n; + int m = 0; // mesh number + for (int i = 0; i < n; ++i) { // LUT index + if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) { + m++; + } + m_LUT[j].push_back(m); + } + } + m_roff = m_mesh[2].size(); // index offset for incrementing r by 1 + m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1 } template<class T> diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h index b82c6bd275b0..1d78885d5312 100644 --- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h +++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMeshZR.h @@ -40,8 +40,8 @@ public: BFieldCacheZR& cache, double scaleFactor = 1.0) const; // accessors - double min(size_t i); - double max(size_t i); + double min(size_t i) const; + double max(size_t i) const; double zmin() const; double zmax() const; double rmin() const; diff --git a/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx b/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx index a3927c4553d7..d9e66122f27f 100644 --- a/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx +++ b/MagneticField/MagFieldElements/src/AtlasFieldCache.cxx @@ -11,6 +11,9 @@ // #include "MagFieldElements/AtlasFieldCache.h" +#include <cmath> +#include <iostream> + void MagField::AtlasFieldCache::getField(const double* ATH_RESTRICT xyz, double* ATH_RESTRICT bxyz, @@ -89,7 +92,7 @@ MagField::AtlasFieldCache::getFieldZR(const double* ATH_RESTRICT xyz, const double x = xyz[0]; const double y = xyz[1]; const double z = xyz[2]; - const double r = sqrt(x * x + y * y); + const double r = std::sqrt(x * x + y * y); // test if the cache was initialized and the ZR cache is valid for current // position diff --git a/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx b/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx index 1956e5f35d9b..d1f928528cd0 100644 --- a/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx +++ b/MagneticField/MagFieldElements/src/AtlasFieldMap.cxx @@ -89,46 +89,30 @@ MagField::AtlasFieldMap::initializeMap(TFile* rootfile, // } int id; double zmin; - double zmax; - double rmin; - double rmax; - double phimin; - double phimax; double bscale; int ncond; bool* finite; double* p1x; - double* p1y; - double* p1z; - double* p2x; - double* p2y; - double* p2z; double* curr; int nmeshz; - int nmeshr; - int nmeshphi; double* meshz; - double* meshr; - double* meshphi; int nfield; short* fieldz; - short* fieldr; - short* fieldphi; // define the fixed-sized branches first tree->SetBranchAddress("id", &id); -- GitLab