STEER/AliMagF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18
  19 Revision 1.4  2000/03/28 12:40:24  fca
  20 Introduce factor for magnetic field
  21
  22
  23 Revision 1.3  1999/09/29 09:24:29  fca
  24 Introduction of the Copyright and cvs Log
  25
  26 */
  27
  28
  29 #include "AliMagF.h"
  30 #include "TSystem.h"
  31
  32 #include <stdlib.h>
  33 #include <stdio.h>
  34
  35 //ZDC part -------------------------------------------------------------------
  36
  37   static const Float_t G1=20.03;
  38   static const Float_t FDIP=-37.34;
  39   static const Float_t FDIMU=6.;
  40   static const Float_t FCORN=11.72;
  41 //
  42 // ZBEG       Beginning of the inner triplet
  43 // D1BEG      Beginning of separator dipole 1
  44 // D2BEG      Beginning of separator dipole 2
  45 // CORBEG     Corrector dipole beginning (because of dimuon arm)
  46 //
  47   static const Float_t CORBEG=1920,COREND=CORBEG+190, CORRA2=4.5*4.5;
  48 //
  49   static const Float_t ZBEG=2300;
  50   static const Float_t Z1BEG=ZBEG+   0,Z1END=Z1BEG+630,Z1RA2=3.5*3.5;
  51   static const Float_t Z2BEG=ZBEG+ 880,Z2END=Z2BEG+550,Z2RA2=3.5*3.5;
  52   static const Float_t Z3BEG=ZBEG+1530,Z3END=Z3BEG+550,Z3RA2=3.5*3.5;
  53   static const Float_t Z4BEG=ZBEG+2430,Z4END=Z4BEG+630,Z4RA2=3.5*3.5;
  54   static const Float_t D1BEG=5843.5   ,D1END=D1BEG+945,D1RA2=4.5*4.5;
  55   static const Float_t D2BEG=12113.2  ,D2END=D2BEG+945,D2RA2=4.5*.5;
  56
  57 //ZDC part -------------------------------------------------------------------
  58
  59 ClassImp(AliMagF)
  60
  61 //________________________________________
  62 AliMagF::AliMagF(const char *name, const char *title, const Int_t integ, const Int_t map,
  63                  const Float_t factor, const Float_t fmax)
  64   : TNamed(name,title)
  65 {
  66   fMap = map;
  67   fType = Undef;
  68   fInteg = integ;
  69   fFactor = factor;
  70   fMax = fmax;
  71 }
  72
  73 //________________________________________
  74 void AliMagF::Field(Float_t*, Float_t *b)
  75 {
  76   printf("Undefined MagF Field called, returning 0\n");
  77   b[0]=b[1]=b[2]=0;
  78 }
  79
  80 ClassImp(AliMagFC)
  81
  82 //________________________________________
  83 AliMagFC::AliMagFC(const char *name, const char *title, const Int_t integ, const Int_t map,
  84                  const Float_t factor, const Float_t fmax)
  85   : AliMagF(name,title,integ,map,factor,fmax)
  86 {
  87   printf("Constant Field %s created: map= %d, factor= %f\n",fName.Data(),map,factor);
  88   fType = Const;
  89 }
  90
  91 //________________________________________
  92 void AliMagFC::Field(Float_t *x, Float_t *b)
  93 {
  94   b[0]=b[1]=b[2]=0;
  95   if(fMap==1) {
  96     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
  97       b[2]=2;
  98     } else {
  99       if ( 725 <= x[2] && x[2] <= 1225 ) {
 100         Float_t dz = TMath::Abs(975-x[2])*0.01;
 101         b[0]=(1-0.1*dz*dz)*7;
 102       }
 103       else {
 104 //This is the ZDC part
 105         Float_t rad2=x[0]*x[0]+x[1]*x[1];
 106         if(rad2<D2RA2) {
 107           if(x[2]>D2BEG) {
 108
 109 //    Separator Dipole D2
 110             if(x[2]<D2END) b[1]=FDIP;
 111           } else if(x[2]>D1BEG) {
 112
 113 //    Separator Dipole D1
 114             if(x[2]<D1END) b[1]=-FDIP;
 115           }
 116           if(rad2<CORRA2) {
 117
 118 //    First quadrupole of inner triplet de-focussing in x-direction
 119 //    Inner triplet
 120             if(x[2]>Z4BEG) {
 121               if(x[2]<Z4END) {
 122
 123 //    2430 <-> 3060
 124                 b[0]=-G1*x[1];
 125                 b[1]=-G1*x[0];
 126               }
 127             } else if(x[2]>Z3BEG) {
 128               if(x[2]<Z3END) {
 129
 130 //    1530 <-> 2080
 131                 b[0]=G1*x[1];
 132                 b[1]=G1*x[0];
 133               }
 134             } else if(x[2]>Z2BEG) {
 135               if(x[2]<Z2END) {
 136
 137 //    890 <-> 1430
 138                 b[0]=G1*x[1];
 139                 b[1]=G1*x[0];
 140               }
 141             } else if(x[2]>Z1BEG) {
 142               if(x[2]<Z1END) {
 143
 144 //    0 <->  630
 145                 b[0]=-G1*x[1];
 146                 b[1]=-G1*x[0];
 147               }
 148             } else if(x[2]>CORBEG) {
 149               if(x[2]<COREND) {
 150 //    Corrector dipole (because of dimuon arm)
 151                 b[0]=FCORN;
 152               }
 153             }
 154           }
 155         }
 156       }
 157     }
 158     if(fFactor!=1) {
 159       b[0]*=fFactor;
 160       b[1]*=fFactor;
 161       b[2]*=fFactor;
 162     }
 163   } else {
 164     printf("Invalid field map for constant field %d\n",fMap);
 165     exit(1);
 166   }
 167 }
 168
 169 ClassImp(AliMagFCM)
 170
 171 //________________________________________
 172 AliMagFCM::AliMagFCM(const char *name, const char *title, const Int_t integ, const Int_t map,
 173                  const Float_t factor, const Float_t fmax)
 174   : AliMagF(name,title,integ,map,factor,fmax)
 175 {
 176   fType = ConMesh;
 177   printf("Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",fName.Data(),map,factor,fTitle.Data());
 178 }
 179
 180 //________________________________________
 181 void AliMagFCM::Field(Float_t *x, Float_t *b)
 182 {
 183   Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
 184     bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
 185   const Double_t one=1;
 186   Int_t ix, iy, iz;
 187
 188   // --- find the position in the grid ---
 189
 190   b[0]=b[1]=b[2]=0;
 191   if(-700<x[2] && x[2]<fZbeg && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
 192     b[2]=2;
 193   } else  {
 194     Bool_t infield=(fZbeg<=x[2] && x[2]<fZbeg+fZdel*(fZn-1)
 195        &&  ( fXbeg <= TMath::Abs(x[0]) && TMath::Abs(x[0]) < fXbeg+fXdel*(fXn-1) )
 196        &&  ( fYbeg <= TMath::Abs(x[1]) && TMath::Abs(x[1]) < fYbeg+fYdel*(fYn-1) ));
 197       if(infield) {
 198       xl[0]=TMath::Abs(x[0])-fXbeg;
 199       xl[1]=TMath::Abs(x[1])-fYbeg;
 200       xl[2]=x[2]-fZbeg;
 201
 202     // --- start with x
 203
 204       hix=xl[0]*fXdeli;
 205       ratx=hix-int(hix);
 206       ix=int(hix);
 207
 208       hiy=xl[1]*fYdeli;
 209       raty=hiy-int(hiy);
 210       iy=int(hiy);
 211
 212       hiz=xl[2]*fZdeli;
 213       ratz=hiz-int(hiz);
 214       iz=int(hiz);
 215
 216       if(fMap==2) {
 217       // ... simple interpolation
 218         ratx1=one-ratx;
 219         raty1=one-raty;
 220         ratz1=one-ratz;
 221         bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
 222         bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
 223         blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
 224         blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
 225         bhz   = blyhz             *raty1+bhyhz             *raty;
 226         blz   = blylz             *raty1+bhylz             *raty;
 227         b[0]  = blz               *ratz1+bhz               *ratz;
 228         //
 229         bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
 230         bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
 231         blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
 232         blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
 233         bhz   = blyhz             *raty1+bhyhz             *raty;
 234         blz   = blylz             *raty1+bhylz             *raty;
 235         b[1]  = blz               *ratz1+bhz               *ratz;
 236         //
 237         bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
 238         bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
 239         blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
 240         blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
 241         bhz   = blyhz             *raty1+bhyhz             *raty;
 242         blz   = blylz             *raty1+bhylz             *raty;
 243         b[2]  = blz               *ratz1+bhz               *ratz;
 244       //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
 245       //ratx,raty,ratz,b[0],b[1],b[2]);
 246       //
 247     // ... use the dipole symmetry
 248         if (x[0]*x[1] < 0) b[1]=-b[1];
 249         if (x[0]<0) b[2]=-b[2];
 250       } else {
 251         printf("Invalid field map for constant mesh %d\n",fMap);
 252       }
 253     } else {
 254 //This is the ZDC part
 255       Float_t rad2=x[0]*x[0]+x[1]*x[1];
 256       if(rad2<D2RA2) {
 257         if(x[2]>D2BEG) {
 258
 259 //    Separator Dipole D2
 260           if(x[2]<D2END) b[1]=FDIP;
 261         } else if(x[2]>D1BEG) {
 262
 263 //    Separator Dipole D1
 264           if(x[2]<D1END) b[1]=-FDIP;
 265         }
 266         if(rad2<CORRA2) {
 267
 268 //    First quadrupole of inner triplet de-focussing in x-direction
 269 //    Inner triplet
 270           if(x[2]>Z4BEG) {
 271             if(x[2]<Z4END) {
 272
 273 //    2430 <-> 3060
 274               b[0]=-G1*x[1];
 275               b[1]=-G1*x[0];
 276             }
 277           } else if(x[2]>Z3BEG) {
 278             if(x[2]<Z3END) {
 279
 280 //    1530 <-> 2080
 281               b[0]=G1*x[1];
 282               b[1]=G1*x[0];
 283             }
 284           } else if(x[2]>Z2BEG) {
 285             if(x[2]<Z2END) {
 286
 287 //    890 <-> 1430
 288               b[0]=G1*x[1];
 289               b[1]=G1*x[0];
 290             }
 291           } else if(x[2]>Z1BEG) {
 292             if(x[2]<Z1END) {
 293
 294 //    0 <->  630
 295               b[0]=-G1*x[1];
 296               b[1]=-G1*x[0];
 297             }
 298           } else if(x[2]>CORBEG) {
 299             if(x[2]<COREND) {
 300 //    Corrector dipole (because of dimuon arm)
 301               b[0]=FCORN;
 302             }
 303           }
 304         }
 305       }
 306     }
 307   }
 308   if(fFactor!=1) {
 309     b[0]*=fFactor;
 310     b[1]*=fFactor;
 311     b[2]*=fFactor;
 312   }
 313 }
 314
 315 //________________________________________
 316 void AliMagFCM::ReadField()
 317 {
 318   FILE *magfile;
 319   Int_t ix, iy, iz, ipx, ipy, ipz;
 320   Float_t bx, by, bz;
 321   char *fname;
 322   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
 323   fname = gSystem->ExpandPathName(fTitle.Data());
 324   magfile=fopen(fname,"r");
 325   delete [] fname;
 326   if (magfile) {
 327     fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
 328            &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
 329     printf("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
 330            fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
 331     fXdeli=1./fXdel;
 332     fYdeli=1./fYdel;
 333     fZdeli=1./fZdel;
 334     fB = new TVector(3*fXn*fYn*fZn);
 335     for (iz=0; iz<fZn; iz++) {
 336       ipz=iz*3*(fXn*fYn);
 337       for (iy=0; iy<fYn; iy++) {
 338         ipy=ipz+iy*3*fXn;
 339         for (ix=0; ix<fXn; ix++) {
 340           ipx=ipy+ix*3;
 341           fscanf(magfile,"%f %f %f",&bz,&by,&bx);
 342           (*fB)(ipx+2)=bz;
 343           (*fB)(ipx+1)=by;
 344           (*fB)(ipx  )=bx;
 345         }
 346       }
 347     }
 348   } else {
 349     printf("File %s not found !\n",fTitle.Data());
 350     exit(1);
 351   }
 352 }
 353 // -------------------------------------------------------
 354
 355 ClassImp(AliMagFDM)
 356
 357 //________________________________________
 358 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
 359 const Int_t map, const Float_t factor, const Float_t fmax)
 360   : AliMagF(name,title,integ,map,factor,fmax)
 361
 362 {
 363   fType = DipoMap;
 364
 365   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());
 366
 367 }
 368
 369 //________________________________________
 370
 371 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
 372 {
 373   static  const Double_t eps=0.1E-06;
 374   static  const Double_t pi2=.6283185E+01;
 375   static  const Double_t one=1;
 376   static  const Double_t fdYaxi = 0.3;
 377
 378   static  const    Int_t  kiip=33;
 379   static  const    Int_t  miip=0;
 380   static  const    Int_t  liip=0;
 381
 382   static  const    Int_t  kiic=0;
 383   static  const    Int_t  miic=0;
 384   static  const    Int_t  liic=0;
 385
 386   static  const Double_t    fdZbg=502.92;  // Start of Map using in z
 387   static  const Double_t    fdZL3=600;  // Beginning of L3 door in z
 388
 389   Double_t   x[3];
 390   Double_t   xL3[3];
 391   Double_t   bint[3];
 392
 393   Double_t r0;
 394
 395   Double_t bbj;
 396   Int_t Kvar,jb;
 397
 398   Double_t Zp1, Zp2,Xp1,Xp2,Yp1,Yp2;
 399   Double_t Zz1, Zz2,Yy1,Yy2,X2,X1;
 400
 401 // --- start the map fiel from z = 502.92 cm ---
 402
 403   x[0] = xfi[0];
 404   x[1] = xfi[1];
 405   x[2] = xfi[2];
 406   b[0]=b[1]=b[2]=0;
 407   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]);
 408
 409   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 410            r0=rr/100;
 411   Double_t Rpmax;
 412   Rpmax=fdRmax;
 413   if ( (-700<x[2] && x[2]<=fdZbg &&
 414         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
 415        || (fdZbg<x[2] && x[2]<=fdZL3 && rr>=Rpmax*100) )
 416        {
 417         b[2]=2;
 418        }
 419
 420   xL3[0]=x[0]/100;
 421   xL3[1]=(x[1]+30)/100;
 422   xL3[2]=x[2]/100;
 423
 424   Double_t xminn=xL3[2]*fdAx1+fdCx1;
 425   Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
 426   Double_t Zcmin,Zcmax,Ycmin,Ycmax;
 427
 428   Zcmin=fdZmin;
 429   Zcmax=fdZmax;
 430   Ycmin=fdYmin;
 431   Ycmax=fdYmax;
 432
 433 if ((fdZbg/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)))
 434     {
 435      if(fMap==3)
 436       {
 437        if (xL3[2]<Zcmin && r0<Rpmax)
 438        {
 439        //---------------------  Polar part ----------------------
 440
 441        Double_t yyp,ph0;
 442        Int_t  kp0, lp0, mp0;
 443        Int_t kpi,lpi,mpi;
 444        Double_t alp1,alp2,alp3;
 445        Double_t zpz,rp,fip,cphi;
 446
 447        kpi=kiip;
 448        lpi=liip;
 449        mpi=miip;
 450
 451        zpz=xL3[2];
 452
 453        FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&Zp1 ,&Zp2,&fdZpl) ;
 454
 455        yyp=xL3[1]- 0.3;
 456        cphi=yyp/r0;
 457        ph0=TMath::ACos(cphi);
 458        if (xL3[0]< 0) {ph0=pi2 - ph0;}
 459
 460        fip=ph0;
 461        FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &Xp1,&Xp2,&fdPhin);
 462
 463        Double_t Rdel;
 464        Rdel=fdRdel;
 465
 466        if (r0<= fdRdel)
 467         {
 468
 469          if(r0< eps)
 470          {
 471
 472           bint[0]=(Zp1*fdBpx[kp0][0][0] + Zp2*fdBpx[kp0+1][0][0])*10;
 473           bint[1]=(Zp1*fdBpy[kp0][0][0] + Zp2*fdBpy[kp0+1][0][0])*10;
 474           bint[2]=(Zp1*fdBpz[kp0][0][0] + Zp2*fdBpz[kp0+1][0][0])*10;
 475
 476          }
 477
 478         alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0];
 479         alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];
 480         alp1= one - alp2 - alp3;
 481
 482         for (jb=0; jb<3 ; jb++)
 483         {
 484          Kvar=jb;
 485          FRfuncBi(&Kvar,&Zp1,&Zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
 486          bint[jb] = bbj*10 ;
 487         }
 488        }
 489         else
 490        {
 491         rp=r0;
 492
 493         FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&Yp1,&Yp2,&fdRn);
 494
 495         for (jb=0; jb<3 ; jb++)
 496         {
 497          Kvar=jb;
 498          FGfuncBi(&Zp1,&Zp2,&Yp1,&Yp2,&Xp1,&Xp2,&Kvar,&kp0,&lp0,&mp0,&bbj);
 499
 500          bint[jb] = bbj*10 ;
 501         }
 502        }
 503
 504        b[0]=bint[0];
 505        b[1]=bint[1];
 506        b[2]=bint[2];
 507
 508 //    fprintf(fitest,"-------------   Freg2 run -------------\n");
 509
 510    }
 511    else
 512    {
 513    //-------------- Cartensian part ------------------
 514
 515    Double_t zzc,yyc;
 516    Int_t k0, l0,m0;
 517    Double_t  xx1, xx2,dx, xxx ,xXl;
 518    Int_t kci,mci,lci;
 519
 520    kci=kiic;
 521    lci=liic;
 522    mci=miic;
 523
 524    xx1 = fdAx1*xL3[2] + fdCx1;
 525    xx2 = fdAx2*xL3[2] + fdCx2;
 526
 527    zzc=xL3[2];
 528    FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &Zz1, &Zz2, &fdZl);
 529
 530    yyc=xL3[1];
 531    FZ(&yyc, fdY , &fdYdel,&lci, &l0, &Yy1, &Yy2,&fdYl);
 532
 533    xXl = fdXl-one;
 534    dx = (xx2-xx1)/xXl;
 535    xxx= xL3[0]-xx1;
 536    //     xm = xxx/dx;
 537    m0 = int(xxx/dx);
 538
 539    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
 540     {
 541       m0=m0+1;
 542       printf(" m0 %d, m0+1 %d\n",m0,m0+1);
 543     }
 544
 545    X2=(xL3[0]-( xx1+m0*dx))/dx;
 546    X1=one-X2;
 547    m0=m0-1;
 548    for (jb=3; jb<6; jb++)
 549      {
 550        Kvar=jb;
 551        FGfuncBi(&Zz1,&Zz2,&Yy1,&Yy2,&X1,&X2,&Kvar,&k0, &l0, &m0, &bbj);
 552        bint[jb-3] = bbj*10 ;
 553      }
 554
 555    b[0]=bint[0];
 556    b[1]=bint[1];
 557    b[2]=bint[2];
 558
 559 //   fprintf(fitest,"------------   Freg1 run -----------------\n");
 560    }
 561
 562   } else {
 563         printf("Unknown map of Dipole region %d\n",fMap);
 564  }
 565
 566 } else {
 567
 568 //This is the ZDC part
 569       Float_t rad2=x[0]*x[0]+x[1]*x[1];
 570       if(rad2<D2RA2) {
 571         if(x[2]>D2BEG) {
 572
 573 //    Separator Dipole D2
 574           if(x[2]<D2END) b[1]=FDIP;
 575         } else if(x[2]>D1BEG) {
 576
 577 //    Separator Dipole D1
 578           if(x[2]<D1END) b[1]=-FDIP;
 579         }
 580         if(rad2<CORRA2) {
 581
 582 //    First quadrupole of inner triplet de-focussing in x-direction
 583 //    Inner triplet
 584           if(x[2]>Z4BEG) {
 585             if(x[2]<Z4END) {
 586
 587 //    2430 <-> 3060
 588               b[0]=-G1*x[1];
 589               b[1]=-G1*x[0];
 590             }
 591           } else if(x[2]>Z3BEG) {
 592             if(x[2]<Z3END) {
 593
 594 //    1530 <-> 2080
 595               b[0]=G1*x[1];
 596               b[1]=G1*x[0];
 597             }
 598           } else if(x[2]>Z2BEG) {
 599             if(x[2]<Z2END) {
 600
 601 //    890 <-> 1430
 602               b[0]=G1*x[1];
 603               b[1]=G1*x[0];
 604             }
 605           } else if(x[2]>Z1BEG) {
 606             if(x[2]<Z1END) {
 607
 608 //    0 <->  630
 609               b[0]=-G1*x[1];
 610               b[1]=-G1*x[0];
 611             }
 612           } else if(x[2]>CORBEG) {
 613             if(x[2]<COREND) {
 614 //    Corrector dipole (because of dimuon arm)
 615 //            b[0]=FCORN;
 616               b[0]=-FCORN;
 617             }
 618           }
 619         }
 620       }
 621     }
 622
 623   if(fFactor!=1) {
 624     b[0]*=fFactor;
 625     b[1]*=fFactor;
 626     b[2]*=fFactor;
 627   }
 628 }
 629
 630 //_________________________________________
 631
 632 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)
 633
 634  {
 635   static  const Double_t one=1;
 636   Int_t l,ik,ikj;
 637   Double_t temp;
 638   Double_t ddu,delu,ar;
 639
 640   Int_t nk,ku;
 641   temp=*u;
 642   nk=*nu;
 643   ik=*ki;
 644   delu=*du;
 645
 646   ar=Ar[ik];
 647   ddu=temp-ar;
 648
 649   ku=int(ddu/delu);
 650   ikj=ik+ku;
 651   if (ddu<=0) ikj=0;
 652
 653    for(l=ikj; l<nk; l++)
 654    {
 655
 656     if(temp < Ar[l])
 657      {
 658       *kf=l;
 659       *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
 660       *a1= one - *a2;
 661       break;
 662      }
 663     }
 664   }
 665
 666 /*-------------FRfuncBi----------------*/
 667
 668 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)
 669
 670 {
 671 Double_t fa11,fa12,fa13;
 672 Double_t fa21,fa22,fa23;
 673 Double_t faY1,faY2;
 674 Double_t bba;
 675
 676 Double_t zaa1,zaa2,alf1,alf2,alf3;
 677 Int_t kaai,kaa,maa;
 678 kaai=*kai;
 679 kaa=*ka;
 680 maa=*ma;
 681 zaa1=*za1;
 682 zaa2=*za2;
 683 alf1=*al1;
 684 alf2=*al2;
 685 alf3=*al3;
 686
 687  if (kaai==0 ) {
 688               fa11 = fdBpx[kaa][0][0];
 689               fa12 = fdBpx[kaa][0][maa];
 690               fa13 = fdBpx[kaa][0][maa+1];
 691               fa21 = fdBpx[kaa+1][0][0];
 692               fa22 = fdBpx[kaa+1][0][maa];
 693               fa23 = fdBpx[kaa+1][0][maa+1];
 694               }
 695  if (kaai==1 ) {
 696               fa11 = fdBpy[kaa][0][0];
 697               fa12 = fdBpy[kaa][0][maa];
 698               fa13 = fdBpy[kaa][0][maa+1];
 699               fa21 = fdBpy[kaa+1][0][0];
 700               fa22 = fdBpy[kaa+1][0][maa];
 701               fa23 = fdBpy[kaa+1][0][maa+1];
 702               }
 703  if (kaai==2 ) {
 704               fa11 = fdBpz[kaa][0][0];
 705               fa12 = fdBpz[kaa][0][maa];
 706               fa13 = fdBpz[kaa][0][maa+1];
 707               fa21 = fdBpz[kaa+1][0][0];
 708               fa22 = fdBpz[kaa+1][0][maa];
 709               fa23 = fdBpz[kaa+1][0][maa+1];
 710               }
 711     faY1=alf1*fa11+alf2*fa12+alf3*fa13;
 712     faY2=alf1*fa21+alf2*fa22+alf3*fa23;
 713     bba =  zaa1*faY1+zaa2*faY2;
 714     *ba=bba;
 715
 716 }
 717
 718
 719 /*----------- FGfuncBi------------*/
 720
 721 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)
 722
 723 {
 724 Double_t fy1, fy2, ffy;
 725 Double_t gy1,gy2,ggy;
 726 Double_t z1,z2,y1,y2,x1,x2;
 727
 728 Int_t k,l,m,kv;
 729 Double_t bbi;
 730
 731 Double_t bf11,bf12,bf21,bf22;
 732 Double_t bg11,bg12,bg21,bg22;
 733 k=*kk;
 734 l=*ll;
 735 m=*mm;
 736
 737 kv=*kvr;
 738
 739 z1=*zz1;
 740 z2=*zz2;
 741 y1=*yy1;
 742 y2=*yy2;
 743 x1=*xx1;
 744 x2=*xx2;
 745
 746 /*-----------------Polar part ------------------*/
 747
 748 if(kv==0) {
 749            bf11=fdBpx[k][l][m];
 750            bf12=fdBpx[k+1][l][m];
 751            bf21=fdBpx[k+1][l+1][m];
 752            bf22=fdBpx[k][l+1][m];
 753
 754            bg11=fdBpx[k][l][m+1];
 755            bg12=fdBpx[k+1][l][m+1];
 756            bg21=fdBpx[k+1][l+1][m+1];
 757            bg22=fdBpx[k][l+1][m+1];
 758            }
 759  if(kv==1) {
 760            bf11=fdBpy[k][l][m];
 761            bf12=fdBpy[k+1][l][m];
 762            bf21=fdBpy[k+1][l+1][m];
 763            bf22=fdBpy[k][l+1][m];
 764
 765            bg11=fdBpy[k][l][m+1];
 766            bg12=fdBpy[k+1][l][m+1];
 767            bg21=fdBpy[k+1][l+1][m+1];
 768            bg22=fdBpy[k][l+1][m+1];
 769            }
 770
 771  if(kv==2) {
 772            bf11=fdBpz[k][l][m];
 773            bf12=fdBpz[k+1][l][m];
 774            bf21=fdBpz[k+1][l+1][m];
 775            bf22=fdBpz[k][l+1][m];
 776
 777            bg11=fdBpz[k][l][m+1];
 778            bg12=fdBpz[k+1][l][m+1];
 779            bg21=fdBpz[k+1][l+1][m+1];
 780            bg22=fdBpz[k][l+1][m+1];
 781            }
 782 /*-----------------Cartensian part ---------------*/
 783
 784  if(kv==3) {
 785            bf11=fdBcx[k][l][m];
 786            bf12=fdBcx[k+1][l][m];
 787            bf21=fdBcx[k+1][l+1][m];
 788            bf22=fdBcx[k][l+1][m];
 789
 790            bg11=fdBcx[k][l][m+1];
 791            bg12=fdBcx[k+1][l][m+1];
 792            bg21=fdBcx[k+1][l+1][m+1];
 793            bg22=fdBcx[k][l+1][m+1];
 794            }
 795
 796  if(kv==4) {
 797            bf11=fdBcy[k][l][m];
 798            bf12=fdBcy[k+1][l][m];
 799            bf21=fdBcy[k+1][l+1][m];
 800            bf22=fdBcy[k][l+1][m];
 801
 802            bg11=fdBcy[k][l][m+1];
 803            bg12=fdBcy[k+1][l][m+1];
 804            bg21=fdBcy[k+1][l+1][m+1];
 805            bg22=fdBcy[k][l+1][m+1];
 806            }
 807  if(kv==5) {
 808            bf11=fdBcz[k][l][m];
 809            bf12=fdBcz[k+1][l][m];
 810            bf21=fdBcz[k+1][l+1][m];
 811            bf22=fdBcz[k][l+1][m];
 812
 813            bg11=fdBcz[k][l][m+1];
 814            bg12=fdBcz[k+1][l][m+1];
 815            bg21=fdBcz[k+1][l+1][m+1];
 816            bg22=fdBcz[k][l+1][m+1];
 817            }
 818
 819
 820
 821      fy1=z1*bf11+z2*bf12;
 822      fy2=z2*bf21+z1* bf22;
 823      ffy=y1*fy1+ y2*fy2;
 824
 825
 826      gy1 = z1*bg11+z2*bg12;
 827      gy2 = z2*bg21+z1*bg22;
 828      ggy= y1*gy1 +  y2*gy2;
 829
 830      bbi =  x1*ffy+x2*ggy;
 831
 832      *bb=bbi;
 833
 834 }
 835 //____________________________________________
 836
 837 void AliMagFDM::ReadField()
 838 {
 839   FILE *magfile;
 840
 841   Int_t ik, il, im;
 842   Float_t zzp, rr,phii;
 843   Float_t zz, yy, bx,by,bz,bb;
 844
 845   char *fname;
 846   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
 847   fname = gSystem->ExpandPathName(fTitle.Data());
 848   magfile=fopen(fname,"r");
 849   delete [] fname;
 850
 851   printf("Cartensian part\n");
 852
 853   if (magfile) {
 854
 855 //  Cartensian part
 856
 857     fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl);
 858
 859     printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);
 860
 861     for (ik=0; ik<fdZl; ik++)
 862     {
 863
 864       fscanf(magfile, " %e ", &zz);
 865       fdZc[ik]=zz;
 866
 867     }
 868
 869     for (ik=0; ik<fdYl; ik++)
 870     {
 871        fscanf(magfile, " %e ", &yy);
 872        fdY[ik]=yy;
 873
 874     }
 875     for (ik=0; ik<81; ik++)
 876     {
 877            printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]);
 878     }
 879
 880     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2);
 881
 882 printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
 883 printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
 884 printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
 885
 886     for (il=0; il<44; il++)  {
 887      for (im=0; im<81; im++)  {
 888       for (ik=0; ik<81; ik++)  {
 889
 890       fscanf(magfile, " %e ", &by);
 891       fdBcy[ik][im][il]=by;
 892       }
 893      }
 894     }
 895
 896     for (il=0; il<44; il++)  {
 897      for (im=0; im<81; im++)  {
 898       for (ik=0; ik<81; ik++)  {
 899
 900       fscanf(magfile, " %e ", &bx);
 901       fdBcx[ik][im][il]=bx;
 902       }
 903      }
 904     }
 905
 906    for (il=0; il<44; il++)  {
 907      for (im=0; im<81; im++)  {
 908       for (ik=0; ik<81; ik++)  {
 909
 910       fscanf(magfile, " %e ", &bz);
 911       fdBcz[ik][im][il]=bz;
 912       }
 913      }
 914     }
 915 //----------------------   Polar part ---------------------------------
 916
 917     printf("Polar part\n");
 918     fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin);
 919     printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);
 920
 921     printf(" fdZp array\n");
 922
 923     for (ik=0; ik<51; ik++)
 924     {
 925      fscanf(magfile, " %e ", &zzp);
 926      fdZp[ik]=zzp;
 927      printf(" %e\n",fdZp[ik]);
 928     }
 929
 930     printf(" fdR array\n");
 931
 932     for (ik=0; ik<10; ik++)
 933     {
 934      fscanf(magfile, " %e ", &rr);
 935      fdR[ik]=rr;
 936      printf(" %e\n",fdR[ik]);
 937     }
 938
 939 //    printf("fdPhi array\n");
 940
 941      for (il=0; il<33; il++)
 942      {
 943        fscanf(magfile, " %e ", &phii);
 944        fdPhi[il]=phii;
 945 //        printf(" %e\n",fdPhi[il]);
 946      }
 947
 948     fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin);
 949
 950 printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
 951
 952
 953     for (il=0; il<33; il++)  {
 954      for (im=0; im<10; im++)  {
 955       for (ik=0; ik<51; ik++)  {
 956       fscanf(magfile, " %e ", &by);
 957         fdBpy[ik][im][il]=by;
 958       }
 959      }
 960     }
 961
 962     for (il=0; il<33; il++)  {
 963      for (im=0; im<10; im++)  {
 964       for (ik=0; ik<51; ik++)  {
 965       fscanf(magfile, " %e ", &bx);
 966         fdBpx[ik][im][il]=bx;
 967       }
 968      }
 969     }
 970
 971
 972     for (il=0; il<33; il++)  {
 973      for (im=0; im<10; im++)  {
 974       for (ik=0; ik<51; ik++)  {
 975       fscanf(magfile, " %e ", &bz);
 976         fdBpz[ik][im][il]=bz;
 977       }
 978      }
 979     }
 980
 981
 982     for (il=0; il<32; il++) {
 983      for (im=0; im<2; im++)  {
 984       for (ik=0; ik<2; ik++)  {
 985       fscanf(magfile, " %e ", &bb);
 986         fdB[ik][im][il]=bb;
 987       }
 988      }
 989     }
 990 //
 991   } else {
 992     printf("File %s not found !\n",fTitle.Data());
 993     exit(1);
 994   }
 995 }
 996 //________________________________
 997
 998
 999
1000
1001
1002
1003