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