X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMagFDM.cxx;h=ecee137df538bf64698ab2508b6ae6c3b8fd0a43;hb=a83e70793facb29c0370dc6e30f3bc2be94fa43b;hp=42e50002eab653eaecf451fa10b8a4b44b393cbb;hpb=aee8290b9e98b230f8f38596bffd47ca1abfbf9d;p=u%2Fmrichter%2FAliRoot.git

diff --git a/STEER/AliMagFDM.cxx b/STEER/AliMagFDM.cxx
index 42e50002eab..ecee137df53 100644
--- a/STEER/AliMagFDM.cxx
+++ b/STEER/AliMagFDM.cxx
@@ -13,57 +13,127 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-*/
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//   Field with Magnetic Field map
+//   Used by AliRun class
+//   Author:
+//-------------------------------------------------------------------------
+
+#include <stdlib.h>
 
-#include "AliMagFDM.h"
 #include "TSystem.h"
+
+#include "AliMagFDM.h"
     
 
 ClassImp(AliMagFDM)
 
-//________________________________________
+//_______________________________________________________________________
+AliMagFDM::AliMagFDM():
+  fSolenoid(0),
+  fInd(0),
+  fZmin(0),
+  fZmax(0),
+  fYmax(0),
+  fYmin(0),
+  fZpmx(0),
+  fZpmn(0),
+  fRmax(0),
+  fRmin(0),
+  fXdel(0),
+  fYdel(0),
+  fZdel(0),
+  fRdel(0),
+  fPhid(0),
+  fZpdl(0),
+  fCx1(0),
+  fCx2(0),
+  fAx1(0),
+  fAx2(0),
+  fXl(0),
+  fYl(0),
+  fZl(0),    
+  fRn(0),
+  fPhin(0),
+  fZpl(0)  
+{
+  //
+  // Default constructor for the Dipole field
+  //
+}
+
+//_______________________________________________________________________
 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
-const Int_t map, const Float_t factor, const Float_t fmax)
-  : AliMagF(name,title,integ,map,factor,fmax)
-  
+                     const Float_t factor, const Float_t fmax):
+  AliMagF(name,title,integ,factor,fmax),
+  fSolenoid(0),
+  fInd(0),
+  fZmin(0),
+  fZmax(0),
+  fYmax(0),
+  fYmin(0),
+  fZpmx(0),
+  fZpmn(0),
+  fRmax(0),
+  fRmin(0),
+  fXdel(0),
+  fYdel(0),
+  fZdel(0),
+  fRdel(0),
+  fPhid(0),
+  fZpdl(0),
+  fCx1(0),
+  fCx2(0),
+  fAx1(0),
+  fAx2(0),
+  fXl(0),
+  fYl(0),
+  fZl(0),    
+  fRn(0),
+  fPhin(0),
+  fZpl(0)  
 {
+  //
+  // Standard constructor for the Dipole field
+  //
   fType = kDipoMap;
-
-  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());
+  fMap  = 3;
+  SetSolenoidField();
   
+  Info("ctor",
+       "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());
+ 
 }
 
-//________________________________________
-
+//_______________________________________________________________________
 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
 {
   //
   // Main routine to compute the field in a point
   //
-  static  const Double_t keps=0.1E-06;
-  static  const Double_t kpi2=.6283185E+01;
-  static  const Double_t kone=1;
-
-  static  const    Int_t  kiip=33; 
-  static  const    Int_t  kmiip=0;    
-  static  const    Int_t  kliip=0;
-
-  static  const    Int_t  kiic=0;
-  static  const    Int_t  kmiic=0;    
-  static  const    Int_t  kliic=0;       
-
-  static  const Double_t    kfdZbg=502.92;  // Start of Map using in z
-  static  const Double_t    kfdZL3=600;  // Beginning of L3 door in z
+  const Double_t keps=0.1E-06;
+  const Double_t kPI2=2.*TMath::Pi();
+  const Double_t kone=1;
+  
+  const    Int_t  kiip=33; 
+  const    Int_t  kmiip=0;    
+  const    Int_t  kliip=0;
+  
+  const    Int_t  kiic=0;
+  const    Int_t  kmiic=0;    
+  const    Int_t  kliic=0;       
+  
+  const Double_t    kfZbg=502.92;  // Start of Map using in z
+  const Double_t    kfZL3=600;  // Beginning of L3 door in z
 
   Double_t   x[3];   
   Double_t   xL3[3]; 
   Double_t   bint[3]; 
   
   Double_t r0;
-  
-  Double_t bbj;
   Int_t iKvar,jb;
 
   Double_t zp1, zp2,xp1,xp2,yp1,yp2; 
@@ -80,32 +150,36 @@ void AliMagFDM::Field(Float_t *xfi, Float_t *b)
   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
            r0=rr/100;
   Double_t rPmax;
-  rPmax=fdRmax; 
-  if ( (-700<x[2] && x[2]<=kfdZbg && 
+  rPmax=fRmax; 
+  if ( (-700<x[2] && x[2]<=kfZbg && 
         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
-       || (kfdZbg<x[2] && x[2]<=kfdZL3 && rr>=rPmax*100) )
+       || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
        {
-        b[2]=2;
+        b[0]=b[1]=0;
+        b[2]=fSolenoid;
        }
 
   xL3[0]=x[0]/100;
   xL3[1]=(x[1]+30)/100;
   xL3[2]=x[2]/100; 
+
  
-  Double_t xminn=xL3[2]*fdAx1+fdCx1;
-  Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
+  if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
+  if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
+
+  Double_t xminn=xL3[2]*fAx1+fCx1;
+  Double_t xmaxx=xL3[2]*fAx2+fCx2;
   Double_t zCmin,zCmax,yCmin,yCmax;
 
-  zCmin=fdZmin;
-  zCmax=fdZmax;
-  yCmin=fdYmin;
-  yCmax=fdYmax;
-          
-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)))
+  zCmin=fZmin;
+  zCmax=fZmax;
+  yCmin=fYmin;
+  yCmax=fYmax;
+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)))
     {
      if(fMap==3) 
       { 
-       if (xL3[2]<zCmin && r0<rPmax)   
+        if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
        {
        //---------------------  Polar part ----------------------
        
@@ -120,55 +194,68 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
        mpi=kmiip;
    
        zpz=xL3[2];
-
-       FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&zp1 ,&zp2,&fdZpl) ;
-
+       kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
+       zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
+       zp1= kone-zp2;
        yyp=xL3[1]- 0.3;
-       cphi=yyp/r0;
+       cphi=TMath::Abs(yyp/r0);
+       Int_t kcphi=0;
+       if (cphi > kone) {
+        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);
+        cphi =kone;
+        kcphi=777;
+       } 
        ph0=TMath::ACos(cphi);
-       if (xL3[0]< 0) {ph0=kpi2 - ph0;}  
-                                  
-       fip=ph0;
-       FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &xp1,&xp2,&fdPhin); 
+       if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}  
+       if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;} 
+       if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}  
+       if (ph0 > kPI2) {       ph0=ph0 - kPI2;}
+       if (kcphi==777) {
+        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);
+       }  
+       fip=ph0; 
+       mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
+       xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
+       xp1=kone-xp2;
 
        Double_t rDel;
-       rDel=fdRdel;
+       rDel=fRdel;
 
-       if (r0<= fdRdel)
+       if (r0<= fRdel)
         {
 
          if(r0< keps) 
          {
 
-          bint[0]=(zp1*fdBpx[kp0][0][0] + zp2*fdBpx[kp0+1][0][0])*10;     
-          bint[1]=(zp1*fdBpy[kp0][0][0] + zp2*fdBpy[kp0+1][0][0])*10;      
-          bint[2]=(zp1*fdBpz[kp0][0][0] + zp2*fdBpz[kp0+1][0][0])*10; 
+          bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;     
+          bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;      
+          bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10; 
 
-         }  
+         } else { 
       
-        alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0]; 
-        alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];      
+        alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0]; 
+        alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];      
         alp1= kone - alp2 - alp3;
      
         for (jb=0; jb<3 ; jb++) 
-        {
-         iKvar=jb;     
-         FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
-         bint[jb] = bbj*10 ;
-        }  
+         {
+          iKvar=jb;     
+          bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
+         } 
+        } // end   of keps <=r0 
        }
         else
        {
         rp=r0;
 
-        FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&yp1,&yp2,&fdRn);
+        lp0=FZ(rp,fR ,fRdel,lpi,fRn);
+        yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
+	yp1=kone-yp2;
 
         for (jb=0; jb<3 ; jb++) 
         {
          iKvar=jb;
-         FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj); 
-
-         bint[jb] = bbj*10 ;
+         bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
         }
        }
 
@@ -176,8 +263,6 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
        b[1]=bint[1];
        b[2]=bint[2];
 
-//    fprintf(fitest,"-------------   Freg2 run -------------\n");       
-  
    }  
    else 
    {
@@ -192,19 +277,21 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
    lci=kliic;
    mci=kmiic;
 
-   xx1 = fdAx1*xL3[2] + fdCx1;
-   xx2 = fdAx2*xL3[2] + fdCx2;
+   xx1 = fAx1*xL3[2] + fCx1;
+   xx2 = fAx2*xL3[2] + fCx2;
 
    zzc=xL3[2];
-   FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &zz1, &zz2, &fdZl);     
+   k0=FZ(zzc, fZc ,fZdel, kci, fZl); 
+   zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
+   zz1=kone - zz2;;    
 
    yyc=xL3[1];
-   FZ(&yyc, fdY , &fdYdel,&lci, &l0, &yy1, &yy2,&fdYl);  
-
-   xXl = fdXl-kone; 
+   l0=FZ(yyc, fY , fYdel,lci,fYl);  
+   yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
+   yy1=kone - yy2;    
+   xXl = fXl-kone; 
    dx = (xx2-xx1)/xXl;
    xxx= xL3[0]-xx1;
-   //     xm = xxx/dx; 
    m0 = int(xxx/dx);   
 
    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx)) 
@@ -219,17 +306,16 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
    for (jb=3; jb<6; jb++) 
      {
        iKvar=jb;     
-       FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj); 
-       bint[jb-3] = bbj*10 ; 
+       bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ; 
      }    
  
    b[0]=bint[0];
    b[1]=bint[1];
    b[2]=bint[2]; 
 
-//   fprintf(fitest,"------------   Freg1 run -----------------\n");            
    } 
 
+
   } else {
         printf("Unknown map of Dipole region %d\n",fMap);
  }
@@ -237,59 +323,49 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
 } else {
 
 //This is the ZDC part
-      Float_t rad2=x[0]*x[0]+x[1]*x[1];
-      if(rad2<kD2RA2) {
-        if(x[2]>kD2BEG) {
-          
-//    Separator Dipole D2
-          if(x[2]<kD2END) b[1]=kFDIP;
-        } else if(x[2]>kD1BEG) {
-
-//    Separator Dipole D1
-          if(x[2]<kD1END) b[1]=-kFDIP;
-        }
-        if(rad2<kCORRA2) {
-
-//    First quadrupole of inner triplet de-focussing in x-direction
-//    Inner triplet
-          if(x[2]>kZ4BEG) {
-            if(x[2]<kZ4END) {
-              
-//    2430 <-> 3060
-              b[0]=-kG1*x[1];
-              b[1]=-kG1*x[0];
-            }
-          } else if(x[2]>kZ3BEG) {
-            if(x[2]<kZ3END) {
-
-//    1530 <-> 2080
-              b[0]=kG1*x[1];
-              b[1]=kG1*x[0];
-            }
-          } else if(x[2]>kZ2BEG) {
-            if(x[2]<kZ2END) {
-
-//    890 <-> 1430
-              b[0]=kG1*x[1];
-              b[1]=kG1*x[0];
-            }
-          } else if(x[2]>kZ1BEG) {
-            if(x[2]<kZ1END) {
-
-//    0 <->  630
-              b[0]=-kG1*x[1];
-              b[1]=-kG1*x[0];
-            }
-          } else if(x[2]>kCORBEG) {
-            if(x[2]<kCOREND) {
-//    Corrector dipole (because of dimuon arm)
-//            b[0]=kFCORN;
-              b[0]=-kFCORN;
-            }
-          }
-        }
+    Float_t rad2=x[0]*x[0]+x[1]*x[1];
+    if(x[2]>kCORBEG2 && x[2]<kCOREND2){
+      if(rad2<kCOR2RA2){
+        b[0] = kFCORN2;
       }
     }
+    else if(x[2]>kZ1BEG && x[2]<kZ1END){  
+      if(rad2<kZ1RA2){
+        b[0] = -kG1*x[1];
+        b[1] = -kG1*x[0];
+      }
+    }
+    else if(x[2]>kZ2BEG && x[2]<kZ2END){  
+      if(rad2<kZ2RA2){
+        b[0] = kG1*x[1];
+        b[1] = kG1*x[0];
+      }
+    }
+    else if(x[2]>kZ3BEG && x[2]<kZ3END){  
+      if(rad2<kZ3RA2){
+        b[0] = kG1*x[1];
+        b[1] = kG1*x[0];
+      }
+    }
+    else if(x[2]>kZ4BEG && x[2]<kZ4END){  
+      if(rad2<kZ4RA2){
+        b[0] = -kG1*x[1];
+        b[1] = -kG1*x[0];
+      }
+    }
+    else if(x[2]>kD1BEG && x[2]<kD1END){ 
+      if(rad2<kD1RA2){
+        b[1] = -kFDIP;
+      }
+    }
+    else if(x[2]>kD2BEG && x[2]<kD2END){
+      if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
+        || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
+	b[1] = kFDIP;
+      }
+    }
+      }
+    
 
   if(fFactor!=1) {
     b[0]*=fFactor;
@@ -298,200 +374,180 @@ if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2]
   }
 }
 
-//_________________________________________
-
-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)
-  
+//_______________________________________________________________________
+Int_t AliMagFDM::FZ(Double_t temp, Float_t *Ar, 
+                    Float_t delu, Int_t ik,Int_t nk) const
 {
   //
-  // Z component of the field
+  // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
   //
-  static  const Double_t kone=1;
-  Int_t l,ik,ikj;
-  Double_t temp;
-  Double_t ddu,delu,ar;
-
-  Int_t nk,ku;
-  temp=*u;
-  nk=*nu;
-  ik=*ki;
-  delu=*du; 
+  Int_t l,ikj;
+  Double_t ddu,ar;
 
+  Int_t ku,kf;
+  kf=0;
   ar=Ar[ik];
   ddu=temp-ar;
 
   ku=int(ddu/delu);
-  ikj=ik+ku;
+  ikj=ik+ku+1;
   if (ddu<=0) ikj=0;
 
    for(l=ikj; l<nk; l++)
    {
 
-    if(temp < Ar[l])
+    if(temp <= Ar[l])
      {
-      *kf=l;
-      *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
-      *a1= kone - *a2;
-      break;     
+       
+      kf=l-1;
+      break;
      }
     }
+    return kf;
   }
 
-/*-------------FRfuncBi----------------*/
-
-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)
+//_______________________________________________________________________
+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)
 {
   //
-  // This method needs to be commented
+  // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
   //
   Double_t fa11,fa12,fa13;
   Double_t fa21,fa22,fa23;
   Double_t faY1,faY2;
   Double_t bba;
-  
-  Double_t zaa1,zaa2,alf1,alf2,alf3;
-  Int_t kaai,kaa,maa;
-  kaai=*kai;
-  kaa=*ka;
-  maa=*ma;
-  zaa1=*za1;
-  zaa2=*za2;
-  alf1=*al1;  
-  alf2=*al2;
-  alf3=*al3; 
-  
-  if (kaai==0 ) {
-    fa11 = fdBpx[kaa][0][0];
-    fa12 = fdBpx[kaa][0][maa];
-    fa13 = fdBpx[kaa][0][maa+1];
-    fa21 = fdBpx[kaa+1][0][0];
-    fa22 = fdBpx[kaa+1][0][maa];               
-    fa23 = fdBpx[kaa+1][0][maa+1]; 
-  }
-  if (kaai==1 ) {
-    fa11 = fdBpy[kaa][0][0];
-    fa12 = fdBpy[kaa][0][maa];
-    fa13 = fdBpy[kaa][0][maa+1];
-    fa21 = fdBpy[kaa+1][0][0];
-    fa22 = fdBpy[kaa+1][0][maa];               
-    fa23 = fdBpy[kaa+1][0][maa+1]; 
-  }
-  if (kaai==2 ) {
-    fa11 = fdBpz[kaa][0][0];
-    fa12 = fdBpz[kaa][0][maa];
-    fa13 = fdBpz[kaa][0][maa+1];
-    fa21 = fdBpz[kaa+1][0][0];
-    fa22 = fdBpz[kaa+1][0][maa];               
-    fa23 = fdBpz[kaa+1][0][maa+1]; 
+
+  switch (kaai) {
+  case 0:
+    fa11 = fBpx[kaa][0][0];
+    fa12 = fBpx[kaa][0][maa];
+    fa13 = fBpx[kaa][0][maa+1];
+    fa21 = fBpx[kaa+1][0][0];
+    fa22 = fBpx[kaa+1][0][maa];               
+    fa23 = fBpx[kaa+1][0][maa+1]; 
+    break;
+  case 1:
+    fa11 = fBpy[kaa][0][0];
+    fa12 = fBpy[kaa][0][maa];
+    fa13 = fBpy[kaa][0][maa+1];
+    fa21 = fBpy[kaa+1][0][0];
+    fa22 = fBpy[kaa+1][0][maa];               
+    fa23 = fBpy[kaa+1][0][maa+1]; 
+    break;
+  case 2:
+    fa11 = fBpz[kaa][0][0];
+    fa12 = fBpz[kaa][0][maa];
+    fa13 = fBpz[kaa][0][maa+1];
+    fa21 = fBpz[kaa+1][0][0];
+    fa22 = fBpz[kaa+1][0][maa];               
+    fa23 = fBpz[kaa+1][0][maa+1]; 
+    break;
+  default:
+    Fatal("Ba","Invalid value of kaai %d\n",kaai);
+    exit(1);
   }                            
   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
   bba =  zaa1*faY1+zaa2*faY2;    
-  *ba=bba;
-  
+  return bba;  
 }
 
   
-/*----------- FGfuncBi------------*/
-
-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)
+//_______________________________________________________________________
+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)
 {  
   //
-  // This method needs to be commented
+  // Calculation of field componet at a given axis (general case)
   //
   Double_t fy1, fy2, ffy;
   Double_t gy1,gy2,ggy;
-  Double_t z1,z2,y1,y2,x1,x2;
-  
-  Int_t k,l,m,kv;
   Double_t bbi;
-  
   Double_t bf11,bf12,bf21,bf22;
   Double_t bg11,bg12,bg21,bg22;
-  k=*kk;
-  l=*ll;
-  m=*mm;
-  
-  kv=*kvr;
-  
-  z1=*zz1;
-  z2=*zz2;
-  y1=*yy1;
-  y2=*yy2;
-  x1=*xx1;
-  x2=*xx2; 
+
   
   /*-----------------Polar part ------------------*/
   
-  if(kv==0) {
-    bf11=fdBpx[k][l][m];
-    bf12=fdBpx[k+1][l][m];
-    bf21=fdBpx[k+1][l+1][m];
-    bf22=fdBpx[k][l+1][m];
+  switch (kv) {
+  case 0:
+    bf11=fBpx[k][l][m];
+    bf12=fBpx[k+1][l][m];
+    bf21=fBpx[k+1][l+1][m];
+    bf22=fBpx[k][l+1][m];
     
-    bg11=fdBpx[k][l][m+1];
-    bg12=fdBpx[k+1][l][m+1];
-    bg21=fdBpx[k+1][l+1][m+1];
-    bg22=fdBpx[k][l+1][m+1];
-  }
-  if(kv==1) {
-    bf11=fdBpy[k][l][m];
-    bf12=fdBpy[k+1][l][m];
-    bf21=fdBpy[k+1][l+1][m];
-    bf22=fdBpy[k][l+1][m];
+    bg11=fBpx[k][l][m+1];
+    bg12=fBpx[k+1][l][m+1];
+    bg21=fBpx[k+1][l+1][m+1];
+    bg22=fBpx[k][l+1][m+1];
+    break;
+
+  case 1:
+    bf11=fBpy[k][l][m];
+    bf12=fBpy[k+1][l][m];
+    bf21=fBpy[k+1][l+1][m];
+    bf22=fBpy[k][l+1][m];
     
-    bg11=fdBpy[k][l][m+1];
-    bg12=fdBpy[k+1][l][m+1];
-    bg21=fdBpy[k+1][l+1][m+1];
-    bg22=fdBpy[k][l+1][m+1];
-  }  
-  
-  if(kv==2) {
-    bf11=fdBpz[k][l][m];
-    bf12=fdBpz[k+1][l][m];
-    bf21=fdBpz[k+1][l+1][m];
-    bf22=fdBpz[k][l+1][m];
+    bg11=fBpy[k][l][m+1];
+    bg12=fBpy[k+1][l][m+1];
+    bg21=fBpy[k+1][l+1][m+1];
+    bg22=fBpy[k][l+1][m+1];
+    break;
+
+  case 2:
+    bf11=fBpz[k][l][m];
+    bf12=fBpz[k+1][l][m];
+    bf21=fBpz[k+1][l+1][m];
+    bf22=fBpz[k][l+1][m];
     
-    bg11=fdBpz[k][l][m+1];
-    bg12=fdBpz[k+1][l][m+1];
-    bg21=fdBpz[k+1][l+1][m+1];
-    bg22=fdBpz[k][l+1][m+1];           
-  } 
+    bg11=fBpz[k][l][m+1];
+    bg12=fBpz[k+1][l][m+1];
+    bg21=fBpz[k+1][l+1][m+1];
+    bg22=fBpz[k][l+1][m+1]; 
+    break;
   /*-----------------Cartensian part ---------------*/ 
   
-  if(kv==3) {
-    bf11=fdBcx[k][l][m];
-    bf12=fdBcx[k+1][l][m];
-    bf21=fdBcx[k+1][l+1][m];
-    bf22=fdBcx[k][l+1][m];
+  case 3:
+    bf11=fBcx[k][l][m];
+    bf12=fBcx[k+1][l][m];
+    bf21=fBcx[k+1][l+1][m];
+    bf22=fBcx[k][l+1][m];
     
-    bg11=fdBcx[k][l][m+1];
-    bg12=fdBcx[k+1][l][m+1];
-    bg21=fdBcx[k+1][l+1][m+1];
-    bg22=fdBcx[k][l+1][m+1];
-  }          
-  
-  if(kv==4) {
-    bf11=fdBcy[k][l][m];
-    bf12=fdBcy[k+1][l][m];
-    bf21=fdBcy[k+1][l+1][m];
-    bf22=fdBcy[k][l+1][m];
+    bg11=fBcx[k][l][m+1];
+    bg12=fBcx[k+1][l][m+1];
+    bg21=fBcx[k+1][l+1][m+1];
+    bg22=fBcx[k][l+1][m+1];
+    break;
+
+  case 4:
+    bf11=fBcy[k][l][m];
+    bf12=fBcy[k+1][l][m];
+    bf21=fBcy[k+1][l+1][m];
+    bf22=fBcy[k][l+1][m];
     
-    bg11=fdBcy[k][l][m+1];
-    bg12=fdBcy[k+1][l][m+1];
-    bg21=fdBcy[k+1][l+1][m+1];
-    bg22=fdBcy[k][l+1][m+1];
-  }          
-  if(kv==5) {
-    bf11=fdBcz[k][l][m];
-    bf12=fdBcz[k+1][l][m];
-    bf21=fdBcz[k+1][l+1][m];
-    bf22=fdBcz[k][l+1][m];
+    bg11=fBcy[k][l][m+1];
+    bg12=fBcy[k+1][l][m+1];
+    bg21=fBcy[k+1][l+1][m+1];
+    bg22=fBcy[k][l+1][m+1];
+    break;
+
+  case 5:
+    bf11=fBcz[k][l][m];
+    bf12=fBcz[k+1][l][m];
+    bf21=fBcz[k+1][l+1][m];
+    bf22=fBcz[k][l+1][m];
     
-    bg11=fdBcz[k][l][m+1];
-    bg12=fdBcz[k+1][l][m+1];
-    bg21=fdBcz[k+1][l+1][m+1];
-    bg22=fdBcz[k][l+1][m+1];
+    bg11=fBcz[k][l][m+1];
+    bg12=fBcz[k+1][l][m+1];
+    bg21=fBcz[k+1][l+1][m+1];
+    bg22=fBcz[k][l+1][m+1];
+    break;
+
+  default:
+    Fatal("Bb","Invalid value of kv %d\n",kv);
+    exit(1);
   }  
   
   
@@ -506,12 +562,12 @@ void AliMagFDM::FGfuncBi(Double_t *zz1,Double_t *zz2, Double_t *yy1,Double_t *yy
   ggy= y1*gy1 +  y2*gy2;
   
   bbi =  x1*ffy+x2*ggy;
-  
-  *bb=bbi; 
+
+  return bbi;
   
 }     
-//____________________________________________
 
+//_______________________________________________________________________
 void AliMagFDM::ReadField()
 {
   //
@@ -535,41 +591,41 @@ void AliMagFDM::ReadField()
   
 //  Cartensian part 
  
-    fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl); 
+    fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl); 
     
-    printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);     
+    printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);     
     
-    for (ik=0; ik<fdZl; ik++)
+    for (ik=0; ik<fZl; ik++)
     { 
    
       fscanf(magfile, " %e ", &zz);
-      fdZc[ik]=zz; 
+      fZc[ik]=zz; 
 
     } 
    
-    for (ik=0; ik<fdYl; ik++)
+    for (ik=0; ik<fYl; ik++)
     {    
        fscanf(magfile, " %e ", &yy); 
-       fdY[ik]=yy;
+       fY[ik]=yy;
  
     } 
     for (ik=0; ik<81; ik++)
     {    
-           printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]); 
+           printf("fZc %e,fY %e\n", fZc[ik],fY[ik]); 
     }   
              
-    fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2); 
+    fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2); 
 
-printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
-printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
-printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
+printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
+printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
+printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
 
     for (il=0; il<44; il++)  { 
      for (im=0; im<81; im++)  {      
       for (ik=0; ik<81; ik++)  {      
       
       fscanf(magfile, " %e ", &by); 
-      fdBcy[ik][im][il]=by;        
+      fBcy[ik][im][il]=by;        
       }
      }
     } 
@@ -579,7 +635,7 @@ printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
       for (ik=0; ik<81; ik++)  {      
       
       fscanf(magfile, " %e ", &bx); 
-      fdBcx[ik][im][il]=bx;        
+      fBcx[ik][im][il]=bx;        
       }    
      }     
     }
@@ -589,53 +645,53 @@ printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
       for (ik=0; ik<81; ik++)  {      
       
       fscanf(magfile, " %e ", &bz); 
-      fdBcz[ik][im][il]=bz;          
+      fBcz[ik][im][il]=bz;          
       }              
      }     
     } 
 //----------------------   Polar part ---------------------------------
 
     printf("Polar part\n");
-    fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin); 
-    printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);   
+    fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin); 
+    printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);   
 
-    printf(" fdZp array\n"); 
+    printf(" fZp array\n"); 
      
     for (ik=0; ik<51; ik++) 
     {    
      fscanf(magfile, " %e ", &zzp);
-     fdZp[ik]=zzp; 
-     printf(" %e\n",fdZp[ik]);      
+     fZp[ik]=zzp; 
+     printf(" %e\n",fZp[ik]);      
     } 
   
-    printf(" fdR array\n"); 
+    printf(" fR array\n"); 
          
     for (ik=0; ik<10; ik++) 
     {    
      fscanf(magfile, " %e ", &rr); 
-     fdR[ik]=rr;
-     printf(" %e\n",fdR[ik]);
+     fR[ik]=rr;
+     printf(" %e\n",fR[ik]);
     } 
     
-//    printf("fdPhi array\n"); 
+//    printf("fPhi array\n"); 
      
      for (il=0; il<33; il++)  
      {
        fscanf(magfile, " %e ", &phii); 
-       fdPhi[il]=phii; 
-//        printf(" %e\n",fdPhi[il]);          
+       fPhi[il]=phii; 
+//        printf(" %e\n",fPhi[il]);          
      }
 
-    fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin); 
+    fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin); 
 
-printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
+printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
 
                       
     for (il=0; il<33; il++)  { 
      for (im=0; im<10; im++)  {      
       for (ik=0; ik<51; ik++)  {
       fscanf(magfile, " %e ", &by); 
-        fdBpy[ik][im][il]=by;        
+        fBpy[ik][im][il]=by;        
       }
      }
     } 
@@ -644,7 +700,7 @@ printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %
      for (im=0; im<10; im++)  {      
       for (ik=0; ik<51; ik++)  {
       fscanf(magfile, " %e ", &bx); 
-        fdBpx[ik][im][il]=bx;                     
+        fBpx[ik][im][il]=bx;                     
       }    
      }     
     }      
@@ -654,7 +710,7 @@ printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %
      for (im=0; im<10; im++)  {      
       for (ik=0; ik<51; ik++)  {
       fscanf(magfile, " %e ", &bz); 
-        fdBpz[ik][im][il]=bz;                      
+        fBpz[ik][im][il]=bz;                      
       }              
      }     
     } 
@@ -664,7 +720,7 @@ printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %
      for (im=0; im<2; im++)  {      
       for (ik=0; ik<2; ik++)  {
       fscanf(magfile, " %e ", &bb);    
-        fdB[ik][im][il]=bb;  
+        fB[ik][im][il]=bb;  
       }              
      } 
     }
@@ -674,4 +730,3 @@ printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %
     exit(1);
   }
 }
-//________________________________