ParametrizedEngine/src/TkBfield.cc

0001 #include "TkBfield.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 // #include <iostream>
0006 // #include <iomanip>
0007 #include <cmath>
0008 #include <algorithm>
0009
0010 using namespace magfieldparam;
0011
0012 namespace {
0013
0014   const BCylParam<float> fpar1{
0015       4.90541f, 17.8768f, 2.02355f, 0.0210538f, 0.000321885f, 2.37511f, 0.00326725f, 2.07656f, 1.71879f};  // 2.0T-2G
0016   const BCylParam<float> fpar2{
0017       4.41982f, 15.7732f, 3.02621f, 0.0197814f, 0.000515759f, 2.43385f, 0.00584258f, 2.11333f, 1.76079f};  // 3.0T-2G
0018   const BCylParam<float> fpar3{
0019       4.30161f, 15.2586f, 3.51926f, 0.0183494f, 0.000606773f, 2.45110f, 0.00709986f, 2.12161f, 1.77038f};  // 3.5T-2G
0020   const BCylParam<float> fpar4{
0021       4.24326f, 15.0201f, 3.81492f, 0.0178712f, 0.000656527f, 2.45818f, 0.00778695f, 2.12500f, 1.77436f};  // 3.8T-2G
0022   const BCylParam<float> fpar5{
0023       4.21136f, 14.8824f, 4.01683f, 0.0175932f, 0.000695541f, 2.45311f, 0.00813447f, 2.11688f, 1.76076f};  // 4.0T-2G
0024   std::string const flds[] = {"2_0T", "3_0T", "3_5T", "3_8T", "4_0T"};
0025   constexpr float flds_f[] = {2.0, 3.0, 3.5, 3.8, 4.0};
0026   BCylParam<float> const fpars[]{fpar1, fpar2, fpar3, fpar4, fpar5};
0027
0028   BCylParam<float> const& findPar(float fld) {
0029     for (int i = 0; i < 5; ++i) {
0030       if (fabs(fld - flds_f[i]) < 0.001) {
0031         return fpars[i];
0032       }
0033     }
0034     throw cms::Exception("BadParameters") << "Undefined field value " << fld;
0035   }
0036
0037   BCylParam<float> const& findPar(std::string fld) {
0038     auto f = std::find(flds, flds + 5, fld);
0039     if (f - flds > 4)
0040       throw cms::Exception("BadParameters") << "Undefined key - "  // abort!\n";
0041                                             << "Defined keys are: \"2_0T\" \"3_0T\" \"3_5T\" \"3_8T\" and \"4_0T\"\n";
0042     return fpars[f - flds];
0043   }
0044
0045 }  // namespace
0046
0047 TkBfield::TkBfield(float fld) : bcyl(findPar(fld)) {}
0048
0049 TkBfield::TkBfield(std::string fld) : bcyl(findPar(fld)) {}
0050
0051 void TkBfield::getBrfz(float const* __restrict__ x, float* __restrict__ Brfz) const {
0052   float br;
0053   float bz;
0054   float r2 = x[0] * x[0] + x[1] * x[1];
0055   bcyl(x[0] * x[0] + x[1] * x[1], x[2], br, bz);
0056   Brfz[0] = std::sqrt(r2) * br;
0057   Brfz[1] = 0;
0058   Brfz[2] = bz;
0059 }
0060
0061 void TkBfield::getBxyz(float const* __restrict__ x, float* __restrict__ Bxyz) const {
0062   float br;
0063   float bz;
0064   float r2 = x[0] * x[0] + x[1] * x[1];
0065   bcyl(r2, x[2], br, bz);
0066   Bxyz[0] = br * x[0];
0067   Bxyz[1] = br * x[1];
0068   Bxyz[2] = bz;
0069 }