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