virtual void Field(double*,double*) version.
Some auxiliary functions with (float*) arguments are converted
to templates.
~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);
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];}
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
};
//__________________________________________________________________________________________
-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;
}
//__________________________________________________________________________________________
-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);
}
//__________________________________________________________________________________________
-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);
}
//__________________________________________________________________________________________
-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;
**************************************************************************/
#include <cstdlib>
-#include "AliCheb3DCalc.h"
#include <TSystem.h>
+#include "AliCheb3DCalc.h"
ClassImp(AliCheb3DCalc)
//
}
-//__________________________________________________________________________________________
-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--;) {
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);
//
}
{
// 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;
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);
//
}
//
return b0 - x*b1 - ddcf0/2;
}
-
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;
//
//
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
}
//_______________________________________________________________________
-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
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));
}
}
//_______________________________________________________________________
-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
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;}
}
//________________________________________
-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
}
//___________________________________________________
-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) {
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;}
}
//_______________________________________________________________________
-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
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 {
}
//__________________________________________________________________________________________
-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
//
}
+//__________________________________________________________________________________________
+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
{
}
//__________________________________________________________________________________________
-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);
//
}
{
// 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);
//
}
}
#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),
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);
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:
//
//__________________________________________________________________________________________
-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;}
}
//__________________________________________________________________________________________________
-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];
}
//__________________________________________________________________________________________________
-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];
}
//__________________________________________________________________________________________________
-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);
}
//__________________________________________________________________________________________________
-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];
}
//__________________________________________________________________________________________________
-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]);
}
//__________________________________________________________________________________________________
-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
}
//_______________________________________________________________________
-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
}
}
+//_______________________________________________________________________
+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
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 {
}
//_______________________________________________________________________
-void AliMagFMaps::Field(Float_t *x, Float_t *b) const
+void AliMagFMaps::Field(float *x, float *b) const
{
//
// Method to calculate the magnetic field
}
}
+//_______________________________________________________________________
+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
{
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;
}
//_______________________________________________________________________
-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
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;}
}
//_______________________________________________________________________
-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
//
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;
//
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff