00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <iostream>
00014
00015
00016
00017 #include "TrkReco/Bfield.h"
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 Bfield *
00040 Bfield::_field[200] = {0};
00041
00042 Bfield *
00043 Bfield::getBfield(int imap) {
00044 if (! _field[imap]) _field[imap] = new Bfield(imap);
00045 return _field[imap];
00046 }
00047
00048
00049 Bfield::Bfield(int imap) {
00050 std::cout << std::endl;
00051 std::cout << "***********************************************" << std::endl;
00052 std::cout << " Bfield class MAP ID = " << imap << std::endl;
00053 std::cout << " #### R < 174 cm, -152 < Z < 246 cm #### " << std::endl;
00054 std::cout << " C++ version 1.00 " << std::endl;
00055 std::cout << "***********************************************" << std::endl;
00056
00057 const float uniformBz[10] = {0, 15, 14.5, 14, 13, 12.5, 12, 11, 10, 15.5};
00058
00059
00060 for(int i=0; i<175; i++)
00061 for(int j=0; j<399; j++) {
00062 _Bz[i][j] = 0;
00063 _Br[i][j] = 0;
00064 _Bz[i][j] = 0;
00065 if(i<101 && j<163) {
00066 _BzQR[i][j] = 0;
00067 _BrQR[i][j] = 0;
00068 _BphiQR[i][j] = 0;
00069 }
00070 }
00071 for(int i=0; i<17; i++)
00072 for(int j=0; j<51; j++)
00073 for(int k=0; k<52; k++) {
00074 _BzQL[i][j][k] = 0;
00075 _BrQL[i][j][k] = 0;
00076 _BphiQL[i][j][k] =0;
00077 }
00078
00079
00080 _fieldID = imap;
00081
00082
00083
00084 if(imap<10){
00085
00086
00087
00088 m_Bx = 0.;
00089 m_By = 0.;
00090 m_Bz = uniformBz[imap];
00091 m_Bfld.setX((double) m_Bx);
00092 m_Bfld.setY((double) m_By);
00093 m_Bfld.setZ((double) m_Bz);
00094 std::cout << "Bfield class >> creating uniform B field with id = " << imap;
00095 std::cout << ", (Bx,By,Bz) = "<<m_Bfld<<std::endl;
00096 } else {
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 }
00111 std::cout << std::endl;
00112
00113 }
00114
00115
00116
00117 const Hep3Vector &
00118 Bfield::fieldMap(float x, float y, float z) const {
00119
00120 if(_fieldID > 10){
00121 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00122 }
00123
00124 return m_Bfld;
00125 }
00126
00127 const Hep3Vector &
00128 Bfield::fieldMap(const HepPoint3D & xyz) const{
00129
00130 if(_fieldID > 10){
00131 float x = xyz.x();
00132 float y = xyz.y();
00133 float z = xyz.z();
00134 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00135 }
00136
00137 return m_Bfld;
00138 }
00139
00140 void
00141 Bfield::fieldMap(float *position, float *field) {
00142 if(_fieldID > 10){
00143 float x = position[0];
00144 float y = position[1];
00145 float z = position[2];
00146 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00147 }
00148 field[0] = m_Bx;
00149 field[1] = m_By;
00150 field[2] = m_Bz;
00151 return;
00152 }
00153
00154 float
00155 Bfield::bx(float x, float y, float z) const {
00156
00157 if(_fieldID > 10){
00158 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00159 }
00160
00161 return m_Bx;
00162 }
00163
00164 float
00165 Bfield::by(float x, float y, float z) const {
00166
00167 if(_fieldID > 10){
00168 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00169 }
00170
00171 return m_By;
00172 }
00173
00174 float
00175 Bfield::bz(float x, float y, float z) const {
00176
00177 if(_fieldID > 10){
00178 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00179 }
00180
00181 return m_Bz;
00182 }
00183
00184 float
00185 Bfield::bx(const HepPoint3D & xyz) const{
00186
00187 if(_fieldID > 10){
00188 float x = xyz.x();
00189 float y = xyz.y();
00190 float z = xyz.z();
00191 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00192 }
00193
00194 return m_Bx;
00195 }
00196
00197 float
00198 Bfield::by(const HepPoint3D & xyz) const{
00199
00200 if(_fieldID > 10){
00201 float x = xyz.x();
00202 float y = xyz.y();
00203 float z = xyz.z();
00204 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00205 }
00206
00207 return m_By;
00208 }
00209
00210 float
00211 Bfield::bz(const HepPoint3D & xyz) const{
00212
00213 if(_fieldID > 10){
00214 float x = xyz.x();
00215 float y = xyz.y();
00216 float z = xyz.z();
00217 if(x != m_x || y != m_y || z != m_z) updateCache(x, y, z);
00218 }
00219 return m_Bz;
00220 }
00221
00222 void
00223 Bfield::updateCache(float x, float y, float z) const{
00224
00225
00226
00227 if( _fieldID <= 10 ) return;
00228
00229 float PI = 3.14159;
00230
00231
00232 float r = (float)sqrt((double)x*(double)x + (double)y*(double)y);
00233 float phi = (float)atan2((double)y, (double)x);
00234
00235
00236 float zmm = z * 10.;
00237 float rmm = r * 10.;
00238
00239 int tz = (int) (( zmm + 1520.)/10.);
00240 int tr = (int) (rmm/10.);
00241
00242
00243 float bz = 0., br = 0., bphi = 0.;
00244
00245 if(zmm >= -1520. && zmm < 2460. && rmm >= 0. && rmm < 1740.){
00246 if(_Bz[tr][tz] && _Bz[tr][tz+1]){
00247 float pz = (zmm + 1520.)/10.- (float)tz;
00248 float pr = rmm/10.- (float)tr;
00249
00250
00251 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
00252 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
00253
00254 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
00255 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
00256
00257 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
00258 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
00259
00260 if(_fieldID == 21) {
00261
00262
00263
00264 if( zmm>=800. && zmm < 2420. && rmm < 1000. ) {
00265 int tqz = (int)( (zmm-800.)/10. );
00266 bz += (((_BzQR[tr][tqz]*(1.-pz)+_BzQR[tr][tqz+1]*pz)*(1.-pr)
00267 +(_BzQR[tr+1][tqz]*(1.-pz)+_BzQR[tr+1][tqz+1]*pz)*pr)
00268 *(float)sin((double)(2.*phi+2./180.*(double)PI)));
00269 br += (((_BrQR[tr][tqz]*(1.-pz)+_BrQR[tr][tqz+1]*pz)*(1.-pr)
00270 +(_BrQR[tr+1][tqz]*(1.-pz)+_BrQR[tr+1][tqz+1]*pz)*pr)
00271 *(float)sin((double)(2.*phi+2./180.*(double)PI)));
00272 bphi += (((_BphiQR[tr][tqz]*(1.-pz)
00273 +_BphiQR[tr][tqz+1]*pz)*(1.-pr)
00274 +(_BphiQR[tr+1][tqz]*(1.-pz)
00275 +_BphiQR[tr+1][tqz+1]*pz)*pr)
00276 *(float)cos((double)(2.*phi+2./180.*(double)PI)));
00277 }
00278
00279
00280
00281 if(zmm<=-500. && zmm>-1520. && rmm<1000.) {
00282 int tqz = (int)((-zmm-500.)/20.);
00283 int tqr = (int)(tr/2.);
00284 pz = (pz+(float)( tz-(int)(tz/2.)*2. ))/2.;
00285 pr = ( pr + (float)(tr-tqr*2) )/2.;
00286 float f = 1.;
00287
00288
00289
00290
00291
00292
00293
00294 float phi_tmp = phi;
00295 if( phi_tmp < (PI/2.) && phi_tmp >= (-PI/2.) ) {
00296 phi_tmp = PI-phi_tmp;
00297 f =-1.;
00298 } else if( phi_tmp < -PI/2. ) {
00299 phi_tmp = 2.*PI+phi_tmp;
00300 }
00301 float pphi = ( phi_tmp-PI/2. )/(11.25/180.*PI);
00302 int tphi = (int)pphi;
00303 pphi -= (float)tphi;
00304 if (tphi >= 16) tphi -= 16;
00305
00306 bz += f*
00307 (((_BzQL[tphi][tqr][tqz]*(1.-pz)+_BzQL[tphi][tqr][tqz+1]*pz)
00308 *(1-pr)+(_BzQL[tphi][tqr+1][tqz]*(1.-pz)
00309 +_BzQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
00310 +((_BzQL[tphi+1][tqr][tqz]*(1.-pz)
00311 +_BzQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
00312 +(_BzQL[tphi+1][tqr+1][tqz]*(1.-pz)
00313 +_BzQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
00314 br += f*
00315 (((_BrQL[tphi][tqr][tqz]*(1.- pz)
00316 +_BrQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
00317 +(_BrQL[tphi][tqr+1][tqz]*(1.-pz)
00318 +_BrQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1-pphi)
00319 +((_BrQL[tphi+1][tqr][tqz]*(1.- pz)
00320 +_BrQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
00321 +(_BrQL[tphi+1][tqr+1][tqz]*(1.-pz)
00322 +_BrQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
00323 bphi += f*
00324 (((_BphiQL[tphi][tqr][tqz]*(1.- pz)
00325 +_BphiQL[tphi][tqr][tqz+1]*pz)*(1.-pr)
00326 +(_BphiQL[tphi][tqr+1][tqz]*(1.-pz)
00327 +_BphiQL[tphi][tqr+1][tqz+1]*pz)*pr)*(1.-pphi)
00328 +((_BphiQL[tphi+1][tqr][tqz]*(1.- pz)
00329 +_BphiQL[tphi+1][tqr][tqz+1]*pz)*(1.-pr)
00330 +(_BphiQL[tphi+1][tqr+1][tqz]*(1.-pz)
00331 +_BphiQL[tphi+1][tqr+1][tqz+1]*pz)*pr)*pphi);
00332 }
00333 }
00334 } else {
00335 bz=0.;
00336 br=0.;
00337 bphi=0.;
00338 }
00339 } else if(zmm>=-1520. && zmm<=2460. && rmm<=0. && rmm<=1740.){
00340 if(tz == 246) tz=tz-1;
00341 if(tr == 174) tr=tr-1;
00342 float pz= (zmm + 1520.)/10.- (float)tz;
00343 float pr= rmm/10.- (float)tr;
00344 bz = (_Bz[tr][tz]*(1.- pz)+_Bz[tr][tz+1]*pz)*(1.-pr)+
00345 (_Bz[tr+1][tz]*(1.-pz)+_Bz[tr+1][tz+1]*pz)*pr;
00346
00347 br = (_Br[tr][tz]*(1.-pz)+_Br[tr][tz+1]*pz)*(1.-pr)+
00348 (_Br[tr+1][tz]*(1.-pz)+_Br[tr+1][tz+1]*pz)*pr;
00349
00350 bphi = (_Bphi[tr][tz]*(1.-pz)+_Bphi[tr][tz+1]*pz)*(1.-pr)+
00351 (_Bphi[tr+1][tz]*(1.-pz)+_Bphi[tr+1][tz+1]*pz)*pr;
00352 } else {
00353 bz =0.;
00354 br =0.;
00355 bphi =0.;
00356 }
00357
00358
00359
00360 float Bmag_xy = (float)sqrt(br*br + bphi*bphi);
00361 double Bphi_rp = (float)atan2( (double)bphi, (double)br );
00362 m_Bx = Bmag_xy * (float)cos((double)phi + Bphi_rp)/1000.;
00363 m_By = Bmag_xy * (float)sin((double)phi + Bphi_rp)/1000.;
00364
00365
00366 m_Bz = bz/1000.;
00367 m_x = x;
00368 m_y = y;
00369 m_z = z;
00370 m_Bfld.setX((double) m_Bx);
00371 m_Bfld.setY((double) m_By);
00372 m_Bfld.setZ((double) m_Bz);
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405