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