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