Overloaded virtual void AliMagXXX:Field(float*,float*) with
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Dec 2008 15:42:04 +0000 (15:42 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Dec 2008 15:42:04 +0000 (15:42 +0000)
virtual void Field(double*,double*) version.
Some auxiliary functions with (float*) arguments are converted
to templates.

21 files changed:
STEER/AliCheb3D.h
STEER/AliCheb3DCalc.cxx
STEER/AliCheb3DCalc.h
STEER/AliFieldMap.cxx
STEER/AliFieldMap.h
STEER/AliMagF.cxx
STEER/AliMagF.h
STEER/AliMagFC.cxx
STEER/AliMagFC.h
STEER/AliMagFCM.cxx
STEER/AliMagFCM.h
STEER/AliMagFCheb.cxx
STEER/AliMagFCheb.h
STEER/AliMagFDM.cxx
STEER/AliMagFDM.h
STEER/AliMagFMaps.cxx
STEER/AliMagFMaps.h
STEER/AliMagFMapsV1.cxx
STEER/AliMagFMapsV1.h
STEER/AliMagWrapCheb.cxx
STEER/AliMagWrapCheb.h

index a3d5c09..f147129 100644 (file)
@@ -97,8 +97,10 @@ class AliCheb3D: public TNamed
   ~AliCheb3D()                                                                 {Clear();}
   //
   AliCheb3D&   operator=(const AliCheb3D& rhs);
-  void         Eval(const Float_t  *par,Float_t  *res);
-  Float_t      Eval(const Float_t  *par,int idim);
+  template <class T>
+    void       Eval(const T  *par, T *res);
+  template <class T>
+    T          Eval(const T  *par,int idim);
   //
   void         EvalDeriv(int dimd, const Float_t  *par, Float_t  *res);
   void         EvalDeriv2(int dimd1, int dimd2, const Float_t  *par,Float_t  *res);
@@ -107,8 +109,9 @@ class AliCheb3D: public TNamed
   void         EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3]); 
   void         EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3]); 
   void         Print(const Option_t* opt="")                             const;
-  Bool_t       IsInside(const Float_t  *par)                             const;
-  Bool_t       IsInside(const Double_t *par)                             const;
+  template <class T>
+    Bool_t       IsInside(const T  *par)                                 const;
+  //
   AliCheb3DCalc*  GetChebCalc(int i)                                     const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
   Float_t      GetBoundMin(int i)                                        const {return fBMin[i];}
   Float_t      GetBoundMax(int i)                                        const {return fBMax[i];}
@@ -145,8 +148,10 @@ class AliCheb3D: public TNamed
   Int_t        ChebFit(int dmOut);
 #endif
   //
-  Float_t      MapToInternal(Float_t  x,Int_t d) const; // map x to [-1:1]
-  Float_t      MapToExternal(Float_t  x,Int_t d) const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
+  template <class T>
+    T          MapToInternal(T  x,Int_t d)       const; // map x to [-1:1]
+  template <class T>
+    T          MapToExternal(T  x,Int_t d)       const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
   //  
  protected:
   Int_t        fDimOut;            // dimension of the ouput array
@@ -171,17 +176,8 @@ class AliCheb3D: public TNamed
 };
 
 //__________________________________________________________________________________________
-inline Bool_t  AliCheb3D::IsInside(const Float_t  *par) const 
-{
-  // check if the point is inside of the fitted box
-  const float kTol = 1.e-4; 
-  for (int i=3;i--;) if (fBMin[i]-par[i]>kTol || par[i]-fBMax[i]>kTol) return kFALSE;
-  //if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
-  return kTRUE;
-}
-
-//__________________________________________________________________________________________
-inline Bool_t  AliCheb3D::IsInside(const Double_t  *par) const 
+template <class T>
+inline Bool_t  AliCheb3D::IsInside(const T  *par) const 
 {
   // check if the point is inside of the fitted box
   const float kTol = 1.e-4; 
@@ -191,7 +187,8 @@ inline Bool_t  AliCheb3D::IsInside(const Double_t  *par) const
 }
 
 //__________________________________________________________________________________________
-inline void AliCheb3D::Eval(const Float_t  *par, Float_t  *res)
+template <class T>
+inline void AliCheb3D::Eval(const T  *par, T  *res)
 {
   // evaluate Chebyshev parameterization for 3d->DimOut function
   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
@@ -200,7 +197,8 @@ inline void AliCheb3D::Eval(const Float_t  *par, Float_t  *res)
 }
 
 //__________________________________________________________________________________________
-inline Float_t AliCheb3D::Eval(const Float_t  *par, int idim)
+template <class T>
+inline T AliCheb3D::Eval(const T  *par, int idim)
 {
   // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
@@ -262,11 +260,12 @@ inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t  *par, i
 }
 
 //__________________________________________________________________________________________
-inline Float_t AliCheb3D::MapToInternal(Float_t  x,Int_t d) const
+template <class T>
+inline T AliCheb3D::MapToInternal(T  x,Int_t d) const
 {
   // map x to [-1:1]
 #ifdef _BRING_TO_BOUNDARY_
-  float res = (x-fBOffset[d])*fBScale[d];
+  T res = (x-fBOffset[d])*fBScale[d];
   if (res<-1) return -1;
   if (res> 1) return 1;
   return res;
index 7c08e1c..43159b5 100644 (file)
@@ -14,8 +14,8 @@
  **************************************************************************/
 
 #include <cstdlib>
-#include "AliCheb3DCalc.h"
 #include <TSystem.h>
+#include "AliCheb3DCalc.h"
 
 ClassImp(AliCheb3DCalc)
 
@@ -158,36 +158,10 @@ void AliCheb3DCalc::Print(const Option_t* ) const
 }
 
 //__________________________________________________________________________________________
-Float_t  AliCheb3DCalc::Eval(const Float_t  *par) const
-{
-  // evaluate Chebyshev parameterization for 3D function.
-  // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
-  const Float_t  &z = par[2];
-  const Float_t  &y = par[1];
-  const Float_t  &x = par[0];
-  //
-  int ncfRC;
-  for (int id0=fNRows;id0--;) {
-    int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
-    int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
-    for (int id1=nCLoc;id1--;) {
-      int id = id1+col0;
-      fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC) : 0.0;
-    }
-    fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(y,fTmpCf1,nCLoc):0.0;
-  }
-  return ChebEval1D(x,fTmpCf0,fNRows);
-  //
-}
-
-//__________________________________________________________________________________________
 Float_t  AliCheb3DCalc::EvalDeriv(int dim, const Float_t  *par) const
 {
   // evaluate Chebyshev parameterization derivative in given dimension  for 3D function.
   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
-  const Float_t  &z = par[2];
-  const Float_t  &y = par[1];
-  const Float_t  &x = par[0];
   //
   int ncfRC;
   for (int id0=fNRows;id0--;) {
@@ -198,13 +172,13 @@ Float_t  AliCheb3DCalc::EvalDeriv(int dim, const Float_t  *par) const
     for (int id1=nCLoc;id1--;) {
       int id = id1+col0;
       if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
-      if (dim==2) fTmpCf1[id1] = ChebEval1Deriv(z,fCoefs + fCoefBound2D1[id], ncfRC);
-      else        fTmpCf1[id1] = ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC);
+      if (dim==2) fTmpCf1[id1] = ChebEval1Deriv(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
+      else        fTmpCf1[id1] = ChebEval1D(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
     }
-    if (dim==1)   fTmpCf0[id0] = ChebEval1Deriv(y,fTmpCf1,nCLoc);
-    else          fTmpCf0[id0] = ChebEval1D(y,fTmpCf1,nCLoc);
+    if (dim==1)   fTmpCf0[id0] = ChebEval1Deriv(par[1],fTmpCf1,nCLoc);
+    else          fTmpCf0[id0] = ChebEval1D(par[1],fTmpCf1,nCLoc);
   }
-  return (dim==0) ? ChebEval1Deriv(x,fTmpCf0,fNRows) : ChebEval1D(x,fTmpCf0,fNRows);
+  return (dim==0) ? ChebEval1Deriv(par[0],fTmpCf0,fNRows) : ChebEval1D(par[0],fTmpCf0,fNRows);
   //
 }
 
@@ -213,9 +187,6 @@ Float_t  AliCheb3DCalc::EvalDeriv2(int dim1,int dim2, const Float_t  *par) const
 {
   // evaluate Chebyshev parameterization 2n derivative in given dimensions  for 3D function.
   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
-  const Float_t  &z = par[2];
-  const Float_t  &y = par[1];
-  const Float_t  &x = par[0];
   //
   Bool_t same = dim1==dim2;
   int ncfRC;
@@ -227,14 +198,15 @@ Float_t  AliCheb3DCalc::EvalDeriv2(int dim1,int dim2, const Float_t  *par) const
     for (int id1=nCLoc;id1--;) {
       int id = id1+col0;
       if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
-      if (dim1==2||dim2==2) fTmpCf1[id1] = same ? ChebEval1Deriv2(z,fCoefs + fCoefBound2D1[id], ncfRC) 
-                             :                   ChebEval1Deriv(z,fCoefs + fCoefBound2D1[id], ncfRC);
-      else        fTmpCf1[id1] = ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC);
+      if (dim1==2||dim2==2) fTmpCf1[id1] = same ? ChebEval1Deriv2(par[2],fCoefs + fCoefBound2D1[id], ncfRC) 
+                             :                   ChebEval1Deriv(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
+      else        fTmpCf1[id1] = ChebEval1D(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
     }
-    if (dim1==1||dim2==1) fTmpCf0[id0] = same ? ChebEval1Deriv2(y,fTmpCf1,nCLoc):ChebEval1Deriv(y,fTmpCf1,nCLoc);
-    else                  fTmpCf0[id0] = ChebEval1D(y,fTmpCf1,nCLoc);
+    if (dim1==1||dim2==1) fTmpCf0[id0] = same ? ChebEval1Deriv2(par[1],fTmpCf1,nCLoc):ChebEval1Deriv(par[1],fTmpCf1,nCLoc);
+    else                  fTmpCf0[id0] = ChebEval1D(par[1],fTmpCf1,nCLoc);
   }
-  return (dim1==0||dim2==0) ? (same ? ChebEval1Deriv2(x,fTmpCf0,fNRows):ChebEval1Deriv(x,fTmpCf0,fNRows)) : ChebEval1D(x,fTmpCf0,fNRows);
+  return (dim1==0||dim2==0) ? (same ? ChebEval1Deriv2(par[0],fTmpCf0,fNRows):ChebEval1Deriv(par[0],fTmpCf0,fNRows)) : 
+    ChebEval1D(par[0],fTmpCf0,fNRows);
   //
 }
 
@@ -445,4 +417,3 @@ Float_t AliCheb3DCalc::ChebEval1Deriv2(Float_t  x, const Float_t * array, int nc
   //
   return b0 - x*b1 - ddcf0/2;
 }
-
index 9f42b9e..8cb0293 100644 (file)
@@ -32,7 +32,6 @@ class AliCheb3DCalc: public TNamed
   AliCheb3DCalc& operator=(const AliCheb3DCalc& rhs);
   void       Print(const Option_t* opt="")                              const;
   void       LoadData(FILE* stream);
-  Float_t    Eval(const Float_t  *par)                                  const;
   Float_t    EvalDeriv(int dim, const Float_t  *par)                    const;
   Float_t    EvalDeriv2(int dim1,int dim2, const Float_t  *par)         const;
   //
@@ -57,6 +56,23 @@ class AliCheb3DCalc: public TNamed
   //
   static void ReadLine(TString& str,FILE* stream);
   //
+  template <class T>
+    T        Eval(const T  *par)                                        const {
+    // evaluate Chebyshev parameterization for 3D function.
+    // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
+    int ncfRC;
+    for (int id0=fNRows;id0--;) {
+      int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
+      int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
+      for (int id1=nCLoc;id1--;) {
+       int id = id1+col0;
+       fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(par[2],fCoefs + fCoefBound2D1[id], ncfRC) : 0.0;
+      }
+      fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(par[1],fTmpCf1,nCLoc):0.0;
+    }
+    return ChebEval1D(par[0],fTmpCf0,fNRows);
+    }
+  //
  protected:
   Int_t      fNCoefs;            // total number of coeeficients
   Int_t      fNRows;             // number of significant rows in the 3D coeffs matrix
index 97a37b6..5b5e8bd 100644 (file)
@@ -211,7 +211,67 @@ void AliFieldMap::ReadField()
 }
 
 //_______________________________________________________________________
-void AliFieldMap::Field(Float_t *x, Float_t *b) const
+void AliFieldMap::Field(float *x, float *b) const
+{
+  //
+  // Use simple interpolation to obtain field at point x
+  //
+    float ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
+      bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
+    const float kone=1;
+    Int_t ix, iy, iz;
+    b[0]=b[1]=b[2]=0;
+    //
+    
+    xl[0] = x[0] - fXbeg;
+    xl[1] = x[1] - fYbeg;
+    xl[2] = x[2] - fZbeg;
+    
+    hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
+    ratx=hix-int(hix);
+    ix=int(hix);
+    
+    hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
+    raty=hiy-int(hiy);
+    iy=int(hiy);
+    
+    hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
+    ratz=hiz-int(hiz);
+    iz=int(hiz);
+
+    ratx1=kone-ratx;
+    raty1=kone-raty;
+    ratz1=kone-ratz;
+
+    if (!fB) return;
+
+    bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
+    bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
+    blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
+    blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[0]  = blz               *ratz1+bhz               *ratz;
+    //
+    bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
+    bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
+    blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
+    blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[1]  = blz               *ratz1+bhz               *ratz;
+    //
+    bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
+    bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
+    blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
+    blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
+    bhz   = blyhz             *raty1+bhyhz             *raty;
+    blz   = blylz             *raty1+bhylz             *raty;
+    b[2]  = blz               *ratz1+bhz               *ratz;
+}
+
+//_______________________________________________________________________
+void AliFieldMap::Field(double *x, double *b) const
 {
   //
   // Use simple interpolation to obtain field at point x
index ee3e6fe..37bb1ac 100644 (file)
@@ -25,8 +25,8 @@ public:
     virtual ~AliFieldMap();
     void Copy(TObject &map) const;
     virtual AliFieldMap & operator=(const AliFieldMap &map);
-
-    virtual void Field(Float_t *x, Float_t *b) const;
+    virtual void Field(float *x, float *b) const;
+    virtual void Field(double *x, double *b) const;
     Float_t Bx(Int_t ix, Int_t iy, Int_t iz) const{
        return (*fB)(3*((ix*fYn+iy)*fZn+iz));
     }
index 3c18881..a3cf3e2 100644 (file)
@@ -86,7 +86,17 @@ AliMagF::AliMagF(const AliMagF &src):
 }
 
 //_______________________________________________________________________
-void AliMagF::Field(Float_t*, Float_t *b) const
+void AliMagF::Field(float*, float *b) const
+{
+  //
+  // Method to return the field in one point -- dummy in this case
+  //
+  AliWarning("Undefined MagF Field called, returning 0");
+  b[0]=b[1]=b[2]=0;
+}
+
+//_______________________________________________________________________
+void AliMagF::Field(double*, double *b) const
 {
   //
   // Method to return the field in one point -- dummy in this case
index 3344d3a..c360d53 100644 (file)
@@ -24,7 +24,8 @@ public:
   AliMagF(const AliMagF& maps);
   virtual ~AliMagF() {}
   AliMagF& operator=(const AliMagF& rhs);
-  virtual void    Field(Float_t *x, Float_t *b) const;
+  virtual void    Field(float *x, float *b)                  const;
+  virtual void    Field(double *x, double *b)                const;
   virtual void    GetTPCInt(Float_t *xyz, Float_t *b)        const;
   virtual void    GetTPCIntCyl(Float_t *rphiz, Float_t *b)   const;
   virtual Int_t   Type() const {return fType;}
index ac6b90f..7f28979 100644 (file)
@@ -68,7 +68,38 @@ AliMagFC::AliMagFC(const char *name, const char *title, Int_t integ,
 }
 
 //________________________________________
-void AliMagFC::Field(Float_t *x, Float_t *b) const
+void AliMagFC::Field(float *x, float *b) const
+{
+  //
+  // Method to return the field in a point
+  //
+  b[0]=b[1]=b[2]=0;
+  if(fMap==1) {
+    if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
+      b[2]=2;
+    } 
+    else {
+      if(-725 >= x[2] && x[2] >= -1225 ){
+       Float_t dz = TMath::Abs(-975-x[2])*0.01;
+       b[0] = - (1-0.1*dz*dz)*7;
+       if(fFactor!=1) {
+           b[0]*=fFactor;
+           b[1]*=fFactor;
+           b[2]*=fFactor;
+       }
+      }
+      else {
+         ZDCField(x, b);
+      }
+    }
+
+  } 
+  else {
+      AliFatal(Form("Invalid field map for constant field %d",fMap));
+  }
+}
+//________________________________________
+void AliMagFC::Field(double *x, double *b) const
 {
   //
   // Method to return the field in a point
@@ -100,11 +131,155 @@ void AliMagFC::Field(Float_t *x, Float_t *b) const
 }
 
 //___________________________________________________
-void AliMagFC::ZDCField(Float_t *x, Float_t *b) const
+void AliMagFC::ZDCField(float *x, float *b) const
+{
+  // ---- This is the ZDC part
+  
+  float rad2 = x[0] * x[0] + x[1] * x[1];
+  static Bool_t init = kFALSE;
+
+  if (! init) {
+      init = kTRUE;
+      //////////////////////////////////////////////////////////////////////
+      // ---- Magnetic field values (according to beam type and energy) ----
+      if(fBeamType==kBeamTypepp && fBeamEnergy == 5000.){
+         // p-p @ 5+5 TeV
+         fQuadGradient = 15.7145;
+         fDipoleField  = 27.0558;
+         // SIDE C
+         fCCorrField   = 9.7017;
+         // SIDE A
+         fACorr1Field  = -13.2143;
+         fACorr2Field  = -11.9909;
+      } else if (fBeamType == kBeamTypepp && fBeamEnergy == 450.) {
+         // p-p 0.45+0.45 TeV
+         Float_t const kEnergyRatio = fBeamEnergy / 7000.;
+         
+         fQuadGradient = 22.0002 * kEnergyRatio;
+         fDipoleField  = 37.8781 * kEnergyRatio;
+         // SIDE C
+         fCCorrField   =  9.6908;
+         // SIDE A
+         fACorr1Field  = -13.2014;
+         fACorr2Field  = -9.6908;
+      } else if ((fBeamType == kBeamTypepp && fBeamEnergy == 7000.) ||
+                (fBeamType == kBeamTypeAA))
+      {
+         // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
+         fQuadGradient = 22.0002;
+         fDipoleField  = 37.8781;
+         // SIDE C
+         fCCorrField   = 9.6908;
+         // SIDE A
+         fACorr1Field  = -13.2014;
+         fACorr2Field  = -9.6908;
+      }
+  }
+  
+  
+  // SIDE C **************************************************
+  if(x[2]<0.){  
+    if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
+       if (fFactor != 0.) {
+           b[0] = fCCorrField;
+           b[1] = 0.;
+           b[2] = 0.;
+       } 
+    }
+    else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
+       b[0] = fQuadGradient*x[1];
+       b[1] = fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
+       b[0] = -fQuadGradient*x[1];
+       b[1] = -fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
+       b[0] = -fQuadGradient*x[1];
+       b[1] = -fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
+       b[0] = fQuadGradient*x[1];
+       b[1] = fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
+       b[1] = fDipoleField;
+       b[2] = 0.;
+       b[2] = 0.;
+    }
+    else if(x[2] < kCD2Begin && x[2] > kCD2End){
+       if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
+          || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
+         b[1] = -fDipoleField;
+         b[2] = 0.;
+         b[2] = 0.;
+       }
+    }
+  }
+  
+  // SIDE A **************************************************
+  else{        
+    if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
+      // Compensator magnet at z = 1075 m 
+       if (fFactor != 0.) {
+           b[0] = fACorr1Field;
+           b[1] = 0.;
+           b[2] = 0.;
+       }
+       return;
+    }
+    
+    if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
+       if (fFactor != 0.) {
+           b[0] = fACorr2Field;
+           b[1] = 0.;
+           b[2] = 0.;
+       }
+    }          
+    else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
+       // First quadrupole of inner triplet de-focussing in x-direction
+       b[0] = -fQuadGradient*x[1];
+       b[1] = -fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
+       b[0] = fQuadGradient*x[1];
+       b[1] = fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
+       b[0] = fQuadGradient*x[1];
+       b[1] = fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
+       b[0] = -fQuadGradient*x[1];
+       b[1] = -fQuadGradient*x[0];
+       b[2] = 0.;
+    }
+    else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
+       b[0] = 0.;
+       b[1] = -fDipoleField;
+       b[2] = 0.;
+    }
+    else if(x[2] > kAD2Begin && x[2] < kAD2End){
+       if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
+          || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
+           b[1] = fDipoleField;
+       }
+    }
+  }
+}
+
+void AliMagFC::ZDCField(double *x, double *b) const
 {
   // ---- This is the ZDC part
   
-  Float_t rad2 = x[0] * x[0] + x[1] * x[1];
+  double rad2 = x[0] * x[0] + x[1] * x[1];
   static Bool_t init = kFALSE;
 
   if (! init) {
index 680ea4a..55cc14b 100644 (file)
@@ -23,9 +23,11 @@ public:
   AliMagFC(const char *name, const char *title, Int_t integ, 
           Float_t factor, Float_t fmax);
   virtual ~AliMagFC(){}
-  virtual void Field(Float_t *x, Float_t *b) const;
+  virtual void Field(float *x, float *b)      const;
+  virtual void Field(double *x, double *b)    const;
   virtual void ReadField() {}
-  virtual void ZDCField(Float_t *x, Float_t *b) const;
+  virtual void ZDCField(float *x, float *b)   const;
+  virtual void ZDCField(double *x, double *b) const;
   virtual void SetBeamType(BeamType_t type)      {fBeamType    = type;}
   virtual void SetBeamEnergy(Float_t energy)     {fBeamEnergy  = energy;}
   virtual void SetCompensatorMagnet(Bool_t flag) {fCompensator = flag;}
index 8d808d4..cd11c36 100644 (file)
@@ -111,7 +111,107 @@ AliMagFCM::AliMagFCM(const AliMagFCM &magf):
 }
 
 //_______________________________________________________________________
-void AliMagFCM::Field(Float_t *x, Float_t *b) const
+void AliMagFCM::Field(float *x, float *b) const
+{
+  //
+  // Method to calculate the magnetic field
+  //
+  Float_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
+    bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
+  const Float_t kone=1;
+  Int_t ix, iy, iz;
+  
+  // --- find the position in the grid ---
+  
+  b[0]=b[1]=b[2]=0;
+
+  
+  if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)*(x[1]+30) < 560*560) {
+    b[2]=  fSolenoid;
+  } else  {
+      // The field map used here was calculated in a coordinate system where the muon arm is at z > 0
+      // Transfom x -> -x and z -> -z 
+      Float_t xm = - x[0];
+      Float_t ym =   x[1];
+      Float_t zm = - x[2];
+      Bool_t infield=(fZbeg <= zm  && zm < fZbeg+fZdel*(fZn-1)
+                     &&  ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
+                     &&  ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
+    if(infield) {
+      xl[0]=TMath::Abs(xm)-fXbeg;
+      xl[1]=TMath::Abs(ym)-fYbeg;
+      xl[2]=zm-fZbeg;
+      
+      // --- start with x
+      
+      hix=xl[0]*fXdeli;
+      ratx=hix-int(hix);
+      ix=int(hix);
+      
+      hiy=xl[1]*fYdeli;
+      raty=hiy-int(hiy);
+      iy=int(hiy);
+      
+      hiz=xl[2]*fZdeli;
+      ratz=hiz-int(hiz);
+      iz=int(hiz);
+      
+      if(fMap==2) {
+       // ... simple interpolation
+       ratx1=kone-ratx;
+       raty1=kone-raty;
+       ratz1=kone-ratz;
+       bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
+       bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
+       blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
+       blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
+       bhz   = blyhz             *raty1+bhyhz             *raty;
+       blz   = blylz             *raty1+bhylz             *raty;
+       b[0]  = blz               *ratz1+bhz               *ratz;
+       //
+       bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
+       bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
+       blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
+       blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
+       bhz   = blyhz             *raty1+bhyhz             *raty;
+       blz   = blylz             *raty1+bhylz             *raty;
+       b[1]  = blz               *ratz1+bhz               *ratz;
+       //
+       bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
+       bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
+       blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
+       blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
+       bhz   = blyhz             *raty1+bhyhz             *raty;
+       blz   = blylz             *raty1+bhylz             *raty;
+       b[2]  = blz               *ratz1+bhz               *ratz;
+       //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
+       //ratx,raty,ratz,b[0],b[1],b[2]);
+       //
+       // ... use the dipole symmetry
+       if (xm*ym < 0) b[1]=-b[1];
+       if (xm<0)      b[2]=-b[2];
+       b[0] = -b[0];
+       b[2] = -b[2];
+       
+      } else {
+       AliError(Form("Invalid field map for constant mesh %d",fMap));
+      }
+    } else {
+//This is the ZDC part
+       ZDCField(x,b);
+    }
+    
+    if(fFactor!=1) {
+       b[0]*=fFactor;
+       b[1]*=fFactor;
+       b[2]*=fFactor;
+    }
+  }
+}
+
+//_______________________________________________________________________
+void AliMagFCM::Field(double *x, double *b) const
 {
   //
   // Method to calculate the magnetic field
index 31a5068..77e3531 100644 (file)
@@ -24,7 +24,8 @@ public:
           Float_t factor, Float_t fmax);
   AliMagFCM(const AliMagFCM &mag);
   virtual ~AliMagFCM() {delete fB;}
-  virtual void Field(Float_t *x, Float_t *b) const;
+  virtual void Field(float *x, float *b) const;
+  virtual void Field(double *x, double *b) const;
   virtual void ReadField();
   virtual void    SetSolenoidField(Float_t field = 2.) {fSolenoid = field;}
   virtual Float_t SolenoidField() const {
index 044b574..058a877 100644 (file)
@@ -237,7 +237,7 @@ void AliMagFCheb::Clear(const Option_t *)
 }
 
 //__________________________________________________________________________________________
-void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
+void AliMagFCheb::Field(float *xyz, float *b) const
 {
   // compute field in cartesian coordinates. If point is outside of the parameterized region
   // get it at closest valid point
@@ -272,6 +272,41 @@ void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
 }
 
 //__________________________________________________________________________________________
+void AliMagFCheb::Field(double *xyz, double *b) const
+{
+  // compute field in cartesian coordinates. If point is outside of the parameterized region
+  // get it at closest valid point
+  static double rphiz[3];
+  //
+#ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
+  if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) && 
+       !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip())  ) {for (int i=3;i--;) b[i]=0; return;}
+#endif
+  //
+  if (xyz[2]<fMaxZDip) {    // dipole part?
+#ifndef _BRING_TO_BOUNDARY_
+    AliCheb3D* par = GetParamDip(FindDipSegment(xyz));
+    if (par->IsInside(xyz)) {par->Eval(xyz,b); return;}
+    for (int i=3;i--;) b[i]=0; return;
+#else
+    GetParamDip(FindDipSegment(xyz))->Eval(xyz,b); return;  
+#endif
+  }
+  //
+  // Sol region: convert coordinates to cyl system
+  CartToCyl(xyz,rphiz);
+#ifndef _BRING_TO_BOUNDARY_
+  if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
+#endif
+  //
+  FieldCylSol(rphiz,b);
+  //
+  // convert field to cartesian system
+  CylToCartCylB(rphiz, b,b);
+  //
+}
+
+//__________________________________________________________________________________________
 void AliMagFCheb::GetTPCInt(Float_t *xyz, Float_t *b) const
 {
   // compute TPC region field integral in cartesian coordinates.
@@ -294,17 +329,29 @@ void AliMagFCheb::GetTPCInt(Float_t *xyz, Float_t *b) const
 }
 
 //__________________________________________________________________________________________
-void AliMagFCheb::FieldCylSol(const Float_t *rphiz, Float_t *b) const
+void AliMagFCheb::FieldCylSol(const float *rphiz, float *b) const
 {
   // compute Solenoid field in Cylindircal coordinates
   // note: if the point is outside the volume get the field in closest parameterized point
-  const float &r = rphiz[0];
-  const float &z = rphiz[2];
   int SolZId = 0;
-  while (z>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
+  while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
   int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
   int SolRMax = SolRId + fNSegRSol[SolZId];
-  while (r>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
+  while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
+  GetParamSol( SolRId )->Eval(rphiz,b);
+  //
+}
+
+//__________________________________________________________________________________________
+void AliMagFCheb::FieldCylSol(const double *rphiz, double *b) const
+{
+  // compute Solenoid field in Cylindircal coordinates
+  // note: if the point is outside the volume get the field in closest parameterized point
+  int SolZId = 0;
+  while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
+  int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
+  int SolRMax = SolRId + fNSegRSol[SolZId];
+  while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
   GetParamSol( SolRId )->Eval(rphiz,b);
   //
 }
@@ -314,13 +361,11 @@ void AliMagFCheb::GetTPCIntCyl(Float_t *rphiz, Float_t *b) const
 {
   // compute field integral in TPC region in Cylindircal coordinates
   // note: the check for the point being inside the parameterized region is done outside
-  const float &r = rphiz[0];
-  const float &z = rphiz[2];
   int tpcIntZId = 0;
-  while (z>fSegZTPCInt[tpcIntZId] && tpcIntZId<fNSegZTPCInt) ++tpcIntZId;    // find Z segment
+  while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZId<fNSegZTPCInt) ++tpcIntZId;    // find Z segment
   int tpcIntRId = fSegZIdTPCInt[tpcIntZId];        // first R segment for this Z
   int tpcIntRIdMax = tpcIntRId + fNSegRTPCInt[tpcIntZId];
-  while (r>fSegRTPCInt[tpcIntRId] && tpcIntRId<tpcIntRIdMax) ++tpcIntRId;    // find R segment
+  while (rphiz[0]>fSegRTPCInt[tpcIntRId] && tpcIntRId<tpcIntRIdMax) ++tpcIntRId;    // find R segment
   GetParamTPCInt( tpcIntRId )->Eval(rphiz,b);
   //
 }
@@ -468,29 +513,9 @@ void AliMagFCheb::LoadData(const char* inpfile)
 }
 #endif
 
-//_________________________________________________________________________
-Int_t AliMagFCheb::FindDipSegment(const float *xyz) const
-{
-  // find the segment containing point xyz. If it is outside find the closest segment 
-  int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,xyz[2]); // find zsegment
-  int ysegBeg = fBegSegYDip[zid];
-  //
-  for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
-  if ( --yid < 0 ) yid = 0;
-  yid +=  ysegBeg;
-  //
-  int xsegBeg = fBegSegXDip[yid];
-  for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
-  if ( --xid < 0) xid = 0;
-  xid +=  xsegBeg;
-  //
-  return fSegIDDip[xid];
-}
-
 //_______________________________________________
 #ifdef  _INC_CREATION_ALICHEB3D_
 
-
 //__________________________________________________________________________________________
 AliMagFCheb::AliMagFCheb(const char* inputFile) : 
   fNParamsSol(0),
index d6402b1..d5f34de 100644 (file)
@@ -76,25 +76,36 @@ class AliMagFCheb: public TNamed
   Float_t    GetMaxZTPCInt()                              const {return fMaxZTPCInt;}
   Float_t    GetMaxRTPCInt()                              const {return fMaxRTPCInt;}
   //
-  Int_t      FindDipSegment(const float *xyz)             const;
   AliCheb3D* GetParamSol(Int_t ipar)                      const {return (AliCheb3D*)fParamsSol->UncheckedAt(ipar);}
   AliCheb3D* GetParamTPCInt(Int_t ipar)                   const {return (AliCheb3D*)fParamsTPCInt->UncheckedAt(ipar);}
   AliCheb3D* GetParamDip(Int_t ipar)                      const {return (AliCheb3D*)fParamsDip->UncheckedAt(ipar);}
   //
   virtual void Print(Option_t * = "")                     const;
   //
-  virtual void Field(Float_t *xyz, Float_t *b)            const;
-  virtual void FieldCyl(const Float_t *rphiz, Float_t *b) const;
+  virtual void Field(float *xyz, float *b)                const;
+  virtual void Field(double *xyz, double *b)              const;
+  //
+  virtual void FieldCyl(const float *rphiz, float *b)     const;
+  virtual void FieldCyl(const double *rphiz, double *b)   const;
   //
   virtual void GetTPCInt(Float_t *xyz, Float_t *b)        const;
   virtual void GetTPCIntCyl(Float_t *rphiz, Float_t *b)   const;
   //
-  static void CylToCartCylB(const float *rphiz, const float *brphiz,float *bxyz);
-  static void CylToCartCartB(const float *xyz,  const float *brphiz,float *bxyz);
-  static void CartToCylCartB(const float *xyz,  const float *bxyz,  float *brphiz);
-  static void CartToCylCylB(const float *rphiz, const float *bxyz,  float *brphiz);
-  static void CartToCyl(const float *xyz,  float *rphiz);
-  static void CylToCart(const float *rphiz,float *xyz);
+  template <class T>
+    Int_t      FindDipSegment(const T *xyz)               const; 
+  //
+  template <class T>
+    static void CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz);
+  template <class T>
+    static void CylToCartCartB(const T *xyz,  const T *brphiz,T *bxyz);
+  template <class T>
+    static void CartToCylCartB(const T *xyz,  const T *bxyz,  T *brphiz);
+  template <class T>  
+    static void CartToCylCylB(const T *rphiz, const T *bxyz,  T *brphiz);
+  template <class T>
+    static void CartToCyl(const T *xyz,  T *rphiz);
+  template <class T>
+    static void CylToCart(const T *rphiz,T *xyz);
   //
 #ifdef  _INC_CREATION_ALICHEB3D_                          // see AliCheb3D.h for explanation
   void         LoadData(const char* inpfile);
@@ -111,11 +122,13 @@ class AliMagFCheb: public TNamed
   void       BuildTableSol();
   void       BuildTableTPCInt();
   void       ResetTPCInt();
-
+  //
+  //
 #endif
   //
  protected:
-  virtual void FieldCylSol(const Float_t *rphiz, Float_t *b)    const;
+    virtual void FieldCylSol(const float *rphiz, float *b)      const;
+    virtual void FieldCylSol(const double *rphiz, double *b)    const;
   //
  protected:
   //
@@ -173,7 +186,16 @@ class AliMagFCheb: public TNamed
 
 
 //__________________________________________________________________________________________
-inline void AliMagFCheb::FieldCyl(const Float_t *rphiz, Float_t *b) const
+inline void AliMagFCheb::FieldCyl(const float *rphiz, float *b) const
+{
+  // compute field in Cylindircal coordinates
+  //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
+  FieldCylSol(rphiz,b);
+}
+
+
+//__________________________________________________________________________________________
+inline void AliMagFCheb::FieldCyl(const double *rphiz, double *b) const
 {
   // compute field in Cylindircal coordinates
   //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
@@ -181,11 +203,12 @@ inline void AliMagFCheb::FieldCyl(const Float_t *rphiz, Float_t *b) const
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CylToCartCylB(const float *rphiz, const float *brphiz,float *bxyz)
+template <class T>
+inline void AliMagFCheb::CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz)
 {
   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
-  float btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
-  float psiPLUSphi = TMath::ATan2(brphiz[1],brphiz[0]) + rphiz[1];
+  T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
+  T psiPLUSphi = TMath::ATan2(brphiz[1],brphiz[0]) + rphiz[1];
   bxyz[0] = btr*TMath::Cos(psiPLUSphi);
   bxyz[1] = btr*TMath::Sin(psiPLUSphi);
   bxyz[2] = brphiz[2];
@@ -193,11 +216,12 @@ inline void AliMagFCheb::CylToCartCylB(const float *rphiz, const float *brphiz,f
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CylToCartCartB(const float *xyz, const float *brphiz,float *bxyz)
+template <class T>
+inline void AliMagFCheb::CylToCartCartB(const T *xyz, const T *brphiz, T *bxyz)
 {
   // convert field in cylindrical coordinates to cartesian system, point is in cart.system
-  float btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
-  float phiPLUSpsi = TMath::ATan2(xyz[1],xyz[0]) +  TMath::ATan2(brphiz[1],brphiz[0]);
+  T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
+  T phiPLUSpsi = TMath::ATan2(xyz[1],xyz[0]) +  TMath::ATan2(brphiz[1],brphiz[0]);
   bxyz[0] = btr*TMath::Cos(phiPLUSpsi);
   bxyz[1] = btr*TMath::Sin(phiPLUSpsi);
   bxyz[2] = brphiz[2];
@@ -205,11 +229,12 @@ inline void AliMagFCheb::CylToCartCartB(const float *xyz, const float *brphiz,fl
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CartToCylCartB(const float *xyz, const float *bxyz, float *brphiz)
+template <class T>
+inline void AliMagFCheb::CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz)
 {
   // convert field in cylindrical coordinates to cartesian system, poin is in cart.system
-  float btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
-  float psiMINphi = TMath::ATan2(bxyz[1],bxyz[0]) - TMath::ATan2(xyz[1],xyz[0]);
+  T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+  T psiMINphi = TMath::ATan2(bxyz[1],bxyz[0]) - TMath::ATan2(xyz[1],xyz[0]);
   //
   brphiz[0] = btr*TMath::Cos(psiMINphi);
   brphiz[1] = btr*TMath::Sin(psiMINphi);
@@ -218,11 +243,12 @@ inline void AliMagFCheb::CartToCylCartB(const float *xyz, const float *bxyz, flo
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CartToCylCylB(const float *rphiz, const float *bxyz, float *brphiz)
+template <class T>
+inline void AliMagFCheb::CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz)
 {
   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
-  float btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
-  float psiMINphi =  TMath::ATan2(bxyz[1],bxyz[0]) - rphiz[1];
+  T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+  T psiMINphi =  TMath::ATan2(bxyz[1],bxyz[0]) - rphiz[1];
   brphiz[0] = btr*TMath::Cos(psiMINphi);
   brphiz[1] = btr*TMath::Sin(psiMINphi);
   brphiz[2] = bxyz[2];
@@ -230,7 +256,8 @@ inline void AliMagFCheb::CartToCylCylB(const float *rphiz, const float *bxyz, fl
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CartToCyl(const float *xyz,float *rphiz)
+template <class T>
+inline void AliMagFCheb::CartToCyl(const T *xyz,T *rphiz)
 {
   rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
   rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
@@ -238,11 +265,32 @@ inline void AliMagFCheb::CartToCyl(const float *xyz,float *rphiz)
 }
 
 //__________________________________________________________________________________________________
-inline void AliMagFCheb::CylToCart(const float *rphiz, float *xyz)
+template <class T>
+inline void AliMagFCheb::CylToCart(const T *rphiz, T *xyz)
 {
   xyz[0] = rphiz[0]*TMath::Cos(rphiz[1]);
   xyz[1] = rphiz[0]*TMath::Sin(rphiz[1]);
   xyz[2] = rphiz[2];
 }
 
+//__________________________________________________________________________________________________
+template <class T>
+Int_t    AliMagFCheb::FindDipSegment(const T *xyz) const 
+{
+  // find the segment containing point xyz. If it is outside find the closest segment 
+  int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(float)xyz[2]); // find zsegment
+  int ysegBeg = fBegSegYDip[zid];
+  //
+  for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
+  if ( --yid < 0 ) yid = 0;
+  yid +=  ysegBeg;
+  //
+  int xsegBeg = fBegSegXDip[yid];
+  for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
+  if ( --xid < 0) xid = 0;
+  xid +=  xsegBeg;
+  //
+  return fSegIDDip[xid];
+}
+
 #endif
index e4cccc9..23ff5af 100644 (file)
@@ -108,7 +108,7 @@ AliMagFDM::AliMagFDM(const char *name, const char *title, Int_t integ,
 }
 
 //_______________________________________________________________________
-void AliMagFDM::Field(Float_t *xfi, Float_t *b) const
+void AliMagFDM::Field(float *xfi, float *b) const
 {
   //
   // Main routine to compute the field in a point
@@ -336,6 +336,235 @@ if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2]
 }
 
 //_______________________________________________________________________
+void AliMagFDM::Field(double *xfi, double *b) const
+{
+  //
+  // Main routine to compute the field in a point
+  //
+  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;
+  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];
+  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]);
+           r0=rr/100;
+  Double_t rPmax;
+  rPmax=fRmax; 
+  if ( (-700<x[2] && x[2]<=kfZbg && 
+        (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
+       || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
+       {
+        b[0]=b[1]=0;
+        b[2]=fSolenoid;
+       }
+
+  xL3[0]=x[0]/100;
+  xL3[1]=(x[1]+30)/100;
+  xL3[2]=x[2]/100; 
+
+  if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
+  if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
+
+  Double_t xminn=xL3[2]*fAx1+fCx1;
+  Double_t xmaxx=xL3[2]*fAx2+fCx2;
+  Double_t zCmin,zCmax,yCmin,yCmax;
+
+  zCmin=fZmin;
+  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(fMap==3) 
+      { 
+        if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
+       {
+       //---------------------  Polar part ----------------------
+       
+       Double_t yyp,ph0;
+       Int_t  kp0, lp0, mp0;
+       Int_t kpi,lpi,mpi;
+       Double_t alp1,alp2,alp3;
+       Double_t zpz,rp,fip,cphi; 
+
+       kpi=kiip; 
+       lpi=kliip;
+       mpi=kmiip;
+   
+       zpz=xL3[2];
+       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=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 && 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;
+
+       if (r0<= fRdel)
+        {
+
+         if(r0< keps) 
+         {
+
+          bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;     
+          bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;      
+          bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10; 
+
+         } else { 
+      
+        alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0]; 
+        alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];      
+        alp1= kone - alp2 - alp3;
+     
+        for (jb=0; jb<3 ; jb++) 
+         {
+          iKvar=jb;     
+          bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
+         } 
+        } // end   of keps <=r0 
+       }
+        else
+       {
+        rp=r0;
+
+        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;
+         bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
+        }
+       }
+
+       b[0]=-bint[0];
+       b[1]=bint[1];
+       b[2]=-bint[2];
+
+   }  
+   else 
+   {
+   //-------------- Cartensian part ------------------
+          
+   Double_t zzc,yyc; 
+   Int_t k0, l0,m0;
+   Double_t  xx1, xx2,dx, xxx ,xXl; 
+   Int_t kci,mci,lci;
+
+   kci=kiic;
+   lci=kliic;
+   mci=kmiic;
+
+   xx1 = fAx1*xL3[2] + fCx1;
+   xx2 = fAx2*xL3[2] + fCx2;
+
+   zzc=xL3[2];
+   k0=FZ(zzc, fZc ,fZdel, kci, fZl); 
+   zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
+   zz1=kone - zz2;;    
+
+   yyc=xL3[1];
+   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;
+   m0 = int(xxx/dx);   
+
+   if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx)) 
+    {
+      m0=m0+1;   
+      AliDebug(2,Form(" m0 %d, m0+1 %d\n",m0,m0+1));  
+    }
+
+   x2=(xL3[0]-( xx1+m0*dx))/dx;
+   x1=kone-x2;
+   m0=m0-1;
+   for (jb=3; jb<6; jb++) 
+     {
+       iKvar=jb;     
+       bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ; 
+     }    
+   b[0]=-bint[0];
+   b[1]=bint[1];
+   b[2]=-bint[2]; 
+
+   } 
+
+
+  } else {
+        AliError(Form("Unknown map of Dipole region %d",fMap));
+ }
+           
+} else {
+    ZDCField(xfi,b);
+
+    }
+
+  if(fFactor!=1) {
+    b[0]*=fFactor;
+    b[1]*=fFactor;
+    b[2]*=fFactor;
+  }
+}
+
+
+//_______________________________________________________________________
 Int_t AliMagFDM::FZ(Double_t temp, const Float_t *Ar, 
                     Float_t delu, Int_t ik, Int_t nk) const
 {
index ee102b5..031e1dd 100644 (file)
@@ -22,7 +22,8 @@ public:
   AliMagFDM(const char *name, const char *title, Int_t integ,
            Float_t factor, Float_t fmax);
   virtual ~AliMagFDM(){} 
-  virtual void Field(Float_t *x, Float_t *b) const;
+  virtual void Field(float *x, float *b) const;
+  virtual void Field(double *x, double *b) const;
   virtual void ReadField(); 
   virtual void SetSolenoidField(Float_t field = 2.) {fSolenoid = field;}
   virtual Float_t SolenoidField() const {
index 5ac11cb..9051b92 100644 (file)
@@ -187,7 +187,7 @@ Float_t AliMagFMaps::SolenoidField() const
 }
 
 //_______________________________________________________________________
-void AliMagFMaps::Field(Float_t *x, Float_t *b) const
+void AliMagFMaps::Field(float *x, float *b) const
 {
   //
   // Method to calculate the magnetic field
@@ -245,6 +245,62 @@ void AliMagFMaps::Field(Float_t *x, Float_t *b) const
 }
 
 //_______________________________________________________________________
+void AliMagFMaps::Field(double *x, double *b) const
+{
+  //
+  // Method to calculate the magnetic field
+  //
+  // --- find the position in the grid ---
+  
+ //    if (!fFieldRead) ReadField();
+
+  //
+  // Field Maps have been calculated for the coordinate system in which 
+  // the Muon Spectrometer is placed at z > 0
+  // Transform coordinates corresponingly 
+  //
+  
+  b[0]=b[1]=b[2]=0;
+  double xm[3];
+  xm[0] = - x[0];
+  xm[1] =   x[1];
+  xm[2] = - x[2];
+  
+  AliFieldMap* map = 0;
+  if (fFieldMap[0]->Inside(xm[0], xm[1], xm[2])) {
+      map = fFieldMap[0];
+      Float_t r = TMath::Sqrt(xm[0] * xm[0] + xm[1] * xm[1]);
+      
+      if (!fL3Option && TMath::Abs(xm[2]) < 370. && r < 550.) {
+      //
+      //     Constant L3 field , if this option was selected
+      //
+       b[2] = (- fSolenoid)*fFactor;
+         return;
+      } 
+  } else if (fFieldMap[1]->Inside(xm[0], xm[1], xm[2])) {
+    map = fFieldMap[1];
+  } else if (fFieldMap[2]->Inside(xm[0], xm[1], xm[2])) {
+    map = fFieldMap[2];
+  }
+  
+  if(map){
+      map->Field(xm,b);
+      b[0] = - b[0];
+      b[2] = - b[2];
+      
+      if(fFactor!=1) {
+         b[0]*=fFactor;
+         b[1]*=fFactor;
+         b[2]*=fFactor;
+      }
+  } else {
+      //This is the ZDC part
+      ZDCField(x, b);
+  }
+}
+
+//_______________________________________________________________________
 void AliMagFMaps::Copy(TObject & /* magf */) const
 {
   //
index 2aa458f..2d5ced5 100644 (file)
@@ -26,7 +26,8 @@ public:
                Int_t l3 = 1);
     AliMagFMaps(const AliMagFMaps &mag);
     virtual ~AliMagFMaps();
-    virtual void    Field(Float_t *x, Float_t *b) const;
+    virtual void    Field(float  *x, float  *b) const;
+    virtual void    Field(double *x, double *b) const;
     AliFieldMap* FieldMap(Int_t i) {return fFieldMap[i];}
     virtual void ReadField();
     virtual Float_t SolenoidField() const;
index 3c9f61e..0eae1a9 100644 (file)
@@ -94,18 +94,49 @@ AliMagFMapsV1::~AliMagFMapsV1()
 }
 
 //_______________________________________________________________________
-void AliMagFMapsV1::Field(Float_t *x, Float_t *b) const
+void AliMagFMapsV1::Field(float *x, float *b) const
 {
   //
   // Method to calculate the magnetic field at position x
   //
-    const Float_t kRmax2 = 500. * 500.;
-    const Float_t kZmax  = 550.; 
-    const Float_t kTeslaTokG = 10.;
-    const Float_t kScale = 0.98838; // matching factor
+    const float kRmax2 = 500. * 500.;
+    const float kZmax  = 550.; 
+    const float kTeslaTokG = 10.;
+    const float kScale = 0.98838; // matching factor
     
     // Check if position inside measured map
-    Float_t r2 = x[0] * x[0] + x[1] * x[1];
+    float r2 = x[0] * x[0] + x[1] * x[1];
+    if (fMeasuredMap              &&
+       r2 < kRmax2               && 
+       TMath::Abs(x[2]) < kZmax
+       ) 
+    {
+       fMeasuredMap->Field(x, b);
+       b[0] *= kTeslaTokG;
+       b[1] *= kTeslaTokG;
+       b[2] *= kTeslaTokG;
+    } else {
+       AliMagFMaps::Field(x, b);
+       // Match to measure map
+       b[0] = - b[0] * kScale;
+       b[2] = - b[2] * kScale;
+       b[1] = - b[1] * kScale;
+    }
+}
+
+//_______________________________________________________________________
+void AliMagFMapsV1::Field(double *x, double *b) const
+{
+  //
+  // Method to calculate the magnetic field at position x
+  //
+    const double kRmax2 = 500. * 500.;
+    const double kZmax  = 550.; 
+    const double kTeslaTokG = 10.;
+    const double kScale = 0.98838; // matching factor
+    
+    // Check if position inside measured map
+    double r2 = x[0] * x[0] + x[1] * x[1];
     if (fMeasuredMap              &&
        r2 < kRmax2               && 
        TMath::Abs(x[2]) < kZmax
index 7ee82d2..b28be52 100644 (file)
@@ -25,7 +25,8 @@ public:
     AliMagFMapsV1(const AliMagFMapsV1& maps);             
     AliMagFMapsV1& operator=(const AliMagFMapsV1& maps) {maps.Copy(*this); return *this;}
     virtual ~AliMagFMapsV1();
-    virtual void    Field(Float_t *x, Float_t *b) const;
+    virtual void    Field(float *x, float *b) const;
+    virtual void    Field(double *x, double *b) const;
     virtual Float_t SolenoidField() const;
     AliMagFCheb* GetMeasuredMap()                   const {return fMeasuredMap;}
     void SetMeasuredMap(AliMagFCheb* parm)               {if (parm) delete parm; fMeasuredMap = parm;}
index 30cfa4e..d0e397b 100644 (file)
@@ -108,7 +108,23 @@ void AliMagWrapCheb::GetTPCIntCyl(Float_t *rphiz, Float_t *b) const
 }
 
 //_______________________________________________________________________
-void AliMagWrapCheb::Field(Float_t *xyz, Float_t *b) const
+void AliMagWrapCheb::Field(float *xyz, float *b) const
+{
+  // Method to calculate the field at point  xyz
+  //
+  b[0]=b[1]=b[2]=0.0;
+    if (xyz[2] > 919. || xyz[2] < -1972.) {
+       ZDCField(xyz, b);
+    } else {
+       if (fMeasuredMap && fFactor !=0.) {
+           fMeasuredMap->Field(xyz,b);
+           for (int i=3;i--;) b[i] *= fFactor;
+       }
+    }
+}
+
+//_______________________________________________________________________
+void AliMagWrapCheb::Field(double *xyz, double *b) const
 {
   // Method to calculate the field at point  xyz
   //
index 4d72000..d623cf9 100644 (file)
@@ -26,7 +26,8 @@ public:
   AliMagWrapCheb& operator=(const AliMagWrapCheb& maps);
   virtual ~AliMagWrapCheb();
   //
-  virtual void Field(Float_t *x, Float_t *b)                    const;
+  virtual void Field(float *x, float *b)                        const;
+  virtual void Field(double *x, double *b)                      const;
   virtual void GetTPCInt(Float_t *xyz, Float_t *b)              const;
   virtual void GetTPCIntCyl(Float_t *rphiz, Float_t *b)         const;
   //