* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.7 2000/12/01 11:20:27 alibrary
-Corrector dipole removed from ZDC
+/* $Id$ */
-Revision 1.6 2000/11/10 18:09:55 fca
-New field map for the ZDC
+//-------------------------------------------------------------------------
+// Field with Magnetic Field map
+// Used by AliRun class
+// Author:
+//-------------------------------------------------------------------------
-Revision 1.5 2000/10/27 14:17:04 morsch
-- Bug causing segmentation violation during muon reconstruction corrected
-- Coding rule violations corrected.
-(Galina Chabratova)
-
-Revision 1.4 2000/10/02 21:28:14 fca
-Removal of useless dependecies via forward declarations
-
-Revision 1.3 2000/07/13 16:19:09 fca
-Mainly coding conventions + some small bug fixes
-
-Revision 1.2 2000/07/12 08:56:25 fca
-Coding convention correction and warning removal
-
-Revision 1.1 2000/07/11 18:24:59 fca
-Coding convention corrections + few minor bug fixes
-
-*/
-
-#include <stdlib.h>
+#include "TSystem.h"
+#include "AliLog.h"
#include "AliMagFDM.h"
-#include "TSystem.h"
-
ClassImp(AliMagFDM)
-//________________________________________
-AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
- const Float_t factor, const Float_t fmax)
- : AliMagF(name,title,integ,factor,fmax)
-
+//_______________________________________________________________________
+AliMagFDM::AliMagFDM():
+ fSolenoid(0),
+ fInd(0),
+ fZmin(0),
+ fZmax(0),
+ fYmax(0),
+ fYmin(0),
+ fZpmx(0),
+ fZpmn(0),
+ fRmax(0),
+ fRmin(0),
+ fXdel(0),
+ fYdel(0),
+ fZdel(0),
+ fRdel(0),
+ fPhid(0),
+ fZpdl(0),
+ fCx1(0),
+ fCx2(0),
+ fAx1(0),
+ fAx2(0),
+ fXl(0),
+ fYl(0),
+ fZl(0),
+ fRn(0),
+ fPhin(0),
+ fZpl(0)
+{
+ //
+ // Default constructor for the Dipole field
+ //
+}
+
+//_______________________________________________________________________
+AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
+ Float_t factor, Float_t fmax):
+ AliMagFC(name,title,integ,factor,fmax),
+ fSolenoid(0),
+ fInd(0),
+ fZmin(0),
+ fZmax(0),
+ fYmax(0),
+ fYmin(0),
+ fZpmx(0),
+ fZpmn(0),
+ fRmax(0),
+ fRmin(0),
+ fXdel(0),
+ fYdel(0),
+ fZdel(0),
+ fRdel(0),
+ fPhid(0),
+ fZpdl(0),
+ fCx1(0),
+ fCx2(0),
+ fAx1(0),
+ fAx2(0),
+ fXl(0),
+ fYl(0),
+ fZl(0),
+ fRn(0),
+ fPhin(0),
+ fZpl(0)
{
//
// Standard constructor for the Dipole field
//
fType = kDipoMap;
fMap = 3;
+ SetSolenoidField();
- printf("Field Map for Muon Arm from IP till muon filter %s created: map= %d, factor= %f, file=%s\n",fName.Data(), fMap ,factor,fTitle.Data());
-
+ AliDebug(1, Form(
+ "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s",
+ fName.Data(), fMap ,integ,factor,fTitle.Data()));
+
}
-//________________________________________
-
-void AliMagFDM::Field(Float_t *xfi, Float_t *b)
+//_______________________________________________________________________
+void AliMagFDM::Field(Float_t *xfi, Float_t *b) const
{
//
// Main routine to compute the field in a point
//
- static const Double_t keps=0.1E-06;
- static const Double_t kpi2=.6283185E+01;
- static const Double_t kone=1;
-
- static const Int_t kiip=33;
- static const Int_t kmiip=0;
- static const Int_t kliip=0;
-
- static const Int_t kiic=0;
- static const Int_t kmiic=0;
- static const Int_t kliic=0;
-
- static const Double_t kfZbg=502.92; // Start of Map using in z
- static const Double_t kfZL3=600; // Beginning of L3 door in z
+ const Double_t keps=0.1E-06;
+ const Double_t kPI2=2.*TMath::Pi();
+ const Double_t kone=1;
+
+ const Int_t kiip=33;
+ const Int_t kmiip=0;
+ const Int_t kliip=0;
+
+ const Int_t kiic=0;
+ const Int_t kmiic=0;
+ const Int_t kliic=0;
+
+ const Double_t kfZbg=502.92; // Start of Map using in z
+ const Double_t kfZL3=600; // Beginning of L3 door in z
Double_t x[3];
Double_t xL3[3];
Double_t bint[3];
Double_t r0;
-
- Double_t bbj;
Int_t iKvar,jb;
Double_t zp1, zp2,xp1,xp2,yp1,yp2;
Double_t zz1, zz2,yy1,yy2,x2,x1;
// --- start the map fiel from z = 502.92 cm ---
+//
+// This map has been calculated in a coordinate system in which the muon spectrometer sits at z > 0
+// Transfor correspondingly.
- x[0] = xfi[0];
- x[1] = xfi[1];
- x[2] = xfi[2];
-
- /*
- if (x[0]==0) x[0]=0.1E-06;
- if (x[1]==0) x[1]=0.1E-06;
- */
- b[0]=b[1]=b[2]=0;
+ x[0] = - xfi[0];
+ x[1] = xfi[1];
+ x[2] = - xfi[2];
+ b[0]=b[1]=b[2]=0;
+//
// printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
|| (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
{
b[0]=b[1]=0;
- b[2]=2;
+ b[2]=fSolenoid;
}
xL3[0]=x[0]/100;
zCmax=fZmax;
yCmin=fYmin;
yCmax=fYmax;
- /*
-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)))
- */
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)))
{
if(fMap==3)
{
if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
- // if (xL3[2]<zCmin && r0<rPmax)
{
//--------------------- Polar part ----------------------
mpi=kmiip;
zpz=xL3[2];
-
- FZ(&zpz, fZp ,&fZpdl,&kpi,&kp0,&zp1 ,&zp2,&fZpl) ;
-
+ kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
+ zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
+ zp1= kone-zp2;
yyp=xL3[1]- 0.3;
- cphi=yyp/r0;
+ cphi=TMath::Abs(yyp/r0);
+ Int_t kcphi=0;
+ if (cphi > kone) {
+ 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));
+ cphi =kone;
+ kcphi=777;
+ }
ph0=TMath::ACos(cphi);
- if (xL3[0]< 0) {ph0=kpi2 - ph0;}
-
- fip=ph0;
- FZ(&fip,fPhi,&fPhid ,&mpi,&mp0, &xp1,&xp2,&fPhin);
+ if (xL3[0] < 0 && yyp > 0 ) {ph0=kPI2/2 - ph0;}
+ if (xL3[0] < 0 && yyp < 0 ) {ph0=kPI2/2 + ph0;}
+ if (xL3[0] > 0 && yyp < 0 ) {ph0=kPI2 - ph0;}
+ if (ph0 > kPI2) { ph0=ph0 - kPI2;}
+ if (kcphi==777) {
+ 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));
+ }
+ fip=ph0;
+ mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
+ xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
+ xp1=kone-xp2;
Double_t rDel;
rDel=fRdel;
for (jb=0; jb<3 ; jb++)
{
iKvar=jb;
- FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
- bint[jb] = bbj*10 ;
+ bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
}
} // end of keps <=r0
}
{
rp=r0;
- FZ(&rp,fR ,&fRdel,&lpi,&lp0,&yp1,&yp2,&fRn);
+ lp0=FZ(rp,fR ,fRdel,lpi,fRn);
+ yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
+ yp1=kone-yp2;
for (jb=0; jb<3 ; jb++)
{
iKvar=jb;
- FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj);
-
- bint[jb] = bbj*10 ;
+ bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
}
}
- b[0]=bint[0];
+ b[0]=-bint[0];
b[1]=bint[1];
- b[2]=bint[2];
+ b[2]=-bint[2];
-// fprintf(fitest,"------------- Freg2 run -------------\n");
-
}
else
{
xx2 = fAx2*xL3[2] + fCx2;
zzc=xL3[2];
- FZ(&zzc, fZc ,&fZdel, &kci,&k0, &zz1, &zz2, &fZl);
+ k0=FZ(zzc, fZc ,fZdel, kci, fZl);
+ zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
+ zz1=kone - zz2;;
yyc=xL3[1];
- FZ(&yyc, fY , &fYdel,&lci, &l0, &yy1, &yy2,&fYl);
-
+ l0=FZ(yyc, fY , fYdel,lci,fYl);
+ yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
+ yy1=kone - yy2;
xXl = fXl-kone;
dx = (xx2-xx1)/xXl;
xxx= xL3[0]-xx1;
- // xm = xxx/dx;
m0 = int(xxx/dx);
if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
{
m0=m0+1;
- printf(" m0 %d, m0+1 %d\n",m0,m0+1);
+ AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));
}
x2=(xL3[0]-( xx1+m0*dx))/dx;
for (jb=3; jb<6; jb++)
{
iKvar=jb;
- FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj);
- bint[jb-3] = bbj*10 ;
+ bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
}
- b[0]=bint[0];
+ b[0]=-bint[0];
b[1]=bint[1];
- b[2]=bint[2];
+ b[2]=-bint[2];
-// fprintf(fitest,"------------ Freg1 run -----------------\n");
}
+
} else {
- printf("Unknown map of Dipole region %d\n",fMap);
+ AliError(Form("Unknown map of Dipole region %d",fMap));
}
} else {
+ ZDCField(xfi,b);
-//This is the ZDC part
- Float_t rad2=x[0]*x[0]+x[1]*x[1];
- if(x[2]>kCORBEG2 && x[2]<kCOREND2){
- if(rad2<kCOR2RA2){
- b[0] = kFCORN2;
- }
}
- else if(x[2]>kZ1BEG && x[2]<kZ1END){
- if(rad2<kZ1RA2){
- b[0] = -kG1*x[1];
- b[1] = -kG1*x[0];
- }
- }
- else if(x[2]>kZ2BEG && x[2]<kZ2END){
- if(rad2<kZ2RA2){
- b[0] = kG1*x[1];
- b[1] = kG1*x[0];
- }
- }
- else if(x[2]>kZ3BEG && x[2]<kZ3END){
- if(rad2<kZ3RA2){
- b[0] = kG1*x[1];
- b[1] = kG1*x[0];
- }
- }
- else if(x[2]>kZ4BEG && x[2]<kZ4END){
- if(rad2<kZ4RA2){
- b[0] = -kG1*x[1];
- b[1] = -kG1*x[0];
- }
- }
- else if(x[2]>kD1BEG && x[2]<kD1END){
- if(rad2<kD1RA2){
- b[1] = -kFDIP;
- }
- }
- else if(x[2]>kD2BEG && x[2]<kD2END){
- if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
- || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
- b[1] = kFDIP;
- }
- }
- }
-
if(fFactor!=1) {
b[0]*=fFactor;
}
}
-//_________________________________________
-
-void 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)
+//_______________________________________________________________________
+Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar,
+ Float_t delu, Int_t ik, Int_t nk) const
{
//
- // Z component of the field
+ // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
//
- static const Double_t kone=1;
- Int_t l,ik,ikj;
- Double_t temp;
- Double_t ddu,delu,ar;
-
- Int_t nk,ku;
- temp=*u;
- nk=*nu;
- ik=*ki;
- delu=*du;
+ Int_t l,ikj;
+ Double_t ddu,ar;
+ Int_t ku,kf;
+ kf=0;
ar=Ar[ik];
ddu=temp-ar;
ku=int(ddu/delu);
- // ikj=ik+ku; //falsh
- ikj=ik+ku+1; //correct
+ ikj=ik+ku+1;
if (ddu<=0) ikj=0;
for(l=ikj; l<nk; l++)
{
- // if(temp < Ar[l])
+
if(temp <= Ar[l])
{
- *kf=l-1;
- *a2=(temp-Ar[l-1])/(Ar[l]-Ar[l-1]);
-
- /*
- *kf=l;
- *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
- */
- break;
+ kf=l-1;
+ break;
}
}
- // *a2=(temp-Ar[*kf])/(Ar[*kf+1]-Ar[*kf]);
- *a1= kone - *a2;
+ return kf;
}
-/*-------------FRfuncBi----------------*/
-
-void 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)
+//_______________________________________________________________________
+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) const
{
//
- // This method needs to be commented
+ // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
//
- Double_t fa11,fa12,fa13;
- Double_t fa21,fa22,fa23;
- Double_t faY1,faY2;
- Double_t bba;
-
- Double_t zaa1,zaa2,alf1,alf2,alf3;
- Int_t kaai,kaa,maa;
- kaai=*kai;
- kaa=*ka;
- maa=*ma;
- zaa1=*za1;
- zaa2=*za2;
- alf1=*al1;
- alf2=*al2;
- alf3=*al3;
-
+ Double_t fa11=0,fa12=0,fa13=0;
+ Double_t fa21=0,fa22=0,fa23=0;
+ Double_t faY1=0,faY2=0;
+ Double_t bba=0;
+
switch (kaai) {
case 0:
fa11 = fBpx[kaa][0][0];
fa23 = fBpz[kaa+1][0][maa+1];
break;
default:
- Fatal("FRfuncBi","Invalid value of kaai %d\n",kaai);
- exit(1);
+ AliFatal(Form("Invalid value of kaai %d",kaai));
}
faY1=alf1*fa11+alf2*fa12+alf3*fa13;
faY2=alf1*fa21+alf2*fa22+alf3*fa23;
bba = zaa1*faY1+zaa2*faY2;
- *ba=bba;
-
+ return bba;
}
-/*----------- FGfuncBi------------*/
-
-void 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)
+//_______________________________________________________________________
+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) const
{
//
- // This method needs to be commented
+ // Calculation of field componet at a given axis (general case)
//
- Double_t fy1, fy2, ffy;
- Double_t gy1,gy2,ggy;
- Double_t z1,z2,y1,y2,x1,x2;
-
- Int_t k,l,m,kv;
- Double_t bbi;
-
- Double_t bf11,bf12,bf21,bf22;
- Double_t bg11,bg12,bg21,bg22;
- k=*kk;
- l=*ll;
- m=*mm;
-
- kv=*kvr;
-
- z1=*zz1;
- z2=*zz2;
- y1=*yy1;
- y2=*yy2;
- x1=*xx1;
- x2=*xx2;
+ Double_t fy1=0,fy2=0,ffy=0;
+ Double_t gy1=0,gy2=0,ggy=0;
+ Double_t bbi=0;
+ Double_t bf11=0,bf12=0,bf21=0,bf22=0;
+ Double_t bg11=0,bg12=0,bg21=0,bg22=0;
+
/*-----------------Polar part ------------------*/
break;
default:
- Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
- exit(1);
+ AliFatal(Form("Invalid value of kv %d",kv));
}
ggy= y1*gy1 + y2*gy2;
bbi = x1*ffy+x2*ggy;
-
- *bb=bbi;
+
+ return bbi;
}
-//____________________________________________
+//_______________________________________________________________________
void AliMagFDM::ReadField()
{
//
Float_t zz, yy, bx,by,bz,bb;
char *fname;
- printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
+ AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
fname = gSystem->ExpandPathName(fTitle.Data());
magfile=fopen(fname,"r");
delete [] fname;
- printf("Cartensian part\n");
+ AliDebug(2,"Cartensian part");
if (magfile) {
fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
- printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
+ AliDebug(3,Form("fYl %d, fXl %d, fZl %d",fYl, fXl, fZl));
for (ik=0; ik<fZl; ik++)
{
}
for (ik=0; ik<81; ik++)
{
- printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
+ AliDebug(4,Form("fZc %e,fY %e", fZc[ik],fY[ik]));
}
fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
-printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
-printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
-printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
+AliDebug(3,Form("fYdel %e, fXdel %e, fZdel %e",fYdel,fXdel,fZdel));
+AliDebug(3,Form("fZmax %e, fZmin %e, fYmax %e,fYmin %e",fZmax,fZmin,fYmax,fYmin));
+AliDebug(3,Form("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e",fAx1,fCx1,fAx2,fCx2));
for (il=0; il<44; il++) {
for (im=0; im<81; im++) {
}
//---------------------- Polar part ---------------------------------
- printf("Polar part\n");
+ AliDebug(2,"Polar part");
fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
- printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
+ AliDebug(3,Form("fZpl %d, fRn %d, fPhin %d",fZpl,fRn,fPhin));
- printf(" fZp array\n");
+ AliDebug(4," fZp array");
for (ik=0; ik<51; ik++)
{
fscanf(magfile, " %e ", &zzp);
fZp[ik]=zzp;
- printf(" %e\n",fZp[ik]);
+ AliDebug(4,Form(" %e",fZp[ik]));
}
- printf(" fR array\n");
+ AliDebug(4," fR array");
for (ik=0; ik<10; ik++)
{
fscanf(magfile, " %e ", &rr);
fR[ik]=rr;
- printf(" %e\n",fR[ik]);
+ AliDebug(4,Form(" %e",fR[ik]));
}
-// printf("fPhi array\n");
+// AliDebug(4,"fPhi array");
for (il=0; il<33; il++)
{
fscanf(magfile, " %e ", &phii);
fPhi[il]=phii;
-// printf(" %e\n",fPhi[il]);
+// AliDebug(4,Form(" %e",fPhi[il]));
}
fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
-printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
+AliDebug(3,Form("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin));
for (il=0; il<33; il++) {
}
//
} else {
- printf("File %s not found !\n",fTitle.Data());
- exit(1);
+ AliFatal(Form("File %s not found !",fTitle.Data()));
}
}
-//________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff