CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/MagneticField/ParametrizedEngine/plugins/TkBfield.cc

Go to the documentation of this file.
00001 #include "TkBfield.h"
00002 
00003 #include "FWCore/Utilities/interface/Exception.h"
00004 
00005 // #include <iostream>
00006 // #include <iomanip>
00007 #include <cmath>
00008 #include "DataFormats/Math/interface/SSEVec.h"
00009 
00010 using namespace magfieldparam;
00011 
00012 
00013 TkBfield::TkBfield(std::string fld) {
00014   double p1[]={4.90541,17.8768,2.02355,0.0210538,0.000321885,2.37511,0.00326725,2.07656,1.71879}; // 2.0T-2G
00015   double p2[]={4.41982,15.7732,3.02621,0.0197814,0.000515759,2.43385,0.00584258,2.11333,1.76079}; // 3.0T-2G
00016   double p3[]={4.30161,15.2586,3.51926,0.0183494,0.000606773,2.45110,0.00709986,2.12161,1.77038}; // 3.5T-2G
00017   double p4[]={4.24326,15.0201,3.81492,0.0178712,0.000656527,2.45818,0.00778695,2.12500,1.77436}; // 3.8T-2G
00018   double p5[]={4.21136,14.8824,4.01683,0.0175932,0.000695541,2.45311,0.00813447,2.11688,1.76076}; // 4.0T-2G
00019   prm[0]=0;
00020   if (fld=="2_0T") for (int i=0; i<9; i++) prm[i]=p1[i];
00021   if (fld=="3_0T") for (int i=0; i<9; i++) prm[i]=p2[i];
00022   if (fld=="3_5T") for (int i=0; i<9; i++) prm[i]=p3[i];
00023   if (fld=="3_8T") for (int i=0; i<9; i++) prm[i]=p4[i];
00024   if (fld=="4_0T") for (int i=0; i<9; i++) prm[i]=p5[i];
00025   //  cout<<std::endl<<"Instantiation of TkBfield with key "<<fld<<endl;
00026   if (!prm[0]) {
00027     throw cms::Exception("BadParameters") << "Undefined key - " // abort!\n";
00028     <<"Defined keys are: \"2_0T\" \"3_0T\" \"3_5T\" \"3_8T\" and \"4_0T\"\n";
00029     // exit(1);
00030   }
00031   ap2=4*prm[0]*prm[0]/(prm[1]*prm[1]);  
00032   hb0=0.5*prm[2]*std::sqrt(1.0+ap2);
00033   hlova=1/std::sqrt(ap2);
00034   ainv=2*hlova/prm[1];
00035   coeff=1/(prm[8]*prm[8]);
00036 }
00037 
00038 namespace {
00039 
00040   template<typename T>
00041   inline void ffunkti(T u, T * __restrict__ ff) __attribute__((always_inline)) __attribute__ ((pure));
00042  
00043   template<typename T>
00044   inline void ffunkti(T u, T * __restrict__ ff) {
00045     // Function and its 3 derivatives
00046     T a,b,a2,u2;
00047     u2=u*u; 
00048     a=T(1)/(T(1)+u2);
00049     a2=-3*a*a;
00050     b=mathSSE::sqrt(a);
00051     ff[0]=u*b;
00052     ff[1]=a*b;
00053     ff[2]=a2*ff[0];
00054     ff[3]=a2*ff[1]*(T(1)-4*u2);
00055   }
00056 }
00057 
00058 #if (defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ > 4)) || defined(__clang__)
00059 
00060 void  TkBfield::Bcyl(double r, double z, double * __restrict__ Bw)  const{
00061   z-=prm[3];                    // max Bz point is shifted in z
00062   double az=std::abs(z);
00063   double zainv=z*ainv;
00064 
00065   mathSSE::Vec2D uv(hlova-zainv,hlova+zainv);
00066   mathSSE::Vec2D fuv[4];
00067   ffunkti(uv,fuv);
00068 
00069   double rat=0.5*r*ainv;
00070   double rat2=rat*rat;
00071   Bw[0]=hb0*rat*(fuv[1][0]-fuv[1][1]-(fuv[3][0]-fuv[3][1])*rat2*0.5);
00072   Bw[1]=0;
00073   Bw[2]=hb0*(fuv[0][0]+fuv[0][1]-(fuv[2][0]+fuv[2][1])*rat2);
00074 
00075   double corBr= prm[4]*r*z*(az-prm[5])*(az-prm[5]);
00076   double corBz=-prm[6]*(exp(-(z-prm[7])*(z-prm[7])*coeff) +
00077                         exp(-(z+prm[7])*(z+prm[7])*coeff)
00078                         ); // double Gaussian
00079   Bw[0]+=corBr;
00080   Bw[2]+=corBz;
00081 }
00082 
00083 #else
00084 
00085 void  TkBfield::Bcyl(double r, double z, double * __restrict__ Bw)  const{
00086   //  if (r<1.15&&fabs(z)<2.8) // NOTE: check omitted, is done already by the wrapper! (NA)
00087   z-=prm[3];                    // max Bz point is shifted in z
00088   double az=std::abs(z);
00089   double zainv=z*ainv;
00090   double u=hlova-zainv;
00091   double v=hlova+zainv;
00092   double fu[4],gv[4];
00093   ffunkti(u,fu);
00094   ffunkti(v,gv);
00095   double rat=0.5*r*ainv;
00096   double rat2=rat*rat;
00097   Bw[0]=hb0*rat*(fu[1]-gv[1]-(fu[3]-gv[3])*rat2*0.5);
00098   Bw[1]=0;
00099   Bw[2]=hb0*(fu[0]+gv[0]-(fu[2]+gv[2])*rat2);
00100   double corBr= prm[4]*r*z*(az-prm[5])*(az-prm[5]);
00101   //    double corBz=-prm[6]*exp(coeff*(az-prm[7])*(az-prm[7]));
00102   //    double corBz=-prm[6]/(1+coeff*(az-prm[7])*(az-prm[7]));
00103   //  double corBz=-prm[6]*(exp(-(z-prm[7])*(z-prm[7])/(prm[8]*prm[8]))
00104   //                    + exp(-(z+prm[7])*(z+prm[7])/(prm[8]*prm[8]))); // double Gaussian
00105   double corBz=-prm[6]*(exp(-(z-prm[7])*(z-prm[7])*coeff) +
00106                         exp(-(z+prm[7])*(z+prm[7])*coeff)
00107                         ); // double Gaussian
00108   Bw[0]+=corBr;
00109   Bw[2]+=corBz;
00110   //   } else {
00111   //     cout <<"TkBfield: The point is outside the region r<1.15m && |z|<2.8m"<<endl;
00112 }
00113 #endif
00114 
00115 
00116 void TkBfield::getBrfz(double const  * __restrict__ x, double * __restrict__ Brfz)  const {
00117   Bcyl(sqrt(x[0]*x[0]+x[1]*x[1]),x[2], Brfz);
00118 }
00119 
00120 void TkBfield::getBxyz(double const  * __restrict__ x, double * __restrict__ Bxyz)  const {
00121   double Bw[3];
00122   double r=sqrt(x[0]*x[0]+x[1]*x[1]);
00123   Bcyl(r, x[2], Bw);
00124   double rinv=(r>0) ? 1/r:0;
00125   Bxyz[0]=Bw[0]*x[0]*rinv;
00126   Bxyz[1]=Bw[0]*x[1]*rinv;
00127   Bxyz[2]=Bw[2];
00128 }
00129