STEER/AliMagFDM.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //   Field with Magnetic Field map
  20 //   Used by AliRun class
  21 //   Author:
  22 //-------------------------------------------------------------------------
  23
  24 #include <TMath.h>
  25 #include <TSystem.h>
  26
  27 #include "AliLog.h"
  28 #include "AliMagFDM.h"
  29
  30 ClassImp(AliMagFDM)
  31
  32 //_______________________________________________________________________
  33 AliMagFDM::AliMagFDM():
  34   fSolenoid(0),
  35   fInd(0),
  36   fZmin(0),
  37   fZmax(0),
  38   fYmax(0),
  39   fYmin(0),
  40   fZpmx(0),
  41   fZpmn(0),
  42   fRmax(0),
  43   fRmin(0),
  44   fXdel(0),
  45   fYdel(0),
  46   fZdel(0),
  47   fRdel(0),
  48   fPhid(0),
  49   fZpdl(0),
  50   fCx1(0),
  51   fCx2(0),
  52   fAx1(0),
  53   fAx2(0),
  54   fXl(0),
  55   fYl(0),
  56   fZl(0),
  57   fRn(0),
  58   fPhin(0),
  59   fZpl(0)
  60 {
  61   //
  62   // Default constructor for the Dipole field
  63   //
  64 }
  65
  66 //_______________________________________________________________________
  67 AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
  68                      Float_t factor, Float_t fmax):
  69   AliMagFC(name,title,integ,factor,fmax),
  70   fSolenoid(0),
  71   fInd(0),
  72   fZmin(0),
  73   fZmax(0),
  74   fYmax(0),
  75   fYmin(0),
  76   fZpmx(0),
  77   fZpmn(0),
  78   fRmax(0),
  79   fRmin(0),
  80   fXdel(0),
  81   fYdel(0),
  82   fZdel(0),
  83   fRdel(0),
  84   fPhid(0),
  85   fZpdl(0),
  86   fCx1(0),
  87   fCx2(0),
  88   fAx1(0),
  89   fAx2(0),
  90   fXl(0),
  91   fYl(0),
  92   fZl(0),
  93   fRn(0),
  94   fPhin(0),
  95   fZpl(0)
  96 {
  97   //
  98   // Standard constructor for the Dipole field
  99   //
 100   fType = kDipoMap;
 101   fMap  = 3;
 102   SetSolenoidField();
 103
 104   AliDebug(1, Form(
 105        "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s",
 106        fName.Data(), fMap ,integ,factor,fTitle.Data()));
 107
 108 }
 109
 110 //_______________________________________________________________________
 111 void AliMagFDM::Field(const float *xfi, float *b) const
 112 {
 113   //
 114   // Main routine to compute the field in a point
 115   //
 116   const Double_t keps=0.1E-06;
 117   const Double_t kPI2=2.*TMath::Pi();
 118   const Double_t kone=1;
 119
 120   const    Int_t  kiip=33;
 121   const    Int_t  kmiip=0;
 122   const    Int_t  kliip=0;
 123
 124   const    Int_t  kiic=0;
 125   const    Int_t  kmiic=0;
 126   const    Int_t  kliic=0;
 127
 128   const Double_t    kfZbg=502.92;  // Start of Map using in z
 129   const Double_t    kfZL3=600;  // Beginning of L3 door in z
 130
 131   Double_t   x[3];
 132   Double_t   xL3[3];
 133   Double_t   bint[3];
 134
 135   Double_t r0;
 136   Int_t iKvar,jb;
 137
 138   Double_t zp1, zp2,xp1,xp2,yp1,yp2;
 139   Double_t zz1, zz2,yy1,yy2,x2,x1;
 140
 141 // --- start the map fiel from z = 502.92 cm ---
 142 //
 143 // This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
 144 // Transfor correspondingly.
 145
 146   x[0] = - xfi[0];
 147   x[1] =   xfi[1];
 148   x[2] = - xfi[2];
 149   b[0]=b[1]=b[2]=0;
 150 //
 151   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]);
 152
 153   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 154            r0=rr/100;
 155   Double_t rPmax;
 156   rPmax=fRmax;
 157   if ( (-700<x[2] && x[2]<=kfZbg &&
 158         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
 159        || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
 160        {
 161         b[0]=b[1]=0;
 162         b[2]=fSolenoid;
 163        }
 164
 165   xL3[0]=x[0]/100;
 166   xL3[1]=(x[1]+30)/100;
 167   xL3[2]=x[2]/100;
 168
 169
 170   if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
 171   if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
 172
 173   Double_t xminn=xL3[2]*fAx1+fCx1;
 174   Double_t xmaxx=xL3[2]*fAx2+fCx2;
 175   Double_t zCmin,zCmax,yCmin,yCmax;
 176
 177   zCmin=fZmin;
 178   zCmax=fZmax;
 179   yCmin=fYmin;
 180   yCmax=fYmax;
 181 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)))
 182     {
 183      if(fMap==3)
 184       {
 185         if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
 186        {
 187        //---------------------  Polar part ----------------------
 188
 189        Double_t yyp,ph0;
 190        Int_t  kp0, lp0, mp0;
 191        Int_t kpi,lpi,mpi;
 192        Double_t alp1,alp2,alp3;
 193        Double_t zpz,rp,fip,cphi;
 194
 195        kpi=kiip;
 196        lpi=kliip;
 197        mpi=kmiip;
 198
 199        zpz=xL3[2];
 200        kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
 201        zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
 202        zp1= kone-zp2;
 203        yyp=xL3[1]- 0.3;
 204        cphi=TMath::Abs(yyp/r0);
 205        Int_t kcphi=0;
 206        if (cphi > kone) {
 207         AliDebug(2,Form("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, cphi %e",xL3[0],xL3[1],xL3[2],yyp,r0,cphi));
 208         cphi =kone;
 209         kcphi=777;
 210        }
 211        ph0=TMath::ACos(cphi);
 212        if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
 213        if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
 214        if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
 215        if (ph0 > kPI2) {       ph0=ph0 - kPI2;}
 216        if (kcphi==777) {
 217         AliDebug(2,Form("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, ph0 %e",xL3[0],xL3[1],xL3[2],yyp,r0,ph0));
 218        }
 219        fip=ph0;
 220        mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
 221        xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
 222        xp1=kone-xp2;
 223
 224        Double_t rDel;
 225        rDel=fRdel;
 226
 227        if (r0<= fRdel)
 228         {
 229
 230          if(r0< keps)
 231          {
 232
 233           bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
 234           bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
 235           bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
 236
 237          } else {
 238
 239         alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
 240         alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
 241         alp1= kone - alp2 - alp3;
 242
 243         for (jb=0; jb<3 ; jb++)
 244          {
 245           iKvar=jb;
 246           bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
 247          }
 248         } // end   of keps <=r0
 249        }
 250         else
 251        {
 252         rp=r0;
 253
 254         lp0=FZ(rp,fR ,fRdel,lpi,fRn);
 255         yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
 256         yp1=kone-yp2;
 257
 258         for (jb=0; jb<3 ; jb++)
 259         {
 260          iKvar=jb;
 261          bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
 262         }
 263        }
 264
 265        b[0]=-bint[0];
 266        b[1]=bint[1];
 267        b[2]=-bint[2];
 268
 269    }
 270    else
 271    {
 272    //-------------- Cartensian part ------------------
 273
 274    Double_t zzc,yyc;
 275    Int_t k0, l0,m0;
 276    Double_t  xx1, xx2,dx, xxx ,xXl;
 277    Int_t kci,mci,lci;
 278
 279    kci=kiic;
 280    lci=kliic;
 281    mci=kmiic;
 282
 283    xx1 = fAx1*xL3[2] + fCx1;
 284    xx2 = fAx2*xL3[2] + fCx2;
 285
 286    zzc=xL3[2];
 287    k0=FZ(zzc, fZc ,fZdel, kci, fZl);
 288    zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
 289    zz1=kone - zz2;;
 290
 291    yyc=xL3[1];
 292    l0=FZ(yyc, fY , fYdel,lci,fYl);
 293    yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
 294    yy1=kone - yy2;
 295    xXl = fXl-kone;
 296    dx = (xx2-xx1)/xXl;
 297    xxx= xL3[0]-xx1;
 298    m0 = int(xxx/dx);
 299
 300    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
 301     {
 302       m0=m0+1;
 303       AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
 304     }
 305
 306    x2=(xL3[0]-( xx1+m0*dx))/dx;
 307    x1=kone-x2;
 308    m0=m0-1;
 309    for (jb=3; jb<6; jb++)
 310      {
 311        iKvar=jb;
 312        bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
 313      }
 314
 315    b[0]=-bint[0];
 316    b[1]=bint[1];
 317    b[2]=-bint[2];
 318
 319    }
 320
 321
 322   } else {
 323         AliError(Form("Unknown map of Dipole region %d",fMap));
 324  }
 325
 326 } else {
 327     ZDCField(xfi,b);
 328
 329     }
 330
 331   if(fFactor!=1) {
 332     b[0]*=fFactor;
 333     b[1]*=fFactor;
 334     b[2]*=fFactor;
 335   }
 336 }
 337
 338 //_______________________________________________________________________
 339 void AliMagFDM::Field(const double *xfi, double *b) const
 340 {
 341   //
 342   // Main routine to compute the field in a point
 343   //
 344   const Double_t keps=0.1E-06;
 345   const Double_t kPI2=2.*TMath::Pi();
 346   const Double_t kone=1;
 347
 348   const    Int_t  kiip=33;
 349   const    Int_t  kmiip=0;
 350   const    Int_t  kliip=0;
 351
 352   const    Int_t  kiic=0;
 353   const    Int_t  kmiic=0;
 354   const    Int_t  kliic=0;
 355
 356   const Double_t    kfZbg=502.92;  // Start of Map using in z
 357   const Double_t    kfZL3=600;  // Beginning of L3 door in z
 358
 359   Double_t   x[3];
 360   Double_t   xL3[3];
 361   Double_t   bint[3];
 362
 363   Double_t r0;
 364   Int_t iKvar,jb;
 365
 366   Double_t zp1, zp2,xp1,xp2,yp1,yp2;
 367   Double_t zz1, zz2,yy1,yy2,x2,x1;
 368
 369 // --- start the map fiel from z = 502.92 cm ---
 370 //
 371 // This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
 372 // Transfor correspondingly.
 373
 374   x[0] = - xfi[0];
 375   x[1] =   xfi[1];
 376   x[2] = - xfi[2];
 377   b[0]=b[1]=b[2]=0;
 378 //
 379   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]);
 380
 381   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 382            r0=rr/100;
 383   Double_t rPmax;
 384   rPmax=fRmax;
 385   if ( (-700<x[2] && x[2]<=kfZbg &&
 386         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
 387        || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
 388        {
 389         b[0]=b[1]=0;
 390         b[2]=fSolenoid;
 391        }
 392
 393   xL3[0]=x[0]/100;
 394   xL3[1]=(x[1]+30)/100;
 395   xL3[2]=x[2]/100;
 396
 397
 398   if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
 399   if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
 400
 401   Double_t xminn=xL3[2]*fAx1+fCx1;
 402   Double_t xmaxx=xL3[2]*fAx2+fCx2;
 403   Double_t zCmin,zCmax,yCmin,yCmax;
 404
 405   zCmin=fZmin;
 406   zCmax=fZmax;
 407   yCmin=fYmin;
 408   yCmax=fYmax;
 409 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)))
 410     {
 411      if(fMap==3)
 412       {
 413         if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
 414        {
 415        //---------------------  Polar part ----------------------
 416
 417        Double_t yyp,ph0;
 418        Int_t  kp0, lp0, mp0;
 419        Int_t kpi,lpi,mpi;
 420        Double_t alp1,alp2,alp3;
 421        Double_t zpz,rp,fip,cphi;
 422
 423        kpi=kiip;
 424        lpi=kliip;
 425        mpi=kmiip;
 426
 427        zpz=xL3[2];
 428        kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
 429        zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
 430        zp1= kone-zp2;
 431        yyp=xL3[1]- 0.3;
 432        cphi=TMath::Abs(yyp/r0);
 433        Int_t kcphi=0;
 434        if (cphi > kone) {
 435         AliDebug(2,Form("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, cphi %e",xL3[0],xL3[1],xL3[2],yyp,r0,cphi));
 436         cphi =kone;
 437         kcphi=777;
 438        }
 439        ph0=TMath::ACos(cphi);
 440        if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
 441        if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
 442        if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
 443        if (ph0 > kPI2) {       ph0=ph0 - kPI2;}
 444        if (kcphi==777) {
 445         AliDebug(2,Form("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, ph0 %e",xL3[0],xL3[1],xL3[2],yyp,r0,ph0));
 446        }
 447        fip=ph0;
 448        mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
 449        xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
 450        xp1=kone-xp2;
 451
 452        Double_t rDel;
 453        rDel=fRdel;
 454
 455        if (r0<= fRdel)
 456         {
 457
 458          if(r0< keps)
 459          {
 460
 461           bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
 462           bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
 463           bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
 464
 465          } else {
 466
 467         alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
 468         alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
 469         alp1= kone - alp2 - alp3;
 470
 471         for (jb=0; jb<3 ; jb++)
 472          {
 473           iKvar=jb;
 474           bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
 475          }
 476         } // end   of keps <=r0
 477        }
 478         else
 479        {
 480         rp=r0;
 481
 482         lp0=FZ(rp,fR ,fRdel,lpi,fRn);
 483         yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
 484         yp1=kone-yp2;
 485
 486         for (jb=0; jb<3 ; jb++)
 487         {
 488          iKvar=jb;
 489          bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
 490         }
 491        }
 492
 493        b[0]=-bint[0];
 494        b[1]=bint[1];
 495        b[2]=-bint[2];
 496
 497    }
 498    else
 499    {
 500    //-------------- Cartensian part ------------------
 501
 502    Double_t zzc,yyc;
 503    Int_t k0, l0,m0;
 504    Double_t  xx1, xx2,dx, xxx ,xXl;
 505    Int_t kci,mci,lci;
 506
 507    kci=kiic;
 508    lci=kliic;
 509    mci=kmiic;
 510
 511    xx1 = fAx1*xL3[2] + fCx1;
 512    xx2 = fAx2*xL3[2] + fCx2;
 513
 514    zzc=xL3[2];
 515    k0=FZ(zzc, fZc ,fZdel, kci, fZl);
 516    zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
 517    zz1=kone - zz2;;
 518
 519    yyc=xL3[1];
 520    l0=FZ(yyc, fY , fYdel,lci,fYl);
 521    yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
 522    yy1=kone - yy2;
 523    xXl = fXl-kone;
 524    dx = (xx2-xx1)/xXl;
 525    xxx= xL3[0]-xx1;
 526    m0 = int(xxx/dx);
 527
 528    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
 529     {
 530       m0=m0+1;
 531       AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
 532     }
 533
 534    x2=(xL3[0]-( xx1+m0*dx))/dx;
 535    x1=kone-x2;
 536    m0=m0-1;
 537    for (jb=3; jb<6; jb++)
 538      {
 539        iKvar=jb;
 540        bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
 541      }
 542
 543    b[0]=-bint[0];
 544    b[1]=bint[1];
 545    b[2]=-bint[2];
 546
 547    }
 548
 549
 550   } else {
 551         AliError(Form("Unknown map of Dipole region %d",fMap));
 552  }
 553
 554 } else {
 555     ZDCField(xfi,b);
 556
 557     }
 558
 559   if(fFactor!=1) {
 560     b[0]*=fFactor;
 561     b[1]*=fFactor;
 562     b[2]*=fFactor;
 563   }
 564 }
 565
 566
 567 //_______________________________________________________________________
 568 Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar,
 569                     Float_t delu, Int_t ik, Int_t nk) const
 570 {
 571   //
 572   // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
 573   //
 574   Int_t l,ikj;
 575   Double_t ddu,ar;
 576
 577   Int_t ku,kf;
 578   kf=0;
 579   ar=Ar[ik];
 580   ddu=temp-ar;
 581
 582   ku=int(ddu/delu);
 583   ikj=ik+ku+1;
 584   if (ddu<=0) ikj=0;
 585
 586    for(l=ikj; l<nk; l++)
 587    {
 588
 589     if(temp <= Ar[l])
 590      {
 591
 592       kf=l-1;
 593       break;
 594      }
 595     }
 596     return kf;
 597   }
 598
 599 //_______________________________________________________________________
 600 Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2,
 601                        Double_t alf1, Double_t alf2, Double_t alf3,
 602                        Int_t kaa, Int_t maa) const
 603 {
 604   //
 605   // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
 606   //
 607   Double_t fa11=0,fa12=0,fa13=0;
 608   Double_t fa21=0,fa22=0,fa23=0;
 609   Double_t faY1=0,faY2=0;
 610   Double_t bba=0;
 611
 612   switch (kaai) {
 613   case 0:
 614     fa11 = fBpx[kaa][0][0];
 615     fa12 = fBpx[kaa][0][maa];
 616     fa13 = fBpx[kaa][0][maa+1];
 617     fa21 = fBpx[kaa+1][0][0];
 618     fa22 = fBpx[kaa+1][0][maa];
 619     fa23 = fBpx[kaa+1][0][maa+1];
 620     break;
 621   case 1:
 622     fa11 = fBpy[kaa][0][0];
 623     fa12 = fBpy[kaa][0][maa];
 624     fa13 = fBpy[kaa][0][maa+1];
 625     fa21 = fBpy[kaa+1][0][0];
 626     fa22 = fBpy[kaa+1][0][maa];
 627     fa23 = fBpy[kaa+1][0][maa+1];
 628     break;
 629   case 2:
 630     fa11 = fBpz[kaa][0][0];
 631     fa12 = fBpz[kaa][0][maa];
 632     fa13 = fBpz[kaa][0][maa+1];
 633     fa21 = fBpz[kaa+1][0][0];
 634     fa22 = fBpz[kaa+1][0][maa];
 635     fa23 = fBpz[kaa+1][0][maa+1];
 636     break;
 637   default:
 638     AliFatal(Form("Invalid value of kaai %d",kaai));
 639   }
 640   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
 641   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
 642   bba =  zaa1*faY1+zaa2*faY2;
 643   return bba;
 644 }
 645
 646
 647 //_______________________________________________________________________
 648 Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2,
 649                        Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m) const
 650 {
 651   //
 652   // Calculation of field componet at a given axis (general case)
 653   //
 654   Double_t fy1=0,fy2=0,ffy=0;
 655   Double_t gy1=0,gy2=0,ggy=0;
 656   Double_t bbi=0;
 657   Double_t bf11=0,bf12=0,bf21=0,bf22=0;
 658   Double_t bg11=0,bg12=0,bg21=0,bg22=0;
 659
 660
 661   /*-----------------Polar part ------------------*/
 662
 663   switch (kv) {
 664   case 0:
 665     bf11=fBpx[k][l][m];
 666     bf12=fBpx[k+1][l][m];
 667     bf21=fBpx[k+1][l+1][m];
 668     bf22=fBpx[k][l+1][m];
 669
 670     bg11=fBpx[k][l][m+1];
 671     bg12=fBpx[k+1][l][m+1];
 672     bg21=fBpx[k+1][l+1][m+1];
 673     bg22=fBpx[k][l+1][m+1];
 674     break;
 675
 676   case 1:
 677     bf11=fBpy[k][l][m];
 678     bf12=fBpy[k+1][l][m];
 679     bf21=fBpy[k+1][l+1][m];
 680     bf22=fBpy[k][l+1][m];
 681
 682     bg11=fBpy[k][l][m+1];
 683     bg12=fBpy[k+1][l][m+1];
 684     bg21=fBpy[k+1][l+1][m+1];
 685     bg22=fBpy[k][l+1][m+1];
 686     break;
 687
 688   case 2:
 689     bf11=fBpz[k][l][m];
 690     bf12=fBpz[k+1][l][m];
 691     bf21=fBpz[k+1][l+1][m];
 692     bf22=fBpz[k][l+1][m];
 693
 694     bg11=fBpz[k][l][m+1];
 695     bg12=fBpz[k+1][l][m+1];
 696     bg21=fBpz[k+1][l+1][m+1];
 697     bg22=fBpz[k][l+1][m+1];
 698     break;
 699   /*-----------------Cartensian part ---------------*/
 700
 701   case 3:
 702     bf11=fBcx[k][l][m];
 703     bf12=fBcx[k+1][l][m];
 704     bf21=fBcx[k+1][l+1][m];
 705     bf22=fBcx[k][l+1][m];
 706
 707     bg11=fBcx[k][l][m+1];
 708     bg12=fBcx[k+1][l][m+1];
 709     bg21=fBcx[k+1][l+1][m+1];
 710     bg22=fBcx[k][l+1][m+1];
 711     break;
 712
 713   case 4:
 714     bf11=fBcy[k][l][m];
 715     bf12=fBcy[k+1][l][m];
 716     bf21=fBcy[k+1][l+1][m];
 717     bf22=fBcy[k][l+1][m];
 718
 719     bg11=fBcy[k][l][m+1];
 720     bg12=fBcy[k+1][l][m+1];
 721     bg21=fBcy[k+1][l+1][m+1];
 722     bg22=fBcy[k][l+1][m+1];
 723     break;
 724
 725   case 5:
 726     bf11=fBcz[k][l][m];
 727     bf12=fBcz[k+1][l][m];
 728     bf21=fBcz[k+1][l+1][m];
 729     bf22=fBcz[k][l+1][m];
 730
 731     bg11=fBcz[k][l][m+1];
 732     bg12=fBcz[k+1][l][m+1];
 733     bg21=fBcz[k+1][l+1][m+1];
 734     bg22=fBcz[k][l+1][m+1];
 735     break;
 736
 737   default:
 738     AliFatal(Form("Invalid value of kv %d",kv));
 739   }
 740
 741
 742
 743   fy1=z1*bf11+z2*bf12;
 744   fy2=z2*bf21+z1* bf22;
 745   ffy=y1*fy1+ y2*fy2;
 746
 747
 748   gy1 = z1*bg11+z2*bg12;
 749   gy2 = z2*bg21+z1*bg22;
 750   ggy= y1*gy1 +  y2*gy2;
 751
 752   bbi =  x1*ffy+x2*ggy;
 753
 754   return bbi;
 755
 756 }
 757
 758 //_______________________________________________________________________
 759 void AliMagFDM::ReadField()
 760 {
 761   //
 762   // Method to read the magnetic field from file
 763   //
 764   FILE *magfile;
 765
 766   Int_t ik, il, im;
 767   Float_t zzp, rr,phii;
 768   Float_t zz, yy, bx,by,bz,bb;
 769
 770   char *fname;
 771   AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
 772   fname = gSystem->ExpandPathName(fTitle.Data());
 773   magfile=fopen(fname,"r");
 774   delete [] fname;
 775
 776   AliDebug(2,"Cartensian part");
 777
 778   if (magfile) {
 779
 780 //  Cartensian part
 781
 782     fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
 783
 784     AliDebug(3,Form("fYl %d, fXl %d, fZl %d",fYl, fXl, fZl));
 785
 786     for (ik=0; ik<fZl; ik++)
 787     {
 788
 789       fscanf(magfile, " %e ", &zz);
 790       fZc[ik]=zz;
 791
 792     }
 793
 794     for (ik=0; ik<fYl; ik++)
 795     {
 796        fscanf(magfile, " %e ", &yy);
 797        fY[ik]=yy;
 798
 799     }
 800     for (ik=0; ik<81; ik++)
 801     {
 802            AliDebug(4,Form("fZc %e,fY %e", fZc[ik],fY[ik]));
 803     }
 804
 805     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
 806
 807 AliDebug(3,Form("fYdel %e, fXdel %e, fZdel %e",fYdel,fXdel,fZdel));
 808 AliDebug(3,Form("fZmax %e, fZmin %e, fYmax %e,fYmin %e",fZmax,fZmin,fYmax,fYmin));
 809 AliDebug(3,Form("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e",fAx1,fCx1,fAx2,fCx2));
 810
 811     for (il=0; il<44; il++)  {
 812      for (im=0; im<81; im++)  {
 813       for (ik=0; ik<81; ik++)  {
 814
 815       fscanf(magfile, " %e ", &by);
 816       fBcy[ik][im][il]=by;
 817       }
 818      }
 819     }
 820
 821     for (il=0; il<44; il++)  {
 822      for (im=0; im<81; im++)  {
 823       for (ik=0; ik<81; ik++)  {
 824
 825       fscanf(magfile, " %e ", &bx);
 826       fBcx[ik][im][il]=bx;
 827       }
 828      }
 829     }
 830
 831    for (il=0; il<44; il++)  {
 832      for (im=0; im<81; im++)  {
 833       for (ik=0; ik<81; ik++)  {
 834
 835       fscanf(magfile, " %e ", &bz);
 836       fBcz[ik][im][il]=bz;
 837       }
 838      }
 839     }
 840 //----------------------   Polar part ---------------------------------
 841
 842     AliDebug(2,"Polar part");
 843     fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
 844     AliDebug(3,Form("fZpl %d, fRn %d, fPhin %d",fZpl,fRn,fPhin));
 845
 846     AliDebug(4," fZp array");
 847
 848     for (ik=0; ik<51; ik++)
 849     {
 850      fscanf(magfile, " %e ", &zzp);
 851      fZp[ik]=zzp;
 852      AliDebug(4,Form(" %e",fZp[ik]));
 853     }
 854
 855     AliDebug(4," fR array");
 856
 857     for (ik=0; ik<10; ik++)
 858     {
 859      fscanf(magfile, " %e ", &rr);
 860      fR[ik]=rr;
 861      AliDebug(4,Form(" %e",fR[ik]));
 862     }
 863
 864 //    AliDebug(4,"fPhi array");
 865
 866      for (il=0; il<33; il++)
 867      {
 868        fscanf(magfile, " %e ", &phii);
 869        fPhi[il]=phii;
 870 //        AliDebug(4,Form(" %e",fPhi[il]));
 871      }
 872
 873     fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
 874
 875 AliDebug(3,Form("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin));
 876
 877
 878     for (il=0; il<33; il++)  {
 879      for (im=0; im<10; im++)  {
 880       for (ik=0; ik<51; ik++)  {
 881       fscanf(magfile, " %e ", &by);
 882         fBpy[ik][im][il]=by;
 883       }
 884      }
 885     }
 886
 887     for (il=0; il<33; il++)  {
 888      for (im=0; im<10; im++)  {
 889       for (ik=0; ik<51; ik++)  {
 890       fscanf(magfile, " %e ", &bx);
 891         fBpx[ik][im][il]=bx;
 892       }
 893      }
 894     }
 895
 896
 897     for (il=0; il<33; il++)  {
 898      for (im=0; im<10; im++)  {
 899       for (ik=0; ik<51; ik++)  {
 900       fscanf(magfile, " %e ", &bz);
 901         fBpz[ik][im][il]=bz;
 902       }
 903      }
 904     }
 905
 906
 907     for (il=0; il<32; il++) {
 908      for (im=0; im<2; im++)  {
 909       for (ik=0; ik<2; ik++)  {
 910       fscanf(magfile, " %e ", &bb);
 911         fB[ik][im][il]=bb;
 912       }
 913      }
 914     }
 915 //
 916   } else {
 917     AliFatal(Form("File %s not found !",fTitle.Data()));
 918   }
 919 }