* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.4 2000/10/02 21:28:14 fca
-Removal of useless dependecies via forward declarations
+/* $Id$ */
-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
-
-*/
+//-------------------------------------------------------------------------
+// Field with Magnetic Field map
+// Used by AliRun class
+// Author:
+//-------------------------------------------------------------------------
#include <stdlib.h>
-#include "AliMagFDM.h"
#include "TSystem.h"
-
+
+#include "AliMagFDM.h"
ClassImp(AliMagFDM)
-//________________________________________
-AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
-const Int_t map, const Float_t factor, const Float_t fmax)
- : AliMagF(name,title,integ,map,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;
-
- 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());
+ fMap = 3;
+ SetSolenoidField();
+ Info("ctor",
+ "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());
+
}
-//________________________________________
-
+//_______________________________________________________________________
void AliMagFDM::Field(Float_t *xfi, Float_t *b)
{
//
// 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) {
+ 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);
+ 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) {
+ 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);
+ }
+ 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))
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);
}
} else {
+ ZDCField(xfi,b);
-//This is the ZDC part
- Float_t rad2=x[0]*x[0]+x[1]*x[1];
- if(rad2<kD2RA2) {
- if(x[2]>kD2BEG) {
-
-// Separator Dipole D2
- if(x[2]<kD2END) b[1]=kFDIP;
- } else if(x[2]>kD1BEG) {
-
-// Separator Dipole D1
- if(x[2]<kD1END) b[1]=-kFDIP;
- }
- if(rad2<kCORRA2) {
-
-// First quadrupole of inner triplet de-focussing in x-direction
-// Inner triplet
- if(x[2]>kZ4BEG) {
- if(x[2]<kZ4END) {
-
-// 2430 <-> 3060
- b[0]=-kG1*x[1];
- b[1]=-kG1*x[0];
- }
- } else if(x[2]>kZ3BEG) {
- if(x[2]<kZ3END) {
-
-// 1530 <-> 2080
- b[0]=kG1*x[1];
- b[1]=kG1*x[0];
- }
- } else if(x[2]>kZ2BEG) {
- if(x[2]<kZ2END) {
-
-// 890 <-> 1430
- b[0]=kG1*x[1];
- b[1]=kG1*x[0];
- }
- } else if(x[2]>kZ1BEG) {
- if(x[2]<kZ1END) {
-
-// 0 <-> 630
- b[0]=-kG1*x[1];
- b[1]=-kG1*x[0];
- }
- } else if(x[2]>kCORBEG) {
- if(x[2]<kCOREND) {
-// Corrector dipole (because of dimuon arm)
-// b[0]=kFCORN;
- b[0]=-kFCORN;
- }
- }
- }
- }
}
if(fFactor!=1) {
}
}
-//_________________________________________
-
-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, 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)
{
//
- // 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;
-
+
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);
+ Fatal("Ba","Invalid value of kaai %d\n",kaai);
exit(1);
}
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)
{
//
- // 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;
+
/*-----------------Polar part ------------------*/
break;
default:
- Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
+ Fatal("Bb","Invalid value of kv %d\n",kv);
exit(1);
}
ggy= y1*gy1 + y2*gy2;
bbi = x1*ffy+x2*ggy;
-
- *bb=bbi;
+
+ return bbi;
}
-//____________________________________________
+//_______________________________________________________________________
void AliMagFDM::ReadField()
{
//
exit(1);
}
}
-//________________________________