ParametrizedEngine/src/BCyl.h

0001 #ifndef ParametrizedEngine_BCyl_h
0002 #define ParametrizedEngine_BCyl_h
0003 /**
0004  *
0005  *
0006  *    B-field in Tracker volume - based on the TOSCA computation version 1103l
0007  *    (tuned on MTCC measured field (fall 2006))
0008  *
0009  *     In:   x[3]: coordinates (m)
0010  *    Out:   B[3]: Bx,By,Bz    (T)    (getBxyz)
0011  *    Out:   B[3]: Br,Bf,Bz    (T)    (getBrfz)
0012  *
0013  *    Valid for r<1.15m and |z|<2.80m
0014  *
0015  *  \author V.Karimaki 080228, 080407
0016  *  new float version V.I. October 2012
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     //  constexpr BCylParam(std::initializer_list<T> init) :
0027     template <typename... Args>
0028     constexpr BCylParam(Args... init)
0029         : prm{std::forward<Args>(init)...},
0030           //prm(std::forward<Args>(init)...),
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       // Function and its 3 derivatives
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   }  // namespace bcylDetails
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     // in meters and T  (Br needs to be multiplied by r)
0073     void compute(T r2, T z, T& Br, T& Bz) const {
0074       using namespace bcylDetails;
0075       //  if (r<1.15&&fabs(z)<2.8) // NOTE: check omitted, is done already by the wrapper! (NA)
0076       z -= pars.prm[3];  // max Bz point is shifted in z
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));  // double Gaussian
0092       Br += corBr;
0093       Bz += corBz;
0094     }
0095
0096   private:
0097     BCylParam<T> pars;
0098   };
0099 }  // namespace magfieldparam
0100
0101 #endif