[u/mrichter/AliRoot.git] / STEER / AliMagFDM.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
Revision 1.5  2000/10/27 14:17:04  morsch
- Bug causing segmentation violation during muon reconstruction corrected
- Coding rule violations corrected.
(Galina Chabratova)

Revision 1.4  2000/10/02 21:28:14  fca
Removal of useless dependecies via forward declarations

Revision 1.3  2000/07/13 16:19:09  fca
Mainly coding conventions + some small bug fixes

Revision 1.2  2000/07/12 08:56:25  fca
Coding convention correction and warning removal

Revision 1.1  2000/07/11 18:24:59  fca
Coding convention corrections + few minor bug fixes

*/

#include <stdlib.h>

#include "AliMagFDM.h"
#include "TSystem.h"
    

ClassImp(AliMagFDM)

//________________________________________
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)
  
{
  //
  // 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());
  
}

//________________________________________

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    kfZbg=502.92;  // Start of Map using in z
  static  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; 
  Double_t zz1, zz2,yy1,yy2,x2,x1; 
     
// --- start the map fiel from z = 502.92 cm ---

  x[0] = xfi[0];
  x[1] = xfi[1];
  x[2] = xfi[2];
    
  /*  
  if (x[0]==0) x[0]=0.1E-06;
  if (x[1]==0) x[1]=0.1E-06;
  */ 
 b[0]=b[1]=b[2]=0;
  //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]); 

  Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
           r0=rr/100;
  Double_t rPmax;
  rPmax=fRmax; 
  if ( (-700<x[2] && x[2]<=kfZbg && 
        (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
       || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
       {
        b[0]=b[1]=0;
        b[2]=2;
       }

  xL3[0]=x[0]/100;
  xL3[1]=(x[1]+30)/100;
  xL3[2]=x[2]/100; 

 
  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=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 ((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 (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
	//       if (xL3[2]<zCmin && r0<rPmax)   
       {
       //---------------------  Polar part ----------------------
       
       Double_t yyp,ph0;
       Int_t  kp0, lp0, mp0;
       Int_t kpi,lpi,mpi;
       Double_t alp1,alp2,alp3;
       Double_t zpz,rp,fip,cphi; 

       kpi=kiip; 
       lpi=kliip;
       mpi=kmiip;
   
       zpz=xL3[2];

       FZ(&zpz, fZp ,&fZpdl,&kpi,&kp0,&zp1 ,&zp2,&fZpl) ;

       yyp=xL3[1]- 0.3;
       cphi=yyp/r0;
       ph0=TMath::ACos(cphi);
       if (xL3[0]< 0) {ph0=kpi2 - ph0;}  
                                  
       fip=ph0;
       FZ(&fip,fPhi,&fPhid ,&mpi,&mp0, &xp1,&xp2,&fPhin); 

       Double_t rDel;
       rDel=fRdel;

       if (r0<= fRdel)
        {

         if(r0< keps) 
         {

          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= 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 ;
         } 
        } // end   of keps <=r0 
       }
        else
       {
        rp=r0;

        FZ(&rp,fR ,&fRdel,&lpi,&lp0,&yp1,&yp2,&fRn);

        for (jb=0; jb<3 ; jb++) 
        {
         iKvar=jb;
         FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj); 

         bint[jb] = bbj*10 ;
        }
       }

       b[0]=bint[0];
       b[1]=bint[1];
       b[2]=bint[2];

//    fprintf(fitest,"-------------   Freg2 run -------------\n");       
  
   }  
   else 
   {
   //-------------- Cartensian part ------------------
          
   Double_t zzc,yyc; 
   Int_t k0, l0,m0;
   Double_t  xx1, xx2,dx, xxx ,xXl; 
   Int_t kci,mci,lci;

   kci=kiic;
   lci=kliic;
   mci=kmiic;

   xx1 = fAx1*xL3[2] + fCx1;
   xx2 = fAx2*xL3[2] + fCx2;

   zzc=xL3[2];
   FZ(&zzc, fZc ,&fZdel, &kci,&k0, &zz1, &zz2, &fZl);     

   yyc=xL3[1];
   FZ(&yyc, fY , &fYdel,&lci, &l0, &yy1, &yy2,&fYl);  

   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)) 
    {
      m0=m0+1;   
      printf(" m0 %d, m0+1 %d\n",m0,m0+1);  
    }

   x2=(xL3[0]-( xx1+m0*dx))/dx;
   x1=kone-x2;
   m0=m0-1;
   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 ; 
     }    
 
   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);
 }
           
} else {

//This is the ZDC part
    Float_t rad2=x[0]*x[0]+x[1]*x[1];
    if(x[2]>kCORBEG1 && x[2]<kCOREND1){
      if(rad2<kCOR1RA2){
        b[0] = kFCORN1;
      }
    }
    else 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;
    b[1]*=fFactor;
    b[2]*=fFactor;
  }
}

//_________________________________________

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)
{
  //
  // Z component of the field
  //
  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; 

  ar=Ar[ik];
  ddu=temp-ar;

  ku=int(ddu/delu);
  //  ikj=ik+ku; //falsh
    ikj=ik+ku+1; //correct
  if (ddu<=0) ikj=0;

   for(l=ikj; l<nk; l++)
   {
     //    if(temp < Ar[l])
    if(temp <= Ar[l])
     {
       
      *kf=l-1;
      *a2=(temp-Ar[l-1])/(Ar[l]-Ar[l-1]);
      
       /*
      *kf=l;
      *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
      */
       break;
     }
    }
   //    *a2=(temp-Ar[*kf])/(Ar[*kf+1]-Ar[*kf]);
    *a1= kone - *a2;
  }

/*-------------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)
{
  //
  // This method needs to be commented
  //
  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; 
  
  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("FRfuncBi","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;
  
}

  
/*----------- 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)
{  
  //
  // This method needs to be commented
  //
  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 ------------------*/
  
  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=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=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=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 ---------------*/ 
  
  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=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=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=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("FGfuncBi","Invalid value of kv %d\n",kv);
    exit(1);
  }  
  
  
  fy1=z1*bf11+z2*bf12;
  fy2=z2*bf21+z1* bf22;
  ffy=y1*fy1+ y2*fy2; 
  
  
  gy1 = z1*bg11+z2*bg12;
  gy2 = z2*bg21+z1*bg22;
  ggy= y1*gy1 +  y2*gy2;
  
  bbi =  x1*ffy+x2*ggy;
  
  *bb=bbi; 
  
}     
//____________________________________________

void AliMagFDM::ReadField()
{
  //
  // Method to read the magnetic field from file
  //
  FILE *magfile;

  Int_t ik, il, im;
  Float_t zzp, rr,phii;
  Float_t zz, yy, bx,by,bz,bb;

  char *fname;
  printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
  fname = gSystem->ExpandPathName(fTitle.Data());
  magfile=fopen(fname,"r");
  delete [] fname;

  printf("Cartensian part\n");
 
  if (magfile) {
  
//  Cartensian part 
 
    fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl); 
    
    printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);     
    
    for (ik=0; ik<fZl; ik++)
    { 
   
      fscanf(magfile, " %e ", &zz);
      fZc[ik]=zz; 

    } 
   
    for (ik=0; ik<fYl; ik++)
    {    
       fscanf(magfile, " %e ", &yy); 
       fY[ik]=yy;
 
    } 
    for (ik=0; ik<81; ik++)
    {    
           printf("fZc %e,fY %e\n", fZc[ik],fY[ik]); 
    }   
             
    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("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); 
      fBcy[ik][im][il]=by;        
      }
     }
    } 

    for (il=0; il<44; il++)  { 
     for (im=0; im<81; im++)  {      
      for (ik=0; ik<81; ik++)  {      
      
      fscanf(magfile, " %e ", &bx); 
      fBcx[ik][im][il]=bx;        
      }    
     }     
    }
    
   for (il=0; il<44; il++)  { 
     for (im=0; im<81; im++)  {      
      for (ik=0; ik<81; ik++)  {      
      
      fscanf(magfile, " %e ", &bz); 
      fBcz[ik][im][il]=bz;          
      }              
     }     
    } 
//----------------------   Polar part ---------------------------------

    printf("Polar part\n");
    fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin); 
    printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);   

    printf(" fZp array\n"); 
     
    for (ik=0; ik<51; ik++) 
    {    
     fscanf(magfile, " %e ", &zzp);
     fZp[ik]=zzp; 
     printf(" %e\n",fZp[ik]);      
    } 
  
    printf(" fR array\n"); 
         
    for (ik=0; ik<10; ik++) 
    {    
     fscanf(magfile, " %e ", &rr); 
     fR[ik]=rr;
     printf(" %e\n",fR[ik]);
    } 
    
//    printf("fPhi array\n"); 
     
     for (il=0; il<33; il++)  
     {
       fscanf(magfile, " %e ", &phii); 
       fPhi[il]=phii; 
//        printf(" %e\n",fPhi[il]);          
     }

    fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin); 

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); 
        fBpy[ik][im][il]=by;        
      }
     }
    } 
    
    for (il=0; il<33; il++)  { 
     for (im=0; im<10; im++)  {      
      for (ik=0; ik<51; ik++)  {
      fscanf(magfile, " %e ", &bx); 
        fBpx[ik][im][il]=bx;                     
      }    
     }     
    }      


    for (il=0; il<33; il++)  { 
     for (im=0; im<10; im++)  {      
      for (ik=0; ik<51; ik++)  {
      fscanf(magfile, " %e ", &bz); 
        fBpz[ik][im][il]=bz;                      
      }              
     }     
    } 
    
    
    for (il=0; il<32; il++) { 
     for (im=0; im<2; im++)  {      
      for (ik=0; ik<2; ik++)  {
      fscanf(magfile, " %e ", &bb);    
        fB[ik][im][il]=bb;  
      }              
     } 
    }
//
  } else { 
    printf("File %s not found !\n",fTitle.Data());
    exit(1);
  }
}
//________________________________
Commit	Line	Data
aee8290b	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$
d95ea23f	18	Revision 1.5 2000/10/27 14:17:04 morsch
	19	- Bug causing segmentation violation during muon reconstruction corrected
	20	- Coding rule violations corrected.
	21	(Galina Chabratova)
	22
803c9752	23	Revision 1.4 2000/10/02 21:28:14 fca
	24	Removal of useless dependecies via forward declarations
	25
94de3818	26	Revision 1.3 2000/07/13 16:19:09 fca
	27	Mainly coding conventions + some small bug fixes
	28
ef42d733	29	Revision 1.2 2000/07/12 08:56:25 fca
	30	Coding convention correction and warning removal
	31
8918e700	32	Revision 1.1 2000/07/11 18:24:59 fca
	33	Coding convention corrections + few minor bug fixes
	34
aee8290b	35	*/
aee8290b	36
8918e700	37	#include <stdlib.h>
8918e700	38
aee8290b	39	#include "AliMagFDM.h"
	40	#include "TSystem.h"
	41
	42
	43	ClassImp(AliMagFDM)
	44
	45	//________________________________________
	46	AliMagFDM::AliMagFDM(const char name, const char title, const Int_t integ,
	47	const Int_t map, const Float_t factor, const Float_t fmax)
	48	: AliMagF(name,title,integ,map,factor,fmax)
	49
	50	{
ef42d733	51	//
	52	// Standard constructor for the Dipole field
	53	//
aee8290b	54	fType = kDipoMap;
	55
	56	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());
	57
	58	}
	59
	60	//________________________________________
	61
	62	void AliMagFDM::Field(Float_t xfi, Float_t b)
	63	{
	64	//
	65	// Main routine to compute the field in a point
	66	//
	67	static const Double_t keps=0.1E-06;
	68	static const Double_t kpi2=.6283185E+01;
	69	static const Double_t kone=1;
	70
	71	static const Int_t kiip=33;
	72	static const Int_t kmiip=0;
	73	static const Int_t kliip=0;
	74
	75	static const Int_t kiic=0;
	76	static const Int_t kmiic=0;
	77	static const Int_t kliic=0;
	78
803c9752	79	static const Double_t kfZbg=502.92; // Start of Map using in z
803c9752	80	static const Double_t kfZL3=600; // Beginning of L3 door in z
aee8290b	81
	82	Double_t x[3];
	83	Double_t xL3[3];
	84	Double_t bint[3];
	85
	86	Double_t r0;
	87
	88	Double_t bbj;
	89	Int_t iKvar,jb;
	90
	91	Double_t zp1, zp2,xp1,xp2,yp1,yp2;
	92	Double_t zz1, zz2,yy1,yy2,x2,x1;
	93
	94	// --- start the map fiel from z = 502.92 cm ---
	95
	96	x[0] = xfi[0];
	97	x[1] = xfi[1];
	98	x[2] = xfi[2];
803c9752	99
	100	/*
	101	if (x[0]==0) x[0]=0.1E-06;
	102	if (x[1]==0) x[1]=0.1E-06;
	103	*/
	104	b[0]=b[1]=b[2]=0;
aee8290b	105	// printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
	106
	107	Double_t rr=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
	108	r0=rr/100;
	109	Double_t rPmax;
803c9752	110	rPmax=fRmax;
803c9752	111	if ( (-700<x[2] && x[2]<=kfZbg &&
aee8290b	112	(x[0]x[0]+(x[1]+30)(x[1]+30))< 560*560)
803c9752	113	\|\| (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
aee8290b	114	{
803c9752	115	b[0]=b[1]=0;
aee8290b	116	b[2]=2;
	117	}
	118
	119	xL3[0]=x[0]/100;
	120	xL3[1]=(x[1]+30)/100;
	121	xL3[2]=x[2]/100;
803c9752	122
aee8290b	123
803c9752	124	if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
	125	if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
	126
	127	Double_t xminn=xL3[2]*fAx1+fCx1;
	128	Double_t xmaxx=xL3[2]*fAx2+fCx2;
aee8290b	129	Double_t zCmin,zCmax,yCmin,yCmax;
aee8290b	130
803c9752	131	zCmin=fZmin;
	132	zCmax=fZmax;
	133	yCmin=fYmin;
	134	yCmax=fYmax;
	135	/*
	136	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)))
	137	*/
	138	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)))
aee8290b	139	{
	140	if(fMap==3)
	141	{
803c9752	142	if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
803c9752	143	// if (xL3[2]<zCmin && r0<rPmax)
aee8290b	144	{
	145	//--------------------- Polar part ----------------------
	146
	147	Double_t yyp,ph0;
	148	Int_t kp0, lp0, mp0;
	149	Int_t kpi,lpi,mpi;
	150	Double_t alp1,alp2,alp3;
	151	Double_t zpz,rp,fip,cphi;
	152
	153	kpi=kiip;
	154	lpi=kliip;
	155	mpi=kmiip;
	156
	157	zpz=xL3[2];
	158
803c9752	159	FZ(&zpz, fZp ,&fZpdl,&kpi,&kp0,&zp1 ,&zp2,&fZpl) ;
aee8290b	160
	161	yyp=xL3[1]- 0.3;
	162	cphi=yyp/r0;
	163	ph0=TMath::ACos(cphi);
	164	if (xL3[0]< 0) {ph0=kpi2 - ph0;}
	165
	166	fip=ph0;
803c9752	167	FZ(&fip,fPhi,&fPhid ,&mpi,&mp0, &xp1,&xp2,&fPhin);
aee8290b	168
aee8290b	169	Double_t rDel;
803c9752	170	rDel=fRdel;
aee8290b	171
803c9752	172	if (r0<= fRdel)
aee8290b	173	{
	174
	175	if(r0< keps)
	176	{
	177
803c9752	178	bint[0]=(zp1fBpx[kp0][0][0] + zp2fBpx[kp0+1][0][0])*10;
	179	bint[1]=(zp1fBpy[kp0][0][0] + zp2fBpy[kp0+1][0][0])*10;
	180	bint[2]=(zp1fBpz[kp0][0][0] + zp2fBpz[kp0+1][0][0])*10;
aee8290b	181
803c9752	182	} else {
aee8290b	183
803c9752	184	alp2= fB[0][0][mp0]yyp + fB[0][1][mp0]xL3[0];
803c9752	185	alp3= fB[1][0][mp0]yyp + fB[1][1][mp0]xL3[0];
aee8290b	186	alp1= kone - alp2 - alp3;
	187
	188	for (jb=0; jb<3 ; jb++)
803c9752	189	{
	190	iKvar=jb;
	191	FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
	192	bint[jb] = bbj*10 ;
	193	}
	194	} // end of keps <=r0
aee8290b	195	}
	196	else
	197	{
	198	rp=r0;
	199
803c9752	200	FZ(&rp,fR ,&fRdel,&lpi,&lp0,&yp1,&yp2,&fRn);
aee8290b	201
	202	for (jb=0; jb<3 ; jb++)
	203	{
	204	iKvar=jb;
	205	FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj);
	206
	207	bint[jb] = bbj*10 ;
	208	}
	209	}
	210
	211	b[0]=bint[0];
	212	b[1]=bint[1];
	213	b[2]=bint[2];
	214
	215	// fprintf(fitest,"------------- Freg2 run -------------\n");
	216
	217	}
	218	else
	219	{
	220	//-------------- Cartensian part ------------------
	221
	222	Double_t zzc,yyc;
	223	Int_t k0, l0,m0;
	224	Double_t xx1, xx2,dx, xxx ,xXl;
	225	Int_t kci,mci,lci;
	226
	227	kci=kiic;
	228	lci=kliic;
	229	mci=kmiic;
	230
803c9752	231	xx1 = fAx1*xL3[2] + fCx1;
803c9752	232	xx2 = fAx2*xL3[2] + fCx2;
aee8290b	233
aee8290b	234	zzc=xL3[2];
803c9752	235	FZ(&zzc, fZc ,&fZdel, &kci,&k0, &zz1, &zz2, &fZl);
aee8290b	236
aee8290b	237	yyc=xL3[1];
803c9752	238	FZ(&yyc, fY , &fYdel,&lci, &l0, &yy1, &yy2,&fYl);
aee8290b	239
803c9752	240	xXl = fXl-kone;
aee8290b	241	dx = (xx2-xx1)/xXl;
	242	xxx= xL3[0]-xx1;
	243	// xm = xxx/dx;
	244	m0 = int(xxx/dx);
	245
	246	if(xL3[0]<(xx1+m0dx) \|\| xL3[0] >(xx1+(m0+1)dx))
	247	{
	248	m0=m0+1;
	249	printf(" m0 %d, m0+1 %d\n",m0,m0+1);
	250	}
	251
	252	x2=(xL3[0]-( xx1+m0*dx))/dx;
	253	x1=kone-x2;
	254	m0=m0-1;
	255	for (jb=3; jb<6; jb++)
	256	{
	257	iKvar=jb;
	258	FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj);
	259	bint[jb-3] = bbj*10 ;
	260	}
	261
	262	b[0]=bint[0];
	263	b[1]=bint[1];
	264	b[2]=bint[2];
	265
	266	// fprintf(fitest,"------------ Freg1 run -----------------\n");
	267	}
	268
	269	} else {
	270	printf("Unknown map of Dipole region %d\n",fMap);
	271	}
	272
	273	} else {
	274
	275	//This is the ZDC part
d95ea23f	276	Float_t rad2=x[0]x[0]+x[1]x[1];
	277	if(x[2]>kCORBEG1 && x[2]<kCOREND1){
	278	if(rad2<kCOR1RA2){
	279	b[0] = kFCORN1;
	280	}
	281	}
	282	else if(x[2]>kCORBEG2 && x[2]<kCOREND2){
	283	if(rad2<kCOR2RA2){
	284	b[0] = kFCORN2;
	285	}
	286	}
	287	else if(x[2]>kZ1BEG && x[2]<kZ1END){
	288	if(rad2<kZ1RA2){
	289	b[0] = -kG1*x[1];
	290	b[1] = -kG1*x[0];
	291	}
	292	}
	293	else if(x[2]>kZ2BEG && x[2]<kZ2END){
	294	if(rad2<kZ2RA2){
	295	b[0] = kG1*x[1];
	296	b[1] = kG1*x[0];
aee8290b	297	}
aee8290b	298	}
d95ea23f	299	else if(x[2]>kZ3BEG && x[2]<kZ3END){
	300	if(rad2<kZ3RA2){
	301	b[0] = kG1*x[1];
	302	b[1] = kG1*x[0];
	303	}
	304	}
	305	else if(x[2]>kZ4BEG && x[2]<kZ4END){
	306	if(rad2<kZ4RA2){
	307	b[0] = -kG1*x[1];
	308	b[1] = -kG1*x[0];
	309	}
	310	}
	311	else if(x[2]>kD1BEG && x[2]<kD1END){
	312	if(rad2<kD1RA2){
	313	b[1] = -kFDIP;
	314	}
	315	}
	316	else if(x[2]>kD2BEG && x[2]<kD2END){
	317	if(((x[0]-kXCEN1D2)(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)(x[1]-kYCEN1D2))<kD2RA2
	318	\|\| ((x[0]-kXCEN2D2)(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)(x[1]-kYCEN2D2))<kD2RA2){
	319	b[1] = kFDIP;
	320	}
	321	}
	322	}
	323
aee8290b	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)
aee8290b	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);
803c9752	354	// ikj=ik+ku; //falsh
803c9752	355	ikj=ik+ku+1; //correct
aee8290b	356	if (ddu<=0) ikj=0;
	357
	358	for(l=ikj; l<nk; l++)
	359	{
803c9752	360	// if(temp < Ar[l])
803c9752	361	if(temp <= Ar[l])
aee8290b	362	{
803c9752	363
	364	*kf=l-1;
	365	*a2=(temp-Ar[l-1])/(Ar[l]-Ar[l-1]);
	366
	367	/*
aee8290b	368	*kf=l;
aee8290b	369	*a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
803c9752	370	*/
803c9752	371	break;
aee8290b	372	}
aee8290b	373	}
803c9752	374	// a2=(temp-Ar[kf])/(Ar[kf+1]-Ar[kf]);
803c9752	375	a1= kone - a2;
aee8290b	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
8918e700	401	switch (kaai) {
8918e700	402	case 0:
803c9752	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];
8918e700	409	break;
8918e700	410	case 1:
803c9752	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];
8918e700	417	break;
8918e700	418	case 2:
803c9752	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];
8918e700	425	break;
	426	default:
	427	Fatal("FRfuncBi","Invalid value of kaai %d\n",kaai);
	428	exit(1);
aee8290b	429	}
	430	faY1=alf1fa11+alf2fa12+alf3*fa13;
	431	faY2=alf1fa21+alf2fa22+alf3*fa23;
	432	bba = zaa1faY1+zaa2faY2;
	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
8918e700	469	switch (kv) {
8918e700	470	case 0:
803c9752	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];
aee8290b	475
803c9752	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];
8918e700	480	break;
	481
	482	case 1:
803c9752	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];
aee8290b	487
803c9752	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];
8918e700	492	break;
	493
	494	case 2:
803c9752	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];
aee8290b	499
803c9752	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];
8918e700	504	break;
aee8290b	505	/-----------------Cartensian part ---------------/
aee8290b	506
8918e700	507	case 3:
803c9752	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];
aee8290b	512
803c9752	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];
8918e700	517	break;
	518
	519	case 4:
803c9752	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];
aee8290b	524
803c9752	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];
8918e700	529	break;
	530
	531	case 5:
803c9752	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];
aee8290b	536
803c9752	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];
8918e700	541	break;
	542
	543	default:
	544	Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
	545	exit(1);
aee8290b	546	}
	547
	548
	549
	550	fy1=z1bf11+z2bf12;
	551	fy2=z2bf21+z1 bf22;
	552	ffy=y1fy1+ y2fy2;
	553
	554
	555	gy1 = z1bg11+z2bg12;
	556	gy2 = z2bg21+z1bg22;
	557	ggy= y1gy1 + y2gy2;
	558
	559	bbi = x1ffy+x2ggy;
	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
803c9752	589	fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
aee8290b	590
803c9752	591	printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
aee8290b	592
803c9752	593	for (ik=0; ik<fZl; ik++)
aee8290b	594	{
	595
	596	fscanf(magfile, " %e ", &zz);
803c9752	597	fZc[ik]=zz;
aee8290b	598
	599	}
	600
803c9752	601	for (ik=0; ik<fYl; ik++)
aee8290b	602	{
aee8290b	603	fscanf(magfile, " %e ", &yy);
803c9752	604	fY[ik]=yy;
aee8290b	605
	606	}
	607	for (ik=0; ik<81; ik++)
	608	{
803c9752	609	printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
aee8290b	610	}
aee8290b	611
803c9752	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);
aee8290b	613
803c9752	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);
aee8290b	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);
803c9752	623	fBcy[ik][im][il]=by;
aee8290b	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);
803c9752	633	fBcx[ik][im][il]=bx;
aee8290b	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);
803c9752	643	fBcz[ik][im][il]=bz;
aee8290b	644	}
	645	}
	646	}
	647	//---------------------- Polar part ---------------------------------
	648
	649	printf("Polar part\n");
803c9752	650	fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
803c9752	651	printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
aee8290b	652
803c9752	653	printf(" fZp array\n");
aee8290b	654
	655	for (ik=0; ik<51; ik++)
	656	{
	657	fscanf(magfile, " %e ", &zzp);
803c9752	658	fZp[ik]=zzp;
803c9752	659	printf(" %e\n",fZp[ik]);
aee8290b	660	}
aee8290b	661
803c9752	662	printf(" fR array\n");
aee8290b	663
	664	for (ik=0; ik<10; ik++)
	665	{
	666	fscanf(magfile, " %e ", &rr);
803c9752	667	fR[ik]=rr;
803c9752	668	printf(" %e\n",fR[ik]);
aee8290b	669	}
aee8290b	670
803c9752	671	// printf("fPhi array\n");
aee8290b	672
	673	for (il=0; il<33; il++)
	674	{
	675	fscanf(magfile, " %e ", &phii);
803c9752	676	fPhi[il]=phii;
803c9752	677	// printf(" %e\n",fPhi[il]);
aee8290b	678	}
aee8290b	679
803c9752	680	fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
aee8290b	681
803c9752	682	printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
aee8290b	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);
803c9752	689	fBpy[ik][im][il]=by;
aee8290b	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);
803c9752	698	fBpx[ik][im][il]=bx;
aee8290b	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);
803c9752	708	fBpz[ik][im][il]=bz;
aee8290b	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);
803c9752	718	fB[ik][im][il]=bb;
aee8290b	719	}
	720	}
	721	}
	722	//
	723	} else {
	724	printf("File %s not found !\n",fTitle.Data());
	725	exit(1);
	726	}
	727	}
	728	//________________________________