git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Introduce new class AliMagFDM - Galina Chabratova
[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx
CommitLineData
4c039060 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$
7a15f6b8 18
19Revision 1.4 2000/03/28 12:40:24 fca
20Introduce factor for magnetic field
21
22
dd045843 23Revision 1.3 1999/09/29 09:24:29 fca
24Introduction of the Copyright and cvs Log
25
4c039060 26*/
27
fe4da5cc 28
29#include "AliMagF.h"
30#include "TSystem.h"
7a15f6b8 31
fe4da5cc 32#include <stdlib.h>
33#include <stdio.h>
34
35//ZDC part -------------------------------------------------------------------
36
37 static const Float_t G1=20.03;
38 static const Float_t FDIP=-37.34;
39 static const Float_t FDIMU=6.;
40 static const Float_t FCORN=11.72;
41//
42// ZBEG Beginning of the inner triplet
43// D1BEG Beginning of separator dipole 1
44// D2BEG Beginning of separator dipole 2
45// CORBEG Corrector dipole beginning (because of dimuon arm)
46//
47 static const Float_t CORBEG=1920,COREND=CORBEG+190, CORRA2=4.5*4.5;
48//
49 static const Float_t ZBEG=2300;
50 static const Float_t Z1BEG=ZBEG+ 0,Z1END=Z1BEG+630,Z1RA2=3.5*3.5;
51 static const Float_t Z2BEG=ZBEG+ 880,Z2END=Z2BEG+550,Z2RA2=3.5*3.5;
52 static const Float_t Z3BEG=ZBEG+1530,Z3END=Z3BEG+550,Z3RA2=3.5*3.5;
53 static const Float_t Z4BEG=ZBEG+2430,Z4END=Z4BEG+630,Z4RA2=3.5*3.5;
54 static const Float_t D1BEG=5843.5 ,D1END=D1BEG+945,D1RA2=4.5*4.5;
55 static const Float_t D2BEG=12113.2 ,D2END=D2BEG+945,D2RA2=4.5*.5;
56
57//ZDC part -------------------------------------------------------------------
58
59ClassImp(AliMagF)
60
61//________________________________________
62AliMagF::AliMagF(const char *name, const char *title, const Int_t integ, const Int_t map,
63 const Float_t factor, const Float_t fmax)
64 : TNamed(name,title)
65{
66 fMap = map;
67 fType = Undef;
68 fInteg = integ;
69 fFactor = factor;
70 fMax = fmax;
71}
72
73//________________________________________
74void AliMagF::Field(Float_t*, Float_t *b)
75{
76 printf("Undefined MagF Field called, returning 0\n");
77 b[0]=b[1]=b[2]=0;
78}
79
80ClassImp(AliMagFC)
81
82//________________________________________
83AliMagFC::AliMagFC(const char *name, const char *title, const Int_t integ, const Int_t map,
84 const Float_t factor, const Float_t fmax)
85 : AliMagF(name,title,integ,map,factor,fmax)
86{
87 printf("Constant Field %s created: map= %d, factor= %f\n",fName.Data(),map,factor);
88 fType = Const;
89}
90
91//________________________________________
92void AliMagFC::Field(Float_t *x, Float_t *b)
93{
94 b[0]=b[1]=b[2]=0;
95 if(fMap==1) {
96 if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
97 b[2]=2;
98 } else {
99 if ( 725 <= x[2] && x[2] <= 1225 ) {
100 Float_t dz = TMath::Abs(975-x[2])*0.01;
101 b[0]=(1-0.1*dz*dz)*7;
102 }
103 else {
104//This is the ZDC part
105 Float_t rad2=x[0]*x[0]+x[1]*x[1];
106 if(rad2<D2RA2) {
107 if(x[2]>D2BEG) {
108
109// Separator Dipole D2
110 if(x[2]<D2END) b[1]=FDIP;
111 } else if(x[2]>D1BEG) {
112
113// Separator Dipole D1
114 if(x[2]<D1END) b[1]=-FDIP;
115 }
116 if(rad2<CORRA2) {
117
118// First quadrupole of inner triplet de-focussing in x-direction
119// Inner triplet
120 if(x[2]>Z4BEG) {
121 if(x[2]<Z4END) {
122
123// 2430 <-> 3060
124 b[0]=-G1*x[1];
125 b[1]=-G1*x[0];
126 }
127 } else if(x[2]>Z3BEG) {
128 if(x[2]<Z3END) {
129
130// 1530 <-> 2080
131 b[0]=G1*x[1];
132 b[1]=G1*x[0];
133 }
134 } else if(x[2]>Z2BEG) {
135 if(x[2]<Z2END) {
136
137// 890 <-> 1430
138 b[0]=G1*x[1];
139 b[1]=G1*x[0];
140 }
141 } else if(x[2]>Z1BEG) {
142 if(x[2]<Z1END) {
143
144// 0 <-> 630
145 b[0]=-G1*x[1];
146 b[1]=-G1*x[0];
147 }
148 } else if(x[2]>CORBEG) {
149 if(x[2]<COREND) {
150// Corrector dipole (because of dimuon arm)
151 b[0]=FCORN;
152 }
153 }
154 }
155 }
156 }
157 }
dd045843 158 if(fFactor!=1) {
159 b[0]*=fFactor;
160 b[1]*=fFactor;
161 b[2]*=fFactor;
162 }
fe4da5cc 163 } else {
164 printf("Invalid field map for constant field %d\n",fMap);
165 exit(1);
166 }
167}
168
169ClassImp(AliMagFCM)
170
171//________________________________________
172AliMagFCM::AliMagFCM(const char *name, const char *title, const Int_t integ, const Int_t map,
173 const Float_t factor, const Float_t fmax)
174 : AliMagF(name,title,integ,map,factor,fmax)
175{
176 fType = ConMesh;
177 printf("Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",fName.Data(),map,factor,fTitle.Data());
178}
179
180//________________________________________
181void AliMagFCM::Field(Float_t *x, Float_t *b)
182{
183 Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
184 bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
185 const Double_t one=1;
186 Int_t ix, iy, iz;
187
188 // --- find the position in the grid ---
189
190 b[0]=b[1]=b[2]=0;
191 if(-700<x[2] && x[2]<fZbeg && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
192 b[2]=2;
193 } else {
11f6d36c 194 Bool_t infield=(fZbeg<=x[2] && x[2]<fZbeg+fZdel*(fZn-1)
fe4da5cc 195 && ( fXbeg <= TMath::Abs(x[0]) && TMath::Abs(x[0]) < fXbeg+fXdel*(fXn-1) )
11f6d36c 196 && ( fYbeg <= TMath::Abs(x[1]) && TMath::Abs(x[1]) < fYbeg+fYdel*(fYn-1) ));
197 if(infield) {
fe4da5cc 198 xl[0]=TMath::Abs(x[0])-fXbeg;
199 xl[1]=TMath::Abs(x[1])-fYbeg;
200 xl[2]=x[2]-fZbeg;
201
202 // --- start with x
203
204 hix=xl[0]*fXdeli;
205 ratx=hix-int(hix);
206 ix=int(hix);
207
208 hiy=xl[1]*fYdeli;
209 raty=hiy-int(hiy);
210 iy=int(hiy);
211
212 hiz=xl[2]*fZdeli;
213 ratz=hiz-int(hiz);
214 iz=int(hiz);
215
216 if(fMap==2) {
217 // ... simple interpolation
218 ratx1=one-ratx;
219 raty1=one-raty;
220 ratz1=one-ratz;
221 bhyhz = Bx(ix ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
222 bhylz = Bx(ix ,iy+1,iz )*ratx1+Bx(ix+1,iy+1,iz )*ratx;
223 blyhz = Bx(ix ,iy ,iz+1)*ratx1+Bx(ix+1,iy ,iz+1)*ratx;
224 blylz = Bx(ix ,iy ,iz )*ratx1+Bx(ix+1,iy ,iz )*ratx;
225 bhz = blyhz *raty1+bhyhz *raty;
226 blz = blylz *raty1+bhylz *raty;
227 b[0] = blz *ratz1+bhz *ratz;
228 //
229 bhyhz = By(ix ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
230 bhylz = By(ix ,iy+1,iz )*ratx1+By(ix+1,iy+1,iz )*ratx;
231 blyhz = By(ix ,iy ,iz+1)*ratx1+By(ix+1,iy ,iz+1)*ratx;
232 blylz = By(ix ,iy ,iz )*ratx1+By(ix+1,iy ,iz )*ratx;
233 bhz = blyhz *raty1+bhyhz *raty;
234 blz = blylz *raty1+bhylz *raty;
235 b[1] = blz *ratz1+bhz *ratz;
236 //
237 bhyhz = Bz(ix ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
238 bhylz = Bz(ix ,iy+1,iz )*ratx1+Bz(ix+1,iy+1,iz )*ratx;
239 blyhz = Bz(ix ,iy ,iz+1)*ratx1+Bz(ix+1,iy ,iz+1)*ratx;
240 blylz = Bz(ix ,iy ,iz )*ratx1+Bz(ix+1,iy ,iz )*ratx;
241 bhz = blyhz *raty1+bhyhz *raty;
242 blz = blylz *raty1+bhylz *raty;
243 b[2] = blz *ratz1+bhz *ratz;
244 //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
245 //ratx,raty,ratz,b[0],b[1],b[2]);
246 //
247 // ... use the dipole symmetry
248 if (x[0]*x[1] < 0) b[1]=-b[1];
249 if (x[0]<0) b[2]=-b[2];
250 } else {
251 printf("Invalid field map for constant mesh %d\n",fMap);
252 }
253 } else {
254//This is the ZDC part
255 Float_t rad2=x[0]*x[0]+x[1]*x[1];
256 if(rad2<D2RA2) {
257 if(x[2]>D2BEG) {
258
259// Separator Dipole D2
260 if(x[2]<D2END) b[1]=FDIP;
261 } else if(x[2]>D1BEG) {
262
263// Separator Dipole D1
264 if(x[2]<D1END) b[1]=-FDIP;
265 }
266 if(rad2<CORRA2) {
267
268// First quadrupole of inner triplet de-focussing in x-direction
269// Inner triplet
270 if(x[2]>Z4BEG) {
271 if(x[2]<Z4END) {
272
273// 2430 <-> 3060
274 b[0]=-G1*x[1];
275 b[1]=-G1*x[0];
276 }
277 } else if(x[2]>Z3BEG) {
278 if(x[2]<Z3END) {
279
280// 1530 <-> 2080
281 b[0]=G1*x[1];
282 b[1]=G1*x[0];
283 }
284 } else if(x[2]>Z2BEG) {
285 if(x[2]<Z2END) {
286
287// 890 <-> 1430
288 b[0]=G1*x[1];
289 b[1]=G1*x[0];
290 }
291 } else if(x[2]>Z1BEG) {
292 if(x[2]<Z1END) {
293
294// 0 <-> 630
295 b[0]=-G1*x[1];
296 b[1]=-G1*x[0];
297 }
298 } else if(x[2]>CORBEG) {
299 if(x[2]<COREND) {
300// Corrector dipole (because of dimuon arm)
301 b[0]=FCORN;
302 }
303 }
304 }
305 }
306 }
307 }
dd045843 308 if(fFactor!=1) {
309 b[0]*=fFactor;
310 b[1]*=fFactor;
311 b[2]*=fFactor;
312 }
fe4da5cc 313}
314
315//________________________________________
316void AliMagFCM::ReadField()
317{
318 FILE *magfile;
319 Int_t ix, iy, iz, ipx, ipy, ipz;
320 Float_t bx, by, bz;
321 char *fname;
322 printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
323 fname = gSystem->ExpandPathName(fTitle.Data());
324 magfile=fopen(fname,"r");
325 delete [] fname;
326 if (magfile) {
327 fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
328 &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
329 printf("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
330 fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
331 fXdeli=1./fXdel;
332 fYdeli=1./fYdel;
333 fZdeli=1./fZdel;
334 fB = new TVector(3*fXn*fYn*fZn);
335 for (iz=0; iz<fZn; iz++) {
336 ipz=iz*3*(fXn*fYn);
337 for (iy=0; iy<fYn; iy++) {
338 ipy=ipz+iy*3*fXn;
339 for (ix=0; ix<fXn; ix++) {
340 ipx=ipy+ix*3;
341 fscanf(magfile,"%f %f %f",&bz,&by,&bx);
342 (*fB)(ipx+2)=bz;
343 (*fB)(ipx+1)=by;
344 (*fB)(ipx )=bx;
345 }
346 }
347 }
348 } else {
349 printf("File %s not found !\n",fTitle.Data());
350 exit(1);
351 }
352}
7a15f6b8 353// -------------------------------------------------------
354
355ClassImp(AliMagFDM)
fe4da5cc 356
7a15f6b8 357//________________________________________
358AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
359const Int_t map, const Float_t factor, const Float_t fmax)
360 : AliMagF(name,title,integ,map,factor,fmax)
361
362{
363 fType = DipoMap;
364
365 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());
366
367}
368
369//________________________________________
370
371void AliMagFDM::Field(Float_t *xfi, Float_t *b)
372{
373 static const Double_t eps=0.1E-06;
374 static const Double_t pi2=.6283185E+01;
375 static const Double_t one=1;
376 static const Double_t fdYaxi = 0.3;
377
378 static const Int_t kiip=33;
379 static const Int_t miip=0;
380 static const Int_t liip=0;
381
382 static const Int_t kiic=0;
383 static const Int_t miic=0;
384 static const Int_t liic=0;
385
386 static const Double_t fdZbg=502.92; // Start of Map using in z
387 static const Double_t fdZL3=600; // Beginning of L3 door in z
388
389 Double_t x[3];
390 Double_t xL3[3];
391 Double_t bint[3];
392
393 Double_t r0;
fe4da5cc 394
7a15f6b8 395 Double_t bbj;
396 Int_t Kvar,jb;
397
398 Double_t Zp1, Zp2,Xp1,Xp2,Yp1,Yp2;
399 Double_t Zz1, Zz2,Yy1,Yy2,X2,X1;
400
401// --- start the map fiel from z = 502.92 cm ---
402
403 x[0] = xfi[0];
404 x[1] = xfi[1];
405 x[2] = xfi[2];
406 b[0]=b[1]=b[2]=0;
407 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
408
409 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
410 r0=rr/100;
411 Double_t Rpmax;
412 Rpmax=fdRmax;
413 if ( (-700<x[2] && x[2]<=fdZbg &&
414 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
415 || (fdZbg<x[2] && x[2]<=fdZL3 && rr>=Rpmax*100) )
416 {
417 b[2]=2;
418 }
419
420 xL3[0]=x[0]/100;
421 xL3[1]=(x[1]+30)/100;
422 xL3[2]=x[2]/100;
423
424 Double_t xminn=xL3[2]*fdAx1+fdCx1;
425 Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
426 Double_t Zcmin,Zcmax,Ycmin,Ycmax;
427
428 Zcmin=fdZmin;
429 Zcmax=fdZmax;
430 Ycmin=fdYmin;
431 Ycmax=fdYmax;
432
433if ((fdZbg/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)))
434 {
435 if(fMap==3)
436 {
437 if (xL3[2]<Zcmin && r0<Rpmax)
438 {
439 //--------------------- Polar part ----------------------
440
441 Double_t yyp,ph0;
442 Int_t kp0, lp0, mp0;
443 Int_t kpi,lpi,mpi;
444 Double_t alp1,alp2,alp3;
445 Double_t zpz,rp,fip,cphi;
446
447 kpi=kiip;
448 lpi=liip;
449 mpi=miip;
450
451 zpz=xL3[2];
452
453 FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&Zp1 ,&Zp2,&fdZpl) ;
454
455 yyp=xL3[1]- 0.3;
456 cphi=yyp/r0;
457 ph0=TMath::ACos(cphi);
458 if (xL3[0]< 0) {ph0=pi2 - ph0;}
459
460 fip=ph0;
461 FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &Xp1,&Xp2,&fdPhin);
462
463 Double_t Rdel;
464 Rdel=fdRdel;
465
466 if (r0<= fdRdel)
467 {
468
469 if(r0< eps)
470 {
471
472 bint[0]=(Zp1*fdBpx[kp0][0][0] + Zp2*fdBpx[kp0+1][0][0])*10;
473 bint[1]=(Zp1*fdBpy[kp0][0][0] + Zp2*fdBpy[kp0+1][0][0])*10;
474 bint[2]=(Zp1*fdBpz[kp0][0][0] + Zp2*fdBpz[kp0+1][0][0])*10;
475
476 }
477
478 alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0];
479 alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];
480 alp1= one - alp2 - alp3;
481
482 for (jb=0; jb<3 ; jb++)
483 {
484 Kvar=jb;
485 FRfuncBi(&Kvar,&Zp1,&Zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
486 bint[jb] = bbj*10 ;
487 }
488 }
489 else
490 {
491 rp=r0;
492
493 FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&Yp1,&Yp2,&fdRn);
494
495 for (jb=0; jb<3 ; jb++)
496 {
497 Kvar=jb;
498 FGfuncBi(&Zp1,&Zp2,&Yp1,&Yp2,&Xp1,&Xp2,&Kvar,&kp0,&lp0,&mp0,&bbj);
499
500 bint[jb] = bbj*10 ;
501 }
502 }
503
504 b[0]=bint[0];
505 b[1]=bint[1];
506 b[2]=bint[2];
507
508// fprintf(fitest,"------------- Freg2 run -------------\n");
509
510 }
511 else
512 {
513 //-------------- Cartensian part ------------------
514
515 Double_t zzc,yyc;
516 Int_t k0, l0,m0;
517 Double_t xx1, xx2,dx, xxx ,xXl;
518 Int_t kci,mci,lci;
519
520 kci=kiic;
521 lci=liic;
522 mci=miic;
523
524 xx1 = fdAx1*xL3[2] + fdCx1;
525 xx2 = fdAx2*xL3[2] + fdCx2;
526
527 zzc=xL3[2];
528 FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &Zz1, &Zz2, &fdZl);
529
530 yyc=xL3[1];
531 FZ(&yyc, fdY , &fdYdel,&lci, &l0, &Yy1, &Yy2,&fdYl);
532
533 xXl = fdXl-one;
534 dx = (xx2-xx1)/xXl;
535 xxx= xL3[0]-xx1;
536 // xm = xxx/dx;
537 m0 = int(xxx/dx);
538
539 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
540 {
541 m0=m0+1;
542 printf(" m0 %d, m0+1 %d\n",m0,m0+1);
543 }
544
545 X2=(xL3[0]-( xx1+m0*dx))/dx;
546 X1=one-X2;
547 m0=m0-1;
548 for (jb=3; jb<6; jb++)
549 {
550 Kvar=jb;
551 FGfuncBi(&Zz1,&Zz2,&Yy1,&Yy2,&X1,&X2,&Kvar,&k0, &l0, &m0, &bbj);
552 bint[jb-3] = bbj*10 ;
553 }
554
555 b[0]=bint[0];
556 b[1]=bint[1];
557 b[2]=bint[2];
558
559// fprintf(fitest,"------------ Freg1 run -----------------\n");
560 }
561
562 } else {
563 printf("Unknown map of Dipole region %d\n",fMap);
564 }
565
566} else {
567
568//This is the ZDC part
569 Float_t rad2=x[0]*x[0]+x[1]*x[1];
570 if(rad2<D2RA2) {
571 if(x[2]>D2BEG) {
572
573// Separator Dipole D2
574 if(x[2]<D2END) b[1]=FDIP;
575 } else if(x[2]>D1BEG) {
576
577// Separator Dipole D1
578 if(x[2]<D1END) b[1]=-FDIP;
579 }
580 if(rad2<CORRA2) {
581
582// First quadrupole of inner triplet de-focussing in x-direction
583// Inner triplet
584 if(x[2]>Z4BEG) {
585 if(x[2]<Z4END) {
586
587// 2430 <-> 3060
588 b[0]=-G1*x[1];
589 b[1]=-G1*x[0];
590 }
591 } else if(x[2]>Z3BEG) {
592 if(x[2]<Z3END) {
593
594// 1530 <-> 2080
595 b[0]=G1*x[1];
596 b[1]=G1*x[0];
597 }
598 } else if(x[2]>Z2BEG) {
599 if(x[2]<Z2END) {
600
601// 890 <-> 1430
602 b[0]=G1*x[1];
603 b[1]=G1*x[0];
604 }
605 } else if(x[2]>Z1BEG) {
606 if(x[2]<Z1END) {
607
608// 0 <-> 630
609 b[0]=-G1*x[1];
610 b[1]=-G1*x[0];
611 }
612 } else if(x[2]>CORBEG) {
613 if(x[2]<COREND) {
614// Corrector dipole (because of dimuon arm)
615// b[0]=FCORN;
616 b[0]=-FCORN;
617 }
618 }
619 }
620 }
621 }
622
623 if(fFactor!=1) {
624 b[0]*=fFactor;
625 b[1]*=fFactor;
626 b[2]*=fFactor;
627 }
628}
629
630//_________________________________________
631
632void 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)
633
634 {
635 static const Double_t one=1;
636 Int_t l,ik,ikj;
637 Double_t temp;
638 Double_t ddu,delu,ar;
639
640 Int_t nk,ku;
641 temp=*u;
642 nk=*nu;
643 ik=*ki;
644 delu=*du;
645
646 ar=Ar[ik];
647 ddu=temp-ar;
648
649 ku=int(ddu/delu);
650 ikj=ik+ku;
651 if (ddu<=0) ikj=0;
652
653 for(l=ikj; l<nk; l++)
654 {
655
656 if(temp < Ar[l])
657 {
658 *kf=l;
659 *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
660 *a1= one - *a2;
661 break;
662 }
663 }
664 }
665
666/*-------------FRfuncBi----------------*/
667
668void 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)
669
670{
671Double_t fa11,fa12,fa13;
672Double_t fa21,fa22,fa23;
673Double_t faY1,faY2;
674Double_t bba;
675
676Double_t zaa1,zaa2,alf1,alf2,alf3;
677Int_t kaai,kaa,maa;
678kaai=*kai;
679kaa=*ka;
680maa=*ma;
681zaa1=*za1;
682zaa2=*za2;
683alf1=*al1;
684alf2=*al2;
685alf3=*al3;
686
687 if (kaai==0 ) {
688 fa11 = fdBpx[kaa][0][0];
689 fa12 = fdBpx[kaa][0][maa];
690 fa13 = fdBpx[kaa][0][maa+1];
691 fa21 = fdBpx[kaa+1][0][0];
692 fa22 = fdBpx[kaa+1][0][maa];
693 fa23 = fdBpx[kaa+1][0][maa+1];
694 }
695 if (kaai==1 ) {
696 fa11 = fdBpy[kaa][0][0];
697 fa12 = fdBpy[kaa][0][maa];
698 fa13 = fdBpy[kaa][0][maa+1];
699 fa21 = fdBpy[kaa+1][0][0];
700 fa22 = fdBpy[kaa+1][0][maa];
701 fa23 = fdBpy[kaa+1][0][maa+1];
702 }
703 if (kaai==2 ) {
704 fa11 = fdBpz[kaa][0][0];
705 fa12 = fdBpz[kaa][0][maa];
706 fa13 = fdBpz[kaa][0][maa+1];
707 fa21 = fdBpz[kaa+1][0][0];
708 fa22 = fdBpz[kaa+1][0][maa];
709 fa23 = fdBpz[kaa+1][0][maa+1];
710 }
711 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
712 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
713 bba = zaa1*faY1+zaa2*faY2;
714 *ba=bba;
715
716}
717
718
719/*----------- FGfuncBi------------*/
720
721void 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)
722
723{
724Double_t fy1, fy2, ffy;
725Double_t gy1,gy2,ggy;
726Double_t z1,z2,y1,y2,x1,x2;
727
728Int_t k,l,m,kv;
729Double_t bbi;
730
731Double_t bf11,bf12,bf21,bf22;
732Double_t bg11,bg12,bg21,bg22;
733k=*kk;
734l=*ll;
735m=*mm;
736
737kv=*kvr;
738
739z1=*zz1;
740z2=*zz2;
741y1=*yy1;
742y2=*yy2;
743x1=*xx1;
744x2=*xx2;
745
746/*-----------------Polar part ------------------*/
747
748if(kv==0) {
749 bf11=fdBpx[k][l][m];
750 bf12=fdBpx[k+1][l][m];
751 bf21=fdBpx[k+1][l+1][m];
752 bf22=fdBpx[k][l+1][m];
753
754 bg11=fdBpx[k][l][m+1];
755 bg12=fdBpx[k+1][l][m+1];
756 bg21=fdBpx[k+1][l+1][m+1];
757 bg22=fdBpx[k][l+1][m+1];
758 }
759 if(kv==1) {
760 bf11=fdBpy[k][l][m];
761 bf12=fdBpy[k+1][l][m];
762 bf21=fdBpy[k+1][l+1][m];
763 bf22=fdBpy[k][l+1][m];
764
765 bg11=fdBpy[k][l][m+1];
766 bg12=fdBpy[k+1][l][m+1];
767 bg21=fdBpy[k+1][l+1][m+1];
768 bg22=fdBpy[k][l+1][m+1];
769 }
770
771 if(kv==2) {
772 bf11=fdBpz[k][l][m];
773 bf12=fdBpz[k+1][l][m];
774 bf21=fdBpz[k+1][l+1][m];
775 bf22=fdBpz[k][l+1][m];
776
777 bg11=fdBpz[k][l][m+1];
778 bg12=fdBpz[k+1][l][m+1];
779 bg21=fdBpz[k+1][l+1][m+1];
780 bg22=fdBpz[k][l+1][m+1];
781 }
782/*-----------------Cartensian part ---------------*/
783
784 if(kv==3) {
785 bf11=fdBcx[k][l][m];
786 bf12=fdBcx[k+1][l][m];
787 bf21=fdBcx[k+1][l+1][m];
788 bf22=fdBcx[k][l+1][m];
789
790 bg11=fdBcx[k][l][m+1];
791 bg12=fdBcx[k+1][l][m+1];
792 bg21=fdBcx[k+1][l+1][m+1];
793 bg22=fdBcx[k][l+1][m+1];
794 }
795
796 if(kv==4) {
797 bf11=fdBcy[k][l][m];
798 bf12=fdBcy[k+1][l][m];
799 bf21=fdBcy[k+1][l+1][m];
800 bf22=fdBcy[k][l+1][m];
801
802 bg11=fdBcy[k][l][m+1];
803 bg12=fdBcy[k+1][l][m+1];
804 bg21=fdBcy[k+1][l+1][m+1];
805 bg22=fdBcy[k][l+1][m+1];
806 }
807 if(kv==5) {
808 bf11=fdBcz[k][l][m];
809 bf12=fdBcz[k+1][l][m];
810 bf21=fdBcz[k+1][l+1][m];
811 bf22=fdBcz[k][l+1][m];
812
813 bg11=fdBcz[k][l][m+1];
814 bg12=fdBcz[k+1][l][m+1];
815 bg21=fdBcz[k+1][l+1][m+1];
816 bg22=fdBcz[k][l+1][m+1];
817 }
818
819
820
821 fy1=z1*bf11+z2*bf12;
822 fy2=z2*bf21+z1* bf22;
823 ffy=y1*fy1+ y2*fy2;
824
825
826 gy1 = z1*bg11+z2*bg12;
827 gy2 = z2*bg21+z1*bg22;
828 ggy= y1*gy1 + y2*gy2;
829
830 bbi = x1*ffy+x2*ggy;
831
832 *bb=bbi;
833
834}
835//____________________________________________
836
837void AliMagFDM::ReadField()
838{
839 FILE *magfile;
840
841 Int_t ik, il, im;
842 Float_t zzp, rr,phii;
843 Float_t zz, yy, bx,by,bz,bb;
844
845 char *fname;
846 printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
847 fname = gSystem->ExpandPathName(fTitle.Data());
848 magfile=fopen(fname,"r");
849 delete [] fname;
850
851 printf("Cartensian part\n");
852
853 if (magfile) {
854
855// Cartensian part
856
857 fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl);
858
859 printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);
860
861 for (ik=0; ik<fdZl; ik++)
862 {
863
864 fscanf(magfile, " %e ", &zz);
865 fdZc[ik]=zz;
866
867 }
868
869 for (ik=0; ik<fdYl; ik++)
870 {
871 fscanf(magfile, " %e ", &yy);
872 fdY[ik]=yy;
873
874 }
875 for (ik=0; ik<81; ik++)
876 {
877 printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]);
878 }
879
880 fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2);
881
882printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
883printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
884printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
885
886 for (il=0; il<44; il++) {
887 for (im=0; im<81; im++) {
888 for (ik=0; ik<81; ik++) {
889
890 fscanf(magfile, " %e ", &by);
891 fdBcy[ik][im][il]=by;
892 }
893 }
894 }
895
896 for (il=0; il<44; il++) {
897 for (im=0; im<81; im++) {
898 for (ik=0; ik<81; ik++) {
899
900 fscanf(magfile, " %e ", &bx);
901 fdBcx[ik][im][il]=bx;
902 }
903 }
904 }
905
906 for (il=0; il<44; il++) {
907 for (im=0; im<81; im++) {
908 for (ik=0; ik<81; ik++) {
909
910 fscanf(magfile, " %e ", &bz);
911 fdBcz[ik][im][il]=bz;
912 }
913 }
914 }
915//---------------------- Polar part ---------------------------------
916
917 printf("Polar part\n");
918 fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin);
919 printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);
920
921 printf(" fdZp array\n");
922
923 for (ik=0; ik<51; ik++)
924 {
925 fscanf(magfile, " %e ", &zzp);
926 fdZp[ik]=zzp;
927 printf(" %e\n",fdZp[ik]);
928 }
929
930 printf(" fdR array\n");
931
932 for (ik=0; ik<10; ik++)
933 {
934 fscanf(magfile, " %e ", &rr);
935 fdR[ik]=rr;
936 printf(" %e\n",fdR[ik]);
937 }
938
939// printf("fdPhi array\n");
940
941 for (il=0; il<33; il++)
942 {
943 fscanf(magfile, " %e ", &phii);
944 fdPhi[il]=phii;
945// printf(" %e\n",fdPhi[il]);
946 }
947
948 fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin);
949
950printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
951
952
953 for (il=0; il<33; il++) {
954 for (im=0; im<10; im++) {
955 for (ik=0; ik<51; ik++) {
956 fscanf(magfile, " %e ", &by);
957 fdBpy[ik][im][il]=by;
958 }
959 }
960 }
961
962 for (il=0; il<33; il++) {
963 for (im=0; im<10; im++) {
964 for (ik=0; ik<51; ik++) {
965 fscanf(magfile, " %e ", &bx);
966 fdBpx[ik][im][il]=bx;
967 }
968 }
969 }
970
971
972 for (il=0; il<33; il++) {
973 for (im=0; im<10; im++) {
974 for (ik=0; ik<51; ik++) {
975 fscanf(magfile, " %e ", &bz);
976 fdBpz[ik][im][il]=bz;
977 }
978 }
979 }
980
981
982 for (il=0; il<32; il++) {
983 for (im=0; im<2; im++) {
984 for (ik=0; ik<2; ik++) {
985 fscanf(magfile, " %e ", &bb);
986 fdB[ik][im][il]=bb;
987 }
988 }
989 }
990//
991 } else {
992 printf("File %s not found !\n",fTitle.Data());
993 exit(1);
994 }
995}
996//________________________________
997
998
999
1000
1001
1002
fe4da5cc 1003