]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - STEER/AliMagFDM.cxx
First Run and Last Run members added, infinite validity of calib obj implemented.
[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/* $Id$ */
17
18//-------------------------------------------------------------------------
19// Field with Magnetic Field map
20// Used by AliRun class
21// Author:
22//-------------------------------------------------------------------------
23
24#include <TMath.h>
25#include <TSystem.h>
26
27#include "AliLog.h"
28#include "AliMagFDM.h"
29
30ClassImp(AliMagFDM)
31
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//_______________________________________________________________________
67AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
68 Float_t factor, Float_t fmax):
69 AliMagFC(name,title,integ,factor,fmax),
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)
96{
97 //
98 // Standard constructor for the Dipole field
99 //
100 fType = kDipoMap;
101 fMap = 3;
102 SetSolenoidField();
103
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()));
107
108}
109
110//_______________________________________________________________________
111void AliMagFDM::Field(Float_t *xfi, Float_t *b) const
112{
113 //
114 // Main routine to compute the field in a point
115 //
116 const Double_t keps=0.1E-06;
117 const Double_t kPI2=2.*TMath::Pi();
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
130
131 Double_t x[3];
132 Double_t xL3[3];
133 Double_t bint[3];
134
135 Double_t r0;
136 Int_t iKvar,jb;
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 ---
142//
143// This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
144// Transfor correspondingly.
145
146 x[0] = - xfi[0];
147 x[1] = xfi[1];
148 x[2] = - xfi[2];
149 b[0]=b[1]=b[2]=0;
150//
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;
156 rPmax=fRmax;
157 if ( (-700<x[2] && x[2]<=kfZbg &&
158 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
159 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
160 {
161 b[0]=b[1]=0;
162 b[2]=fSolenoid;
163 }
164
165 xL3[0]=x[0]/100;
166 xL3[1]=(x[1]+30)/100;
167 xL3[2]=x[2]/100;
168
169
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;
175 Double_t zCmin,zCmax,yCmin,yCmax;
176
177 zCmin=fZmin;
178 zCmax=fZmax;
179 yCmin=fYmin;
180 yCmax=fYmax;
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)))
182 {
183 if(fMap==3)
184 {
185 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
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];
200 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
201 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
202 zp1= kone-zp2;
203 yyp=xL3[1]- 0.3;
204 cphi=TMath::Abs(yyp/r0);
205 Int_t kcphi=0;
206 if (cphi > kone) {
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));
208 cphi =kone;
209 kcphi=777;
210 }
211 ph0=TMath::ACos(cphi);
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;}
216 if (kcphi==777) {
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));
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;
223
224 Double_t rDel;
225 rDel=fRdel;
226
227 if (r0<= fRdel)
228 {
229
230 if(r0< keps)
231 {
232
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;
236
237 } else {
238
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];
241 alp1= kone - alp2 - alp3;
242
243 for (jb=0; jb<3 ; jb++)
244 {
245 iKvar=jb;
246 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
247 }
248 } // end of keps <=r0
249 }
250 else
251 {
252 rp=r0;
253
254 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
255 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
256 yp1=kone-yp2;
257
258 for (jb=0; jb<3 ; jb++)
259 {
260 iKvar=jb;
261 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
262 }
263 }
264
265 b[0]=-bint[0];
266 b[1]=bint[1];
267 b[2]=-bint[2];
268
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
283 xx1 = fAx1*xL3[2] + fCx1;
284 xx2 = fAx2*xL3[2] + fCx2;
285
286 zzc=xL3[2];
287 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
288 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
289 zz1=kone - zz2;;
290
291 yyc=xL3[1];
292 l0=FZ(yyc, fY , fYdel,lci,fYl);
293 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
294 yy1=kone - yy2;
295 xXl = fXl-kone;
296 dx = (xx2-xx1)/xXl;
297 xxx= xL3[0]-xx1;
298 m0 = int(xxx/dx);
299
300 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
301 {
302 m0=m0+1;
303 AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
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;
312 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
313 }
314
315 b[0]=-bint[0];
316 b[1]=bint[1];
317 b[2]=-bint[2];
318
319 }
320
321
322 } else {
323 AliError(Form("Unknown map of Dipole region %d",fMap));
324 }
325
326} else {
327 ZDCField(xfi,b);
328
329 }
330
331 if(fFactor!=1) {
332 b[0]*=fFactor;
333 b[1]*=fFactor;
334 b[2]*=fFactor;
335 }
336}
337
338//_______________________________________________________________________
339Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar,
340 Float_t delu, Int_t ik, Int_t nk) const
341{
342 //
343 // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
344 //
345 Int_t l,ikj;
346 Double_t ddu,ar;
347
348 Int_t ku,kf;
349 kf=0;
350 ar=Ar[ik];
351 ddu=temp-ar;
352
353 ku=int(ddu/delu);
354 ikj=ik+ku+1;
355 if (ddu<=0) ikj=0;
356
357 for(l=ikj; l<nk; l++)
358 {
359
360 if(temp <= Ar[l])
361 {
362
363 kf=l-1;
364 break;
365 }
366 }
367 return kf;
368 }
369
370//_______________________________________________________________________
371Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2,
372 Double_t alf1, Double_t alf2, Double_t alf3,
373 Int_t kaa, Int_t maa) const
374{
375 //
376 // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
377 //
378 Double_t fa11=0,fa12=0,fa13=0;
379 Double_t fa21=0,fa22=0,fa23=0;
380 Double_t faY1=0,faY2=0;
381 Double_t bba=0;
382
383 switch (kaai) {
384 case 0:
385 fa11 = fBpx[kaa][0][0];
386 fa12 = fBpx[kaa][0][maa];
387 fa13 = fBpx[kaa][0][maa+1];
388 fa21 = fBpx[kaa+1][0][0];
389 fa22 = fBpx[kaa+1][0][maa];
390 fa23 = fBpx[kaa+1][0][maa+1];
391 break;
392 case 1:
393 fa11 = fBpy[kaa][0][0];
394 fa12 = fBpy[kaa][0][maa];
395 fa13 = fBpy[kaa][0][maa+1];
396 fa21 = fBpy[kaa+1][0][0];
397 fa22 = fBpy[kaa+1][0][maa];
398 fa23 = fBpy[kaa+1][0][maa+1];
399 break;
400 case 2:
401 fa11 = fBpz[kaa][0][0];
402 fa12 = fBpz[kaa][0][maa];
403 fa13 = fBpz[kaa][0][maa+1];
404 fa21 = fBpz[kaa+1][0][0];
405 fa22 = fBpz[kaa+1][0][maa];
406 fa23 = fBpz[kaa+1][0][maa+1];
407 break;
408 default:
409 AliFatal(Form("Invalid value of kaai %d",kaai));
410 }
411 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
412 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
413 bba = zaa1*faY1+zaa2*faY2;
414 return bba;
415}
416
417
418//_______________________________________________________________________
419Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2,
420 Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m) const
421{
422 //
423 // Calculation of field componet at a given axis (general case)
424 //
425 Double_t fy1=0,fy2=0,ffy=0;
426 Double_t gy1=0,gy2=0,ggy=0;
427 Double_t bbi=0;
428 Double_t bf11=0,bf12=0,bf21=0,bf22=0;
429 Double_t bg11=0,bg12=0,bg21=0,bg22=0;
430
431
432 /*-----------------Polar part ------------------*/
433
434 switch (kv) {
435 case 0:
436 bf11=fBpx[k][l][m];
437 bf12=fBpx[k+1][l][m];
438 bf21=fBpx[k+1][l+1][m];
439 bf22=fBpx[k][l+1][m];
440
441 bg11=fBpx[k][l][m+1];
442 bg12=fBpx[k+1][l][m+1];
443 bg21=fBpx[k+1][l+1][m+1];
444 bg22=fBpx[k][l+1][m+1];
445 break;
446
447 case 1:
448 bf11=fBpy[k][l][m];
449 bf12=fBpy[k+1][l][m];
450 bf21=fBpy[k+1][l+1][m];
451 bf22=fBpy[k][l+1][m];
452
453 bg11=fBpy[k][l][m+1];
454 bg12=fBpy[k+1][l][m+1];
455 bg21=fBpy[k+1][l+1][m+1];
456 bg22=fBpy[k][l+1][m+1];
457 break;
458
459 case 2:
460 bf11=fBpz[k][l][m];
461 bf12=fBpz[k+1][l][m];
462 bf21=fBpz[k+1][l+1][m];
463 bf22=fBpz[k][l+1][m];
464
465 bg11=fBpz[k][l][m+1];
466 bg12=fBpz[k+1][l][m+1];
467 bg21=fBpz[k+1][l+1][m+1];
468 bg22=fBpz[k][l+1][m+1];
469 break;
470 /*-----------------Cartensian part ---------------*/
471
472 case 3:
473 bf11=fBcx[k][l][m];
474 bf12=fBcx[k+1][l][m];
475 bf21=fBcx[k+1][l+1][m];
476 bf22=fBcx[k][l+1][m];
477
478 bg11=fBcx[k][l][m+1];
479 bg12=fBcx[k+1][l][m+1];
480 bg21=fBcx[k+1][l+1][m+1];
481 bg22=fBcx[k][l+1][m+1];
482 break;
483
484 case 4:
485 bf11=fBcy[k][l][m];
486 bf12=fBcy[k+1][l][m];
487 bf21=fBcy[k+1][l+1][m];
488 bf22=fBcy[k][l+1][m];
489
490 bg11=fBcy[k][l][m+1];
491 bg12=fBcy[k+1][l][m+1];
492 bg21=fBcy[k+1][l+1][m+1];
493 bg22=fBcy[k][l+1][m+1];
494 break;
495
496 case 5:
497 bf11=fBcz[k][l][m];
498 bf12=fBcz[k+1][l][m];
499 bf21=fBcz[k+1][l+1][m];
500 bf22=fBcz[k][l+1][m];
501
502 bg11=fBcz[k][l][m+1];
503 bg12=fBcz[k+1][l][m+1];
504 bg21=fBcz[k+1][l+1][m+1];
505 bg22=fBcz[k][l+1][m+1];
506 break;
507
508 default:
509 AliFatal(Form("Invalid value of kv %d",kv));
510 }
511
512
513
514 fy1=z1*bf11+z2*bf12;
515 fy2=z2*bf21+z1* bf22;
516 ffy=y1*fy1+ y2*fy2;
517
518
519 gy1 = z1*bg11+z2*bg12;
520 gy2 = z2*bg21+z1*bg22;
521 ggy= y1*gy1 + y2*gy2;
522
523 bbi = x1*ffy+x2*ggy;
524
525 return bbi;
526
527}
528
529//_______________________________________________________________________
530void AliMagFDM::ReadField()
531{
532 //
533 // Method to read the magnetic field from file
534 //
535 FILE *magfile;
536
537 Int_t ik, il, im;
538 Float_t zzp, rr,phii;
539 Float_t zz, yy, bx,by,bz,bb;
540
541 char *fname;
542 AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
543 fname = gSystem->ExpandPathName(fTitle.Data());
544 magfile=fopen(fname,"r");
545 delete [] fname;
546
547 AliDebug(2,"Cartensian part");
548
549 if (magfile) {
550
551// Cartensian part
552
553 fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
554
555 AliDebug(3,Form("fYl %d, fXl %d, fZl %d",fYl, fXl, fZl));
556
557 for (ik=0; ik<fZl; ik++)
558 {
559
560 fscanf(magfile, " %e ", &zz);
561 fZc[ik]=zz;
562
563 }
564
565 for (ik=0; ik<fYl; ik++)
566 {
567 fscanf(magfile, " %e ", &yy);
568 fY[ik]=yy;
569
570 }
571 for (ik=0; ik<81; ik++)
572 {
573 AliDebug(4,Form("fZc %e,fY %e", fZc[ik],fY[ik]));
574 }
575
576 fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
577
578AliDebug(3,Form("fYdel %e, fXdel %e, fZdel %e",fYdel,fXdel,fZdel));
579AliDebug(3,Form("fZmax %e, fZmin %e, fYmax %e,fYmin %e",fZmax,fZmin,fYmax,fYmin));
580AliDebug(3,Form("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e",fAx1,fCx1,fAx2,fCx2));
581
582 for (il=0; il<44; il++) {
583 for (im=0; im<81; im++) {
584 for (ik=0; ik<81; ik++) {
585
586 fscanf(magfile, " %e ", &by);
587 fBcy[ik][im][il]=by;
588 }
589 }
590 }
591
592 for (il=0; il<44; il++) {
593 for (im=0; im<81; im++) {
594 for (ik=0; ik<81; ik++) {
595
596 fscanf(magfile, " %e ", &bx);
597 fBcx[ik][im][il]=bx;
598 }
599 }
600 }
601
602 for (il=0; il<44; il++) {
603 for (im=0; im<81; im++) {
604 for (ik=0; ik<81; ik++) {
605
606 fscanf(magfile, " %e ", &bz);
607 fBcz[ik][im][il]=bz;
608 }
609 }
610 }
611//---------------------- Polar part ---------------------------------
612
613 AliDebug(2,"Polar part");
614 fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
615 AliDebug(3,Form("fZpl %d, fRn %d, fPhin %d",fZpl,fRn,fPhin));
616
617 AliDebug(4," fZp array");
618
619 for (ik=0; ik<51; ik++)
620 {
621 fscanf(magfile, " %e ", &zzp);
622 fZp[ik]=zzp;
623 AliDebug(4,Form(" %e",fZp[ik]));
624 }
625
626 AliDebug(4," fR array");
627
628 for (ik=0; ik<10; ik++)
629 {
630 fscanf(magfile, " %e ", &rr);
631 fR[ik]=rr;
632 AliDebug(4,Form(" %e",fR[ik]));
633 }
634
635// AliDebug(4,"fPhi array");
636
637 for (il=0; il<33; il++)
638 {
639 fscanf(magfile, " %e ", &phii);
640 fPhi[il]=phii;
641// AliDebug(4,Form(" %e",fPhi[il]));
642 }
643
644 fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
645
646AliDebug(3,Form("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin));
647
648
649 for (il=0; il<33; il++) {
650 for (im=0; im<10; im++) {
651 for (ik=0; ik<51; ik++) {
652 fscanf(magfile, " %e ", &by);
653 fBpy[ik][im][il]=by;
654 }
655 }
656 }
657
658 for (il=0; il<33; il++) {
659 for (im=0; im<10; im++) {
660 for (ik=0; ik<51; ik++) {
661 fscanf(magfile, " %e ", &bx);
662 fBpx[ik][im][il]=bx;
663 }
664 }
665 }
666
667
668 for (il=0; il<33; il++) {
669 for (im=0; im<10; im++) {
670 for (ik=0; ik<51; ik++) {
671 fscanf(magfile, " %e ", &bz);
672 fBpz[ik][im][il]=bz;
673 }
674 }
675 }
676
677
678 for (il=0; il<32; il++) {
679 for (im=0; im<2; im++) {
680 for (ik=0; ik<2; ik++) {
681 fscanf(magfile, " %e ", &bb);
682 fB[ik][im][il]=bb;
683 }
684 }
685 }
686//
687 } else {
688 AliFatal(Form("File %s not found !",fTitle.Data()));
689 }
690}