diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
index f18209c09a259b22ff5419509bc3f93a5bd2ce36..6675121ca1879628975b9250517ce4f53442e3e6 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
@@ -18,6 +18,7 @@
#define BFIELDCACHE_H
#include "CxxUtils/restrict.h"
+#include "CxxUtils/vec.h"
#include "MagFieldElements/BFieldVector.h"
class BFieldCache
@@ -36,8 +37,10 @@ public:
double phimin,
double phimax);
- // set field array, filled externally
- void setField(double field[][8]);
+ // set the 3x8 field array.
+ void setField(const CxxUtils::vec<double, 8>& field1,
+ const CxxUtils::vec<double, 8>& field2,
+ const CxxUtils::vec<double, 8>& field3);
// set the multiplicative factor for the field vectors
void setBscale(double bscale);
@@ -54,7 +57,7 @@ public:
double* ATH_RESTRICT deriv = nullptr) const;
private:
- // bin range in z
+ // bin range in z
double m_zmin = 0.0;
double m_zmax = 0.0;
// bin range in r
@@ -62,12 +65,12 @@ private:
double m_rmax = 0.0;
// bin range in phi
double m_phimin = 0.0;
- double m_phimax = -1.0;
+ double m_phimax = -1.0;
// 1/(bin size) in z, r, phi
double m_invz;
double m_invr;
- double m_invphi;
- double m_scale; // unit of m_field in kT
+ double m_invphi;
+ double m_scale; // unit of m_field in kT
alignas(16) double m_field[3][8]; // (Bz,Br,Bphi) at 8 corners of the bin
};
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
index d9528098145ac978a94c4156e11ae0a231243f16..20dadc075fe577df1f358f806e721d67cedd7ff0 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
@@ -1,6 +1,7 @@
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
+#include "CxxUtils/vec.h"
#include <cmath>
inline void
BFieldCache::invalidate()
@@ -29,9 +30,13 @@ BFieldCache::setRange(double zmin,
// set field array, filled externally
inline void
-BFieldCache::setField(double field[][8])
+BFieldCache::setField(const CxxUtils::vec<double, 8>& field1,
+ const CxxUtils::vec<double, 8>& field2,
+ const CxxUtils::vec<double, 8>& field3)
{
- std::copy(&field[0][0], &field[0][0] + 3 * 8, &m_field[0][0]);
+ CxxUtils::vstore(&m_field[0][0], field1);
+ CxxUtils::vstore(&m_field[1][0], field2);
+ CxxUtils::vstore(&m_field[2][0], field3);
}
inline void
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
index 1abe3474d8d6176af368125686d235dd71f3c59f..b43cfc406b0513d8e4bb90dca6ede70cced3f0de 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
@@ -2,6 +2,7 @@
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
+#include "CxxUtils/vec.h"
//
// Reserve space in the vectors to avoid unnecessary memory re-allocations.
//
@@ -70,37 +71,48 @@ 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
const double sf = scaleFactor;
- double field1[3][8] = { { sf * m_field[im0][0],
- sf * m_field[im0 + m_roff][0],
- sf * m_field[im0 + m_zoff][0],
- sf * m_field[im0 + m_zoff + m_roff][0],
- sf * m_field[im0 + 1][0],
- sf * m_field[im0 + m_roff + 1][0],
- sf * m_field[im0 + m_zoff + 1][0],
- sf * m_field[im0 + m_zoff + m_roff + 1][0] },
- { sf * m_field[im0][1],
- sf * m_field[im0 + m_roff][1],
- sf * m_field[im0 + m_zoff][1],
- sf * m_field[im0 + m_zoff + m_roff][1],
- sf * m_field[im0 + 1][1],
- sf * m_field[im0 + m_roff + 1][1],
- sf * m_field[im0 + m_zoff + 1][1],
- sf * m_field[im0 + m_zoff + m_roff + 1][1] },
- { sf * m_field[im0][2],
- sf * m_field[im0 + m_roff][2],
- sf * m_field[im0 + m_zoff][2],
- sf * m_field[im0 + m_zoff + m_roff][2],
- sf * m_field[im0 + 1][2],
- sf * m_field[im0 + m_roff + 1][2],
- sf * m_field[im0 + m_zoff + 1][2],
- sf * m_field[im0 + m_zoff + m_roff + 1][2] } };
- cache.setField(field1);
+ /*
+ * The typical case is that we convert from short to double
+ * (explicitly or implicitly) in order to multiply with the double sf
+ */
+ CxxUtils::vec<double, 8> field1 = {
+ static_cast<double>(m_field[im0][0]),
+ static_cast<double>(m_field[im0 + m_roff][0]),
+ static_cast<double>(m_field[im0 + m_zoff][0]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff][0]),
+ static_cast<double>(m_field[im0 + 1][0]),
+ static_cast<double>(m_field[im0 + m_roff + 1][0]),
+ static_cast<double>(m_field[im0 + m_zoff + 1][0]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff + 1][0])
+ };
+
+ CxxUtils::vec<double, 8> field2 = {
+ static_cast<double>(m_field[im0][1]),
+ static_cast<double>(m_field[im0 + m_roff][1]),
+ static_cast<double>(m_field[im0 + m_zoff][1]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff][1]),
+ static_cast<double>(m_field[im0 + 1][1]),
+ static_cast<double>(m_field[im0 + m_roff + 1][1]),
+ static_cast<double>(m_field[im0 + m_zoff + 1][1]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff + 1][1])
+ };
+
+ CxxUtils::vec<double, 8> field3 = {
+ static_cast<double>(m_field[im0][2]),
+ static_cast<double>(m_field[im0 + m_roff][2]),
+ static_cast<double>(m_field[im0 + m_zoff][2]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff][2]),
+ static_cast<double>(m_field[im0 + 1][2]),
+ static_cast<double>(m_field[im0 + m_roff + 1][2]),
+ static_cast<double>(m_field[im0 + m_zoff + 1][2]),
+ static_cast<double>(m_field[im0 + m_zoff + m_roff + 1][2])
+ };
+ cache.setField(sf * field1, sf * field2, sf * field3);
// store the B scale
cache.setBscale(m_scale);
}
@@ -157,14 +169,16 @@ BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz,
}
// 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] };
+ 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;
@@ -194,9 +208,9 @@ BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz,
const double sr = m_scale / (mr[ir + 1] - mr[ir]);
const double sphi = m_scale / (mphi[iphi + 1] - mphi[iphi]);
- std::array<double,3> dBdz;
- std::array<double,3> dBdr;
- std::array<double,3> dBdphi;
+ 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])) +
@@ -280,4 +294,3 @@ BFieldMesh<T>::memSize() const
return size;
}
-