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] 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 7161f00aa43..5878c8b0690 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 b3b166f30fb..b6dc1f34ffc 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 7b72c093943..b43cfc406b0 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 b82c6bd275b..1d78885d531 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 a3927c4553d..d9e66122f27 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 1956e5f35d9..d1f928528cd 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