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