CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
OAE85lParametrizedMagneticField.cc
Go to the documentation of this file.
1 
11 
12 #include "TkBfield.h"
13 
14 using namespace std;
15 using namespace magfieldparam;
16 
17 
19  float a_,
20  float l_)
21  : b0(b0_),
22  l(l_),
23  a(a_)
24 {
25  init();
26 }
27 
28 
30  b0 = parameters.getParameter<double>("b0");
31  l = parameters.getParameter<double>("l");
32  a = parameters.getParameter<double>("a");
33  init();
34 }
35 
37  ap2=4.0*a*a/(l*l);
38  hb0=0.5*b0*sqrt(1.0+ap2);
39  hlova=1.0/sqrt(ap2);
40  ainv=2.0*hlova/l;
41 }
42 
43 
44 
46 
47 
50 
51  if (isDefined(gp)) {
52  return inTeslaUnchecked(gp);
53  } else {
54  edm::LogWarning("MagneticField|FieldOutsideValidity") << " Point " << gp << " is outside the validity region of OAE85lParametrizedMagneticField";
55  return GlobalVector();
56  }
57 }
58 
59 
62 
63 // Method formerly named "trackerField"
64 
65 //
66 // B-field in Tracker volume
67 //
68 // In: xyz[3]: coordinates (m)
69 // Out: bxyz[3]: Bx,By,Bz (kG)
70 //
71 // Valid for r<1.2 and |z|<3.0 V.Karimäki 040301
72 // Updated for CMSSW field 070424
73 //
74 
75 // b0=field at centre, l=solenoid length, a=radius (m) (phenomen. parameters)
76 
77  // FIXME: beware of statics...
78 // static float ap2=4.0*a*a/(l*l);
79 // static float hb0=0.5*b0*sqrt(1.0+ap2);
80 // static float hlova=1.0/sqrt(ap2);
81 // static float ainv=2.0*hlova/l;
82  float xyz[3];//, bxyz[3];
83  // convert to m (CMSSW)
84  xyz[0]=0.01*gp.x();
85  xyz[1]=0.01*gp.y();
86  xyz[2]=0.01*gp.z();
87 
88  float r=sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
89  float z=xyz[2];
90  float az = fabs(z);
91  // if (r<1.2&&az<3.0) // NOTE: check omitted, is done already by inTesla (NA)
92  {
93  float zainv=z*ainv;
94  float rinv=(r>0.0) ? 1.0/r : 0.0;
95  float u=hlova-zainv;
96  float v=hlova+zainv;
97  float fu[5],gv[5];
98  ffunkti(u,fu);
99  ffunkti(v,gv);
100  float rat=r*ainv;
101  float corrr=0.00894*r*z*(az-2.2221)*(az-2.2221);
102  float corrz=-0.02996*exp(-0.5*(az-1.9820)*(az-1.9820)/(0.78915*0.78915));
103  float br=hb0*0.5*rat*(fu[1]-gv[1]-0.125*(fu[3]-gv[3])*rat*rat)+corrr;
104  float bz=hb0*(fu[0]+gv[0]-(fu[2]+gv[2])*0.25*rat*rat)+corrz;
105  //bxyz[0]=br*xyz[0]*rinv;
106  //bxyz[1]=br*xyz[1]*rinv;
107  //bxyz[2]=bz;
109  // bvec.x()=0.1*br*xyz[0]*rinv;
110  //bvec.y()=0.1*br*xyz[1]*rinv;
111  //bvec.z()=0.1*bz;
112  return 0.1*GlobalVector(br*xyz[0]*rinv, br*xyz[1]*rinv, bz);
113  }
114  // else {
115  // cout <<"The point is outside the region r<1.2m && |z|<3.0m"<<endl;
116  //}
117  // return GlobalVector();
118 }
119 
120 void OAE85lParametrizedMagneticField::ffunkti(float u, float* ff) const {
121 // Function and its 3 derivatives
122  float a,b;
123  a=1.0/(1.0+u*u);
124  b=sqrt(a);
125  ff[0]=u*b;
126  ff[1]=a*b;
127  ff[2]=-3.0*u*a*ff[1];
128  ff[3]=a*ff[2]*((1.0/u)-4.0*u);
129  return;
130 }
131 
132 
133 bool
135  return (gp.perp()<120. && fabs(gp.z())<300.);
136 }
137 
T getParameter(std::string const &) const
dictionary parameters
Definition: Parameters.py:2
T perp() const
Definition: PV3DBase.h:72
T y() const
Definition: PV3DBase.h:63
float float float z
OAE85lParametrizedMagneticField(float b0_=40.681, float a_=4.6430, float l_=15.284)
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
GlobalVector inTeslaUnchecked(const GlobalPoint &gp) const
GlobalVector inTesla(const GlobalPoint &gp) const
Field value ad specified global point, in Tesla.
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
T x() const
Definition: PV3DBase.h:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
bool isDefined(const GlobalPoint &gp) const
True if the point is within the region where the concrete field.