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