ecee137df538bf64698ab2508b6ae6c3b8fd0a43
[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 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //   Field with Magnetic Field map
20 //   Used by AliRun class
21 //   Author:
22 //-------------------------------------------------------------------------
23
24 #include <stdlib.h>
25
26 #include "TSystem.h"
27
28 #include "AliMagFDM.h"
29     
30
31 ClassImp(AliMagFDM)
32
33 //_______________________________________________________________________
34 AliMagFDM::AliMagFDM():
35   fSolenoid(0),
36   fInd(0),
37   fZmin(0),
38   fZmax(0),
39   fYmax(0),
40   fYmin(0),
41   fZpmx(0),
42   fZpmn(0),
43   fRmax(0),
44   fRmin(0),
45   fXdel(0),
46   fYdel(0),
47   fZdel(0),
48   fRdel(0),
49   fPhid(0),
50   fZpdl(0),
51   fCx1(0),
52   fCx2(0),
53   fAx1(0),
54   fAx2(0),
55   fXl(0),
56   fYl(0),
57   fZl(0),    
58   fRn(0),
59   fPhin(0),
60   fZpl(0)  
61 {
62   //
63   // Default constructor for the Dipole field
64   //
65 }
66
67 //_______________________________________________________________________
68 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
69                      const Float_t factor, const Float_t fmax):
70   AliMagF(name,title,integ,factor,fmax),
71   fSolenoid(0),
72   fInd(0),
73   fZmin(0),
74   fZmax(0),
75   fYmax(0),
76   fYmin(0),
77   fZpmx(0),
78   fZpmn(0),
79   fRmax(0),
80   fRmin(0),
81   fXdel(0),
82   fYdel(0),
83   fZdel(0),
84   fRdel(0),
85   fPhid(0),
86   fZpdl(0),
87   fCx1(0),
88   fCx2(0),
89   fAx1(0),
90   fAx2(0),
91   fXl(0),
92   fYl(0),
93   fZl(0),    
94   fRn(0),
95   fPhin(0),
96   fZpl(0)  
97 {
98   //
99   // Standard constructor for the Dipole field
100   //
101   fType = kDipoMap;
102   fMap  = 3;
103   SetSolenoidField();
104   
105   Info("ctor",
106        "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s\n",
107        fName.Data(), fMap ,integ,factor,fTitle.Data());
108  
109 }
110
111 //_______________________________________________________________________
112 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
113 {
114   //
115   // Main routine to compute the field in a point
116   //
117   const Double_t keps=0.1E-06;
118   const Double_t kPI2=2.*TMath::Pi();
119   const Double_t kone=1;
120   
121   const    Int_t  kiip=33; 
122   const    Int_t  kmiip=0;    
123   const    Int_t  kliip=0;
124   
125   const    Int_t  kiic=0;
126   const    Int_t  kmiic=0;    
127   const    Int_t  kliic=0;       
128   
129   const Double_t    kfZbg=502.92;  // Start of Map using in z
130   const Double_t    kfZL3=600;  // Beginning of L3 door in z
131
132   Double_t   x[3];   
133   Double_t   xL3[3]; 
134   Double_t   bint[3]; 
135   
136   Double_t r0;
137   Int_t iKvar,jb;
138
139   Double_t zp1, zp2,xp1,xp2,yp1,yp2; 
140   Double_t zz1, zz2,yy1,yy2,x2,x1; 
141      
142 // --- start the map fiel from z = 502.92 cm ---
143
144   x[0] = xfi[0];
145   x[1] = xfi[1];
146   x[2] = xfi[2];
147   b[0]=b[1]=b[2]=0;
148   //       printf("x[0]  %f,x[1] %f,x[2]  %f\n",x[0],x[1],x[2]); 
149
150   Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
151            r0=rr/100;
152   Double_t rPmax;
153   rPmax=fRmax; 
154   if ( (-700<x[2] && x[2]<=kfZbg && 
155         (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
156        || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
157        {
158         b[0]=b[1]=0;
159         b[2]=fSolenoid;
160        }
161
162   xL3[0]=x[0]/100;
163   xL3[1]=(x[1]+30)/100;
164   xL3[2]=x[2]/100; 
165
166  
167   if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
168   if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
169
170   Double_t xminn=xL3[2]*fAx1+fCx1;
171   Double_t xmaxx=xL3[2]*fAx2+fCx2;
172   Double_t zCmin,zCmax,yCmin,yCmax;
173
174   zCmin=fZmin;
175   zCmax=fZmax;
176   yCmin=fYmin;
177   yCmax=fYmax;
178 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)))
179     {
180      if(fMap==3) 
181       { 
182         if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
183        {
184        //---------------------  Polar part ----------------------
185        
186        Double_t yyp,ph0;
187        Int_t  kp0, lp0, mp0;
188        Int_t kpi,lpi,mpi;
189        Double_t alp1,alp2,alp3;
190        Double_t zpz,rp,fip,cphi; 
191
192        kpi=kiip; 
193        lpi=kliip;
194        mpi=kmiip;
195    
196        zpz=xL3[2];
197        kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
198        zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
199        zp1= kone-zp2;
200        yyp=xL3[1]- 0.3;
201        cphi=TMath::Abs(yyp/r0);
202        Int_t kcphi=0;
203        if (cphi > kone) {
204         printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, cphi %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,cphi);
205         cphi =kone;
206         kcphi=777;
207        } 
208        ph0=TMath::ACos(cphi);
209        if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}  
210        if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;} 
211        if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}  
212        if (ph0 > kPI2) {       ph0=ph0 - kPI2;}
213        if (kcphi==777) {
214         printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, ph0 %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,ph0);
215        }  
216        fip=ph0; 
217        mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
218        xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
219        xp1=kone-xp2;
220
221        Double_t rDel;
222        rDel=fRdel;
223
224        if (r0<= fRdel)
225         {
226
227          if(r0< keps) 
228          {
229
230           bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;     
231           bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;      
232           bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10; 
233
234          } else { 
235       
236         alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0]; 
237         alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];      
238         alp1= kone - alp2 - alp3;
239      
240         for (jb=0; jb<3 ; jb++) 
241          {
242           iKvar=jb;     
243           bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
244          } 
245         } // end   of keps <=r0 
246        }
247         else
248        {
249         rp=r0;
250
251         lp0=FZ(rp,fR ,fRdel,lpi,fRn);
252         yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
253         yp1=kone-yp2;
254
255         for (jb=0; jb<3 ; jb++) 
256         {
257          iKvar=jb;
258          bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
259         }
260        }
261
262        b[0]=bint[0];
263        b[1]=bint[1];
264        b[2]=bint[2];
265
266    }  
267    else 
268    {
269    //-------------- Cartensian part ------------------
270           
271    Double_t zzc,yyc; 
272    Int_t k0, l0,m0;
273    Double_t  xx1, xx2,dx, xxx ,xXl; 
274    Int_t kci,mci,lci;
275
276    kci=kiic;
277    lci=kliic;
278    mci=kmiic;
279
280    xx1 = fAx1*xL3[2] + fCx1;
281    xx2 = fAx2*xL3[2] + fCx2;
282
283    zzc=xL3[2];
284    k0=FZ(zzc, fZc ,fZdel, kci, fZl); 
285    zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
286    zz1=kone - zz2;;    
287
288    yyc=xL3[1];
289    l0=FZ(yyc, fY , fYdel,lci,fYl);  
290    yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
291    yy1=kone - yy2;    
292    xXl = fXl-kone; 
293    dx = (xx2-xx1)/xXl;
294    xxx= xL3[0]-xx1;
295    m0 = int(xxx/dx);   
296
297    if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx)) 
298     {
299       m0=m0+1;   
300       printf(" m0 %d, m0+1 %d\n",m0,m0+1);  
301     }
302
303    x2=(xL3[0]-( xx1+m0*dx))/dx;
304    x1=kone-x2;
305    m0=m0-1;
306    for (jb=3; jb<6; jb++) 
307      {
308        iKvar=jb;     
309        bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ; 
310      }    
311  
312    b[0]=bint[0];
313    b[1]=bint[1];
314    b[2]=bint[2]; 
315
316    } 
317
318
319   } else {
320         printf("Unknown map of Dipole region %d\n",fMap);
321  }
322            
323 } else {
324
325 //This is the ZDC part
326     Float_t rad2=x[0]*x[0]+x[1]*x[1];
327     if(x[2]>kCORBEG2 && x[2]<kCOREND2){
328       if(rad2<kCOR2RA2){
329         b[0] = kFCORN2;
330       }
331     }
332     else if(x[2]>kZ1BEG && x[2]<kZ1END){  
333       if(rad2<kZ1RA2){
334         b[0] = -kG1*x[1];
335         b[1] = -kG1*x[0];
336       }
337     }
338     else if(x[2]>kZ2BEG && x[2]<kZ2END){  
339       if(rad2<kZ2RA2){
340         b[0] = kG1*x[1];
341         b[1] = kG1*x[0];
342       }
343     }
344     else if(x[2]>kZ3BEG && x[2]<kZ3END){  
345       if(rad2<kZ3RA2){
346         b[0] = kG1*x[1];
347         b[1] = kG1*x[0];
348       }
349     }
350     else if(x[2]>kZ4BEG && x[2]<kZ4END){  
351       if(rad2<kZ4RA2){
352         b[0] = -kG1*x[1];
353         b[1] = -kG1*x[0];
354       }
355     }
356     else if(x[2]>kD1BEG && x[2]<kD1END){ 
357       if(rad2<kD1RA2){
358         b[1] = -kFDIP;
359       }
360     }
361     else if(x[2]>kD2BEG && x[2]<kD2END){
362       if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
363         || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
364         b[1] = kFDIP;
365       }
366     }
367       }
368     
369
370   if(fFactor!=1) {
371     b[0]*=fFactor;
372     b[1]*=fFactor;
373     b[2]*=fFactor;
374   }
375 }
376
377 //_______________________________________________________________________
378 Int_t AliMagFDM::FZ(Double_t temp, Float_t *Ar, 
379                     Float_t delu, Int_t ik,Int_t nk) const
380 {
381   //
382   // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
383   //
384   Int_t l,ikj;
385   Double_t ddu,ar;
386
387   Int_t ku,kf;
388   kf=0;
389   ar=Ar[ik];
390   ddu=temp-ar;
391
392   ku=int(ddu/delu);
393   ikj=ik+ku+1;
394   if (ddu<=0) ikj=0;
395
396    for(l=ikj; l<nk; l++)
397    {
398
399     if(temp <= Ar[l])
400      {
401        
402       kf=l-1;
403       break;
404      }
405     }
406     return kf;
407   }
408
409 //_______________________________________________________________________
410 Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2, 
411                        Double_t alf1, Double_t alf2, Double_t alf3, 
412                        Int_t kaa, Int_t maa)
413 {
414   //
415   // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
416   //
417   Double_t fa11,fa12,fa13;
418   Double_t fa21,fa22,fa23;
419   Double_t faY1,faY2;
420   Double_t bba;
421
422   switch (kaai) {
423   case 0:
424     fa11 = fBpx[kaa][0][0];
425     fa12 = fBpx[kaa][0][maa];
426     fa13 = fBpx[kaa][0][maa+1];
427     fa21 = fBpx[kaa+1][0][0];
428     fa22 = fBpx[kaa+1][0][maa];               
429     fa23 = fBpx[kaa+1][0][maa+1]; 
430     break;
431   case 1:
432     fa11 = fBpy[kaa][0][0];
433     fa12 = fBpy[kaa][0][maa];
434     fa13 = fBpy[kaa][0][maa+1];
435     fa21 = fBpy[kaa+1][0][0];
436     fa22 = fBpy[kaa+1][0][maa];               
437     fa23 = fBpy[kaa+1][0][maa+1]; 
438     break;
439   case 2:
440     fa11 = fBpz[kaa][0][0];
441     fa12 = fBpz[kaa][0][maa];
442     fa13 = fBpz[kaa][0][maa+1];
443     fa21 = fBpz[kaa+1][0][0];
444     fa22 = fBpz[kaa+1][0][maa];               
445     fa23 = fBpz[kaa+1][0][maa+1]; 
446     break;
447   default:
448     Fatal("Ba","Invalid value of kaai %d\n",kaai);
449     exit(1);
450   }                            
451   faY1=alf1*fa11+alf2*fa12+alf3*fa13;
452   faY2=alf1*fa21+alf2*fa22+alf3*fa23;
453   bba =  zaa1*faY1+zaa2*faY2;    
454   return bba;  
455 }
456
457   
458 //_______________________________________________________________________
459 Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2, 
460                        Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m)
461 {  
462   //
463   // Calculation of field componet at a given axis (general case)
464   //
465   Double_t fy1, fy2, ffy;
466   Double_t gy1,gy2,ggy;
467   Double_t bbi;
468   Double_t bf11,bf12,bf21,bf22;
469   Double_t bg11,bg12,bg21,bg22;
470
471   
472   /*-----------------Polar part ------------------*/
473   
474   switch (kv) {
475   case 0:
476     bf11=fBpx[k][l][m];
477     bf12=fBpx[k+1][l][m];
478     bf21=fBpx[k+1][l+1][m];
479     bf22=fBpx[k][l+1][m];
480     
481     bg11=fBpx[k][l][m+1];
482     bg12=fBpx[k+1][l][m+1];
483     bg21=fBpx[k+1][l+1][m+1];
484     bg22=fBpx[k][l+1][m+1];
485     break;
486
487   case 1:
488     bf11=fBpy[k][l][m];
489     bf12=fBpy[k+1][l][m];
490     bf21=fBpy[k+1][l+1][m];
491     bf22=fBpy[k][l+1][m];
492     
493     bg11=fBpy[k][l][m+1];
494     bg12=fBpy[k+1][l][m+1];
495     bg21=fBpy[k+1][l+1][m+1];
496     bg22=fBpy[k][l+1][m+1];
497     break;
498
499   case 2:
500     bf11=fBpz[k][l][m];
501     bf12=fBpz[k+1][l][m];
502     bf21=fBpz[k+1][l+1][m];
503     bf22=fBpz[k][l+1][m];
504     
505     bg11=fBpz[k][l][m+1];
506     bg12=fBpz[k+1][l][m+1];
507     bg21=fBpz[k+1][l+1][m+1];
508     bg22=fBpz[k][l+1][m+1]; 
509     break;
510   /*-----------------Cartensian part ---------------*/ 
511   
512   case 3:
513     bf11=fBcx[k][l][m];
514     bf12=fBcx[k+1][l][m];
515     bf21=fBcx[k+1][l+1][m];
516     bf22=fBcx[k][l+1][m];
517     
518     bg11=fBcx[k][l][m+1];
519     bg12=fBcx[k+1][l][m+1];
520     bg21=fBcx[k+1][l+1][m+1];
521     bg22=fBcx[k][l+1][m+1];
522     break;
523
524   case 4:
525     bf11=fBcy[k][l][m];
526     bf12=fBcy[k+1][l][m];
527     bf21=fBcy[k+1][l+1][m];
528     bf22=fBcy[k][l+1][m];
529     
530     bg11=fBcy[k][l][m+1];
531     bg12=fBcy[k+1][l][m+1];
532     bg21=fBcy[k+1][l+1][m+1];
533     bg22=fBcy[k][l+1][m+1];
534     break;
535
536   case 5:
537     bf11=fBcz[k][l][m];
538     bf12=fBcz[k+1][l][m];
539     bf21=fBcz[k+1][l+1][m];
540     bf22=fBcz[k][l+1][m];
541     
542     bg11=fBcz[k][l][m+1];
543     bg12=fBcz[k+1][l][m+1];
544     bg21=fBcz[k+1][l+1][m+1];
545     bg22=fBcz[k][l+1][m+1];
546     break;
547
548   default:
549     Fatal("Bb","Invalid value of kv %d\n",kv);
550     exit(1);
551   }  
552   
553   
554   
555   fy1=z1*bf11+z2*bf12;
556   fy2=z2*bf21+z1* bf22;
557   ffy=y1*fy1+ y2*fy2; 
558   
559   
560   gy1 = z1*bg11+z2*bg12;
561   gy2 = z2*bg21+z1*bg22;
562   ggy= y1*gy1 +  y2*gy2;
563   
564   bbi =  x1*ffy+x2*ggy;
565
566   return bbi;
567   
568 }     
569
570 //_______________________________________________________________________
571 void AliMagFDM::ReadField()
572 {
573   //
574   // Method to read the magnetic field from file
575   //
576   FILE *magfile;
577
578   Int_t ik, il, im;
579   Float_t zzp, rr,phii;
580   Float_t zz, yy, bx,by,bz,bb;
581
582   char *fname;
583   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
584   fname = gSystem->ExpandPathName(fTitle.Data());
585   magfile=fopen(fname,"r");
586   delete [] fname;
587
588   printf("Cartensian part\n");
589  
590   if (magfile) {
591   
592 //  Cartensian part 
593  
594     fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl); 
595     
596     printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);     
597     
598     for (ik=0; ik<fZl; ik++)
599     { 
600    
601       fscanf(magfile, " %e ", &zz);
602       fZc[ik]=zz; 
603
604     } 
605    
606     for (ik=0; ik<fYl; ik++)
607     {    
608        fscanf(magfile, " %e ", &yy); 
609        fY[ik]=yy;
610  
611     } 
612     for (ik=0; ik<81; ik++)
613     {    
614            printf("fZc %e,fY %e\n", fZc[ik],fY[ik]); 
615     }   
616              
617     fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2); 
618
619 printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
620 printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
621 printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
622
623     for (il=0; il<44; il++)  { 
624      for (im=0; im<81; im++)  {      
625       for (ik=0; ik<81; ik++)  {      
626       
627       fscanf(magfile, " %e ", &by); 
628       fBcy[ik][im][il]=by;        
629       }
630      }
631     } 
632
633     for (il=0; il<44; il++)  { 
634      for (im=0; im<81; im++)  {      
635       for (ik=0; ik<81; ik++)  {      
636       
637       fscanf(magfile, " %e ", &bx); 
638       fBcx[ik][im][il]=bx;        
639       }    
640      }     
641     }
642     
643    for (il=0; il<44; il++)  { 
644      for (im=0; im<81; im++)  {      
645       for (ik=0; ik<81; ik++)  {      
646       
647       fscanf(magfile, " %e ", &bz); 
648       fBcz[ik][im][il]=bz;          
649       }              
650      }     
651     } 
652 //----------------------   Polar part ---------------------------------
653
654     printf("Polar part\n");
655     fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin); 
656     printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);   
657
658     printf(" fZp array\n"); 
659      
660     for (ik=0; ik<51; ik++) 
661     {    
662      fscanf(magfile, " %e ", &zzp);
663      fZp[ik]=zzp; 
664      printf(" %e\n",fZp[ik]);      
665     } 
666   
667     printf(" fR array\n"); 
668          
669     for (ik=0; ik<10; ik++) 
670     {    
671      fscanf(magfile, " %e ", &rr); 
672      fR[ik]=rr;
673      printf(" %e\n",fR[ik]);
674     } 
675     
676 //    printf("fPhi array\n"); 
677      
678      for (il=0; il<33; il++)  
679      {
680        fscanf(magfile, " %e ", &phii); 
681        fPhi[il]=phii; 
682 //        printf(" %e\n",fPhi[il]);          
683      }
684
685     fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin); 
686
687 printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
688
689                       
690     for (il=0; il<33; il++)  { 
691      for (im=0; im<10; im++)  {      
692       for (ik=0; ik<51; ik++)  {
693       fscanf(magfile, " %e ", &by); 
694         fBpy[ik][im][il]=by;        
695       }
696      }
697     } 
698     
699     for (il=0; il<33; il++)  { 
700      for (im=0; im<10; im++)  {      
701       for (ik=0; ik<51; ik++)  {
702       fscanf(magfile, " %e ", &bx); 
703         fBpx[ik][im][il]=bx;                     
704       }    
705      }     
706     }      
707
708
709     for (il=0; il<33; il++)  { 
710      for (im=0; im<10; im++)  {      
711       for (ik=0; ik<51; ik++)  {
712       fscanf(magfile, " %e ", &bz); 
713         fBpz[ik][im][il]=bz;                      
714       }              
715      }     
716     } 
717     
718     
719     for (il=0; il<32; il++) { 
720      for (im=0; im<2; im++)  {      
721       for (ik=0; ik<2; ik++)  {
722       fscanf(magfile, " %e ", &bb);    
723         fB[ik][im][il]=bb;  
724       }              
725      } 
726     }
727 //
728   } else { 
729     printf("File %s not found !\n",fTitle.Data());
730     exit(1);
731   }
732 }