]> 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
Removing the hard-wired particle masses (B. Hippolyte)
[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
acd84897 16/* $Id$ */
aee8290b 17
116cbefd 18//-------------------------------------------------------------------------
116cbefd 19// Field with Magnetic Field map
5d8718b8 20// Used by AliRun class
21// Author:
116cbefd 22//-------------------------------------------------------------------------
23
090026bf 24#include <TMath.h>
25#include <TSystem.h>
e2afb3b6 26
594d8990 27#include "AliLog.h"
e2afb3b6 28#include "AliMagFDM.h"
aee8290b 29
30ClassImp(AliMagFDM)
31
e2afb3b6 32//_______________________________________________________________________
33AliMagFDM::AliMagFDM():
34 fSolenoid(0),
35 fInd(0),
36 fZmin(0),
37 fZmax(0),
38 fYmax(0),
39 fYmin(0),
40 fZpmx(0),
41 fZpmn(0),
42 fRmax(0),
43 fRmin(0),
44 fXdel(0),
45 fYdel(0),
46 fZdel(0),
47 fRdel(0),
48 fPhid(0),
49 fZpdl(0),
50 fCx1(0),
51 fCx2(0),
52 fAx1(0),
53 fAx2(0),
54 fXl(0),
55 fYl(0),
56 fZl(0),
57 fRn(0),
58 fPhin(0),
59 fZpl(0)
60{
61 //
62 // Default constructor for the Dipole field
63 //
64}
65
66//_______________________________________________________________________
d0f1ee3b 67AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
68 Float_t factor, Float_t fmax):
57754f18 69 AliMagFC(name,title,integ,factor,fmax),
e2afb3b6 70 fSolenoid(0),
71 fInd(0),
72 fZmin(0),
73 fZmax(0),
74 fYmax(0),
75 fYmin(0),
76 fZpmx(0),
77 fZpmn(0),
78 fRmax(0),
79 fRmin(0),
80 fXdel(0),
81 fYdel(0),
82 fZdel(0),
83 fRdel(0),
84 fPhid(0),
85 fZpdl(0),
86 fCx1(0),
87 fCx2(0),
88 fAx1(0),
89 fAx2(0),
90 fXl(0),
91 fYl(0),
92 fZl(0),
93 fRn(0),
94 fPhin(0),
95 fZpl(0)
aee8290b 96{
ef42d733 97 //
98 // Standard constructor for the Dipole field
99 //
aee8290b 100 fType = kDipoMap;
d8408e76 101 fMap = 3;
4cc8933f 102 SetSolenoidField();
d8408e76 103
594d8990 104 AliDebug(1, Form(
105 "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s",
106 fName.Data(), fMap ,integ,factor,fTitle.Data()));
e678332c 107
aee8290b 108}
109
e2afb3b6 110//_______________________________________________________________________
611fa94a 111void AliMagFDM::Field(const float *xfi, float *b) const
aee8290b 112{
113 //
114 // Main routine to compute the field in a point
115 //
e2afb3b6 116 const Double_t keps=0.1E-06;
116cbefd 117 const Double_t kPI2=2.*TMath::Pi();
e2afb3b6 118 const Double_t kone=1;
119
120 const Int_t kiip=33;
121 const Int_t kmiip=0;
122 const Int_t kliip=0;
123
124 const Int_t kiic=0;
125 const Int_t kmiic=0;
126 const Int_t kliic=0;
127
128 const Double_t kfZbg=502.92; // Start of Map using in z
129 const Double_t kfZL3=600; // Beginning of L3 door in z
aee8290b 130
131 Double_t x[3];
132 Double_t xL3[3];
133 Double_t bint[3];
134
135 Double_t r0;
e2afb3b6 136 Int_t iKvar,jb;
aee8290b 137
138 Double_t zp1, zp2,xp1,xp2,yp1,yp2;
139 Double_t zz1, zz2,yy1,yy2,x2,x1;
140
141// --- start the map fiel from z = 502.92 cm ---
57754f18 142//
143// This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
144// Transfor correspondingly.
aee8290b 145
57754f18 146 x[0] = - xfi[0];
147 x[1] = xfi[1];
148 x[2] = - xfi[2];
e678332c 149 b[0]=b[1]=b[2]=0;
57754f18 150//
aee8290b 151 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
152
153 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
154 r0=rr/100;
155 Double_t rPmax;
803c9752 156 rPmax=fRmax;
157 if ( (-700<x[2] && x[2]<=kfZbg &&
aee8290b 158 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
803c9752 159 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
aee8290b 160 {
803c9752 161 b[0]=b[1]=0;
4cc8933f 162 b[2]=fSolenoid;
aee8290b 163 }
164
165 xL3[0]=x[0]/100;
166 xL3[1]=(x[1]+30)/100;
167 xL3[2]=x[2]/100;
803c9752 168
aee8290b 169
803c9752 170 if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
171 if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
172
173 Double_t xminn=xL3[2]*fAx1+fCx1;
174 Double_t xmaxx=xL3[2]*fAx2+fCx2;
aee8290b 175 Double_t zCmin,zCmax,yCmin,yCmax;
176
803c9752 177 zCmin=fZmin;
178 zCmax=fZmax;
179 yCmin=fYmin;
180 yCmax=fYmax;
803c9752 181if ((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 182 {
183 if(fMap==3)
184 {
803c9752 185 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
aee8290b 186 {
187 //--------------------- Polar part ----------------------
188
189 Double_t yyp,ph0;
190 Int_t kp0, lp0, mp0;
191 Int_t kpi,lpi,mpi;
192 Double_t alp1,alp2,alp3;
193 Double_t zpz,rp,fip,cphi;
194
195 kpi=kiip;
196 lpi=kliip;
197 mpi=kmiip;
198
199 zpz=xL3[2];
e678332c 200 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
201 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
202 zp1= kone-zp2;
aee8290b 203 yyp=xL3[1]- 0.3;
e678332c 204 cphi=TMath::Abs(yyp/r0);
205 Int_t kcphi=0;
206 if (cphi > kone) {
594d8990 207 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));
e678332c 208 cphi =kone;
209 kcphi=777;
210 }
aee8290b 211 ph0=TMath::ACos(cphi);
116cbefd 212 if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
213 if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
214 if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
215 if (ph0 > kPI2) { ph0=ph0 - kPI2;}
e678332c 216 if (kcphi==777) {
594d8990 217 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));
e678332c 218 }
219 fip=ph0;
220 mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
221 xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
222 xp1=kone-xp2;
aee8290b 223
224 Double_t rDel;
803c9752 225 rDel=fRdel;
aee8290b 226
803c9752 227 if (r0<= fRdel)
aee8290b 228 {
229
230 if(r0< keps)
231 {
232
803c9752 233 bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
234 bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
235 bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
aee8290b 236
803c9752 237 } else {
aee8290b 238
803c9752 239 alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
240 alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
aee8290b 241 alp1= kone - alp2 - alp3;
242
243 for (jb=0; jb<3 ; jb++)
803c9752 244 {
245 iKvar=jb;
e678332c 246 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
803c9752 247 }
248 } // end of keps <=r0
aee8290b 249 }
250 else
251 {
252 rp=r0;
253
e678332c 254 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
255 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
256 yp1=kone-yp2;
aee8290b 257
258 for (jb=0; jb<3 ; jb++)
259 {
260 iKvar=jb;
e678332c 261 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
aee8290b 262 }
263 }
264
57754f18 265 b[0]=-bint[0];
aee8290b 266 b[1]=bint[1];
57754f18 267 b[2]=-bint[2];
aee8290b 268
aee8290b 269 }
270 else
271 {
272 //-------------- Cartensian part ------------------
273
274 Double_t zzc,yyc;
275 Int_t k0, l0,m0;
276 Double_t xx1, xx2,dx, xxx ,xXl;
277 Int_t kci,mci,lci;
278
279 kci=kiic;
280 lci=kliic;
281 mci=kmiic;
282
803c9752 283 xx1 = fAx1*xL3[2] + fCx1;
284 xx2 = fAx2*xL3[2] + fCx2;
aee8290b 285
286 zzc=xL3[2];
e678332c 287 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
288 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
289 zz1=kone - zz2;;
aee8290b 290
291 yyc=xL3[1];
e678332c 292 l0=FZ(yyc, fY , fYdel,lci,fYl);
293 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
294 yy1=kone - yy2;
803c9752 295 xXl = fXl-kone;
aee8290b 296 dx = (xx2-xx1)/xXl;
297 xxx= xL3[0]-xx1;
aee8290b 298 m0 = int(xxx/dx);
299
300 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
301 {
302 m0=m0+1;
594d8990 303 AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
aee8290b 304 }
305
306 x2=(xL3[0]-( xx1+m0*dx))/dx;
307 x1=kone-x2;
308 m0=m0-1;
309 for (jb=3; jb<6; jb++)
310 {
311 iKvar=jb;
e678332c 312 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
aee8290b 313 }
314
57754f18 315 b[0]=-bint[0];
aee8290b 316 b[1]=bint[1];
57754f18 317 b[2]=-bint[2];
aee8290b 318
aee8290b 319 }
e678332c 320
aee8290b 321
322 } else {
594d8990 323 AliError(Form("Unknown map of Dipole region %d",fMap));
aee8290b 324 }
325
326} else {
57754f18 327 ZDCField(xfi,b);
aee8290b 328
d95ea23f 329 }
aee8290b 330
331 if(fFactor!=1) {
332 b[0]*=fFactor;
333 b[1]*=fFactor;
334 b[2]*=fFactor;
335 }
336}
337
ff66b122 338//_______________________________________________________________________
611fa94a 339void AliMagFDM::Field(const double *xfi, double *b) const
ff66b122 340{
341 //
342 // Main routine to compute the field in a point
343 //
344 const Double_t keps=0.1E-06;
345 const Double_t kPI2=2.*TMath::Pi();
346 const Double_t kone=1;
347
348 const Int_t kiip=33;
349 const Int_t kmiip=0;
350 const Int_t kliip=0;
351
352 const Int_t kiic=0;
353 const Int_t kmiic=0;
354 const Int_t kliic=0;
355
356 const Double_t kfZbg=502.92; // Start of Map using in z
357 const Double_t kfZL3=600; // Beginning of L3 door in z
358
359 Double_t x[3];
360 Double_t xL3[3];
361 Double_t bint[3];
362
363 Double_t r0;
364 Int_t iKvar,jb;
365
366 Double_t zp1, zp2,xp1,xp2,yp1,yp2;
367 Double_t zz1, zz2,yy1,yy2,x2,x1;
368
369// --- start the map fiel from z = 502.92 cm ---
370//
371// This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
372// Transfor correspondingly.
373
374 x[0] = - xfi[0];
375 x[1] = xfi[1];
376 x[2] = - xfi[2];
377 b[0]=b[1]=b[2]=0;
378//
379 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
380
381 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
382 r0=rr/100;
383 Double_t rPmax;
384 rPmax=fRmax;
385 if ( (-700<x[2] && x[2]<=kfZbg &&
386 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
387 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
388 {
389 b[0]=b[1]=0;
390 b[2]=fSolenoid;
391 }
392
393 xL3[0]=x[0]/100;
394 xL3[1]=(x[1]+30)/100;
395 xL3[2]=x[2]/100;
396
397
398 if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
399 if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
400
401 Double_t xminn=xL3[2]*fAx1+fCx1;
402 Double_t xmaxx=xL3[2]*fAx2+fCx2;
403 Double_t zCmin,zCmax,yCmin,yCmax;
404
405 zCmin=fZmin;
406 zCmax=fZmax;
407 yCmin=fYmin;
408 yCmax=fYmax;
409if ((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)))
410 {
411 if(fMap==3)
412 {
413 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
414 {
415 //--------------------- Polar part ----------------------
416
417 Double_t yyp,ph0;
418 Int_t kp0, lp0, mp0;
419 Int_t kpi,lpi,mpi;
420 Double_t alp1,alp2,alp3;
421 Double_t zpz,rp,fip,cphi;
422
423 kpi=kiip;
424 lpi=kliip;
425 mpi=kmiip;
426
427 zpz=xL3[2];
428 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
429 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
430 zp1= kone-zp2;
431 yyp=xL3[1]- 0.3;
432 cphi=TMath::Abs(yyp/r0);
433 Int_t kcphi=0;
434 if (cphi > kone) {
435 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));
436 cphi =kone;
437 kcphi=777;
438 }
439 ph0=TMath::ACos(cphi);
440 if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
441 if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
442 if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
443 if (ph0 > kPI2) { ph0=ph0 - kPI2;}
444 if (kcphi==777) {
445 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));
446 }
447 fip=ph0;
448 mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
449 xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
450 xp1=kone-xp2;
451
452 Double_t rDel;
453 rDel=fRdel;
454
455 if (r0<= fRdel)
456 {
457
458 if(r0< keps)
459 {
460
461 bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
462 bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
463 bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
464
465 } else {
466
467 alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
468 alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
469 alp1= kone - alp2 - alp3;
470
471 for (jb=0; jb<3 ; jb++)
472 {
473 iKvar=jb;
474 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
475 }
476 } // end of keps <=r0
477 }
478 else
479 {
480 rp=r0;
481
482 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
483 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
484 yp1=kone-yp2;
485
486 for (jb=0; jb<3 ; jb++)
487 {
488 iKvar=jb;
489 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
490 }
491 }
492
493 b[0]=-bint[0];
494 b[1]=bint[1];
495 b[2]=-bint[2];
496
497 }
498 else
499 {
500 //-------------- Cartensian part ------------------
501
502 Double_t zzc,yyc;
503 Int_t k0, l0,m0;
504 Double_t xx1, xx2,dx, xxx ,xXl;
505 Int_t kci,mci,lci;
506
507 kci=kiic;
508 lci=kliic;
509 mci=kmiic;
510
511 xx1 = fAx1*xL3[2] + fCx1;
512 xx2 = fAx2*xL3[2] + fCx2;
513
514 zzc=xL3[2];
515 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
516 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
517 zz1=kone - zz2;;
518
519 yyc=xL3[1];
520 l0=FZ(yyc, fY , fYdel,lci,fYl);
521 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
522 yy1=kone - yy2;
523 xXl = fXl-kone;
524 dx = (xx2-xx1)/xXl;
525 xxx= xL3[0]-xx1;
526 m0 = int(xxx/dx);
527
528 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
529 {
530 m0=m0+1;
531 AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
532 }
533
534 x2=(xL3[0]-( xx1+m0*dx))/dx;
535 x1=kone-x2;
536 m0=m0-1;
537 for (jb=3; jb<6; jb++)
538 {
539 iKvar=jb;
540 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
541 }
542
543 b[0]=-bint[0];
544 b[1]=bint[1];
545 b[2]=-bint[2];
546
547 }
548
549
550 } else {
551 AliError(Form("Unknown map of Dipole region %d",fMap));
552 }
553
554} else {
555 ZDCField(xfi,b);
556
557 }
558
559 if(fFactor!=1) {
560 b[0]*=fFactor;
561 b[1]*=fFactor;
562 b[2]*=fFactor;
563 }
564}
565
566
e2afb3b6 567//_______________________________________________________________________
6f3038e9 568Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar,
ecd57058 569 Float_t delu, Int_t ik, Int_t nk) const
aee8290b 570{
571 //
e678332c 572 // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
aee8290b 573 //
e678332c 574 Int_t l,ikj;
575 Double_t ddu,ar;
aee8290b 576
e678332c 577 Int_t ku,kf;
578 kf=0;
aee8290b 579 ar=Ar[ik];
580 ddu=temp-ar;
581
582 ku=int(ddu/delu);
e678332c 583 ikj=ik+ku+1;
aee8290b 584 if (ddu<=0) ikj=0;
585
586 for(l=ikj; l<nk; l++)
587 {
e678332c 588
803c9752 589 if(temp <= Ar[l])
aee8290b 590 {
803c9752 591
e678332c 592 kf=l-1;
593 break;
aee8290b 594 }
595 }
e678332c 596 return kf;
aee8290b 597 }
598
e2afb3b6 599//_______________________________________________________________________
600Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2,
601 Double_t alf1, Double_t alf2, Double_t alf3,
6f3038e9 602 Int_t kaa, Int_t maa) const
aee8290b 603{
604 //
e678332c 605 // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
aee8290b 606 //
0742d588 607 Double_t fa11=0,fa12=0,fa13=0;
608 Double_t fa21=0,fa22=0,fa23=0;
609 Double_t faY1=0,faY2=0;
610 Double_t bba=0;
e678332c 611
8918e700 612 switch (kaai) {
613 case 0:
803c9752 614 fa11 = fBpx[kaa][0][0];
615 fa12 = fBpx[kaa][0][maa];
616 fa13 = fBpx[kaa][0][maa+1];
617 fa21 = fBpx[kaa+1][0][0];
618 fa22 = fBpx[kaa+1][0][maa];
619 fa23 = fBpx[kaa+1][0][maa+1];
8918e700 620 break;
621 case 1:
803c9752 622 fa11 = fBpy[kaa][0][0];
623 fa12 = fBpy[kaa][0][maa];
624 fa13 = fBpy[kaa][0][maa+1];
625 fa21 = fBpy[kaa+1][0][0];
626 fa22 = fBpy[kaa+1][0][maa];
627 fa23 = fBpy[kaa+1][0][maa+1];
8918e700 628 break;
629 case 2:
803c9752 630 fa11 = fBpz[kaa][0][0];
631 fa12 = fBpz[kaa][0][maa];
632 fa13 = fBpz[kaa][0][maa+1];
633 fa21 = fBpz[kaa+1][0][0];
634 fa22 = fBpz[kaa+1][0][maa];
635 fa23 = fBpz[kaa+1][0][maa+1];
8918e700 636 break;
637 default:
594d8990 638 AliFatal(Form("Invalid value of kaai %d",kaai));
aee8290b 639 }
640 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
641 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
642 bba = zaa1*faY1+zaa2*faY2;
e678332c 643 return bba;
aee8290b 644}
645
646
e2afb3b6 647//_______________________________________________________________________
648Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2,
6f3038e9 649 Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m) const
aee8290b 650{
651 //
e678332c 652 // Calculation of field componet at a given axis (general case)
aee8290b 653 //
0742d588 654 Double_t fy1=0,fy2=0,ffy=0;
655 Double_t gy1=0,gy2=0,ggy=0;
656 Double_t bbi=0;
657 Double_t bf11=0,bf12=0,bf21=0,bf22=0;
658 Double_t bg11=0,bg12=0,bg21=0,bg22=0;
e678332c 659
aee8290b 660
661 /*-----------------Polar part ------------------*/
662
8918e700 663 switch (kv) {
664 case 0:
803c9752 665 bf11=fBpx[k][l][m];
666 bf12=fBpx[k+1][l][m];
667 bf21=fBpx[k+1][l+1][m];
668 bf22=fBpx[k][l+1][m];
aee8290b 669
803c9752 670 bg11=fBpx[k][l][m+1];
671 bg12=fBpx[k+1][l][m+1];
672 bg21=fBpx[k+1][l+1][m+1];
673 bg22=fBpx[k][l+1][m+1];
8918e700 674 break;
675
676 case 1:
803c9752 677 bf11=fBpy[k][l][m];
678 bf12=fBpy[k+1][l][m];
679 bf21=fBpy[k+1][l+1][m];
680 bf22=fBpy[k][l+1][m];
aee8290b 681
803c9752 682 bg11=fBpy[k][l][m+1];
683 bg12=fBpy[k+1][l][m+1];
684 bg21=fBpy[k+1][l+1][m+1];
685 bg22=fBpy[k][l+1][m+1];
8918e700 686 break;
687
688 case 2:
803c9752 689 bf11=fBpz[k][l][m];
690 bf12=fBpz[k+1][l][m];
691 bf21=fBpz[k+1][l+1][m];
692 bf22=fBpz[k][l+1][m];
aee8290b 693
803c9752 694 bg11=fBpz[k][l][m+1];
695 bg12=fBpz[k+1][l][m+1];
696 bg21=fBpz[k+1][l+1][m+1];
697 bg22=fBpz[k][l+1][m+1];
8918e700 698 break;
aee8290b 699 /*-----------------Cartensian part ---------------*/
700
8918e700 701 case 3:
803c9752 702 bf11=fBcx[k][l][m];
703 bf12=fBcx[k+1][l][m];
704 bf21=fBcx[k+1][l+1][m];
705 bf22=fBcx[k][l+1][m];
aee8290b 706
803c9752 707 bg11=fBcx[k][l][m+1];
708 bg12=fBcx[k+1][l][m+1];
709 bg21=fBcx[k+1][l+1][m+1];
710 bg22=fBcx[k][l+1][m+1];
8918e700 711 break;
712
713 case 4:
803c9752 714 bf11=fBcy[k][l][m];
715 bf12=fBcy[k+1][l][m];
716 bf21=fBcy[k+1][l+1][m];
717 bf22=fBcy[k][l+1][m];
aee8290b 718
803c9752 719 bg11=fBcy[k][l][m+1];
720 bg12=fBcy[k+1][l][m+1];
721 bg21=fBcy[k+1][l+1][m+1];
722 bg22=fBcy[k][l+1][m+1];
8918e700 723 break;
724
725 case 5:
803c9752 726 bf11=fBcz[k][l][m];
727 bf12=fBcz[k+1][l][m];
728 bf21=fBcz[k+1][l+1][m];
729 bf22=fBcz[k][l+1][m];
aee8290b 730
803c9752 731 bg11=fBcz[k][l][m+1];
732 bg12=fBcz[k+1][l][m+1];
733 bg21=fBcz[k+1][l+1][m+1];
734 bg22=fBcz[k][l+1][m+1];
8918e700 735 break;
736
737 default:
594d8990 738 AliFatal(Form("Invalid value of kv %d",kv));
aee8290b 739 }
740
741
742
743 fy1=z1*bf11+z2*bf12;
744 fy2=z2*bf21+z1* bf22;
745 ffy=y1*fy1+ y2*fy2;
746
747
748 gy1 = z1*bg11+z2*bg12;
749 gy2 = z2*bg21+z1*bg22;
750 ggy= y1*gy1 + y2*gy2;
751
752 bbi = x1*ffy+x2*ggy;
e678332c 753
754 return bbi;
aee8290b 755
756}
aee8290b 757
e2afb3b6 758//_______________________________________________________________________
aee8290b 759void AliMagFDM::ReadField()
760{
761 //
762 // Method to read the magnetic field from file
763 //
764 FILE *magfile;
765
766 Int_t ik, il, im;
767 Float_t zzp, rr,phii;
768 Float_t zz, yy, bx,by,bz,bb;
769
770 char *fname;
594d8990 771 AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
aee8290b 772 fname = gSystem->ExpandPathName(fTitle.Data());
773 magfile=fopen(fname,"r");
774 delete [] fname;
775
594d8990 776 AliDebug(2,"Cartensian part");
aee8290b 777
778 if (magfile) {
779
780// Cartensian part
781
803c9752 782 fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
aee8290b 783
594d8990 784 AliDebug(3,Form("fYl %d, fXl %d, fZl %d",fYl, fXl, fZl));
aee8290b 785
803c9752 786 for (ik=0; ik<fZl; ik++)
aee8290b 787 {
788
789 fscanf(magfile, " %e ", &zz);
803c9752 790 fZc[ik]=zz;
aee8290b 791
792 }
793
803c9752 794 for (ik=0; ik<fYl; ik++)
aee8290b 795 {
796 fscanf(magfile, " %e ", &yy);
803c9752 797 fY[ik]=yy;
aee8290b 798
799 }
800 for (ik=0; ik<81; ik++)
801 {
594d8990 802 AliDebug(4,Form("fZc %e,fY %e", fZc[ik],fY[ik]));
aee8290b 803 }
804
803c9752 805 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 806
594d8990 807AliDebug(3,Form("fYdel %e, fXdel %e, fZdel %e",fYdel,fXdel,fZdel));
808AliDebug(3,Form("fZmax %e, fZmin %e, fYmax %e,fYmin %e",fZmax,fZmin,fYmax,fYmin));
809AliDebug(3,Form("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e",fAx1,fCx1,fAx2,fCx2));
aee8290b 810
811 for (il=0; il<44; il++) {
812 for (im=0; im<81; im++) {
813 for (ik=0; ik<81; ik++) {
814
815 fscanf(magfile, " %e ", &by);
803c9752 816 fBcy[ik][im][il]=by;
aee8290b 817 }
818 }
819 }
820
821 for (il=0; il<44; il++) {
822 for (im=0; im<81; im++) {
823 for (ik=0; ik<81; ik++) {
824
825 fscanf(magfile, " %e ", &bx);
803c9752 826 fBcx[ik][im][il]=bx;
aee8290b 827 }
828 }
829 }
830
831 for (il=0; il<44; il++) {
832 for (im=0; im<81; im++) {
833 for (ik=0; ik<81; ik++) {
834
835 fscanf(magfile, " %e ", &bz);
803c9752 836 fBcz[ik][im][il]=bz;
aee8290b 837 }
838 }
839 }
840//---------------------- Polar part ---------------------------------
841
594d8990 842 AliDebug(2,"Polar part");
803c9752 843 fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
594d8990 844 AliDebug(3,Form("fZpl %d, fRn %d, fPhin %d",fZpl,fRn,fPhin));
aee8290b 845
594d8990 846 AliDebug(4," fZp array");
aee8290b 847
848 for (ik=0; ik<51; ik++)
849 {
850 fscanf(magfile, " %e ", &zzp);
803c9752 851 fZp[ik]=zzp;
594d8990 852 AliDebug(4,Form(" %e",fZp[ik]));
aee8290b 853 }
854
594d8990 855 AliDebug(4," fR array");
aee8290b 856
857 for (ik=0; ik<10; ik++)
858 {
859 fscanf(magfile, " %e ", &rr);
803c9752 860 fR[ik]=rr;
594d8990 861 AliDebug(4,Form(" %e",fR[ik]));
aee8290b 862 }
863
594d8990 864// AliDebug(4,"fPhi array");
aee8290b 865
866 for (il=0; il<33; il++)
867 {
868 fscanf(magfile, " %e ", &phii);
803c9752 869 fPhi[il]=phii;
594d8990 870// AliDebug(4,Form(" %e",fPhi[il]));
aee8290b 871 }
872
803c9752 873 fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
aee8290b 874
594d8990 875AliDebug(3,Form("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin));
aee8290b 876
877
878 for (il=0; il<33; il++) {
879 for (im=0; im<10; im++) {
880 for (ik=0; ik<51; ik++) {
881 fscanf(magfile, " %e ", &by);
803c9752 882 fBpy[ik][im][il]=by;
aee8290b 883 }
884 }
885 }
886
887 for (il=0; il<33; il++) {
888 for (im=0; im<10; im++) {
889 for (ik=0; ik<51; ik++) {
890 fscanf(magfile, " %e ", &bx);
803c9752 891 fBpx[ik][im][il]=bx;
aee8290b 892 }
893 }
894 }
895
896
897 for (il=0; il<33; il++) {
898 for (im=0; im<10; im++) {
899 for (ik=0; ik<51; ik++) {
900 fscanf(magfile, " %e ", &bz);
803c9752 901 fBpz[ik][im][il]=bz;
aee8290b 902 }
903 }
904 }
905
906
907 for (il=0; il<32; il++) {
908 for (im=0; im<2; im++) {
909 for (ik=0; ik<2; ik++) {
910 fscanf(magfile, " %e ", &bb);
803c9752 911 fB[ik][im][il]=bb;
aee8290b 912 }
913 }
914 }
915//
916 } else {
594d8990 917 AliFatal(Form("File %s not found !",fTitle.Data()));
aee8290b 918 }
919}