00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00083
00084
00085
00087
00123 #include <string.h>
00124 #include <assert.h>
00125 #include "StarMagField.h"
00126 #include "StarCallf77.h"
00127 #include <string>
00128 #include "TMath.h"
00129 #define myMax(A,B) (((A)>(B))? (A):(B))
00130 #define myMin(A,B) (((A)<(B))? (A):(B))
00131 #define mySign(A,B) (((B)>=0)? fabs(A):-fabs(A))
00132
00133 StarMagField *StarMagField::fgInstance = 0;
00134
00135
00136 #define agufld F77_NAME(agufld,AGUFLD)
00137 #define gufld F77_NAME( gufld, GUFLD)
00138 #define mfldgeo F77_NAME(mfldgeo,MFLDGEO)
00139
00140 R__EXTERN "C" {
00141 Float_t type_of_call gufld(Float_t *x, Float_t *bf);
00142
00143 Float_t type_of_call agufld(Float_t *x, Float_t *bf) {
00144 bf[0] = bf[1] = bf[2] = 0;
00145 if (StarMagField::Instance())
00146 StarMagField::Instance()->BField(x,bf);
00147 else {
00148 printf("agufld:: request for non initialized mag.field, return 0\n");
00149 assert(StarMagField::Instance());
00150 }
00151 return 0;
00152 }
00153
00154 void type_of_call mfldgeo(float &factor) {
00155 if (StarMagField::Instance()) {
00156 printf("StarMagField mfldgeo: The field has been already instantiated.\n");
00157 } else {
00158 printf("StarMagField instantiate starsim field=%g\n",factor);
00159 (new StarMagField(StarMagField::kMapped,factor/5.))->SetLock();
00160 }
00161 float x[3]={0},b[3];
00162 gufld(x,b);
00163 printf("StarMagField:mfldgeo(%g) Bz=%g\n",factor,b[2]);
00164 }
00165 }
00166
00167
00168 struct BFLD_t {
00169 Int_t version;
00170 Char_t *code;
00171 Float_t date; Int_t kz; Float_t rmaxx, zmaxx, rrm, zz1, zz2;
00172 Float_t RmaxInn, ZmaxInn;
00173 Int_t nrp, nzp;
00174 };
00175 static const BFLD_t BFLD = {
00176 3 ,
00177 "opt1" ,
00178 22.10 ,
00179 22 ,
00180 270 ,
00181 290 ,
00182 400 ,
00183 270 ,
00184 800 ,
00185 264.265 ,
00186 312.500 ,
00187 200 ,
00188 800
00189 };
00190
00191 struct BDAT_t {
00192 Int_t N;
00193 Float_t Zi, Ri[20], Bzi[20], Bri[20];
00194 };
00195
00196 static const Int_t nZext = 23;
00197 static const BDAT_t BDAT[nZext] = {
00198 { 15 ,
00199 0.0 ,
00200 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00201 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00202 {4704.6,4768.0,4946.6,5148.3,5128.6,4927.3,4844.9,4830.1,
00203 4823.3,4858.0,5110.9,3402.4, -18.1, -13.6, -25.0 } ,
00204 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00205 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } },
00206
00207
00208 { 15 ,
00209 270.0 ,
00210 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00211 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00212 {4704.6,4768.0,4946.6,5148.3,5128.6,4927.3,4844.9,4830.1,
00213 4823.3,4858.0,5110.9,3402.4, -18.1, -13.6, -25.0 } ,
00214 { 0.0, 131.3, 188.4, 74.1,-148.6,-164.8, -53.2, 23.8,
00215 97.4, 213.7, 329.3, 75.3, 18.2, -44.3, -36.5 } },
00216 { 15 ,
00217 280.0 ,
00218 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00219 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00220 {4568.8,4649.8,4898.6,5241.5,5234.5,4883.9,4806.5,4802.4,
00221 4781.7,4771.8,5057.8,3504.2,-144.8, -15.0, -28.0 } ,
00222 { 0.0, 188.6, 297.5, 151.9,-241.2,-242.5, -60.1, 19.5,
00223 92.3, 244.3, 541.5, 396.8, 83.6, -49.9, -40.6 } },
00224 { 15 ,
00225 290.0 ,
00226 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00227 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00228 {4383.8,4478.4,4801.1,5378.7,5431.1,4771.4,4765.5,4778.2,
00229 4741.4,4651.9,4852.9,3684.4, 6.9, -16.9, -31.6 } ,
00230 { 0.0, 260.2, 456.5, 312.9,-414.5,-349.8, -51.7, 14.4,
00231 74.7, 234.0, 858.0, 726.3, 355.0, -56.5, -45.0 } },
00232 { 15 ,
00233 300.0 ,
00234 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00235 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00236 {4142.1,4240.8,4614.8,5546.4,5829.1,4450.0,4737.6,4761.4,
00237 4711.3,4534.1,4231.0,4067.5,-880.0, -19.3, -36.2 } ,
00238 { 0.0, 341.1, 669.5, 661.0,-766.7,-480.9, -24.5, 8.8,
00239 43.5, 149.9,1333.6, 999.3, 53.6, -64.2, -49.8 } },
00240 { 15 ,
00241 310.0 ,
00242 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
00243 200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } ,
00244 {3842.7,3930.2,4292.5,5589.1,6643.0,3236.8,4733.0,4755.1,
00245 4699.4,4485.0,1931.8,4782.0, 50.2, -22.8, -42.0 } ,
00246 { 0.0, 421.2, 915.6,1382.6,-1482.8,-1019.7, 1.2, 2.0,
00247 1.9, -2.3,2069.4, 791.7, 240.6, -73.6, -54.9 } },
00248 { 8 ,
00249 320.0 ,
00250 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 375.0, 400.0 } ,
00251 {3491.2,3552.1,3807.3,4923.7,7889.6,1983.9, -28.0, -49.4 } ,
00252 { 0.0, 485.7,1133.5,2502.8, -38.8,-174.8, -85.1, -60.0 } },
00253 { 7 ,
00254 330.0 ,
00255 { 0.0, 25.0, 50.0, 75.0, 250.0, 375.0, 400.0 } ,
00256 {3105.3,3127.0,3200.4,3268.9, -3.5, -36.6, -59.0 } ,
00257 { 0.0, 521.1,1246.1,3029.5,9199.2, -99.4, -64.5 } },
00258 { 6 ,
00259 340.0 ,
00260 { 0.0, 25.0, 50.0, 75.0, 375.0, 400.0 } ,
00261 {2706.4,2686.8,2574.5,1826.7, -51.8, -71.0 } ,
00262 { 0.0, 520.6,1218.1,2485.3,-116.9, -67.3 } },
00263 { 6 ,
00264 350.0 ,
00265 { 0.0, 25.0, 50.0, 75.0, 375.0, 400.0 } ,
00266 {2317.7,2264.6,2026.3,1142.6, -80.8, -85.1 } ,
00267 { 0.0, 487.6,1082.3,1787.2,-133.8, -67.0 } },
00268 { 8 ,
00269 360.0 ,
00270 { 0.0, 25.0, 50.0, 75.0, 100.0, 250.0, 375.0, 400.0 } ,
00271 {1958.5,1885.6,1595.6, 829.2,-563.7,4895.8,-127.6, -99.8 } ,
00272 { 0.0, 432.4, 901.7,1265.8, 788.0,9507.4,-134.0, -62.2 } },
00273 { 17 ,
00274 370.0 ,
00275 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
00276 225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } ,
00277 {1637.8,1562.2,1276.7, 678.9, 15.7, 251.6, 384.9, 503.7, 683.3,
00278 1087.1,1868.1,-1320.5,-593.9,-391.5,-345.9,-168.2,-112.9 } ,
00279 { 0.0, 367.3, 720.6, 900.1, 421.6, 60.4, 37.1, 44.5, 79.7,
00280 229.6,2339.4, 654.6, 114.6, 35.9, -30.0,-101.8, -52.4 } },
00281 { 17 ,
00282 380.0 ,
00283 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
00284 225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } ,
00285 {1373.5,1296.6,1045.5, 603.2, 221.4, 278.6, 382.7, 488.2, 638.7,
00286 892.4, 708.6,-709.9,-515.0,-364.7,-293.1,-181.5,-122.1 } ,
00287 { 0.0, 302.3, 563.3, 650.3, 369.7, 120.0, 79.6, 96.2, 169.1,
00288 430.1,1454.7, 860.7, 228.6, 77.5, -10.8, -60.2, -39.4 } },
00289 { 17 ,
00290 390.0 ,
00291 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
00292 225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } ,
00293 {1151.2,1083.6, 877.2, 557.6, 308.9, 305.2, 377.6, 463.3, 573.3,
00294 684.5, 377.5,-376.2,-415.2,-326.0,-258.2,-179.8,-126.9 } ,
00295 { 0.0, 243.7, 437.7, 486.1, 319.9, 155.1, 115.2, 139.4, 232.4,
00296 494.6,1019.1, 751.4, 289.6, 112.2, 19.4, -26.7, -25.0 } },
00297 { 17 ,
00298 400.0 ,
00299 { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
00300 225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } ,
00301 {971.6, 914.8, 751.6, 520.8, 348.0, 323.7, 369.1, 432.0, 500.9,
00302 520.4, 251.2,-214.3,-320.8,-282.3,-230.2,-171.7,-127.7 } ,
00303 { 0.0, 194.5, 341.1, 375.8, 277.5, 171.7, 142.1, 172.1, 269.4,
00304 486.6, 769.0, 624.1, 308.0, 137.2, 44.9, -1.4, -11.2 } },
00305 { 9 ,
00306 450.0 ,
00307 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00308 {481.5, 423.2, 325.1, 283.8, 242.6, 88.2, -85.2,-123.1,-103.9 } ,
00309 { 0.0, 119.3, 157.5, 174.6, 248.4, 314.8, 220.8, 95.4, 32.3 } },
00310 { 9 ,
00311 500.0 ,
00312 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00313 {291.3, 273.2, 234.4, 192.6, 136.7, 53.6, -25.6, -64.6, -72.5 } ,
00314 { 0.0, 60.7, 103.2, 135.9, 168.1, 177.4, 140.6, 84.7, 41.9 } },
00315 { 9 ,
00316 550.0 ,
00317 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00318 {190.4, 181.9, 159.6, 127.7, 86.0, 37.1, -7.3, -35.9, -48.9 } ,
00319 { 0.0, 37.8, 69.7, 94.3, 110.3, 110.8, 92.6, 64.2, 37.3 } },
00320 { 9 ,
00321 600.0 ,
00322 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00323 {126.9, 121.7, 107.1, 84.9, 56.8, 26.4, -1.2, -21.2, -32.9 } ,
00324 { 0.0, 25.0, 46.9, 63.4, 72.6, 72.2, 62.3, 46.3, 29.2 } },
00325 { 9 ,
00326 650.0 ,
00327 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00328 { 84.8, 81.4, 71.7, 56.8, 38.3, 18.7, 0.8, -13.1, -22.2 } ,
00329 { 0.0, 16.7, 31.4, 42.4, 48.2, 48.1, 42.4, 32.8, 21.6 } },
00330 { 9 ,
00331 700.0 ,
00332 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00333 { 56.7, 54.4, 47.9, 38.0, 25.9, 13.1, 1.3, -8.3, -15.0 } ,
00334 { 0.0, 11.2, 21.1, 28.4, 32.4, 32.5, 29.1, 23.0, 15.5 } },
00335 { 9 ,
00336 750.0 ,
00337 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00338 { 37.7, 36.2, 31.9, 25.4, 17.5, 9.1, 1.2, -5.4, -10.1 } ,
00339 { 0.0, 7.6, 14.3, 19.3, 22.0, 22.2, 20.1, 16.1, 11.0 } },
00340 { 9 ,
00341 800.0 ,
00342 { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } ,
00343 { 24.8, 23.8, 21.0, 16.8, 11.6, 6.1, 0.9, -3.5, -6.7 } ,
00344 { 0.0, 5.2, 9.8, 13.3, 15.2, 15.4, 14.1, 11.4, 7.9 } },
00345 };
00346
00347 StarMagField::StarMagField ( EBField map, Float_t factor,
00348 Bool_t lock, Float_t rescale,
00349 Float_t BDipole, Float_t RmaxDip,
00350 Float_t ZminDip, Float_t ZmaxDip) :
00351 fMap(map),
00352 fFactor(factor), fRescale(rescale),
00353 fBDipole(BDipole), fRmaxDip(RmaxDip),
00354 fZminDip(ZminDip), fZmaxDip(ZmaxDip),
00355 fLock(lock)
00356 {
00357 if (fgInstance) {
00358 printf("Cannot initialise twice StarMagField class\n");
00359 assert(0);
00360 }
00361 fgInstance = this;
00362 if (fMap == kUndefined) {
00363 printf("StarMagField is instantiated with predefined factor %f and map %i\n",fFactor, fMap);
00364 } else {
00365 if (fLock) printf("StarMagField is locked, no modification from DB will be accepted\n");
00366 }
00367 ReadField() ;
00368 }
00369
00371 void StarMagField::BField( const Double_t x[], Double_t B[] ) {
00372 Float_t xx[3] = {x[0], x[1], x[2]};
00373 Float_t bb[3];
00374 BField(xx,bb);
00375 B[0] = bb[0]; B[1] = bb[1]; B[2] = bb[2];
00376 }
00377
00378
00379
00380
00381 void StarMagField::BField( const Float_t x[], Float_t B[] )
00382
00383 {
00384
00385 Float_t r, z, Br_value, Bz_value ;
00386 Float_t phi, Bphi_value,phi1;
00387 Bphi_value=0;
00388 Br_value = Bz_value = 0;
00389 B[0] = B[1] = B[2] = 0;
00390 z = x[2] ;
00391 r = sqrt( x[0]*x[0] + x[1]*x[1] ) ;
00392 phi = atan2( x[1], x[0] ) ;
00393 if ( phi < 0 ) phi += 2*TMath::Pi() ;
00394
00395
00396 Float_t za = fabs(z);
00397 if (za > fZminDip && za < fZmaxDip && r < fRmaxDip) {
00398 B[1] = mySign(fBDipole, z);
00399 B[2] = fabs(B[1]/1000.);
00400 return;
00401 }
00402 if (z >= ZList[0] && z <= ZList[nZ-1] && r <= Radius[nR-1]) {
00403 Interpolate2DBfield( r, z, Br_value, Bz_value ) ;
00404 B[2] = Bz_value ;
00405 if ( r != 0.0 ) {
00406 B[0] = Br_value * (x[0]/r) ;
00407 B[1] = Br_value * (x[1]/r) ;
00408 }
00409 return;
00410 }
00411
00412
00413
00414
00415 if (za <=342.20 && r>=303.29 && r <= 364.25) {
00416
00417 phi1=phi*TMath::RadToDeg();
00418 if(phi1>12) phi1=phi1-int(phi1/12)*12;
00419
00420 Interpolate3DBSteelfield( r, za, phi1, Br_value, Bz_value, Bphi_value ) ;
00421 B[0] = Br_value * (x[0]/r) - Bphi_value * (x[1]/r) ;
00422 B[1] = Br_value * (x[1]/r) + Bphi_value * (x[0]/r) ;
00423 B[2] = Bz_value ;
00424
00425 if(z<0) {
00426 B[0]=-B[0];
00427 B[1]=-B[1];
00428 }
00429
00430
00431
00432 return;
00433 }
00434
00435
00436
00437
00438
00439
00440 Interpolate2ExtDBfield( r, z, Br_value, Bz_value ) ;
00441 if (za <= BFLD.zmaxx && r <= BFLD.rmaxx) {
00442 static const Float_t zero = 0;
00443 static const Float_t one = 1;
00444 Float_t wz = (za - ZList[nZ-1] )/(BFLD.zmaxx - ZList[nZ-1]);
00445 Float_t wr = (r - Radius[nR-1])/(BFLD.rmaxx - Radius[nR-1]);
00446 Float_t w = myMin(myMax(zero,myMax(wz,wr)),one);
00447 Float_t rm = myMin(r,Radius[nR-1]);
00448 Float_t zm = mySign(myMin(za,ZList[nZ-1]),z);
00449 Float_t BrI, BzI;
00450 Interpolate2DBfield( rm, zm, BrI, BzI ) ;
00451 Br_value = (1-w)*BrI + w*Br_value;
00452 Bz_value = (1-w)*BzI + w*Bz_value;
00453 }
00454 B[2] = Bz_value ;
00455 if ( r != 0.0 ) {
00456 B[0] = Br_value * (x[0]/r) ;
00457 B[1] = Br_value * (x[1]/r) ;
00458 }
00459
00460
00461 return;
00462
00463
00464
00465
00466
00467 }
00468
00469
00470
00471
00473 void StarMagField::B3DField( const Float_t x[], Float_t B[] )
00474 {
00475 Float_t r, z, phi, Br_value, Bz_value, Bphi_value ;
00476 Bphi_value=0;
00477 Br_value = Bz_value = 0;
00478 B[0] = B[1] = B[2] = 0;
00479 #if 0
00480 Float_t phi1;
00481 #endif
00482 z = x[2] ;
00483 r = sqrt( x[0]*x[0] + x[1]*x[1] ) ;
00484
00485 if ( r != 0.0 )
00486 {
00487 phi = TMath::ATan2( x[1], x[0] ) ;
00488 if ( phi < 0 ) phi += 2*TMath::Pi() ;
00489 #if 0
00490
00491 phi1=phi*TMath::RadToDeg();
00492
00493
00494 Interpolate3DBfield( r, z, phi1, Br_value, Bz_value, Bphi_value ) ;
00495 #else
00496 Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
00497 #endif
00498 B[0] = Br_value * (x[0]/r) - Bphi_value * (x[1]/r) ;
00499 B[1] = Br_value * (x[1]/r) + Bphi_value * (x[0]/r) ;
00500 B[2] = Bz_value ;
00501 }
00502 else
00503 {
00504 phi = 0 ;
00505 Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
00506 B[0] = Br_value ;
00507 B[1] = Bphi_value ;
00508 B[2] = Bz_value ;
00509 }
00510
00511 return ;
00512
00513 }
00514
00515
00516
00517
00518
00519
00520
00522
00523 void StarMagField::BrBzField( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value )
00524
00525 {
00526
00527
00528 Br_value = Bz_value = 0;
00529 if(r>0) Interpolate2DBfield( r, z, Br_value, Bz_value ) ;
00530 return;
00531
00532 }
00533
00534
00535
00536
00537 void StarMagField::BrBz3DField( const Float_t r, const Float_t z, const Float_t phi,
00538 Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
00539
00540 {
00541
00542 Bphi_value=0;
00543 Br_value = Bz_value = 0;
00544
00545 #if 0
00546
00547 Float_t phiprime ;
00548
00549 phiprime = phi ;
00550 if ( phiprime < 0 ) phiprime += 2*TMath::Pi() ;
00551
00552 phiprime=phiprime*TMath::RadToDeg();
00553
00554 #endif
00555
00556 if(r>0) {
00557 #if 0
00558 Interpolate3DBfield( r, z, phiprime, Br_value, Bz_value, Bphi_value ) ;
00559 #else
00560 Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
00561 #endif
00562 }
00563
00564 return;
00565
00566 }
00567
00568
00569
00570
00572
00573 void StarMagField::ReadField( )
00574
00575 {
00576
00577 FILE *magfile, *b3Dfile ;
00578 std::string comment, filename, filename3D ;
00579 std::string MapLocation ;
00580 std::string BaseLocation = getenv("STAR") ;
00581 BaseLocation += "/StarDb/StMagF/" ;
00582
00583 if ( fMap == kMapped )
00584 {
00585 if ( fabs(fFactor) > 0.8 )
00586 {
00587 if ( fFactor > 0 )
00588 {
00589 filename = "bfield_full_positive_2D.dat" ;
00590 filename3D = "bfield_full_positive_3D.dat" ;
00591 comment = "Measured Full Field" ;
00592 fRescale = 1 ;
00593 }
00594 else
00595 {
00596 filename = "bfield_full_negative_2D.dat" ;
00597 filename3D = "bfield_full_negative_3D.dat" ;
00598 comment = "Measured Full Field Reversed" ;
00599 fRescale = -1 ;
00600 }
00601 }
00602 else
00603 {
00604 filename = "bfield_half_positive_2D.dat" ;
00605 filename3D = "bfield_half_positive_3D.dat" ;
00606 comment = "Measured Half Field" ;
00607 fRescale = 2 ;
00608 }
00609 }
00610 else if ( fMap == kConstant )
00611 {
00612 filename = "const_full_positive_2D.dat" ;
00613 comment = "Constant Full Field" ;
00614 fRescale = 1 ;
00615 }
00616 else
00617 {
00618 fprintf(stderr,"StarMagField::ReadField No map available - you must choose a mapped field or a constant field\n");
00619 exit(1) ;
00620 }
00621
00622 printf("StarMagField::ReadField Reading Magnetic Field %s, Scale factor = %f \n",comment.c_str(),fFactor);
00623 printf("StarMagField::ReadField Filename is %s, Adjusted Scale factor = %f \n",filename.c_str(),fFactor*fRescale);
00624
00625 MapLocation = BaseLocation + filename ;
00626 magfile = fopen(MapLocation.c_str(),"r") ;
00627 printf("StarMagField::ReadField Reading 2D Magnetic Field file: %s \n",filename.c_str());
00628
00629 if (magfile)
00630
00631 {
00632 Char_t cname[128] ;
00633 fgets ( cname, sizeof(cname) , magfile ) ;
00634 fgets ( cname, sizeof(cname) , magfile ) ;
00635 fgets ( cname, sizeof(cname) , magfile ) ;
00636 fgets ( cname, sizeof(cname) , magfile ) ;
00637 fgets ( cname, sizeof(cname) , magfile ) ;
00638
00639 for ( Int_t j=0 ; j < nZ ; j++ )
00640 {
00641 for ( Int_t k=0 ; k < nR ; k++ )
00642 {
00643 fgets ( cname, sizeof(cname) , magfile ) ;
00644 sscanf ( cname, " %f %f %f %f ", &Radius[k], &ZList[j], &Br[j][k], &Bz[j][k] ) ;
00645 }
00646 }
00647 }
00648
00649 else
00650
00651 {
00652 fprintf(stderr,"StarMagField::ReadField File %s not found !\n",MapLocation.c_str());
00653 exit(1);
00654 }
00655
00656 fclose(magfile) ;
00657
00658 MapLocation = BaseLocation + filename3D ;
00659 b3Dfile = fopen(MapLocation.c_str(),"r") ;
00660 printf("StarMagField::ReadField Reading 3D Magnetic Field file: %s \n",filename3D.c_str());
00661
00662 if (b3Dfile)
00663
00664 {
00665 Char_t cname[128] ;
00666 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00667 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00668 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00669 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00670 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00671 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00672
00673 for ( Int_t i=0 ; i < nPhi ; i++ )
00674 {
00675 for ( Int_t j=0 ; j < nZ ; j++ )
00676 {
00677 for ( Int_t k=0 ; k < nR ; k++ )
00678 {
00679 fgets ( cname, sizeof(cname) , b3Dfile ) ;
00680 sscanf ( cname, " %f %f %f %f %f %f ",
00681 &R3D[k], &Z3D[j], &Phi3D[i], &Br3D[i][j][k], &Bz3D[i][j][k], &Bphi3D[i][j][k] ) ;
00682 Phi3D[i] *= TMath::Pi() / 180. ;
00683 }
00684 }
00685 }
00686 }
00687
00688 else if ( fMap == kConstant )
00689
00690 {
00691 for ( Int_t i=0 ; i < nPhi ; i++ )
00692 {
00693 for ( Int_t j=0 ; j < nZ ; j++ )
00694 {
00695 for ( Int_t k=0 ; k < nR ; k++ )
00696 {
00697 Br3D[i][j][k] = Br[j][k] ;
00698 Bz3D[i][j][k] = Bz[j][k] ;
00699 Bphi3D[i][j][k] = 0 ;
00700 }
00701 }
00702 }
00703 }
00704
00705 else
00706
00707 {
00708 fprintf(stderr,"StarMagField::ReadField File %s not found !\n",MapLocation.c_str());
00709 exit(1);
00710 }
00711
00712 fclose(b3Dfile) ;
00713 #if 1
00714
00715
00716
00717
00718
00719
00720
00721 MapLocation = BaseLocation + "steel_magfieldmap.dat";
00722 magfile = fopen(MapLocation.c_str(),"r") ;
00723 if (magfile) {
00724 printf("StarMagField::ReadField Reading 3D Magnetic Field file: %s \n",filename.c_str());
00725 Char_t cname[128] ;
00726 for (;;) {
00727 fgets ( cname, sizeof(cname) , magfile ) ;
00728 if (cname[0] == '#') continue;
00729 break;
00730 }
00731 for ( Int_t i=0 ; i < nPhiSteel ; i++ )
00732 {
00733
00734 for ( Int_t k=0 ; k < nRSteel ; k++ )
00735 {
00736 for ( Int_t j=0 ; j < nZSteel ; j++ )
00737 {
00738
00739
00740 fgets ( cname, sizeof(cname) , magfile ) ;
00741 sscanf ( cname, " %f %f %f %f %f %f ",
00742 &R3DSteel[k], &Z3DSteel[j], &Phi3DSteel[i], &Bx3DSteel[i][j][k], &Bz3DSteel[i][j][k], &By3DSteel[i][j][k] ) ;
00743
00744
00745 Br3DSteel[i][j][k]=cos(Phi3DSteel[i]*TMath::DegToRad())*Bx3DSteel[i][j][k]+sin(Phi3DSteel[i]*TMath::DegToRad())*By3DSteel[i][j][k];
00746
00747 Bphi3DSteel[i][j][k]=0-sin(Phi3DSteel[i]*TMath::DegToRad())*Bx3DSteel[i][j][k]+cos(Phi3DSteel[i]*TMath::DegToRad())*By3DSteel[i][j][k];
00748
00749
00750
00751
00752
00753
00754
00755
00756 }
00757 }
00758 }
00759 fclose(magfile);
00760 }
00761 #endif
00762 return ;
00763
00764 }
00765
00766
00767
00768
00770
00771 void StarMagField::Interpolate2DBfield( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value )
00772
00773 {
00774
00775 Float_t fscale ;
00776
00777 fscale = 0.001*fFactor*fRescale ;
00778
00779
00780 const Int_t ORDER = 1 ;
00781 static Int_t jlow=0, klow=0 ;
00782 Float_t save_Br[ORDER+1] ;
00783 Float_t save_Bz[ORDER+1] ;
00784
00785 Search ( nZ, ZList, z, jlow ) ;
00786 Search ( nR, Radius, r, klow ) ;
00787 if ( jlow < 0 ) jlow = 0 ;
00788 if ( klow < 0 ) klow = 0 ;
00789 if ( jlow + ORDER >= nZ - 1 ) jlow = nZ - 1 - ORDER ;
00790 if ( klow + ORDER >= nR - 1 ) klow = nR - 1 - ORDER ;
00791
00792 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
00793 {
00794 save_Br[j-jlow] = Interpolate( &Radius[klow], &Br[j][klow], ORDER, r ) ;
00795 save_Bz[j-jlow] = Interpolate( &Radius[klow], &Bz[j][klow], ORDER, r ) ;
00796 }
00797 Br_value = fscale * Interpolate( &ZList[jlow], save_Br, ORDER, z ) ;
00798 Bz_value = fscale * Interpolate( &ZList[jlow], save_Bz, ORDER, z ) ;
00799
00800 }
00801
00802 void StarMagField::Interpolate2ExtDBfield( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value ) {
00803 static Float_t ZExtList[nZext];
00804 static Bool_t first = kTRUE;
00805 if (first) {
00806 for (Int_t j = 0; j < nZext; j++) ZExtList[j] = BDAT[j].Zi;
00807 first = kFALSE;
00808 }
00809 Float_t za = fabs(z);
00810 if (za > BFLD.zz2 || r > BFLD.rrm) return;
00811 if (za < ZList[nZ-1] && r < Radius[nR-1]) return;
00812
00813
00814 if (za <=342.20 && r>=303.29 && r <= 363.29) return;
00815
00816
00817
00818 Float_t fscale = 0.001*fFactor;
00819
00820 const Int_t ORDER = 1 ;
00821 static Int_t jlow=0, klow=0 ;
00822 Float_t save_Br[ORDER+1] ;
00823 Float_t save_Bz[ORDER+1] ;
00824 Search ( nZext, ZExtList, za, jlow ) ;
00825 if ( jlow < 0 ) jlow = 0 ;
00826 if ( jlow + ORDER >= nZext - 1 ) jlow = nZext - 1 - ORDER ;
00827
00828 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ ) {
00829 Int_t N = BDAT[j].N;
00830 Search ( N, (Float_t *) (&BDAT[j].Ri[0]), r, klow ) ;
00831 if ( klow < 0 ) klow = 0 ;
00832 if ( klow + ORDER >= BDAT[j].N - 1 ) klow = BDAT[j].N - 1 - ORDER ;
00833 save_Br[j-jlow] = Interpolate( &BDAT[j].Ri[klow], &BDAT[j].Bri[klow], ORDER, r ) ;
00834 save_Bz[j-jlow] = Interpolate( &BDAT[j].Ri[klow], &BDAT[j].Bzi[klow], ORDER, r ) ;
00835 }
00836 Br_value = fscale * Interpolate( &ZExtList[jlow], save_Br, ORDER, za ) ;
00837 Bz_value = fscale * Interpolate( &ZExtList[jlow], save_Bz, ORDER, za ) ;
00838 if (z < 0) Br_value = - Br_value;
00839 }
00840
00842
00843 void StarMagField::Interpolate3DBfield( const Float_t r, const Float_t z, const Float_t phi,
00844 Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
00845 {
00846
00847 Float_t fscale ;
00848
00849 fscale = 0.001*fFactor*fRescale ;
00850
00851 const Int_t ORDER = 1 ;
00852 static Int_t ilow=0, jlow=0, klow=0 ;
00853 Float_t save_Br[ORDER+1], saved_Br[ORDER+1] ;
00854 Float_t save_Bz[ORDER+1], saved_Bz[ORDER+1] ;
00855 Float_t save_Bphi[ORDER+1], saved_Bphi[ORDER+1] ;
00856
00857
00858
00859
00860 if(r<0) return;
00861
00862 Search( nPhi, Phi3D, phi, ilow ) ;
00863 Search( nZ, Z3D, z, jlow ) ;
00864 Search( nR, R3D, r, klow ) ;
00865 if ( ilow < 0 ) ilow = 0 ;
00866 if ( jlow < 0 ) jlow = 0 ;
00867 if ( klow < 0 ) klow = 0 ;
00868
00869 if ( ilow + ORDER >= nPhi - 1 ) ilow = nPhi - 1 - ORDER ;
00870 if ( jlow + ORDER >= nZ - 1 ) jlow = nZ - 1 - ORDER ;
00871 if ( klow + ORDER >= nR - 1 ) klow = nR - 1 - ORDER ;
00872
00873 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
00874 {
00875 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
00876 {
00877 save_Br[j-jlow] = Interpolate( &R3D[klow], &Br3D[i][j][klow], ORDER, r ) ;
00878 save_Bz[j-jlow] = Interpolate( &R3D[klow], &Bz3D[i][j][klow], ORDER, r ) ;
00879 save_Bphi[j-jlow] = Interpolate( &R3D[klow], &Bphi3D[i][j][klow], ORDER, r ) ;
00880 }
00881 saved_Br[i-ilow] = Interpolate( &Z3D[jlow], save_Br, ORDER, z ) ;
00882 saved_Bz[i-ilow] = Interpolate( &Z3D[jlow], save_Bz, ORDER, z ) ;
00883 saved_Bphi[i-ilow] = Interpolate( &Z3D[jlow], save_Bphi, ORDER, z ) ;
00884 }
00885 Br_value = fscale * Interpolate( &Phi3D[ilow], saved_Br, ORDER, phi ) ;
00886 Bz_value = fscale * Interpolate( &Phi3D[ilow], saved_Bz, ORDER, phi ) ;
00887 Bphi_value = fscale * Interpolate( &Phi3D[ilow], saved_Bphi, ORDER, phi ) ;
00888
00889 }
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00902
00903 void StarMagField::Interpolate3DBSteelfield( const Float_t r, const Float_t z, const Float_t phi,
00904 Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
00905 {
00906
00907 Float_t fscale ;
00908
00909
00910
00911
00912 fscale = 0.001*fFactor;
00913
00914 const Int_t ORDER = 1 ;
00915 static Int_t ilow=0, jlow=0, klow=0 ;
00916 Float_t save_Br[ORDER+1], saved_Br[ORDER+1] ;
00917 Float_t save_Bz[ORDER+1], saved_Bz[ORDER+1] ;
00918 Float_t save_Bphi[ORDER+1], saved_Bphi[ORDER+1] ;
00919
00920
00921 Search( nPhiSteel, Phi3DSteel, phi, ilow ) ;
00922 Search( nZSteel, Z3DSteel, z, jlow ) ;
00923 Search( nRSteel, R3DSteel, r, klow ) ;
00924 if ( ilow < 0 ) ilow = 0 ;
00925 if ( jlow < 0 ) jlow = 0 ;
00926 if ( klow < 0 ) klow = 0 ;
00927
00928 if ( ilow + ORDER >= nPhiSteel - 1 ) ilow = nPhiSteel - 1 - ORDER ;
00929 if ( jlow + ORDER >= nZSteel - 1 ) jlow = nZSteel - 1 - ORDER ;
00930 if ( klow + ORDER >= nRSteel - 1 ) klow = nRSteel - 1 - ORDER ;
00931
00932 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
00933 {
00934 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
00935 {
00936 save_Br[j-jlow] = Interpolate( &R3DSteel[klow], &Br3DSteel[i][j][klow], ORDER, r ) ;
00937 save_Bz[j-jlow] = Interpolate( &R3DSteel[klow], &Bz3DSteel[i][j][klow], ORDER, r ) ;
00938 save_Bphi[j-jlow] = Interpolate( &R3DSteel[klow], &Bphi3DSteel[i][j][klow], ORDER, r ) ;
00939 }
00940 saved_Br[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Br, ORDER, z ) ;
00941 saved_Bz[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Bz, ORDER, z ) ;
00942 saved_Bphi[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Bphi, ORDER, z ) ;
00943 }
00944 Br_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Br, ORDER, phi ) ;
00945 Bz_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Bz, ORDER, phi ) ;
00946 Bphi_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Bphi, ORDER, phi ) ;
00947
00948 }
00949
00950
00951
00952
00953
00954
00955
00956
00957 #if 0
00959
00960 void StarMagField::Interpolate2DEdistortion( const Float_t r, const Float_t z,
00961 const Float_t Er[neZ][neR], Float_t &Er_value )
00962
00963 {
00964
00965 const Int_t ORDER = 1 ;
00966 static Int_t jlow=0, klow=0 ;
00967 Float_t save_Er[ORDER+1] ;
00968
00969 Search( neZ, eZList, z, jlow ) ;
00970 Search( neR, eRadius, r, klow ) ;
00971 if ( jlow < 0 ) jlow = 0 ;
00972 if ( klow < 0 ) klow = 0 ;
00973 if ( jlow + ORDER >= neZ - 1 ) jlow = neZ - 1 - ORDER ;
00974 if ( klow + ORDER >= neR - 1 ) klow = neR - 1 - ORDER ;
00975
00976 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
00977 {
00978 save_Er[j-jlow] = Interpolate( &eRadius[klow], &Er[j][klow], ORDER, r ) ;
00979 }
00980 Er_value = Interpolate( &eZList[jlow], save_Er, ORDER, z ) ;
00981
00982 }
00983
00985
00986 void StarMagField::Interpolate3DEdistortion( const Float_t r, const Float_t phi, const Float_t z,
00987 const Float_t Er[neZ][nePhi][neR], const Float_t Ephi[neZ][nePhi][neR],
00988 Float_t &Er_value, Float_t &Ephi_value )
00989
00990 {
00991
00992 const Int_t ORDER = 1 ;
00993 static Int_t ilow=0, jlow=0, klow=0 ;
00994 Float_t save_Er[ORDER+1], saved_Er[ORDER+1] ;
00995 Float_t save_Ephi[ORDER+1], saved_Ephi[ORDER+1] ;
00996
00997 Search( neZ, eZList, z, ilow ) ;
00998 Search( nePhi, ePhiList, phi, jlow ) ;
00999 Search( neR, eRadius, r, klow ) ;
01000 if ( ilow < 0 ) ilow = 0 ;
01001 if ( jlow < 0 ) jlow = 0 ;
01002 if ( klow < 0 ) klow = 0 ;
01003
01004 if ( ilow + ORDER >= neZ - 1 ) ilow = neZ - 1 - ORDER ;
01005 if ( jlow + ORDER >= nePhi - 1 ) jlow = nePhi - 1 - ORDER ;
01006 if ( klow + ORDER >= neR - 1 ) klow = neR - 1 - ORDER ;
01007
01008 for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
01009 {
01010 for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
01011 {
01012 save_Er[j-jlow] = Interpolate( &eRadius[klow], &Er[i][j][klow], ORDER, r ) ;
01013 save_Ephi[j-jlow] = Interpolate( &eRadius[klow], &Ephi[i][j][klow], ORDER, r ) ;
01014 }
01015 saved_Er[i-ilow] = Interpolate( &ePhiList[jlow], save_Er, ORDER, phi ) ;
01016 saved_Ephi[i-ilow] = Interpolate( &ePhiList[jlow], save_Ephi, ORDER, phi ) ;
01017 }
01018 Er_value = Interpolate( &eZList[ilow], saved_Er, ORDER, z ) ;
01019 Ephi_value = Interpolate( &eZList[ilow], saved_Ephi, ORDER, z ) ;
01020
01021 }
01022 #endif
01023
01024
01025
01027
01028 Float_t StarMagField::Interpolate( const Float_t Xarray[], const Float_t Yarray[],
01029 const Int_t ORDER, const Float_t x )
01030
01031 {
01032
01033 Float_t y ;
01034
01035
01036 if ( ORDER == 2 )
01037
01038 {
01039 y = (x-Xarray[1]) * (x-Xarray[2]) * Yarray[0] / ( (Xarray[0]-Xarray[1]) * (Xarray[0]-Xarray[2]) ) ;
01040 y += (x-Xarray[2]) * (x-Xarray[0]) * Yarray[1] / ( (Xarray[1]-Xarray[2]) * (Xarray[1]-Xarray[0]) ) ;
01041 y += (x-Xarray[0]) * (x-Xarray[1]) * Yarray[2] / ( (Xarray[2]-Xarray[0]) * (Xarray[2]-Xarray[1]) ) ;
01042
01043 }
01044
01045 else
01046
01047 {
01048 y = Yarray[0] + ( Yarray[1]-Yarray[0] ) * ( x-Xarray[0] ) / ( Xarray[1] - Xarray[0] ) ;
01049 }
01050
01051 return (y) ;
01052
01053 }
01054
01055
01056
01057
01059
01060 void StarMagField::Search( Int_t N, const Float_t Xarray[], Float_t x, Int_t &low )
01061
01062 {
01063 assert(! TMath::IsNaN(x));
01064 Long_t middle, high ;
01065 Int_t ascend = 0, increment = 1 ;
01066
01067 if ( Xarray[N-1] >= Xarray[0] ) ascend = 1 ;
01068
01069 if ( low < 0 || low > N-1 ) { low = -1 ; high = N ; }
01070
01071 else
01072 {
01073 if ( (Int_t)( x >= Xarray[low] ) == ascend )
01074 {
01075 if ( low == N-1 ) return ;
01076 high = low + 1 ;
01077 while ( (Int_t)( x >= Xarray[high] ) == ascend )
01078 {
01079 low = high ;
01080 increment *= 2 ;
01081 high = low + increment ;
01082 if ( high > N-1 ) { high = N ; break ; }
01083 }
01084 }
01085 else
01086 {
01087 if ( low == 0 ) { low = -1 ; return ; }
01088 high = low - 1 ;
01089 while ( (Int_t)( x < Xarray[low] ) == ascend )
01090 {
01091 high = low ;
01092 increment *= 2 ;
01093 if ( increment >= high ) { low = -1 ; break ; }
01094 else low = high - increment ;
01095 }
01096 }
01097 }
01098
01099 while ( (high-low) != 1 )
01100 {
01101 middle = ( high + low ) / 2 ;
01102 if ( (Int_t)( x >= Xarray[middle] ) == ascend )
01103 low = middle ;
01104 else
01105 high = middle ;
01106 }
01107
01108 if ( x == Xarray[N-1] ) low = N-2 ;
01109 if ( x == Xarray[0] ) low = 0 ;
01110
01111 }
01112
01113 #define PPLOCK(A) \
01114 void StarMagField::Set ## A (Float_t m) { \
01115 if (!fLock) f ## A = m; \
01116 else printf("StarMagField::Set"#A"() "#A" is locked at %f; Set to %f is ignored\n", f ## A ,m); \
01117 }
01118 PPLOCK(Factor)
01119 PPLOCK(Rescale)
01120 PPLOCK(BDipole)
01121 PPLOCK(RmaxDip)
01122 PPLOCK(ZminDip)
01123 PPLOCK(ZmaxDip)
01124 #undef PPLOCK
01125
01126 void StarMagField::SetLock () {
01127 if (! fLock) {
01128 fLock = kTRUE;
01129 printf("StarMagField::SetLock lock StarMagField parameters\n");
01130 Print();
01131 }
01132 }
01133
01134 #define PrintPar(A) printf("StarMagField:: "#A"\t%f\n",f ## A)
01135 void StarMagField::Print () {
01136 if (fLock) printf("StarMagField parameters are locked\n");
01137 printf("StarMagField:: Map\t%i\n",fMap );
01138 PrintPar(Factor );
01139 PrintPar(Rescale);
01140 PrintPar(BDipole);
01141 PrintPar(RmaxDip);
01142 PrintPar(ZminDip);
01143 PrintPar(ZmaxDip);
01144 }
01145 #undef PrintPar
