CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
TrackUtil/TrackUtil-00-00-12/src/Bfield.cxx
Go to the documentation of this file.
1//
2// Bfield class
3//
4// 27-Mar-1999 : KUNIYA Toshio
5// Enabled QCS compornent
6//
7// 21-Feb-1999 : KUNIYA Toshio
8// Keeping comatibility, Bfield class is modified.
9// No longer fortran common block is used for bfield map.
10// Access functions are prepared for fortran call.
11//
12
13#include <iostream>
14#include "TrackUtil/Bfield.h"
15
16Bfield *
17Bfield::_field[200] = {0}; // All of 200 elements are initialized.
18
19Bfield *
20Bfield::getBfield(int imap) {
21 if (! _field[imap]) _field[imap] = new Bfield(imap);
22 return _field[imap];
23}
24
25
26Bfield::Bfield(int imap) {
27 std::cout << std::endl;
28 std::cout << "***********************************************" << std::endl;
29 std::cout << " Bfield class MAP ID = " << imap << std::endl;
30 std::cout << " #### R < 174 cm, -152 < Z < 246 cm #### " << std::endl;
31 std::cout << " C++ version 1.00 " << std::endl;
32 std::cout << "***********************************************" << std::endl;
33
34 const float uniformBz[10] = {0, 15, 14.5, 14, 13, 12.5, 12, 11, 10, 15.5};
35
36 //...initialization
37 for(int i=0; i<175; i++)
38 for(int j=0; j<399; j++) {
39 _Bz[i][j] = 0;
40 _Br[i][j] = 0;
41 _Bz[i][j] = 0;
42 if(i<101 && j<163) {
43 _BzQR[i][j] = 0;
44 _BrQR[i][j] = 0;
45 _BphiQR[i][j] = 0;
46 }
47 }
48 for(int i=0; i<17; i++)
49 for(int j=0; j<51; j++)
50 for(int k=0; k<52; k++) {
51 _BzQL[i][j][k] = 0;
52 _BrQL[i][j][k] = 0;
53 _BphiQL[i][j][k] =0;
54 }
55
56 //...
57 _fieldID = imap;
58
59 //...read B field map
60
61 if(imap<10){
62 //
63 // uniform B field map
64 //
65 m_Bx = 0.;
66 m_By = 0.;
67 m_Bz = uniformBz[imap];
68 m_Bfld.setX((double) m_Bx);
69 m_Bfld.setY((double) m_By);
70 m_Bfld.setZ((double) m_Bz);
71 std::cout << "Bfield class >> creating uniform B field with id = " << imap;
72 std::cout << ", (Bx,By,Bz) = "<<m_Bfld<<std::endl;
73 } else {
74 //
75 // non-uniform B field map
76 //
77 /*
78 std::cout << "Bfield class >> loading non-uniform B field map" << std::endl;
79 geo_coil_readmap_(&imap, _Bz, _Br, _Bphi);
80 if( _fieldID == 21 ) {
81 std::cout << "Bfield class >> loading QCS" << std::endl;
82 geo_coil_readqcsrmap_(_BzQR,_BrQR, _BphiQR);
83 geo_coil_readqcslmap_(_BzQL,_BrQL, _BphiQL);
84 }
85 updateCache(0., 0., 0.);
86 */
87 }
88 std::cout << std::endl;
89
90}
91
92//Bfield::~Bfield(){};
93
94const Hep3Vector &
95Bfield::fieldMap(float x, float y, float z) const {
96
97 if(_fieldID > 10){
98 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
99 }
100
101 return m_Bfld;
102}
103
104const Hep3Vector &
105Bfield::fieldMap(const HepPoint3D & xyz) const{
106
107 if(_fieldID > 10){
108 float x = xyz.x();
109 float y = xyz.y();
110 float z = xyz.z();
111 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
112 }
113
114 return m_Bfld;
115}
116
117void
118Bfield::fieldMap(float *position, float *field) {
119 if(_fieldID > 10){
120 float x = position[0];
121 float y = position[1];
122 float z = position[2];
123 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
124 }
125 field[0] = m_Bx;
126 field[1] = m_By;
127 field[2] = m_Bz;
128 return;
129}
130
131float
132Bfield::bx(float x, float y, float z) const {
133
134 if(_fieldID > 10){
135 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
136 }
137
138 return m_Bx;
139}
140
141float
142Bfield::by(float x, float y, float z) const {
143
144 if(_fieldID > 10){
145 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
146 }
147
148 return m_By;
149}
150
151float
152Bfield::bz(float x, float y, float z) const {
153
154 if(_fieldID > 10){
155 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
156 }
157
158 return m_Bz;
159}
160
161float
162Bfield::bx(const HepPoint3D & xyz) const{
163
164 if(_fieldID > 10){
165 float x = xyz.x();
166 float y = xyz.y();
167 float z = xyz.z();
168 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
169 }
170
171 return m_Bx;
172}
173
174float
175Bfield::by(const HepPoint3D & xyz) const{
176
177 if(_fieldID > 10){
178 float x = xyz.x();
179 float y = xyz.y();
180 float z = xyz.z();
181 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
182 }
183
184 return m_By;
185}
186
187float
188Bfield::bz(const HepPoint3D & xyz) const{
189
190 if(_fieldID > 10){
191 float x = xyz.x();
192 float y = xyz.y();
193 float z = xyz.z();
194 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
195 }
196 return m_Bz;
197}
198
199void
200Bfield::updateCache(float x, float y, float z) const{
201
202 // this function is only for non-uniform B field
203
204 if( _fieldID <= 10 ) return;
205
206 float PI = 3.14159;
207
208 //...
209 float r = (float)sqrt((double)x*(double)x + (double)y*(double)y);
210 float phi = (float)atan2((double)y, (double)x);
211
212 //... [cm] --> [mm]
213 float zmm = z * 10.;
214 float rmm = r * 10.;
215 //... make index
216 int tz = (int) (( zmm + 1520.)/10.);
217 int tr = (int) (rmm/10.);
218
219 //...
220 float bz = 0., br = 0., bphi = 0.;
221
222 if(zmm >= -1520. && zmm < 2460. && rmm >= 0. && rmm < 1740.){
223 if(_Bz[tr][tz] && _Bz[tr][tz+1]){
224 float pz = (zmm + 1520.)/10.- (float)tz;
225 float pr = rmm/10.- (float)tr;
226
227 //...
228 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
229 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
230 //...
231 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
232 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
233 //...
234 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
235 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
236
237 if(_fieldID == 21) {
238 //
239 // QCS Right
240 //
241 if( zmm>=800. && zmm < 2420. && rmm < 1000. ) {
242 int tqz = (int)( (zmm-800.)/10. );
243 bz += (((_BzQR[tr][tqz]*(1.-pz)+_BzQR[tr][tqz+1]*pz)*(1.-pr)
244 +(_BzQR[tr+1][tqz]*(1.-pz)+_BzQR[tr+1][tqz+1]*pz)*pr)
245 *(float)sin((double)(2.*phi+2./180.*(double)PI)));
246 br += (((_BrQR[tr][tqz]*(1.-pz)+_BrQR[tr][tqz+1]*pz)*(1.-pr)
247 +(_BrQR[tr+1][tqz]*(1.-pz)+_BrQR[tr+1][tqz+1]*pz)*pr)
248 *(float)sin((double)(2.*phi+2./180.*(double)PI)));
249 bphi += (((_BphiQR[tr][tqz]*(1.-pz)
250 +_BphiQR[tr][tqz+1]*pz)*(1.-pr)
251 +(_BphiQR[tr+1][tqz]*(1.-pz)
252 +_BphiQR[tr+1][tqz+1]*pz)*pr)
253 *(float)cos((double)(2.*phi+2./180.*(double)PI)));
254 }
255 //
256 // QCS Left
257 //
258 if(zmm<=-500. && zmm>-1520. && rmm<1000.) {
259 int tqz = (int)((-zmm-500.)/20.);
260 int tqr = (int)(tr/2.);
261 pz = (pz+(float)( tz-(int)(tz/2.)*2. ))/2.;
262 pr = ( pr + (float)(tr-tqr*2) )/2.;
263 float f = 1.;
264 float phi_tmp = phi;
265 if( phi_tmp < (PI/2.) && phi_tmp >= (-PI/2.) ) {
266 phi_tmp = PI-phi_tmp;
267 f =-1.;
268 } else if( phi_tmp < -PI/2. ) {
269 phi_tmp = 2.*PI+phi_tmp;
270 }
271 float pphi = ( phi_tmp-PI/2. )/(11.25/180.*PI);
272 int tphi = (int)pphi;
273 pphi -= (float)tphi;
274 if (tphi >= 16) tphi -= 16;
275
276 bz += f*
277 (((_BzQL[tphi][tqr][tqz]*(1.-pz)+_BzQL[tphi][tqr][tqz+1]*pz)
278 *(1-pr)+(_BzQL[tphi][tqr+1][tqz]*(1.-pz)
279 +_BzQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
280 +((_BzQL[tphi+1][tqr][tqz]*(1.-pz)
281 +_BzQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
282 +(_BzQL[tphi+1][tqr+1][tqz]*(1.-pz)
283 +_BzQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
284 br += f*
285 (((_BrQL[tphi][tqr][tqz]*(1.- pz)
286 +_BrQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
287 +(_BrQL[tphi][tqr+1][tqz]*(1.-pz)
288 +_BrQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1-pphi)
289 +((_BrQL[tphi+1][tqr][tqz]*(1.- pz)
290 +_BrQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
291 +(_BrQL[tphi+1][tqr+1][tqz]*(1.-pz)
292 +_BrQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
293 bphi += f*
294 (((_BphiQL[tphi][tqr][tqz]*(1.- pz)
295 +_BphiQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
296 +(_BphiQL[tphi][tqr+1][tqz]*(1.-pz)
297 +_BphiQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
298 +((_BphiQL[tphi+1][tqr][tqz]*(1.- pz)
299 +_BphiQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
300 +(_BphiQL[tphi+1][tqr+1][tqz]*(1.-pz)
301 +_BphiQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
302 }
303 }
304 } else {
305 bz=0.;
306 br=0.;
307 bphi=0.;
308 }
309 } else if(zmm>=-1520. && zmm<=2460. && rmm<=0. && rmm<=1740.){
310 if(tz == 246) tz=tz-1;
311 if(tr == 174) tr=tr-1;
312 float pz= (zmm + 1520.)/10.- (float)tz;
313 float pr= rmm/10.- (float)tr;
314 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
315 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
316
317 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
318 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
319
320 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
321 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
322 } else {
323 bz =0.;
324 br =0.;
325 bphi =0.;
326 }
327
328
329 //... Set B field
330 float Bmag_xy = (float)sqrt(br*br + bphi*bphi);
331 double Bphi_rp = (float)atan2( (double)bphi, (double)br );
332 m_Bx = Bmag_xy * (float)cos((double)phi + Bphi_rp)/1000.;
333 m_By = Bmag_xy * (float)sin((double)phi + Bphi_rp)/1000.;
334 m_Bz = bz/1000.;
335 m_x = x;
336 m_y = y;
337 m_z = z;
338 m_Bfld.setX((double) m_Bx);
339 m_Bfld.setY((double) m_By);
340 m_Bfld.setZ((double) m_Bz);
341}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
PI
const double PI
Definition PipiJpsi.cxx:55
Bfield::Bfield
Bfield(int)
Constructor, Destructor.
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:58
Bfield::bz
double bz(double x, double y, double z) const
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:184
Bfield::fieldMap
const Hep3Vector & fieldMap(double x, double y, double z) const
returns B field
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:127
Bfield::by
double by(double x, double y, double z) const
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:174
Bfield::bx
double bx(double x, double y, double z) const
returns an element of B field
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:164
Bfield::getBfield
static Bfield * getBfield(int)
returns Bfield object.
Definition KalFitAlg/KalFitAlg-00-15-20/src/coil/Bfield.cxx:52