BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/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
14#include "KalFitAlg/coil/Bfield.h"
15
16
17#include "CLHEP/Matrix/Vector.h"
18#include "CLHEP/Matrix/Matrix.h"
19#include "CLHEP/Matrix/SymMatrix.h"
20#include "CLHEP/Vector/ThreeVector.h"
21#include "CLHEP/Geometry/Point3D.h"
22#ifndef ENABLE_BACKWARDS_COMPATIBILITY
23 typedef HepGeom::Point3D<double> HepPoint3D;
24#endif
25
26
27
28using CLHEP::HepVector;
29using CLHEP::Hep3Vector;
30using CLHEP::HepMatrix;
31using CLHEP::HepSymMatrix;
32
33
34// prototype of file scoped function to call fortran routine
35// ... read B field map with ID = imap from a file
36/*wangdy
37extern "C" {
38 extern void geo_coil_readmap_
39 (int *imap, double (*bz)[399], double (*br)[399], double (*bphi)[399]);
40 extern void geo_coil_readqcsrmap_
41 (double (*bzqr)[163], double (*brqr)[163], double (*bphiqr)[163]);
42 extern void geo_coil_readqcslmap_
43 (double (*bzrl)[51][52], double (*brql)[51][52], double (*bphiqr)[51][52]);
44}
45 */
46
47
48Bfield *
49Bfield::_field[200] = {0}; // All of 200 elements are initialized.
50
51Bfield *
52Bfield::getBfield(int imap) {
53 if (! _field[imap]) _field[imap] = new Bfield(imap);
54 return _field[imap];
55}
56
57
58Bfield::Bfield(int imap) {
59 std::cout << std::endl;
60 std::cout << "***********************************************" << std::endl;
61 std::cout << " Bfield class MAP ID = " << imap << std::endl;
62 std::cout << " #### R < 174 cm, -152 < Z < 246 cm #### " << std::endl;
63 std::cout << " C++ version 1.00 " << std::endl;
64 std::cout << "***********************************************" << std::endl;
65
66 const double uniformBz[10] = {0, 15, 14.5, 14, 13, 12.5, 12, 11, 10, 15.5};
67
68 //...initialization
69 for(int i=0; i<175; i++)
70 for(int j=0; j<399; j++) {
71 _Bz[i][j] = 0;
72 _Br[i][j] = 0;
73 _Bz[i][j] = 0;
74 if(i<101 && j<163) {
75 _BzQR[i][j] = 0;
76 _BrQR[i][j] = 0;
77 _BphiQR[i][j] = 0;
78 }
79 }
80 for(int i=0; i<17; i++)
81 for(int j=0; j<51; j++)
82 for(int k=0; k<52; k++) {
83 _BzQL[i][j][k] = 0;
84 _BrQL[i][j][k] = 0;
85 _BphiQL[i][j][k] =0;
86 }
87
88 //...
89 _fieldID = imap;
90
91 //...read B field map
92
93 if(imap<10){
94 //
95 // uniform B field map
96 //
97 m_Bx = 0.;
98 m_By = 0.;
99 m_Bz = uniformBz[imap];
100 m_Bfld.setX((double) m_Bx);
101 m_Bfld.setY((double) m_By);
102 m_Bfld.setZ((double) m_Bz);
103 std::cout << "Bfield class >> creating uniform B field with id = " << imap;
104 std::cout << ", (Bx,By,Bz) = "<<m_Bfld<<std::endl;
105 } else {
106 //
107 // non-uniform B field map
108 //
109 std::cout << "Bfield class >> loading non-uniform B field map" << std::endl;
110/*wangdy
111 geo_coil_readmap_(&imap, _Bz, _Br, _Bphi);
112 if( _fieldID == 21 ) {
113 std::cout << "Bfield class >> loading QCS" << std::endl;
114 geo_coil_readqcsrmap_(_BzQR,_BrQR, _BphiQR);
115 geo_coil_readqcslmap_(_BzQL,_BrQL, _BphiQL);
116 }
117 */
118 updateCache(0., 0., 0.);
119 }
120 std::cout << std::endl;
121
122}
123
124//Bfield::~Bfield(){};
125
126const Hep3Vector &
127Bfield::fieldMap(double x, double y, double z) const {
128
129 if(_fieldID > 10){
130 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
131 }
132
133 return m_Bfld;
134}
135
136const Hep3Vector &
137Bfield::fieldMap(const HepPoint3D &xyz) const{
138
139 if(_fieldID > 10){
140 double x = xyz.x();
141 double y = xyz.y();
142 double z = xyz.z();
143 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
144 }
145
146 return m_Bfld;
147}
148
149void
150Bfield::fieldMap(double *position, double *field) {
151 if(_fieldID > 10){
152 double x = position[0];
153 double y = position[1];
154 double z = position[2];
155 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
156 }
157 field[0] = m_Bx;
158 field[1] = m_By;
159 field[2] = m_Bz;
160 return;
161}
162
163double
164Bfield::bx(double x, double y, double z) const {
165
166 if(_fieldID > 10){
167 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
168 }
169
170 return m_Bx;
171}
172
173double
174Bfield::by(double x, double y, double z) const {
175
176 if(_fieldID > 10){
177 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
178 }
179
180 return m_By;
181}
182
183double
184Bfield::bz(double x, double y, double z) const {
185
186 if(_fieldID > 10){
187 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
188 }
189
190 return m_Bz;
191}
192
193double
194Bfield::bx(const HepPoint3D &xyz) const{
195
196 if(_fieldID > 10){
197 double x = xyz.x();
198 double y = xyz.y();
199 double z = xyz.z();
200 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
201 }
202
203 return m_Bx;
204}
205
206double
207Bfield::by(const HepPoint3D &xyz) const{
208
209 if(_fieldID > 10){
210 double x = xyz.x();
211 double y = xyz.y();
212 double z = xyz.z();
213 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
214 }
215
216 return m_By;
217}
218
219double
220Bfield::bz(const HepPoint3D &xyz) const{
221
222 if(_fieldID > 10){
223 double x = xyz.x();
224 double y = xyz.y();
225 double z = xyz.z();
226 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
227 }
228 return m_Bz;
229}
230
231void
232Bfield::updateCache(double x, double y, double z) const{
233
234 // this function is only for non-uniform B field
235
236 if( _fieldID <= 10 ) return;
237
238 double PI = 3.14159;
239
240 //...
241 double r = (double)sqrt((double)x*(double)x + (double)y*(double)y);
242 double phi = (double)atan2((double)y, (double)x);
243
244 //... [cm] --> [mm]
245 double zmm = z * 10.;
246 double rmm = r * 10.;
247 //... make index
248 int tz = (int) (( zmm + 1520.)/10.);
249 int tr = (int) (rmm/10.);
250
251 //...
252 double bz = 0., br = 0., bphi = 0.;
253
254 if(zmm >= -1520. && zmm < 2460. && rmm >= 0. && rmm < 1740.){
255 if(_Bz[tr][tz] && _Bz[tr][tz+1]){
256 double pz = (zmm + 1520.)/10.- (double)tz;
257 double pr = rmm/10.- (double)tr;
258
259 //...
260 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
261 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
262 //...
263 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
264 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
265 //...
266 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
267 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
268
269 if(_fieldID == 21) {
270 //
271 // QCS Right
272 //
273 if( zmm>=800. && zmm < 2420. && rmm < 1000. ) {
274 int tqz = (int)( (zmm-800.)/10. );
275 bz += (((_BzQR[tr][tqz]*(1.-pz)+_BzQR[tr][tqz+1]*pz)*(1.-pr)
276 +(_BzQR[tr+1][tqz]*(1.-pz)+_BzQR[tr+1][tqz+1]*pz)*pr)
277 *(double)sin((double)(2.*phi+2./180.*(double)PI)));
278 br += (((_BrQR[tr][tqz]*(1.-pz)+_BrQR[tr][tqz+1]*pz)*(1.-pr)
279 +(_BrQR[tr+1][tqz]*(1.-pz)+_BrQR[tr+1][tqz+1]*pz)*pr)
280 *(double)sin((double)(2.*phi+2./180.*(double)PI)));
281 bphi += (((_BphiQR[tr][tqz]*(1.-pz)
282 +_BphiQR[tr][tqz+1]*pz)*(1.-pr)
283 +(_BphiQR[tr+1][tqz]*(1.-pz)
284 +_BphiQR[tr+1][tqz+1]*pz)*pr)
285 *(double)cos((double)(2.*phi+2./180.*(double)PI)));
286 }
287 //
288 // QCS Left
289 //
290 if(zmm<=-500. && zmm>-1520. && rmm<1000.) {
291 int tqz = (int)((-zmm-500.)/20.);
292 int tqr = (int)(tr/2.);
293 pz = (pz+(double)( tz-(int)(tz/2.)*2. ))/2.;
294 pr = ( pr + (double)(tr-tqr*2) )/2.;
295 double f = 1.;
296 // if( phi < (PI/2.) && phi >= (-PI/2.) ) {
297 // phi = PI-phi;
298 // f =-1.;
299 // } else if( phi < -PI/2. ) {
300 // phi = 2.*PI+phi;
301 // }
302 // double pphi = ( phi-PI/2. )/(11.25/180.*PI);
303 double phi_tmp = phi;
304 if( phi_tmp < (PI/2.) && phi_tmp >= (-PI/2.) ) {
305 phi_tmp = PI-phi_tmp;
306 f =-1.;
307 } else if( phi_tmp < -PI/2. ) {
308 phi_tmp = 2.*PI+phi_tmp;
309 }
310 double pphi = ( phi_tmp-PI/2. )/(11.25/180.*PI);
311 int tphi = (int)pphi;
312 pphi -= (double)tphi;
313 if (tphi >= 16) tphi -= 16;
314
315 bz += f*
316 (((_BzQL[tphi][tqr][tqz]*(1.-pz)+_BzQL[tphi][tqr][tqz+1]*pz)
317 *(1-pr)+(_BzQL[tphi][tqr+1][tqz]*(1.-pz)
318 +_BzQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
319 +((_BzQL[tphi+1][tqr][tqz]*(1.-pz)
320 +_BzQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
321 +(_BzQL[tphi+1][tqr+1][tqz]*(1.-pz)
322 +_BzQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
323 br += f*
324 (((_BrQL[tphi][tqr][tqz]*(1.- pz)
325 +_BrQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
326 +(_BrQL[tphi][tqr+1][tqz]*(1.-pz)
327 +_BrQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1-pphi)
328 +((_BrQL[tphi+1][tqr][tqz]*(1.- pz)
329 +_BrQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
330 +(_BrQL[tphi+1][tqr+1][tqz]*(1.-pz)
331 +_BrQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
332 bphi += f*
333 (((_BphiQL[tphi][tqr][tqz]*(1.- pz)
334 +_BphiQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
335 +(_BphiQL[tphi][tqr+1][tqz]*(1.-pz)
336 +_BphiQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
337 +((_BphiQL[tphi+1][tqr][tqz]*(1.- pz)
338 +_BphiQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
339 +(_BphiQL[tphi+1][tqr+1][tqz]*(1.-pz)
340 +_BphiQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
341 }
342 }
343 } else {
344 bz=0.;
345 br=0.;
346 bphi=0.;
347 }
348 } else if(zmm>=-1520. && zmm<=2460. && rmm<=0. && rmm<=1740.){
349 if(tz == 246) tz=tz-1;
350 if(tr == 174) tr=tr-1;
351 double pz= (zmm + 1520.)/10.- (double)tz;
352 double pr= rmm/10.- (double)tr;
353 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
354 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
355
356 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
357 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
358
359 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
360 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
361 } else {
362 bz =0.;
363 br =0.;
364 bphi =0.;
365 }
366
367
368 //... Set B field
369 double Bmag_xy = (double)sqrt(br*br + bphi*bphi);
370 double Bphi_rp = (double)atan2( (double)bphi, (double)br );
371 m_Bx = Bmag_xy * (double)cos((double)phi + Bphi_rp)/1000.;
372 m_By = Bmag_xy * (double)sin((double)phi + Bphi_rp)/1000.;
373 //m_Bx = br * (double)cos((double)phi)/1000.;
374 //m_By = br * (double)sin((double)phi)/1000.;
375 m_Bz = bz/1000.;
376 m_x = x;
377 m_y = y;
378 m_z = z;
379 m_Bfld.setX((double) m_Bx);
380 m_Bfld.setY((double) m_By);
381 m_Bfld.setZ((double) m_Bz);
382}
383
384
385//
386//... Fortran callable access functions to Bfield class.
387//
388
389// initialize bfield map
390extern "C"
391void init_bfield_(int *imap) {
392 // create Bfiled map ID = imap
393 // Bfield *thisMap = Bfield::getBfield(*imap);
394 (void) Bfield::getBfield(*imap);
395 // It is OK even though 'thisMap' losts its scope at here.
396 // Because address of field map is not deleted from memory
397 // due to static linkaged Bfield class.
398 return;
399}
400
401// retrieving B field corresponding to the POSition
402extern "C"
403void get_bfield_(int *imap, double *pos, double *field, int *error) {
404 Bfield *thisMap = Bfield::getBfield(*imap);
405 // std::cout << " > accessing Bfield class from fortran routine." << std::endl;
406 if( thisMap != 0 ) {
407 thisMap->fieldMap(pos,field);
408 *error = 0;
409 }
410 else *error = 1;
411 return;
412}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
tr
DOUBLE_PRECISION tr[3]
Definition: EvtPhokharaDef.hh:64
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:23
get_bfield_
void get_bfield_(int *imap, double *pos, double *field, int *error)
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:403
PI
const double PI
Definition: PipiJpsi.cxx:55
Bfield::Bfield
Bfield(int)
Constructor, Destructor.
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:58
Bfield::bz
double bz(double x, double y, double z) const
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:184
Bfield::fieldMap
const Hep3Vector & fieldMap(double x, double y, double z) const
returns B field
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:127
Bfield::by
double by(double x, double y, double z) const
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/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-07-55-p03/src/coil/Bfield.cxx:164
Bfield::getBfield
static Bfield * getBfield(int)
returns Bfield object.
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/coil/Bfield.cxx:52
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
x
double x[1000]
Definition: draw_charge_space_coarse_II.cxx:8
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")