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