]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMagFDM.cxx
replaces AliPHOSAliEnFile: references to AliEn removed
[u/mrichter/AliRoot.git] / STEER / AliMagFDM.cxx
index cae64b0e1de5135308614d167314c886e4e03cb3..5604727cbac68f08273663e47b5fecbd32fec476 100644 (file)
  * 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]);
@@ -108,7 +159,7 @@ void AliMagFDM::Field(Float_t *xfi, Float_t *b)
        || (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;
@@ -127,15 +178,11 @@ void AliMagFDM::Field(Float_t *xfi, Float_t *b)
   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 ----------------------
        
@@ -150,16 +197,29 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
        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;
@@ -183,8 +243,7 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
         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 
        }
@@ -192,23 +251,21 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
        {
         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 
    {
@@ -227,15 +284,17 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
    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)) 
@@ -250,76 +309,23 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
    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) {
@@ -329,75 +335,51 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
   }
 }
 
-//_________________________________________
-
-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];
@@ -424,45 +406,29 @@ void AliMagFDM::FRfuncBi(Int_t *kai,Double_t *za1, Double_t *za2, Double_t *al1,
     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 ------------------*/
   
@@ -541,7 +507,7 @@ void AliMagFDM::FGfuncBi(Double_t *zz1,Double_t *zz2, Double_t *yy1,Double_t *yy
     break;
 
   default:
-    Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
+    Fatal("Bb","Invalid value of kv %d\n",kv);
     exit(1);
   }  
   
@@ -557,12 +523,12 @@ void AliMagFDM::FGfuncBi(Double_t *zz1,Double_t *zz2, Double_t *yy1,Double_t *yy
   ggy= y1*gy1 +  y2*gy2;
   
   bbi =  x1*ffy+x2*ggy;
-  
-  *bb=bbi; 
+
+  return bbi;
   
 }     
-//____________________________________________
 
+//_______________________________________________________________________
 void AliMagFDM::ReadField()
 {
   //
@@ -725,4 +691,3 @@ printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n",
     exit(1);
   }
 }
-//________________________________