File indexing completed on 2023-03-17 11:14:35
0001 #ifndef ParametrizedEngine_BCyl_h
0002 #define ParametrizedEngine_BCyl_h
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <cmath>
0020
0021 #include "DataFormats/Math/interface/approx_exp.h"
0022
0023 namespace magfieldparam {
0024 template <typename T>
0025 struct BCylParam {
0026
0027 template <typename... Args>
0028 constexpr BCylParam(Args... init)
0029 : prm{std::forward<Args>(init)...},
0030
0031 ap2(4 * prm[0] * prm[0] / (prm[1] * prm[1])),
0032 hb0(0.5 * prm[2] * std::sqrt(1.0 + ap2)),
0033 hlova(1 / std::sqrt(ap2)),
0034 ainv(2 * hlova / prm[1]),
0035 coeff(1 / (prm[8] * prm[8])) {}
0036
0037 T prm[9];
0038 T ap2, hb0, hlova, ainv, coeff;
0039 };
0040
0041 namespace bcylDetails {
0042
0043 template <typename T>
0044 inline void ffunkti(T u, T* __restrict__ ff) __attribute__((always_inline));
0045
0046 template <typename T>
0047 inline void ffunkti(T u, T* __restrict__ ff) {
0048
0049 T a, b, a2, u2;
0050 u2 = u * u;
0051 a = T(1) / (T(1) + u2);
0052 a2 = -T(3) * a * a;
0053 b = std::sqrt(a);
0054 ff[0] = u * b;
0055 ff[1] = a * b;
0056 ff[2] = a2 * ff[0];
0057 ff[3] = a2 * ff[1] * (T(1) - 4 * u2);
0058 }
0059
0060 inline double myExp(double x) { return std::exp(x); }
0061 inline float myExp(float x) { return unsafe_expf<3>(x); }
0062
0063 }
0064
0065 template <typename T>
0066 class BCycl {
0067 public:
0068 BCycl(BCylParam<T> const& ipar) : pars(ipar) {}
0069
0070 void operator()(T r2, T z, T& Br, T& Bz) const { compute(r2, z, Br, Bz); }
0071
0072
0073 void compute(T r2, T z, T& Br, T& Bz) const {
0074 using namespace bcylDetails;
0075
0076 z -= pars.prm[3];
0077 T az = std::abs(z);
0078 T zainv = z * pars.ainv;
0079 T u = pars.hlova - zainv;
0080 T v = pars.hlova + zainv;
0081 T fu[4], gv[4];
0082 ffunkti(u, fu);
0083 ffunkti(v, gv);
0084 T rat = T(0.5) * pars.ainv;
0085 T rat2 = rat * rat * r2;
0086 Br = pars.hb0 * rat * (fu[1] - gv[1] - (fu[3] - gv[3]) * rat2 * T(0.5));
0087 Bz = pars.hb0 * (fu[0] + gv[0] - (fu[2] + gv[2]) * rat2);
0088
0089 T corBr = pars.prm[4] * z * (az - pars.prm[5]) * (az - pars.prm[5]);
0090 T corBz = -pars.prm[6] * (myExp(-(z - pars.prm[7]) * (z - pars.prm[7]) * pars.coeff) +
0091 myExp(-(z + pars.prm[7]) * (z + pars.prm[7]) * pars.coeff));
0092 Br += corBr;
0093 Bz += corBz;
0094 }
0095
0096 private:
0097 BCylParam<T> pars;
0098 };
0099 }
0100
0101 #endif