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 /*
  17 $Log$
  18 Revision 1.1  2000/07/11 18:24:59  fca
  19 Coding convention corrections + few minor bug fixes
  20
  21 */
  22
  23 #include <stdlib.h>
  24
  25 #include "AliMagFDM.h"
  26 #include "TSystem.h"
  27
  28
  29 ClassImp(AliMagFDM)
  30
  31 //________________________________________
  32 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
  33 const Int_t map, const Float_t factor, const Float_t fmax)
  34   : AliMagF(name,title,integ,map,factor,fmax)
  35
  36 {
  37   fType = kDipoMap;
  38
  39   printf("Field Map for Muon Arm from IP till muon filter %s created: map= %d, factor= %f, file=%s\n",fName.Data(),map,factor,fTitle.Data());
  40
  41 }
  42
  43 //________________________________________
  44
  45 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
  46 {
  47   //
  48   // Main routine to compute the field in a point
  49   //
  50   static  const Double_t keps=0.1E-06;
  51   static  const Double_t kpi2=.6283185E+01;
  52   static  const Double_t kone=1;
  53
  54   static  const    Int_t  kiip=33;
  55   static  const    Int_t  kmiip=0;
  56   static  const    Int_t  kliip=0;
  57
  58   static  const    Int_t  kiic=0;
  59   static  const    Int_t  kmiic=0;
  60   static  const    Int_t  kliic=0;
  61
  62   static  const Double_t    kfdZbg=502.92;  // Start of Map using in z
  63   static  const Double_t    kfdZL3=600;  // Beginning of L3 door in z
  64
  65   Double_t   x[3];
  66   Double_t   xL3[3];
  67   Double_t   bint[3];
  68
  69   Double_t r0;
  70
  71   Double_t bbj;
  72   Int_t iKvar,jb;
  73
  74   Double_t zp1, zp2,xp1,xp2,yp1,yp2;
  75   Double_t zz1, zz2,yy1,yy2,x2,x1;
  76
  77 // --- start the map fiel from z = 502.92 cm ---
  78
  79   x[0] = xfi[0];
  80   x[1] = xfi[1];
  81   x[2] = xfi[2];
  82   b[0]=b[1]=b[2]=0;
  83   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]);
  84
  85   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
  86            r0=rr/100;
  87   Double_t rPmax;
  88   rPmax=fdRmax;
  89   if ( (-700<x[2] && x[2]<=kfdZbg &&
  90         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
  91        || (kfdZbg<x[2] && x[2]<=kfdZL3 && rr>=rPmax*100) )
  92        {
  93         b[2]=2;
  94        }
  95
  96   xL3[0]=x[0]/100;
  97   xL3[1]=(x[1]+30)/100;
  98   xL3[2]=x[2]/100;
  99
 100   Double_t xminn=xL3[2]*fdAx1+fdCx1;
 101   Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
 102   Double_t zCmin,zCmax,yCmin,yCmax;
 103
 104   zCmin=fdZmin;
 105   zCmax=fdZmax;
 106   yCmin=fdYmin;
 107   yCmax=fdYmax;
 108
 109 if ((kfdZbg/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)))
 110     {
 111      if(fMap==3)
 112       {
 113        if (xL3[2]<zCmin && r0<rPmax)
 114        {
 115        //---------------------  Polar part ----------------------
 116
 117        Double_t yyp,ph0;
 118        Int_t  kp0, lp0, mp0;
 119        Int_t kpi,lpi,mpi;
 120        Double_t alp1,alp2,alp3;
 121        Double_t zpz,rp,fip,cphi;
 122
 123        kpi=kiip;
 124        lpi=kliip;
 125        mpi=kmiip;
 126
 127        zpz=xL3[2];
 128
 129        FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&zp1 ,&zp2,&fdZpl) ;
 130
 131        yyp=xL3[1]- 0.3;
 132        cphi=yyp/r0;
 133        ph0=TMath::ACos(cphi);
 134        if (xL3[0]< 0) {ph0=kpi2 - ph0;}
 135
 136        fip=ph0;
 137        FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &xp1,&xp2,&fdPhin);
 138
 139        Double_t rDel;
 140        rDel=fdRdel;
 141
 142        if (r0<= fdRdel)
 143         {
 144
 145          if(r0< keps)
 146          {
 147
 148           bint[0]=(zp1*fdBpx[kp0][0][0] + zp2*fdBpx[kp0+1][0][0])*10;
 149           bint[1]=(zp1*fdBpy[kp0][0][0] + zp2*fdBpy[kp0+1][0][0])*10;
 150           bint[2]=(zp1*fdBpz[kp0][0][0] + zp2*fdBpz[kp0+1][0][0])*10;
 151
 152          }
 153
 154         alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0];
 155         alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];
 156         alp1= kone - alp2 - alp3;
 157
 158         for (jb=0; jb<3 ; jb++)
 159         {
 160          iKvar=jb;
 161          FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
 162          bint[jb] = bbj*10 ;
 163         }
 164        }
 165         else
 166        {
 167         rp=r0;
 168
 169         FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&yp1,&yp2,&fdRn);
 170
 171         for (jb=0; jb<3 ; jb++)
 172         {
 173          iKvar=jb;
 174          FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj);
 175
 176          bint[jb] = bbj*10 ;
 177         }
 178        }
 179
 180        b[0]=bint[0];
 181        b[1]=bint[1];
 182        b[2]=bint[2];
 183
 184 //    fprintf(fitest,"-------------   Freg2 run -------------\n");
 185
 186    }
 187    else
 188    {
 189    //-------------- Cartensian part ------------------
 190
 191    Double_t zzc,yyc;
 192    Int_t k0, l0,m0;
 193    Double_t  xx1, xx2,dx, xxx ,xXl;
 194    Int_t kci,mci,lci;
 195
 196    kci=kiic;
 197    lci=kliic;
 198    mci=kmiic;
 199
 200    xx1 = fdAx1*xL3[2] + fdCx1;
 201    xx2 = fdAx2*xL3[2] + fdCx2;
 202
 203    zzc=xL3[2];
 204    FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &zz1, &zz2, &fdZl);
 205
 206    yyc=xL3[1];
 207    FZ(&yyc, fdY , &fdYdel,&lci, &l0, &yy1, &yy2,&fdYl);
 208
 209    xXl = fdXl-kone;
 210    dx = (xx2-xx1)/xXl;
 211    xxx= xL3[0]-xx1;
 212    //     xm = xxx/dx;
 213    m0 = int(xxx/dx);
 214
 215    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
 216     {
 217       m0=m0+1;
 218       printf(" m0 %d, m0+1 %d\n",m0,m0+1);
 219     }
 220
 221    x2=(xL3[0]-( xx1+m0*dx))/dx;
 222    x1=kone-x2;
 223    m0=m0-1;
 224    for (jb=3; jb<6; jb++)
 225      {
 226        iKvar=jb;
 227        FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj);
 228        bint[jb-3] = bbj*10 ;
 229      }
 230
 231    b[0]=bint[0];
 232    b[1]=bint[1];
 233    b[2]=bint[2];
 234
 235 //   fprintf(fitest,"------------   Freg1 run -----------------\n");
 236    }
 237
 238   } else {
 239         printf("Unknown map of Dipole region %d\n",fMap);
 240  }
 241
 242 } else {
 243
 244 //This is the ZDC part
 245       Float_t rad2=x[0]*x[0]+x[1]*x[1];
 246       if(rad2<kD2RA2) {
 247         if(x[2]>kD2BEG) {
 248
 249 //    Separator Dipole D2
 250           if(x[2]<kD2END) b[1]=kFDIP;
 251         } else if(x[2]>kD1BEG) {
 252
 253 //    Separator Dipole D1
 254           if(x[2]<kD1END) b[1]=-kFDIP;
 255         }
 256         if(rad2<kCORRA2) {
 257
 258 //    First quadrupole of inner triplet de-focussing in x-direction
 259 //    Inner triplet
 260           if(x[2]>kZ4BEG) {
 261             if(x[2]<kZ4END) {
 262
 263 //    2430 <-> 3060
 264               b[0]=-kG1*x[1];
 265               b[1]=-kG1*x[0];
 266             }
 267           } else if(x[2]>kZ3BEG) {
 268             if(x[2]<kZ3END) {
 269
 270 //    1530 <-> 2080
 271               b[0]=kG1*x[1];
 272               b[1]=kG1*x[0];
 273             }
 274           } else if(x[2]>kZ2BEG) {
 275             if(x[2]<kZ2END) {
 276
 277 //    890 <-> 1430
 278               b[0]=kG1*x[1];
 279               b[1]=kG1*x[0];
 280             }
 281           } else if(x[2]>kZ1BEG) {
 282             if(x[2]<kZ1END) {
 283
 284 //    0 <->  630
 285               b[0]=-kG1*x[1];
 286               b[1]=-kG1*x[0];
 287             }
 288           } else if(x[2]>kCORBEG) {
 289             if(x[2]<kCOREND) {
 290 //    Corrector dipole (because of dimuon arm)
 291 //            b[0]=kFCORN;
 292               b[0]=-kFCORN;
 293             }
 294           }
 295         }
 296       }
 297     }
 298
 299   if(fFactor!=1) {
 300     b[0]*=fFactor;
 301     b[1]*=fFactor;
 302     b[2]*=fFactor;
 303   }
 304 }
 305
 306 //_________________________________________
 307
 308 void AliMagFDM::FZ(Double_t *u, Float_t  *Ar, Float_t *du,Int_t *ki,Int_t *kf,Double_t *a1,Double_t *a2 ,Int_t *nu)
 309
 310 {
 311   //
 312   // Z component of the field
 313   //
 314   static  const Double_t kone=1;
 315   Int_t l,ik,ikj;
 316   Double_t temp;
 317   Double_t ddu,delu,ar;
 318
 319   Int_t nk,ku;
 320   temp=*u;
 321   nk=*nu;
 322   ik=*ki;
 323   delu=*du;
 324
 325   ar=Ar[ik];
 326   ddu=temp-ar;
 327
 328   ku=int(ddu/delu);
 329   ikj=ik+ku;
 330   if (ddu<=0) ikj=0;
 331
 332    for(l=ikj; l<nk; l++)
 333    {
 334
 335     if(temp < Ar[l])
 336      {
 337       *kf=l;
 338       *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
 339       *a1= kone - *a2;
 340       break;
 341      }
 342     }
 343   }
 344
 345 /*-------------FRfuncBi----------------*/
 346
 347 void AliMagFDM::FRfuncBi(Int_t *kai,Double_t *za1, Double_t *za2, Double_t *al1, Double_t *al2, Double_t *al3, Int_t *ka, Int_t *ma, Double_t *ba)
 348 {
 349   //
 350   // This method needs to be commented
 351   //
 352   Double_t fa11,fa12,fa13;
 353   Double_t fa21,fa22,fa23;
 354   Double_t faY1,faY2;
 355   Double_t bba;
 356
 357   Double_t zaa1,zaa2,alf1,alf2,alf3;
 358   Int_t kaai,kaa,maa;
 359   kaai=*kai;
 360   kaa=*ka;
 361   maa=*ma;
 362   zaa1=*za1;
 363   zaa2=*za2;
 364   alf1=*al1;
 365   alf2=*al2;
 366   alf3=*al3;
 367
 368   switch (kaai) {
 369   case 0:
 370     fa11 = fdBpx[kaa][0][0];
 371     fa12 = fdBpx[kaa][0][maa];
 372     fa13 = fdBpx[kaa][0][maa+1];
 373     fa21 = fdBpx[kaa+1][0][0];
 374     fa22 = fdBpx[kaa+1][0][maa];
 375     fa23 = fdBpx[kaa+1][0][maa+1];
 376     break;
 377   case 1:
 378     fa11 = fdBpy[kaa][0][0];
 379     fa12 = fdBpy[kaa][0][maa];
 380     fa13 = fdBpy[kaa][0][maa+1];
 381     fa21 = fdBpy[kaa+1][0][0];
 382     fa22 = fdBpy[kaa+1][0][maa];
 383     fa23 = fdBpy[kaa+1][0][maa+1];
 384     break;
 385   case 2:
 386     fa11 = fdBpz[kaa][0][0];
 387     fa12 = fdBpz[kaa][0][maa];
 388     fa13 = fdBpz[kaa][0][maa+1];
 389     fa21 = fdBpz[kaa+1][0][0];
 390     fa22 = fdBpz[kaa+1][0][maa];
 391     fa23 = fdBpz[kaa+1][0][maa+1];
 392     break;
 393   default:
 394     Fatal("FRfuncBi","Invalid value of kaai %d\n",kaai);
 395     exit(1);
 396   }
 397   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
 398   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
 399   bba =  zaa1*faY1+zaa2*faY2;
 400   *ba=bba;
 401
 402 }
 403
 404
 405 /*----------- FGfuncBi------------*/
 406
 407 void AliMagFDM::FGfuncBi(Double_t *zz1,Double_t *zz2, Double_t *yy1,Double_t *yy2, Double_t *xx1,Double_t *xx2, Int_t *kvr, Int_t *kk, Int_t *ll, Int_t *mm, Double_t *bb)
 408 {
 409   //
 410   // This method needs to be commented
 411   //
 412   Double_t fy1, fy2, ffy;
 413   Double_t gy1,gy2,ggy;
 414   Double_t z1,z2,y1,y2,x1,x2;
 415
 416   Int_t k,l,m,kv;
 417   Double_t bbi;
 418
 419   Double_t bf11,bf12,bf21,bf22;
 420   Double_t bg11,bg12,bg21,bg22;
 421   k=*kk;
 422   l=*ll;
 423   m=*mm;
 424
 425   kv=*kvr;
 426
 427   z1=*zz1;
 428   z2=*zz2;
 429   y1=*yy1;
 430   y2=*yy2;
 431   x1=*xx1;
 432   x2=*xx2;
 433
 434   /*-----------------Polar part ------------------*/
 435
 436   switch (kv) {
 437   case 0:
 438     bf11=fdBpx[k][l][m];
 439     bf12=fdBpx[k+1][l][m];
 440     bf21=fdBpx[k+1][l+1][m];
 441     bf22=fdBpx[k][l+1][m];
 442
 443     bg11=fdBpx[k][l][m+1];
 444     bg12=fdBpx[k+1][l][m+1];
 445     bg21=fdBpx[k+1][l+1][m+1];
 446     bg22=fdBpx[k][l+1][m+1];
 447     break;
 448
 449   case 1:
 450     bf11=fdBpy[k][l][m];
 451     bf12=fdBpy[k+1][l][m];
 452     bf21=fdBpy[k+1][l+1][m];
 453     bf22=fdBpy[k][l+1][m];
 454
 455     bg11=fdBpy[k][l][m+1];
 456     bg12=fdBpy[k+1][l][m+1];
 457     bg21=fdBpy[k+1][l+1][m+1];
 458     bg22=fdBpy[k][l+1][m+1];
 459     break;
 460
 461   case 2:
 462     bf11=fdBpz[k][l][m];
 463     bf12=fdBpz[k+1][l][m];
 464     bf21=fdBpz[k+1][l+1][m];
 465     bf22=fdBpz[k][l+1][m];
 466
 467     bg11=fdBpz[k][l][m+1];
 468     bg12=fdBpz[k+1][l][m+1];
 469     bg21=fdBpz[k+1][l+1][m+1];
 470     bg22=fdBpz[k][l+1][m+1];
 471     break;
 472   /*-----------------Cartensian part ---------------*/
 473
 474   case 3:
 475     bf11=fdBcx[k][l][m];
 476     bf12=fdBcx[k+1][l][m];
 477     bf21=fdBcx[k+1][l+1][m];
 478     bf22=fdBcx[k][l+1][m];
 479
 480     bg11=fdBcx[k][l][m+1];
 481     bg12=fdBcx[k+1][l][m+1];
 482     bg21=fdBcx[k+1][l+1][m+1];
 483     bg22=fdBcx[k][l+1][m+1];
 484     break;
 485
 486   case 4:
 487     bf11=fdBcy[k][l][m];
 488     bf12=fdBcy[k+1][l][m];
 489     bf21=fdBcy[k+1][l+1][m];
 490     bf22=fdBcy[k][l+1][m];
 491
 492     bg11=fdBcy[k][l][m+1];
 493     bg12=fdBcy[k+1][l][m+1];
 494     bg21=fdBcy[k+1][l+1][m+1];
 495     bg22=fdBcy[k][l+1][m+1];
 496     break;
 497
 498   case 5:
 499     bf11=fdBcz[k][l][m];
 500     bf12=fdBcz[k+1][l][m];
 501     bf21=fdBcz[k+1][l+1][m];
 502     bf22=fdBcz[k][l+1][m];
 503
 504     bg11=fdBcz[k][l][m+1];
 505     bg12=fdBcz[k+1][l][m+1];
 506     bg21=fdBcz[k+1][l+1][m+1];
 507     bg22=fdBcz[k][l+1][m+1];
 508     break;
 509
 510   default:
 511     Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
 512     exit(1);
 513   }
 514
 515
 516
 517   fy1=z1*bf11+z2*bf12;
 518   fy2=z2*bf21+z1* bf22;
 519   ffy=y1*fy1+ y2*fy2;
 520
 521
 522   gy1 = z1*bg11+z2*bg12;
 523   gy2 = z2*bg21+z1*bg22;
 524   ggy= y1*gy1 +  y2*gy2;
 525
 526   bbi =  x1*ffy+x2*ggy;
 527
 528   *bb=bbi;
 529
 530 }
 531 //____________________________________________
 532
 533 void AliMagFDM::ReadField()
 534 {
 535   //
 536   // Method to read the magnetic field from file
 537   //
 538   FILE *magfile;
 539
 540   Int_t ik, il, im;
 541   Float_t zzp, rr,phii;
 542   Float_t zz, yy, bx,by,bz,bb;
 543
 544   char *fname;
 545   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
 546   fname = gSystem->ExpandPathName(fTitle.Data());
 547   magfile=fopen(fname,"r");
 548   delete [] fname;
 549
 550   printf("Cartensian part\n");
 551
 552   if (magfile) {
 553
 554 //  Cartensian part
 555
 556     fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl);
 557
 558     printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);
 559
 560     for (ik=0; ik<fdZl; ik++)
 561     {
 562
 563       fscanf(magfile, " %e ", &zz);
 564       fdZc[ik]=zz;
 565
 566     }
 567
 568     for (ik=0; ik<fdYl; ik++)
 569     {
 570        fscanf(magfile, " %e ", &yy);
 571        fdY[ik]=yy;
 572
 573     }
 574     for (ik=0; ik<81; ik++)
 575     {
 576            printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]);
 577     }
 578
 579     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2);
 580
 581 printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
 582 printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
 583 printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
 584
 585     for (il=0; il<44; il++)  {
 586      for (im=0; im<81; im++)  {
 587       for (ik=0; ik<81; ik++)  {
 588
 589       fscanf(magfile, " %e ", &by);
 590       fdBcy[ik][im][il]=by;
 591       }
 592      }
 593     }
 594
 595     for (il=0; il<44; il++)  {
 596      for (im=0; im<81; im++)  {
 597       for (ik=0; ik<81; ik++)  {
 598
 599       fscanf(magfile, " %e ", &bx);
 600       fdBcx[ik][im][il]=bx;
 601       }
 602      }
 603     }
 604
 605    for (il=0; il<44; il++)  {
 606      for (im=0; im<81; im++)  {
 607       for (ik=0; ik<81; ik++)  {
 608
 609       fscanf(magfile, " %e ", &bz);
 610       fdBcz[ik][im][il]=bz;
 611       }
 612      }
 613     }
 614 //----------------------   Polar part ---------------------------------
 615
 616     printf("Polar part\n");
 617     fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin);
 618     printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);
 619
 620     printf(" fdZp array\n");
 621
 622     for (ik=0; ik<51; ik++)
 623     {
 624      fscanf(magfile, " %e ", &zzp);
 625      fdZp[ik]=zzp;
 626      printf(" %e\n",fdZp[ik]);
 627     }
 628
 629     printf(" fdR array\n");
 630
 631     for (ik=0; ik<10; ik++)
 632     {
 633      fscanf(magfile, " %e ", &rr);
 634      fdR[ik]=rr;
 635      printf(" %e\n",fdR[ik]);
 636     }
 637
 638 //    printf("fdPhi array\n");
 639
 640      for (il=0; il<33; il++)
 641      {
 642        fscanf(magfile, " %e ", &phii);
 643        fdPhi[il]=phii;
 644 //        printf(" %e\n",fdPhi[il]);
 645      }
 646
 647     fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin);
 648
 649 printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
 650
 651
 652     for (il=0; il<33; il++)  {
 653      for (im=0; im<10; im++)  {
 654       for (ik=0; ik<51; ik++)  {
 655       fscanf(magfile, " %e ", &by);
 656         fdBpy[ik][im][il]=by;
 657       }
 658      }
 659     }
 660
 661     for (il=0; il<33; il++)  {
 662      for (im=0; im<10; im++)  {
 663       for (ik=0; ik<51; ik++)  {
 664       fscanf(magfile, " %e ", &bx);
 665         fdBpx[ik][im][il]=bx;
 666       }
 667      }
 668     }
 669
 670
 671     for (il=0; il<33; il++)  {
 672      for (im=0; im<10; im++)  {
 673       for (ik=0; ik<51; ik++)  {
 674       fscanf(magfile, " %e ", &bz);
 675         fdBpz[ik][im][il]=bz;
 676       }
 677      }
 678     }
 679
 680
 681     for (il=0; il<32; il++) {
 682      for (im=0; im<2; im++)  {
 683       for (ik=0; ik<2; ik++)  {
 684       fscanf(magfile, " %e ", &bb);
 685         fdB[ik][im][il]=bb;
 686       }
 687      }
 688     }
 689 //
 690   } else {
 691     printf("File %s not found !\n",fTitle.Data());
 692     exit(1);
 693   }
 694 }
 695 //________________________________