]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliMagFDM.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added the HCAL section
[u/mrichter/AliRoot.git] / STEER / AliMagFDM.cxx
CommitLineData
aee8290b 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$
e2afb3b6 18Revision 1.12 2001/05/28 14:10:35 morsch
19SetSolenoidField method to set the L3 field strength. 2 kG is default.
20
4cc8933f 21Revision 1.11 2001/02/08 13:18:00 hristov
22Print removed (J.Gosset)
23
454d82e0 24Revision 1.10 2001/01/18 13:21:30 morsch
25Take pi from TMath.
26
11a3be04 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
e678332c 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
d8408e76 39Revision 1.7 2000/12/01 11:20:27 alibrary
40Corrector dipole removed from ZDC
41
006cb30c 42Revision 1.6 2000/11/10 18:09:55 fca
43New field map for the ZDC
44
d95ea23f 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
803c9752 50Revision 1.4 2000/10/02 21:28:14 fca
51Removal of useless dependecies via forward declarations
52
94de3818 53Revision 1.3 2000/07/13 16:19:09 fca
54Mainly coding conventions + some small bug fixes
55
ef42d733 56Revision 1.2 2000/07/12 08:56:25 fca
57Coding convention correction and warning removal
58
8918e700 59Revision 1.1 2000/07/11 18:24:59 fca
60Coding convention corrections + few minor bug fixes
61
aee8290b 62*/
63
8918e700 64#include <stdlib.h>
65
aee8290b 66#include "TSystem.h"
e2afb3b6 67
68#include "AliMagFDM.h"
aee8290b 69
70
71ClassImp(AliMagFDM)
72
e2afb3b6 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//_______________________________________________________________________
aee8290b 108AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
e2afb3b6 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)
aee8290b 137{
ef42d733 138 //
139 // Standard constructor for the Dipole field
140 //
aee8290b 141 fType = kDipoMap;
d8408e76 142 fMap = 3;
4cc8933f 143 SetSolenoidField();
d8408e76 144
e2afb3b6 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());
e678332c 148
aee8290b 149}
150
e2afb3b6 151//_______________________________________________________________________
aee8290b 152void AliMagFDM::Field(Float_t *xfi, Float_t *b)
153{
154 //
155 // Main routine to compute the field in a point
156 //
e2afb3b6 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
aee8290b 171
172 Double_t x[3];
173 Double_t xL3[3];
174 Double_t bint[3];
175
176 Double_t r0;
e2afb3b6 177 Int_t iKvar,jb;
aee8290b 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];
e678332c 187 b[0]=b[1]=b[2]=0;
aee8290b 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;
803c9752 193 rPmax=fRmax;
194 if ( (-700<x[2] && x[2]<=kfZbg &&
aee8290b 195 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
803c9752 196 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
aee8290b 197 {
803c9752 198 b[0]=b[1]=0;
4cc8933f 199 b[2]=fSolenoid;
aee8290b 200 }
201
202 xL3[0]=x[0]/100;
203 xL3[1]=(x[1]+30)/100;
204 xL3[2]=x[2]/100;
803c9752 205
aee8290b 206
803c9752 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;
aee8290b 212 Double_t zCmin,zCmax,yCmin,yCmax;
213
803c9752 214 zCmin=fZmin;
215 zCmax=fZmax;
216 yCmin=fYmin;
217 yCmax=fYmax;
803c9752 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)))
aee8290b 219 {
220 if(fMap==3)
221 {
803c9752 222 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
aee8290b 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];
e678332c 237 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
238 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
239 zp1= kone-zp2;
aee8290b 240 yyp=xL3[1]- 0.3;
e678332c 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 }
aee8290b 248 ph0=TMath::ACos(cphi);
e678332c 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;
aee8290b 260
261 Double_t rDel;
803c9752 262 rDel=fRdel;
aee8290b 263
803c9752 264 if (r0<= fRdel)
aee8290b 265 {
266
267 if(r0< keps)
268 {
269
803c9752 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;
aee8290b 273
803c9752 274 } else {
aee8290b 275
803c9752 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];
aee8290b 278 alp1= kone - alp2 - alp3;
279
280 for (jb=0; jb<3 ; jb++)
803c9752 281 {
282 iKvar=jb;
e678332c 283 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
803c9752 284 }
285 } // end of keps <=r0
aee8290b 286 }
287 else
288 {
289 rp=r0;
290
e678332c 291 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
292 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
293 yp1=kone-yp2;
aee8290b 294
295 for (jb=0; jb<3 ; jb++)
296 {
297 iKvar=jb;
e678332c 298 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
aee8290b 299 }
300 }
301
302 b[0]=bint[0];
303 b[1]=bint[1];
304 b[2]=bint[2];
305
aee8290b 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
803c9752 320 xx1 = fAx1*xL3[2] + fCx1;
321 xx2 = fAx2*xL3[2] + fCx2;
aee8290b 322
323 zzc=xL3[2];
e678332c 324 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
325 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
326 zz1=kone - zz2;;
aee8290b 327
328 yyc=xL3[1];
e678332c 329 l0=FZ(yyc, fY , fYdel,lci,fYl);
330 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
331 yy1=kone - yy2;
803c9752 332 xXl = fXl-kone;
aee8290b 333 dx = (xx2-xx1)/xXl;
334 xxx= xL3[0]-xx1;
aee8290b 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;
e678332c 349 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
aee8290b 350 }
351
352 b[0]=bint[0];
353 b[1]=bint[1];
354 b[2]=bint[2];
355
aee8290b 356 }
e678332c 357
aee8290b 358
359 } else {
360 printf("Unknown map of Dipole region %d\n",fMap);
361 }
362
363} else {
364
365//This is the ZDC part
d95ea23f 366 Float_t rad2=x[0]*x[0]+x[1]*x[1];
006cb30c 367 if(x[2]>kCORBEG2 && x[2]<kCOREND2){
d95ea23f 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];
aee8290b 382 }
383 }
d95ea23f 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
aee8290b 409
410 if(fFactor!=1) {
411 b[0]*=fFactor;
412 b[1]*=fFactor;
413 b[2]*=fFactor;
414 }
415}
416
e2afb3b6 417//_______________________________________________________________________
418Int_t AliMagFDM::FZ(Double_t temp, Float_t *Ar,
419 Float_t delu, Int_t ik,Int_t nk)
aee8290b 420{
421 //
e678332c 422 // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
aee8290b 423 //
e678332c 424 Int_t l,ikj;
425 Double_t ddu,ar;
aee8290b 426
e678332c 427 Int_t ku,kf;
428 kf=0;
aee8290b 429 ar=Ar[ik];
430 ddu=temp-ar;
431
432 ku=int(ddu/delu);
e678332c 433 ikj=ik+ku+1;
aee8290b 434 if (ddu<=0) ikj=0;
435
436 for(l=ikj; l<nk; l++)
437 {
e678332c 438
803c9752 439 if(temp <= Ar[l])
aee8290b 440 {
803c9752 441
e678332c 442 kf=l-1;
443 break;
aee8290b 444 }
445 }
e678332c 446 return kf;
aee8290b 447 }
448
e2afb3b6 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)
aee8290b 453{
454 //
e678332c 455 // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
aee8290b 456 //
457 Double_t fa11,fa12,fa13;
458 Double_t fa21,fa22,fa23;
459 Double_t faY1,faY2;
460 Double_t bba;
e678332c 461
8918e700 462 switch (kaai) {
463 case 0:
803c9752 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];
8918e700 470 break;
471 case 1:
803c9752 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];
8918e700 478 break;
479 case 2:
803c9752 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];
8918e700 486 break;
487 default:
e678332c 488 Fatal("Ba","Invalid value of kaai %d\n",kaai);
8918e700 489 exit(1);
aee8290b 490 }
491 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
492 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
493 bba = zaa1*faY1+zaa2*faY2;
e678332c 494 return bba;
aee8290b 495}
496
497
e2afb3b6 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)
aee8290b 501{
502 //
e678332c 503 // Calculation of field componet at a given axis (general case)
aee8290b 504 //
505 Double_t fy1, fy2, ffy;
506 Double_t gy1,gy2,ggy;
aee8290b 507 Double_t bbi;
aee8290b 508 Double_t bf11,bf12,bf21,bf22;
509 Double_t bg11,bg12,bg21,bg22;
e678332c 510
aee8290b 511
512 /*-----------------Polar part ------------------*/
513
8918e700 514 switch (kv) {
515 case 0:
803c9752 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];
aee8290b 520
803c9752 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];
8918e700 525 break;
526
527 case 1:
803c9752 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];
aee8290b 532
803c9752 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];
8918e700 537 break;
538
539 case 2:
803c9752 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];
aee8290b 544
803c9752 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];
8918e700 549 break;
aee8290b 550 /*-----------------Cartensian part ---------------*/
551
8918e700 552 case 3:
803c9752 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];
aee8290b 557
803c9752 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];
8918e700 562 break;
563
564 case 4:
803c9752 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];
aee8290b 569
803c9752 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];
8918e700 574 break;
575
576 case 5:
803c9752 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];
aee8290b 581
803c9752 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];
8918e700 586 break;
587
588 default:
e678332c 589 Fatal("Bb","Invalid value of kv %d\n",kv);
8918e700 590 exit(1);
aee8290b 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;
e678332c 605
606 return bbi;
aee8290b 607
608}
aee8290b 609
e2afb3b6 610//_______________________________________________________________________
aee8290b 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
803c9752 634 fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
aee8290b 635
803c9752 636 printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
aee8290b 637
803c9752 638 for (ik=0; ik<fZl; ik++)
aee8290b 639 {
640
641 fscanf(magfile, " %e ", &zz);
803c9752 642 fZc[ik]=zz;
aee8290b 643
644 }
645
803c9752 646 for (ik=0; ik<fYl; ik++)
aee8290b 647 {
648 fscanf(magfile, " %e ", &yy);
803c9752 649 fY[ik]=yy;
aee8290b 650
651 }
652 for (ik=0; ik<81; ik++)
653 {
803c9752 654 printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
aee8290b 655 }
656
803c9752 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);
aee8290b 658
803c9752 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);
aee8290b 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);
803c9752 668 fBcy[ik][im][il]=by;
aee8290b 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);
803c9752 678 fBcx[ik][im][il]=bx;
aee8290b 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);
803c9752 688 fBcz[ik][im][il]=bz;
aee8290b 689 }
690 }
691 }
692//---------------------- Polar part ---------------------------------
693
694 printf("Polar part\n");
803c9752 695 fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
696 printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
aee8290b 697
803c9752 698 printf(" fZp array\n");
aee8290b 699
700 for (ik=0; ik<51; ik++)
701 {
702 fscanf(magfile, " %e ", &zzp);
803c9752 703 fZp[ik]=zzp;
704 printf(" %e\n",fZp[ik]);
aee8290b 705 }
706
803c9752 707 printf(" fR array\n");
aee8290b 708
709 for (ik=0; ik<10; ik++)
710 {
711 fscanf(magfile, " %e ", &rr);
803c9752 712 fR[ik]=rr;
713 printf(" %e\n",fR[ik]);
aee8290b 714 }
715
803c9752 716// printf("fPhi array\n");
aee8290b 717
718 for (il=0; il<33; il++)
719 {
720 fscanf(magfile, " %e ", &phii);
803c9752 721 fPhi[il]=phii;
722// printf(" %e\n",fPhi[il]);
aee8290b 723 }
724
803c9752 725 fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
aee8290b 726
803c9752 727printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
aee8290b 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);
803c9752 734 fBpy[ik][im][il]=by;
aee8290b 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);
803c9752 743 fBpx[ik][im][il]=bx;
aee8290b 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);
803c9752 753 fBpz[ik][im][il]=bz;
aee8290b 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);
803c9752 763 fB[ik][im][il]=bb;
aee8290b 764 }
765 }
766 }
767//
768 } else {
769 printf("File %s not found !\n",fTitle.Data());
770 exit(1);
771 }
772}