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