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