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