Introduce new class AliMagFDM - Galina Chabratova
[u/mrichter/AliRoot.git] / 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