1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.11 2001/02/08 13:18:00 hristov
19 Print removed (J.Gosset)
21 Revision 1.10 2001/01/18 13:21:30 morsch
24 Revision 1.9 2001/01/17 20:02:20 morsch
25 In the AliMagFDM tree call-by-reference functions were changed to
26 call-by-value, what is more adequate for our task. There were added
27 a few comments and put protection to values of cos > 1.000 in
28 AliMagFDM.cxx. (Galina Chabratova)
30 Revision 1.8 2000/12/18 10:44:01 morsch
31 Possibility to set field map by passing pointer to objet of type AliMagF via
34 gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
36 Revision 1.7 2000/12/01 11:20:27 alibrary
37 Corrector dipole removed from ZDC
39 Revision 1.6 2000/11/10 18:09:55 fca
40 New field map for the ZDC
42 Revision 1.5 2000/10/27 14:17:04 morsch
43 - Bug causing segmentation violation during muon reconstruction corrected
44 - Coding rule violations corrected.
47 Revision 1.4 2000/10/02 21:28:14 fca
48 Removal of useless dependecies via forward declarations
50 Revision 1.3 2000/07/13 16:19:09 fca
51 Mainly coding conventions + some small bug fixes
53 Revision 1.2 2000/07/12 08:56:25 fca
54 Coding convention correction and warning removal
56 Revision 1.1 2000/07/11 18:24:59 fca
57 Coding convention corrections + few minor bug fixes
63 #include "AliMagFDM.h"
69 //________________________________________
70 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
71 const Float_t factor, const Float_t fmax)
72 : AliMagF(name,title,integ,factor,fmax)
75 // Standard constructor for the Dipole field
81 printf("Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s\n",fName.Data(), fMap ,integ,factor,fTitle.Data());
85 //________________________________________
87 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
90 // Main routine to compute the field in a point
92 static const Double_t keps=0.1E-06;
93 static const Double_t PI2=2.*TMath::Pi();
94 static const Double_t kone=1;
96 static const Int_t kiip=33;
97 static const Int_t kmiip=0;
98 static const Int_t kliip=0;
100 static const Int_t kiic=0;
101 static const Int_t kmiic=0;
102 static const Int_t kliic=0;
104 static const Double_t kfZbg=502.92; // Start of Map using in z
105 static const Double_t kfZL3=600; // Beginning of L3 door in z
114 Double_t zp1, zp2,xp1,xp2,yp1,yp2;
115 Double_t zz1, zz2,yy1,yy2,x2,x1;
117 // --- start the map fiel from z = 502.92 cm ---
123 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
125 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
129 if ( (-700<x[2] && x[2]<=kfZbg &&
130 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
131 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
138 xL3[1]=(x[1]+30)/100;
142 if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
143 if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
145 Double_t xminn=xL3[2]*fAx1+fCx1;
146 Double_t xmaxx=xL3[2]*fAx2+fCx2;
147 Double_t zCmin,zCmax,yCmin,yCmax;
153 if ((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)))
157 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
159 //--------------------- Polar part ----------------------
164 Double_t alp1,alp2,alp3;
165 Double_t zpz,rp,fip,cphi;
172 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
173 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
176 cphi=TMath::Abs(yyp/r0);
179 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);
183 ph0=TMath::ACos(cphi);
184 if (xL3[0] < 0 && yyp > 0 ) {ph0=PI2/2 - ph0;}
185 if (xL3[0] < 0 && yyp < 0 ) {ph0=PI2/2 + ph0;}
186 if (xL3[0] > 0 && yyp < 0 ) {ph0=PI2 - ph0;}
187 if (ph0 > PI2) { ph0=ph0 - PI2;}
189 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);
192 mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
193 xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
205 bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
206 bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
207 bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
211 alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
212 alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
213 alp1= kone - alp2 - alp3;
215 for (jb=0; jb<3 ; jb++)
218 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
220 } // end of keps <=r0
226 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
227 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
230 for (jb=0; jb<3 ; jb++)
233 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
244 //-------------- Cartensian part ------------------
248 Double_t xx1, xx2,dx, xxx ,xXl;
255 xx1 = fAx1*xL3[2] + fCx1;
256 xx2 = fAx2*xL3[2] + fCx2;
259 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
260 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
264 l0=FZ(yyc, fY , fYdel,lci,fYl);
265 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
272 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
275 printf(" m0 %d, m0+1 %d\n",m0,m0+1);
278 x2=(xL3[0]-( xx1+m0*dx))/dx;
281 for (jb=3; jb<6; jb++)
284 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
295 printf("Unknown map of Dipole region %d\n",fMap);
300 //This is the ZDC part
301 Float_t rad2=x[0]*x[0]+x[1]*x[1];
302 if(x[2]>kCORBEG2 && x[2]<kCOREND2){
307 else if(x[2]>kZ1BEG && x[2]<kZ1END){
313 else if(x[2]>kZ2BEG && x[2]<kZ2END){
319 else if(x[2]>kZ3BEG && x[2]<kZ3END){
325 else if(x[2]>kZ4BEG && x[2]<kZ4END){
331 else if(x[2]>kD1BEG && x[2]<kD1END){
336 else if(x[2]>kD2BEG && x[2]<kD2END){
337 if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
338 || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
352 //_____________________ FZ ____________________
354 Int_t AliMagFDM::FZ(Double_t temp, Float_t *Ar, Float_t delu,Int_t ik,Int_t nk)
357 // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
371 for(l=ikj; l<nk; l++)
384 /*---------------------Ba------------------*/
386 Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2, Double_t alf1, Double_t alf2, Double_t alf3, Int_t kaa, Int_t maa)
389 // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
391 Double_t fa11,fa12,fa13;
392 Double_t fa21,fa22,fa23;
398 fa11 = fBpx[kaa][0][0];
399 fa12 = fBpx[kaa][0][maa];
400 fa13 = fBpx[kaa][0][maa+1];
401 fa21 = fBpx[kaa+1][0][0];
402 fa22 = fBpx[kaa+1][0][maa];
403 fa23 = fBpx[kaa+1][0][maa+1];
406 fa11 = fBpy[kaa][0][0];
407 fa12 = fBpy[kaa][0][maa];
408 fa13 = fBpy[kaa][0][maa+1];
409 fa21 = fBpy[kaa+1][0][0];
410 fa22 = fBpy[kaa+1][0][maa];
411 fa23 = fBpy[kaa+1][0][maa+1];
414 fa11 = fBpz[kaa][0][0];
415 fa12 = fBpz[kaa][0][maa];
416 fa13 = fBpz[kaa][0][maa+1];
417 fa21 = fBpz[kaa+1][0][0];
418 fa22 = fBpz[kaa+1][0][maa];
419 fa23 = fBpz[kaa+1][0][maa+1];
422 Fatal("Ba","Invalid value of kaai %d\n",kaai);
425 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
426 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
427 bba = zaa1*faY1+zaa2*faY2;
432 /*------------------------Bb--------------------------*/
434 Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2, Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m)
437 // Calculation of field componet at a given axis (general case)
439 Double_t fy1, fy2, ffy;
440 Double_t gy1,gy2,ggy;
442 Double_t bf11,bf12,bf21,bf22;
443 Double_t bg11,bg12,bg21,bg22;
446 /*-----------------Polar part ------------------*/
451 bf12=fBpx[k+1][l][m];
452 bf21=fBpx[k+1][l+1][m];
453 bf22=fBpx[k][l+1][m];
455 bg11=fBpx[k][l][m+1];
456 bg12=fBpx[k+1][l][m+1];
457 bg21=fBpx[k+1][l+1][m+1];
458 bg22=fBpx[k][l+1][m+1];
463 bf12=fBpy[k+1][l][m];
464 bf21=fBpy[k+1][l+1][m];
465 bf22=fBpy[k][l+1][m];
467 bg11=fBpy[k][l][m+1];
468 bg12=fBpy[k+1][l][m+1];
469 bg21=fBpy[k+1][l+1][m+1];
470 bg22=fBpy[k][l+1][m+1];
475 bf12=fBpz[k+1][l][m];
476 bf21=fBpz[k+1][l+1][m];
477 bf22=fBpz[k][l+1][m];
479 bg11=fBpz[k][l][m+1];
480 bg12=fBpz[k+1][l][m+1];
481 bg21=fBpz[k+1][l+1][m+1];
482 bg22=fBpz[k][l+1][m+1];
484 /*-----------------Cartensian part ---------------*/
488 bf12=fBcx[k+1][l][m];
489 bf21=fBcx[k+1][l+1][m];
490 bf22=fBcx[k][l+1][m];
492 bg11=fBcx[k][l][m+1];
493 bg12=fBcx[k+1][l][m+1];
494 bg21=fBcx[k+1][l+1][m+1];
495 bg22=fBcx[k][l+1][m+1];
500 bf12=fBcy[k+1][l][m];
501 bf21=fBcy[k+1][l+1][m];
502 bf22=fBcy[k][l+1][m];
504 bg11=fBcy[k][l][m+1];
505 bg12=fBcy[k+1][l][m+1];
506 bg21=fBcy[k+1][l+1][m+1];
507 bg22=fBcy[k][l+1][m+1];
512 bf12=fBcz[k+1][l][m];
513 bf21=fBcz[k+1][l+1][m];
514 bf22=fBcz[k][l+1][m];
516 bg11=fBcz[k][l][m+1];
517 bg12=fBcz[k+1][l][m+1];
518 bg21=fBcz[k+1][l+1][m+1];
519 bg22=fBcz[k][l+1][m+1];
523 Fatal("Bb","Invalid value of kv %d\n",kv);
530 fy2=z2*bf21+z1* bf22;
534 gy1 = z1*bg11+z2*bg12;
535 gy2 = z2*bg21+z1*bg22;
536 ggy= y1*gy1 + y2*gy2;
543 //____________________________________________
545 void AliMagFDM::ReadField()
548 // Method to read the magnetic field from file
553 Float_t zzp, rr,phii;
554 Float_t zz, yy, bx,by,bz,bb;
557 printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
558 fname = gSystem->ExpandPathName(fTitle.Data());
559 magfile=fopen(fname,"r");
562 printf("Cartensian part\n");
568 fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
570 printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
572 for (ik=0; ik<fZl; ik++)
575 fscanf(magfile, " %e ", &zz);
580 for (ik=0; ik<fYl; ik++)
582 fscanf(magfile, " %e ", &yy);
586 for (ik=0; ik<81; ik++)
588 printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
591 fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
593 printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
594 printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
595 printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
597 for (il=0; il<44; il++) {
598 for (im=0; im<81; im++) {
599 for (ik=0; ik<81; ik++) {
601 fscanf(magfile, " %e ", &by);
607 for (il=0; il<44; il++) {
608 for (im=0; im<81; im++) {
609 for (ik=0; ik<81; ik++) {
611 fscanf(magfile, " %e ", &bx);
617 for (il=0; il<44; il++) {
618 for (im=0; im<81; im++) {
619 for (ik=0; ik<81; ik++) {
621 fscanf(magfile, " %e ", &bz);
626 //---------------------- Polar part ---------------------------------
628 printf("Polar part\n");
629 fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
630 printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
632 printf(" fZp array\n");
634 for (ik=0; ik<51; ik++)
636 fscanf(magfile, " %e ", &zzp);
638 printf(" %e\n",fZp[ik]);
641 printf(" fR array\n");
643 for (ik=0; ik<10; ik++)
645 fscanf(magfile, " %e ", &rr);
647 printf(" %e\n",fR[ik]);
650 // printf("fPhi array\n");
652 for (il=0; il<33; il++)
654 fscanf(magfile, " %e ", &phii);
656 // printf(" %e\n",fPhi[il]);
659 fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
661 printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
664 for (il=0; il<33; il++) {
665 for (im=0; im<10; im++) {
666 for (ik=0; ik<51; ik++) {
667 fscanf(magfile, " %e ", &by);
673 for (il=0; il<33; il++) {
674 for (im=0; im<10; im++) {
675 for (ik=0; ik<51; ik++) {
676 fscanf(magfile, " %e ", &bx);
683 for (il=0; il<33; il++) {
684 for (im=0; im<10; im++) {
685 for (ik=0; ik<51; ik++) {
686 fscanf(magfile, " %e ", &bz);
693 for (il=0; il<32; il++) {
694 for (im=0; im<2; im++) {
695 for (ik=0; ik<2; ik++) {
696 fscanf(magfile, " %e ", &bb);
703 printf("File %s not found !\n",fTitle.Data());
707 //________________________________