42e50002eab653eaecf451fa10b8a4b44b393cbb
[u/mrichter/AliRoot.git] / 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 /*
17 $Log$
18 */
19
20 #include "AliMagFDM.h"
21 #include "TSystem.h"
22     
23
24 ClassImp(AliMagFDM)
25
26 //________________________________________
27 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
28 const Int_t map, const Float_t factor, const Float_t fmax)
29   : AliMagF(name,title,integ,map,factor,fmax)
30   
31 {
32   fType = kDipoMap;
33
34   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());
35   
36 }
37
38 //________________________________________
39
40 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
41 {
42   //
43   // Main routine to compute the field in a point
44   //
45   static  const Double_t keps=0.1E-06;
46   static  const Double_t kpi2=.6283185E+01;
47   static  const Double_t kone=1;
48
49   static  const    Int_t  kiip=33; 
50   static  const    Int_t  kmiip=0;    
51   static  const    Int_t  kliip=0;
52
53   static  const    Int_t  kiic=0;
54   static  const    Int_t  kmiic=0;    
55   static  const    Int_t  kliic=0;       
56
57   static  const Double_t    kfdZbg=502.92;  // Start of Map using in z
58   static  const Double_t    kfdZL3=600;  // Beginning of L3 door in z
59
60   Double_t   x[3];   
61   Double_t   xL3[3]; 
62   Double_t   bint[3]; 
63   
64   Double_t r0;
65   
66   Double_t bbj;
67   Int_t iKvar,jb;
68
69   Double_t zp1, zp2,xp1,xp2,yp1,yp2; 
70   Double_t zz1, zz2,yy1,yy2,x2,x1; 
71      
72 // --- start the map fiel from z = 502.92 cm ---
73
74   x[0] = xfi[0];
75   x[1] = xfi[1];
76   x[2] = xfi[2];
77   b[0]=b[1]=b[2]=0;
78   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]); 
79
80   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
81            r0=rr/100;
82   Double_t rPmax;
83   rPmax=fdRmax; 
84   if ( (-700<x[2] && x[2]<=kfdZbg && 
85         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
86        || (kfdZbg<x[2] && x[2]<=kfdZL3 && rr>=rPmax*100) )
87        {
88         b[2]=2;
89        }
90
91   xL3[0]=x[0]/100;
92   xL3[1]=(x[1]+30)/100;
93   xL3[2]=x[2]/100; 
94  
95   Double_t xminn=xL3[2]*fdAx1+fdCx1;
96   Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
97   Double_t zCmin,zCmax,yCmin,yCmax;
98
99   zCmin=fdZmin;
100   zCmax=fdZmax;
101   yCmin=fdYmin;
102   yCmax=fdYmax;
103           
104 if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2] <= zCmax ) && (yCmin<=xL3[1] && xL3[1]<= yCmax) &&  (xminn <= xL3[0] && xL3[0] <= xmaxx)))
105     {
106      if(fMap==3) 
107       { 
108        if (xL3[2]<zCmin && r0<rPmax)   
109        {
110        //---------------------  Polar part ----------------------
111        
112        Double_t yyp,ph0;
113        Int_t  kp0, lp0, mp0;
114        Int_t kpi,lpi,mpi;
115        Double_t alp1,alp2,alp3;
116        Double_t zpz,rp,fip,cphi; 
117
118        kpi=kiip; 
119        lpi=kliip;
120        mpi=kmiip;
121    
122        zpz=xL3[2];
123
124        FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&zp1 ,&zp2,&fdZpl) ;
125
126        yyp=xL3[1]- 0.3;
127        cphi=yyp/r0;
128        ph0=TMath::ACos(cphi);
129        if (xL3[0]< 0) {ph0=kpi2 - ph0;}  
130                                   
131        fip=ph0;
132        FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &xp1,&xp2,&fdPhin); 
133
134        Double_t rDel;
135        rDel=fdRdel;
136
137        if (r0<= fdRdel)
138         {
139
140          if(r0< keps) 
141          {
142
143           bint[0]=(zp1*fdBpx[kp0][0][0] + zp2*fdBpx[kp0+1][0][0])*10;     
144           bint[1]=(zp1*fdBpy[kp0][0][0] + zp2*fdBpy[kp0+1][0][0])*10;      
145           bint[2]=(zp1*fdBpz[kp0][0][0] + zp2*fdBpz[kp0+1][0][0])*10; 
146
147          }  
148       
149         alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0]; 
150         alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];      
151         alp1= kone - alp2 - alp3;
152      
153         for (jb=0; jb<3 ; jb++) 
154         {
155          iKvar=jb;     
156          FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
157          bint[jb] = bbj*10 ;
158         }  
159        }
160         else
161        {
162         rp=r0;
163
164         FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&yp1,&yp2,&fdRn);
165
166         for (jb=0; jb<3 ; jb++) 
167         {
168          iKvar=jb;
169          FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj); 
170
171          bint[jb] = bbj*10 ;
172         }
173        }
174
175        b[0]=bint[0];
176        b[1]=bint[1];
177        b[2]=bint[2];
178
179 //    fprintf(fitest,"-------------   Freg2 run -------------\n");       
180   
181    }  
182    else 
183    {
184    //-------------- Cartensian part ------------------
185           
186    Double_t zzc,yyc; 
187    Int_t k0, l0,m0;
188    Double_t  xx1, xx2,dx, xxx ,xXl; 
189    Int_t kci,mci,lci;
190
191    kci=kiic;
192    lci=kliic;
193    mci=kmiic;
194
195    xx1 = fdAx1*xL3[2] + fdCx1;
196    xx2 = fdAx2*xL3[2] + fdCx2;
197
198    zzc=xL3[2];
199    FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &zz1, &zz2, &fdZl);     
200
201    yyc=xL3[1];
202    FZ(&yyc, fdY , &fdYdel,&lci, &l0, &yy1, &yy2,&fdYl);  
203
204    xXl = fdXl-kone; 
205    dx = (xx2-xx1)/xXl;
206    xxx= xL3[0]-xx1;
207    //     xm = xxx/dx; 
208    m0 = int(xxx/dx);   
209
210    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx)) 
211     {
212       m0=m0+1;   
213       printf(" m0 %d, m0+1 %d\n",m0,m0+1);  
214     }
215
216    x2=(xL3[0]-( xx1+m0*dx))/dx;
217    x1=kone-x2;
218    m0=m0-1;
219    for (jb=3; jb<6; jb++) 
220      {
221        iKvar=jb;     
222        FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj); 
223        bint[jb-3] = bbj*10 ; 
224      }    
225  
226    b[0]=bint[0];
227    b[1]=bint[1];
228    b[2]=bint[2]; 
229
230 //   fprintf(fitest,"------------   Freg1 run -----------------\n");            
231    } 
232
233   } else {
234         printf("Unknown map of Dipole region %d\n",fMap);
235  }
236            
237 } else {
238
239 //This is the ZDC part
240       Float_t rad2=x[0]*x[0]+x[1]*x[1];
241       if(rad2<kD2RA2) {
242         if(x[2]>kD2BEG) {
243           
244 //    Separator Dipole D2
245           if(x[2]<kD2END) b[1]=kFDIP;
246         } else if(x[2]>kD1BEG) {
247
248 //    Separator Dipole D1
249           if(x[2]<kD1END) b[1]=-kFDIP;
250         }
251         if(rad2<kCORRA2) {
252
253 //    First quadrupole of inner triplet de-focussing in x-direction
254 //    Inner triplet
255           if(x[2]>kZ4BEG) {
256             if(x[2]<kZ4END) {
257               
258 //    2430 <-> 3060
259               b[0]=-kG1*x[1];
260               b[1]=-kG1*x[0];
261             }
262           } else if(x[2]>kZ3BEG) {
263             if(x[2]<kZ3END) {
264
265 //    1530 <-> 2080
266               b[0]=kG1*x[1];
267               b[1]=kG1*x[0];
268             }
269           } else if(x[2]>kZ2BEG) {
270             if(x[2]<kZ2END) {
271
272 //    890 <-> 1430
273               b[0]=kG1*x[1];
274               b[1]=kG1*x[0];
275             }
276           } else if(x[2]>kZ1BEG) {
277             if(x[2]<kZ1END) {
278
279 //    0 <->  630
280               b[0]=-kG1*x[1];
281               b[1]=-kG1*x[0];
282             }
283           } else if(x[2]>kCORBEG) {
284             if(x[2]<kCOREND) {
285 //    Corrector dipole (because of dimuon arm)
286 //            b[0]=kFCORN;
287               b[0]=-kFCORN;
288             }
289           }
290         }
291       }
292     }
293
294   if(fFactor!=1) {
295     b[0]*=fFactor;
296     b[1]*=fFactor;
297     b[2]*=fFactor;
298   }
299 }
300
301 //_________________________________________
302
303 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)
304   
305 {
306   //
307   // Z component of the field
308   //
309   static  const Double_t kone=1;
310   Int_t l,ik,ikj;
311   Double_t temp;
312   Double_t ddu,delu,ar;
313
314   Int_t nk,ku;
315   temp=*u;
316   nk=*nu;
317   ik=*ki;
318   delu=*du; 
319
320   ar=Ar[ik];
321   ddu=temp-ar;
322
323   ku=int(ddu/delu);
324   ikj=ik+ku;
325   if (ddu<=0) ikj=0;
326
327    for(l=ikj; l<nk; l++)
328    {
329
330     if(temp < Ar[l])
331      {
332       *kf=l;
333       *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
334       *a1= kone - *a2;
335       break;     
336      }
337     }
338   }
339
340 /*-------------FRfuncBi----------------*/
341
342 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)
343 {
344   //
345   // This method needs to be commented
346   //
347   Double_t fa11,fa12,fa13;
348   Double_t fa21,fa22,fa23;
349   Double_t faY1,faY2;
350   Double_t bba;
351   
352   Double_t zaa1,zaa2,alf1,alf2,alf3;
353   Int_t kaai,kaa,maa;
354   kaai=*kai;
355   kaa=*ka;
356   maa=*ma;
357   zaa1=*za1;
358   zaa2=*za2;
359   alf1=*al1;  
360   alf2=*al2;
361   alf3=*al3; 
362   
363   if (kaai==0 ) {
364     fa11 = fdBpx[kaa][0][0];
365     fa12 = fdBpx[kaa][0][maa];
366     fa13 = fdBpx[kaa][0][maa+1];
367     fa21 = fdBpx[kaa+1][0][0];
368     fa22 = fdBpx[kaa+1][0][maa];               
369     fa23 = fdBpx[kaa+1][0][maa+1]; 
370   }
371   if (kaai==1 ) {
372     fa11 = fdBpy[kaa][0][0];
373     fa12 = fdBpy[kaa][0][maa];
374     fa13 = fdBpy[kaa][0][maa+1];
375     fa21 = fdBpy[kaa+1][0][0];
376     fa22 = fdBpy[kaa+1][0][maa];               
377     fa23 = fdBpy[kaa+1][0][maa+1]; 
378   }
379   if (kaai==2 ) {
380     fa11 = fdBpz[kaa][0][0];
381     fa12 = fdBpz[kaa][0][maa];
382     fa13 = fdBpz[kaa][0][maa+1];
383     fa21 = fdBpz[kaa+1][0][0];
384     fa22 = fdBpz[kaa+1][0][maa];               
385     fa23 = fdBpz[kaa+1][0][maa+1]; 
386   }                            
387   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
388   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
389   bba =  zaa1*faY1+zaa2*faY2;    
390   *ba=bba;
391   
392 }
393
394   
395 /*----------- FGfuncBi------------*/
396
397 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)
398 {  
399   //
400   // This method needs to be commented
401   //
402   Double_t fy1, fy2, ffy;
403   Double_t gy1,gy2,ggy;
404   Double_t z1,z2,y1,y2,x1,x2;
405   
406   Int_t k,l,m,kv;
407   Double_t bbi;
408   
409   Double_t bf11,bf12,bf21,bf22;
410   Double_t bg11,bg12,bg21,bg22;
411   k=*kk;
412   l=*ll;
413   m=*mm;
414   
415   kv=*kvr;
416   
417   z1=*zz1;
418   z2=*zz2;
419   y1=*yy1;
420   y2=*yy2;
421   x1=*xx1;
422   x2=*xx2; 
423   
424   /*-----------------Polar part ------------------*/
425   
426   if(kv==0) {
427     bf11=fdBpx[k][l][m];
428     bf12=fdBpx[k+1][l][m];
429     bf21=fdBpx[k+1][l+1][m];
430     bf22=fdBpx[k][l+1][m];
431     
432     bg11=fdBpx[k][l][m+1];
433     bg12=fdBpx[k+1][l][m+1];
434     bg21=fdBpx[k+1][l+1][m+1];
435     bg22=fdBpx[k][l+1][m+1];
436   }
437   if(kv==1) {
438     bf11=fdBpy[k][l][m];
439     bf12=fdBpy[k+1][l][m];
440     bf21=fdBpy[k+1][l+1][m];
441     bf22=fdBpy[k][l+1][m];
442     
443     bg11=fdBpy[k][l][m+1];
444     bg12=fdBpy[k+1][l][m+1];
445     bg21=fdBpy[k+1][l+1][m+1];
446     bg22=fdBpy[k][l+1][m+1];
447   }  
448   
449   if(kv==2) {
450     bf11=fdBpz[k][l][m];
451     bf12=fdBpz[k+1][l][m];
452     bf21=fdBpz[k+1][l+1][m];
453     bf22=fdBpz[k][l+1][m];
454     
455     bg11=fdBpz[k][l][m+1];
456     bg12=fdBpz[k+1][l][m+1];
457     bg21=fdBpz[k+1][l+1][m+1];
458     bg22=fdBpz[k][l+1][m+1];           
459   } 
460   /*-----------------Cartensian part ---------------*/ 
461   
462   if(kv==3) {
463     bf11=fdBcx[k][l][m];
464     bf12=fdBcx[k+1][l][m];
465     bf21=fdBcx[k+1][l+1][m];
466     bf22=fdBcx[k][l+1][m];
467     
468     bg11=fdBcx[k][l][m+1];
469     bg12=fdBcx[k+1][l][m+1];
470     bg21=fdBcx[k+1][l+1][m+1];
471     bg22=fdBcx[k][l+1][m+1];
472   }          
473   
474   if(kv==4) {
475     bf11=fdBcy[k][l][m];
476     bf12=fdBcy[k+1][l][m];
477     bf21=fdBcy[k+1][l+1][m];
478     bf22=fdBcy[k][l+1][m];
479     
480     bg11=fdBcy[k][l][m+1];
481     bg12=fdBcy[k+1][l][m+1];
482     bg21=fdBcy[k+1][l+1][m+1];
483     bg22=fdBcy[k][l+1][m+1];
484   }          
485   if(kv==5) {
486     bf11=fdBcz[k][l][m];
487     bf12=fdBcz[k+1][l][m];
488     bf21=fdBcz[k+1][l+1][m];
489     bf22=fdBcz[k][l+1][m];
490     
491     bg11=fdBcz[k][l][m+1];
492     bg12=fdBcz[k+1][l][m+1];
493     bg21=fdBcz[k+1][l+1][m+1];
494     bg22=fdBcz[k][l+1][m+1];
495   }  
496   
497   
498   
499   fy1=z1*bf11+z2*bf12;
500   fy2=z2*bf21+z1* bf22;
501   ffy=y1*fy1+ y2*fy2; 
502   
503   
504   gy1 = z1*bg11+z2*bg12;
505   gy2 = z2*bg21+z1*bg22;
506   ggy= y1*gy1 +  y2*gy2;
507   
508   bbi =  x1*ffy+x2*ggy;
509   
510   *bb=bbi; 
511   
512 }     
513 //____________________________________________
514
515 void AliMagFDM::ReadField()
516 {
517   //
518   // Method to read the magnetic field from file
519   //
520   FILE *magfile;
521
522   Int_t ik, il, im;
523   Float_t zzp, rr,phii;
524   Float_t zz, yy, bx,by,bz,bb;
525
526   char *fname;
527   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
528   fname = gSystem->ExpandPathName(fTitle.Data());
529   magfile=fopen(fname,"r");
530   delete [] fname;
531
532   printf("Cartensian part\n");
533  
534   if (magfile) {
535   
536 //  Cartensian part 
537  
538     fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl); 
539     
540     printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);     
541     
542     for (ik=0; ik<fdZl; ik++)
543     { 
544    
545       fscanf(magfile, " %e ", &zz);
546       fdZc[ik]=zz; 
547
548     } 
549    
550     for (ik=0; ik<fdYl; ik++)
551     {    
552        fscanf(magfile, " %e ", &yy); 
553        fdY[ik]=yy;
554  
555     } 
556     for (ik=0; ik<81; ik++)
557     {    
558            printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]); 
559     }   
560              
561     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2); 
562
563 printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
564 printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
565 printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
566
567     for (il=0; il<44; il++)  { 
568      for (im=0; im<81; im++)  {      
569       for (ik=0; ik<81; ik++)  {      
570       
571       fscanf(magfile, " %e ", &by); 
572       fdBcy[ik][im][il]=by;        
573       }
574      }
575     } 
576
577     for (il=0; il<44; il++)  { 
578      for (im=0; im<81; im++)  {      
579       for (ik=0; ik<81; ik++)  {      
580       
581       fscanf(magfile, " %e ", &bx); 
582       fdBcx[ik][im][il]=bx;        
583       }    
584      }     
585     }
586     
587    for (il=0; il<44; il++)  { 
588      for (im=0; im<81; im++)  {      
589       for (ik=0; ik<81; ik++)  {      
590       
591       fscanf(magfile, " %e ", &bz); 
592       fdBcz[ik][im][il]=bz;          
593       }              
594      }     
595     } 
596 //----------------------   Polar part ---------------------------------
597
598     printf("Polar part\n");
599     fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin); 
600     printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);   
601
602     printf(" fdZp array\n"); 
603      
604     for (ik=0; ik<51; ik++) 
605     {    
606      fscanf(magfile, " %e ", &zzp);
607      fdZp[ik]=zzp; 
608      printf(" %e\n",fdZp[ik]);      
609     } 
610   
611     printf(" fdR array\n"); 
612          
613     for (ik=0; ik<10; ik++) 
614     {    
615      fscanf(magfile, " %e ", &rr); 
616      fdR[ik]=rr;
617      printf(" %e\n",fdR[ik]);
618     } 
619     
620 //    printf("fdPhi array\n"); 
621      
622      for (il=0; il<33; il++)  
623      {
624        fscanf(magfile, " %e ", &phii); 
625        fdPhi[il]=phii; 
626 //        printf(" %e\n",fdPhi[il]);          
627      }
628
629     fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin); 
630
631 printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
632
633                       
634     for (il=0; il<33; il++)  { 
635      for (im=0; im<10; im++)  {      
636       for (ik=0; ik<51; ik++)  {
637       fscanf(magfile, " %e ", &by); 
638         fdBpy[ik][im][il]=by;        
639       }
640      }
641     } 
642     
643     for (il=0; il<33; il++)  { 
644      for (im=0; im<10; im++)  {      
645       for (ik=0; ik<51; ik++)  {
646       fscanf(magfile, " %e ", &bx); 
647         fdBpx[ik][im][il]=bx;                     
648       }    
649      }     
650     }      
651
652
653     for (il=0; il<33; il++)  { 
654      for (im=0; im<10; im++)  {      
655       for (ik=0; ik<51; ik++)  {
656       fscanf(magfile, " %e ", &bz); 
657         fdBpz[ik][im][il]=bz;                      
658       }              
659      }     
660     } 
661     
662     
663     for (il=0; il<32; il++) { 
664      for (im=0; im<2; im++)  {      
665       for (ik=0; ik<2; ik++)  {
666       fscanf(magfile, " %e ", &bb);    
667         fdB[ik][im][il]=bb;  
668       }              
669      } 
670     }
671 //
672   } else { 
673     printf("File %s not found !\n",fTitle.Data());
674     exit(1);
675   }
676 }
677 //________________________________