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