Skip to content
Snippets Groups Projects
Commit a8b79151 authored by Christos Anastopoulos's avatar Christos Anastopoulos
Browse files

MagFieldField Elements,use CxxUtils::vec notation also in the getCache

parent 4edc4b5c
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
+ 9
− 6
View file @ a8b79151
......@@ -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
};
......
MagneticField/MagFieldElements/MagFieldElements/BFieldCache.icc
+ 7
− 2
View file @ a8b79151
/*
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
......
MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+ 53
− 40
View file @ a8b79151
......@@ -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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment