STEER/AliMagFDM.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //   Field with Magnetic Field map
  20 //   Used by AliRun class
  21 //   Author:
  22 //-------------------------------------------------------------------------
  23
  24 #include "TSystem.h"
  25
  26 #include "AliMagFDM.h"
  27
  28 ClassImp(AliMagFDM)
  29
  30 //_______________________________________________________________________
  31 AliMagFDM::AliMagFDM():
  32   fSolenoid(0),
  33   fInd(0),
  34   fZmin(0),
  35   fZmax(0),
  36   fYmax(0),
  37   fYmin(0),
  38   fZpmx(0),
  39   fZpmn(0),
  40   fRmax(0),
  41   fRmin(0),
  42   fXdel(0),
  43   fYdel(0),
  44   fZdel(0),
  45   fRdel(0),
  46   fPhid(0),
  47   fZpdl(0),
  48   fCx1(0),
  49   fCx2(0),
  50   fAx1(0),
  51   fAx2(0),
  52   fXl(0),
  53   fYl(0),
  54   fZl(0),
  55   fRn(0),
  56   fPhin(0),
  57   fZpl(0)
  58 {
  59   //
  60   // Default constructor for the Dipole field
  61   //
  62 }
  63
  64 //_______________________________________________________________________
  65 AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
  66                      Float_t factor, Float_t fmax):
  67   AliMagFC(name,title,integ,factor,fmax),
  68   fSolenoid(0),
  69   fInd(0),
  70   fZmin(0),
  71   fZmax(0),
  72   fYmax(0),
  73   fYmin(0),
  74   fZpmx(0),
  75   fZpmn(0),
  76   fRmax(0),
  77   fRmin(0),
  78   fXdel(0),
  79   fYdel(0),
  80   fZdel(0),
  81   fRdel(0),
  82   fPhid(0),
  83   fZpdl(0),
  84   fCx1(0),
  85   fCx2(0),
  86   fAx1(0),
  87   fAx2(0),
  88   fXl(0),
  89   fYl(0),
  90   fZl(0),
  91   fRn(0),
  92   fPhin(0),
  93   fZpl(0)
  94 {
  95   //
  96   // Standard constructor for the Dipole field
  97   //
  98   fType = kDipoMap;
  99   fMap  = 3;
 100   SetSolenoidField();
 101
 102   Info("ctor",
 103        "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s\n",
 104        fName.Data(), fMap ,integ,factor,fTitle.Data());
 105
 106 }
 107
 108 //_______________________________________________________________________
 109 void AliMagFDM::Field(Float_t *xfi, Float_t *b) const
 110 {
 111   //
 112   // Main routine to compute the field in a point
 113   //
 114   const Double_t keps=0.1E-06;
 115   const Double_t kPI2=2.*TMath::Pi();
 116   const Double_t kone=1;
 117
 118   const    Int_t  kiip=33;
 119   const    Int_t  kmiip=0;
 120   const    Int_t  kliip=0;
 121
 122   const    Int_t  kiic=0;
 123   const    Int_t  kmiic=0;
 124   const    Int_t  kliic=0;
 125
 126   const Double_t    kfZbg=502.92;  // Start of Map using in z
 127   const Double_t    kfZL3=600;  // Beginning of L3 door in z
 128
 129   Double_t   x[3];
 130   Double_t   xL3[3];
 131   Double_t   bint[3];
 132
 133   Double_t r0;
 134   Int_t iKvar,jb;
 135
 136   Double_t zp1, zp2,xp1,xp2,yp1,yp2;
 137   Double_t zz1, zz2,yy1,yy2,x2,x1;
 138
 139 // --- start the map fiel from z = 502.92 cm ---
 140 //
 141 // This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
 142 // Transfor correspondingly.
 143
 144   x[0] = - xfi[0];
 145   x[1] =   xfi[1];
 146   x[2] = - xfi[2];
 147   b[0]=b[1]=b[2]=0;
 148 //
 149   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]);
 150
 151   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 152            r0=rr/100;
 153   Double_t rPmax;
 154   rPmax=fRmax;
 155   if ( (-700<x[2] && x[2]<=kfZbg &&
 156         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
 157        || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
 158        {
 159         b[0]=b[1]=0;
 160         b[2]=fSolenoid;
 161        }
 162
 163   xL3[0]=x[0]/100;
 164   xL3[1]=(x[1]+30)/100;
 165   xL3[2]=x[2]/100;
 166
 167
 168   if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
 169   if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
 170
 171   Double_t xminn=xL3[2]*fAx1+fCx1;
 172   Double_t xmaxx=xL3[2]*fAx2+fCx2;
 173   Double_t zCmin,zCmax,yCmin,yCmax;
 174
 175   zCmin=fZmin;
 176   zCmax=fZmax;
 177   yCmin=fYmin;
 178   yCmax=fYmax;
 179 if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2] <= zCmax ) && (yCmin<xL3[1] && xL3[1]<= yCmax) &&  (xminn<xL3[0] && xL3[0] <= xmaxx)))
 180     {
 181      if(fMap==3)
 182       {
 183         if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
 184        {
 185        //---------------------  Polar part ----------------------
 186
 187        Double_t yyp,ph0;
 188        Int_t  kp0, lp0, mp0;
 189        Int_t kpi,lpi,mpi;
 190        Double_t alp1,alp2,alp3;
 191        Double_t zpz,rp,fip,cphi;
 192
 193        kpi=kiip;
 194        lpi=kliip;
 195        mpi=kmiip;
 196
 197        zpz=xL3[2];
 198        kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
 199        zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
 200        zp1= kone-zp2;
 201        yyp=xL3[1]- 0.3;
 202        cphi=TMath::Abs(yyp/r0);
 203        Int_t kcphi=0;
 204        if (cphi > kone) {
 205         printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, cphi %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,cphi);
 206         cphi =kone;
 207         kcphi=777;
 208        }
 209        ph0=TMath::ACos(cphi);
 210        if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
 211        if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
 212        if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
 213        if (ph0 > kPI2) {       ph0=ph0 - kPI2;}
 214        if (kcphi==777) {
 215         printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, ph0 %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,ph0);
 216        }
 217        fip=ph0;
 218        mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
 219        xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
 220        xp1=kone-xp2;
 221
 222        Double_t rDel;
 223        rDel=fRdel;
 224
 225        if (r0<= fRdel)
 226         {
 227
 228          if(r0< keps)
 229          {
 230
 231           bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
 232           bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
 233           bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
 234
 235          } else {
 236
 237         alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
 238         alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
 239         alp1= kone - alp2 - alp3;
 240
 241         for (jb=0; jb<3 ; jb++)
 242          {
 243           iKvar=jb;
 244           bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
 245          }
 246         } // end   of keps <=r0
 247        }
 248         else
 249        {
 250         rp=r0;
 251
 252         lp0=FZ(rp,fR ,fRdel,lpi,fRn);
 253         yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
 254         yp1=kone-yp2;
 255
 256         for (jb=0; jb<3 ; jb++)
 257         {
 258          iKvar=jb;
 259          bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
 260         }
 261        }
 262
 263        b[0]=-bint[0];
 264        b[1]=bint[1];
 265        b[2]=-bint[2];
 266
 267    }
 268    else
 269    {
 270    //-------------- Cartensian part ------------------
 271
 272    Double_t zzc,yyc;
 273    Int_t k0, l0,m0;
 274    Double_t  xx1, xx2,dx, xxx ,xXl;
 275    Int_t kci,mci,lci;
 276
 277    kci=kiic;
 278    lci=kliic;
 279    mci=kmiic;
 280
 281    xx1 = fAx1*xL3[2] + fCx1;
 282    xx2 = fAx2*xL3[2] + fCx2;
 283
 284    zzc=xL3[2];
 285    k0=FZ(zzc, fZc ,fZdel, kci, fZl);
 286    zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
 287    zz1=kone - zz2;;
 288
 289    yyc=xL3[1];
 290    l0=FZ(yyc, fY , fYdel,lci,fYl);
 291    yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
 292    yy1=kone - yy2;
 293    xXl = fXl-kone;
 294    dx = (xx2-xx1)/xXl;
 295    xxx= xL3[0]-xx1;
 296    m0 = int(xxx/dx);
 297
 298    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
 299     {
 300       m0=m0+1;
 301       printf(" m0 %d, m0+1 %d\n",m0,m0+1);
 302     }
 303
 304    x2=(xL3[0]-( xx1+m0*dx))/dx;
 305    x1=kone-x2;
 306    m0=m0-1;
 307    for (jb=3; jb<6; jb++)
 308      {
 309        iKvar=jb;
 310        bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
 311      }
 312
 313    b[0]=-bint[0];
 314    b[1]=bint[1];
 315    b[2]=-bint[2];
 316
 317    }
 318
 319
 320   } else {
 321         printf("Unknown map of Dipole region %d\n",fMap);
 322  }
 323
 324 } else {
 325     ZDCField(xfi,b);
 326
 327     }
 328
 329   if(fFactor!=1) {
 330     b[0]*=fFactor;
 331     b[1]*=fFactor;
 332     b[2]*=fFactor;
 333   }
 334 }
 335
 336 //_______________________________________________________________________
 337 Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar,
 338                     Float_t delu, Int_t ik, Int_t nk) const
 339 {
 340   //
 341   // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
 342   //
 343   Int_t l,ikj;
 344   Double_t ddu,ar;
 345
 346   Int_t ku,kf;
 347   kf=0;
 348   ar=Ar[ik];
 349   ddu=temp-ar;
 350
 351   ku=int(ddu/delu);
 352   ikj=ik+ku+1;
 353   if (ddu<=0) ikj=0;
 354
 355    for(l=ikj; l<nk; l++)
 356    {
 357
 358     if(temp <= Ar[l])
 359      {
 360
 361       kf=l-1;
 362       break;
 363      }
 364     }
 365     return kf;
 366   }
 367
 368 //_______________________________________________________________________
 369 Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2,
 370                        Double_t alf1, Double_t alf2, Double_t alf3,
 371                        Int_t kaa, Int_t maa) const
 372 {
 373   //
 374   // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
 375   //
 376   Double_t fa11=0,fa12=0,fa13=0;
 377   Double_t fa21=0,fa22=0,fa23=0;
 378   Double_t faY1=0,faY2=0;
 379   Double_t bba=0;
 380
 381   switch (kaai) {
 382   case 0:
 383     fa11 = fBpx[kaa][0][0];
 384     fa12 = fBpx[kaa][0][maa];
 385     fa13 = fBpx[kaa][0][maa+1];
 386     fa21 = fBpx[kaa+1][0][0];
 387     fa22 = fBpx[kaa+1][0][maa];
 388     fa23 = fBpx[kaa+1][0][maa+1];
 389     break;
 390   case 1:
 391     fa11 = fBpy[kaa][0][0];
 392     fa12 = fBpy[kaa][0][maa];
 393     fa13 = fBpy[kaa][0][maa+1];
 394     fa21 = fBpy[kaa+1][0][0];
 395     fa22 = fBpy[kaa+1][0][maa];
 396     fa23 = fBpy[kaa+1][0][maa+1];
 397     break;
 398   case 2:
 399     fa11 = fBpz[kaa][0][0];
 400     fa12 = fBpz[kaa][0][maa];
 401     fa13 = fBpz[kaa][0][maa+1];
 402     fa21 = fBpz[kaa+1][0][0];
 403     fa22 = fBpz[kaa+1][0][maa];
 404     fa23 = fBpz[kaa+1][0][maa+1];
 405     break;
 406   default:
 407     Fatal("Ba","Invalid value of kaai %d\n",kaai);
 408   }
 409   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
 410   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
 411   bba =  zaa1*faY1+zaa2*faY2;
 412   return bba;
 413 }
 414
 415
 416 //_______________________________________________________________________
 417 Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2,
 418                        Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m) const
 419 {
 420   //
 421   // Calculation of field componet at a given axis (general case)
 422   //
 423   Double_t fy1=0,fy2=0,ffy=0;
 424   Double_t gy1=0,gy2=0,ggy=0;
 425   Double_t bbi=0;
 426   Double_t bf11=0,bf12=0,bf21=0,bf22=0;
 427   Double_t bg11=0,bg12=0,bg21=0,bg22=0;
 428
 429
 430   /*-----------------Polar part ------------------*/
 431
 432   switch (kv) {
 433   case 0:
 434     bf11=fBpx[k][l][m];
 435     bf12=fBpx[k+1][l][m];
 436     bf21=fBpx[k+1][l+1][m];
 437     bf22=fBpx[k][l+1][m];
 438
 439     bg11=fBpx[k][l][m+1];
 440     bg12=fBpx[k+1][l][m+1];
 441     bg21=fBpx[k+1][l+1][m+1];
 442     bg22=fBpx[k][l+1][m+1];
 443     break;
 444
 445   case 1:
 446     bf11=fBpy[k][l][m];
 447     bf12=fBpy[k+1][l][m];
 448     bf21=fBpy[k+1][l+1][m];
 449     bf22=fBpy[k][l+1][m];
 450
 451     bg11=fBpy[k][l][m+1];
 452     bg12=fBpy[k+1][l][m+1];
 453     bg21=fBpy[k+1][l+1][m+1];
 454     bg22=fBpy[k][l+1][m+1];
 455     break;
 456
 457   case 2:
 458     bf11=fBpz[k][l][m];
 459     bf12=fBpz[k+1][l][m];
 460     bf21=fBpz[k+1][l+1][m];
 461     bf22=fBpz[k][l+1][m];
 462
 463     bg11=fBpz[k][l][m+1];
 464     bg12=fBpz[k+1][l][m+1];
 465     bg21=fBpz[k+1][l+1][m+1];
 466     bg22=fBpz[k][l+1][m+1];
 467     break;
 468   /*-----------------Cartensian part ---------------*/
 469
 470   case 3:
 471     bf11=fBcx[k][l][m];
 472     bf12=fBcx[k+1][l][m];
 473     bf21=fBcx[k+1][l+1][m];
 474     bf22=fBcx[k][l+1][m];
 475
 476     bg11=fBcx[k][l][m+1];
 477     bg12=fBcx[k+1][l][m+1];
 478     bg21=fBcx[k+1][l+1][m+1];
 479     bg22=fBcx[k][l+1][m+1];
 480     break;
 481
 482   case 4:
 483     bf11=fBcy[k][l][m];
 484     bf12=fBcy[k+1][l][m];
 485     bf21=fBcy[k+1][l+1][m];
 486     bf22=fBcy[k][l+1][m];
 487
 488     bg11=fBcy[k][l][m+1];
 489     bg12=fBcy[k+1][l][m+1];
 490     bg21=fBcy[k+1][l+1][m+1];
 491     bg22=fBcy[k][l+1][m+1];
 492     break;
 493
 494   case 5:
 495     bf11=fBcz[k][l][m];
 496     bf12=fBcz[k+1][l][m];
 497     bf21=fBcz[k+1][l+1][m];
 498     bf22=fBcz[k][l+1][m];
 499
 500     bg11=fBcz[k][l][m+1];
 501     bg12=fBcz[k+1][l][m+1];
 502     bg21=fBcz[k+1][l+1][m+1];
 503     bg22=fBcz[k][l+1][m+1];
 504     break;
 505
 506   default:
 507     Fatal("Bb","Invalid value of kv %d\n",kv);
 508   }
 509
 510
 511
 512   fy1=z1*bf11+z2*bf12;
 513   fy2=z2*bf21+z1* bf22;
 514   ffy=y1*fy1+ y2*fy2;
 515
 516
 517   gy1 = z1*bg11+z2*bg12;
 518   gy2 = z2*bg21+z1*bg22;
 519   ggy= y1*gy1 +  y2*gy2;
 520
 521   bbi =  x1*ffy+x2*ggy;
 522
 523   return bbi;
 524
 525 }
 526
 527 //_______________________________________________________________________
 528 void AliMagFDM::ReadField()
 529 {
 530   //
 531   // Method to read the magnetic field from file
 532   //
 533   FILE *magfile;
 534
 535   Int_t ik, il, im;
 536   Float_t zzp, rr,phii;
 537   Float_t zz, yy, bx,by,bz,bb;
 538
 539   char *fname;
 540   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
 541   fname = gSystem->ExpandPathName(fTitle.Data());
 542   magfile=fopen(fname,"r");
 543   delete [] fname;
 544
 545   printf("Cartensian part\n");
 546
 547   if (magfile) {
 548
 549 //  Cartensian part
 550
 551     fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
 552
 553     printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
 554
 555     for (ik=0; ik<fZl; ik++)
 556     {
 557
 558       fscanf(magfile, " %e ", &zz);
 559       fZc[ik]=zz;
 560
 561     }
 562
 563     for (ik=0; ik<fYl; ik++)
 564     {
 565        fscanf(magfile, " %e ", &yy);
 566        fY[ik]=yy;
 567
 568     }
 569     for (ik=0; ik<81; ik++)
 570     {
 571            printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
 572     }
 573
 574     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
 575
 576 printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
 577 printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
 578 printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
 579
 580     for (il=0; il<44; il++)  {
 581      for (im=0; im<81; im++)  {
 582       for (ik=0; ik<81; ik++)  {
 583
 584       fscanf(magfile, " %e ", &by);
 585       fBcy[ik][im][il]=by;
 586       }
 587      }
 588     }
 589
 590     for (il=0; il<44; il++)  {
 591      for (im=0; im<81; im++)  {
 592       for (ik=0; ik<81; ik++)  {
 593
 594       fscanf(magfile, " %e ", &bx);
 595       fBcx[ik][im][il]=bx;
 596       }
 597      }
 598     }
 599
 600    for (il=0; il<44; il++)  {
 601      for (im=0; im<81; im++)  {
 602       for (ik=0; ik<81; ik++)  {
 603
 604       fscanf(magfile, " %e ", &bz);
 605       fBcz[ik][im][il]=bz;
 606       }
 607      }
 608     }
 609 //----------------------   Polar part ---------------------------------
 610
 611     printf("Polar part\n");
 612     fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
 613     printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
 614
 615     printf(" fZp array\n");
 616
 617     for (ik=0; ik<51; ik++)
 618     {
 619      fscanf(magfile, " %e ", &zzp);
 620      fZp[ik]=zzp;
 621      printf(" %e\n",fZp[ik]);
 622     }
 623
 624     printf(" fR array\n");
 625
 626     for (ik=0; ik<10; ik++)
 627     {
 628      fscanf(magfile, " %e ", &rr);
 629      fR[ik]=rr;
 630      printf(" %e\n",fR[ik]);
 631     }
 632
 633 //    printf("fPhi array\n");
 634
 635      for (il=0; il<33; il++)
 636      {
 637        fscanf(magfile, " %e ", &phii);
 638        fPhi[il]=phii;
 639 //        printf(" %e\n",fPhi[il]);
 640      }
 641
 642     fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
 643
 644 printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
 645
 646
 647     for (il=0; il<33; il++)  {
 648      for (im=0; im<10; im++)  {
 649       for (ik=0; ik<51; ik++)  {
 650       fscanf(magfile, " %e ", &by);
 651         fBpy[ik][im][il]=by;
 652       }
 653      }
 654     }
 655
 656     for (il=0; il<33; il++)  {
 657      for (im=0; im<10; im++)  {
 658       for (ik=0; ik<51; ik++)  {
 659       fscanf(magfile, " %e ", &bx);
 660         fBpx[ik][im][il]=bx;
 661       }
 662      }
 663     }
 664
 665
 666     for (il=0; il<33; il++)  {
 667      for (im=0; im<10; im++)  {
 668       for (ik=0; ik<51; ik++)  {
 669       fscanf(magfile, " %e ", &bz);
 670         fBpz[ik][im][il]=bz;
 671       }
 672      }
 673     }
 674
 675
 676     for (il=0; il<32; il++) {
 677      for (im=0; im<2; im++)  {
 678       for (ik=0; ik<2; ik++)  {
 679       fscanf(magfile, " %e ", &bb);
 680         fB[ik][im][il]=bb;
 681       }
 682      }
 683     }
 684 //
 685   } else {
 686     Fatal("ReadField","File %s not found !\n",fTitle.Data());
 687   }
 688 }