/* $Id$ */
//-----------------------------------------------------------------------
-// Class to read the magnetic field map and to calculate
-// the magnetic field in a given point. For the moment uses a simple
-// linear interpolation between the nodes of a rectangular grid
+//
+// Class to handle the field
+// I/O and interpolation
+// of the field map to return the value
+// of the magnetic field in an arbitrary position
// Author: Andreas Morsch <andreas.morsch@cern.ch>
+//
//-----------------------------------------------------------------------
+#include <cstdlib>
+#include <TClass.h>
#include <TSystem.h>
-#include <TVector.h>
+#include "AliLog.h"
#include "AliFieldMap.h"
ClassImp(AliFieldMap)
Int_t ix, iy, iz, ipx, ipy, ipz;
Float_t bx, by, bz;
char *fname = 0;
- printf("%s: Reading Magnetic Field Map %s from file %s\n",
- ClassName(),fName.Data(),fTitle.Data());
+ AliInfo(Form("Reading Magnetic Field Map %s from file %s",
+ fName.Data(),fTitle.Data()));
fname = gSystem->ExpandPathName(fTitle.Data());
magfile = fopen(fname,"r");
}
//_______________________________________________________________________
-// void AliFieldMap::Field(Float_t *x, Float_t *b)
-// {
-// //
-// // Use simple interpolation to obtain field at point x
-// //
-// Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
-// bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
-// const Double_t kone=1;
-// Int_t ix, iy, iz;
-// b[0]=b[1]=b[2]=0;
-// //
+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]=TMath::Abs(x[0])-fXbeg;
-// xl[1]=TMath::Abs(x[1])-fYbeg;
-// xl[2]=x[2]-fZbeg;
+ xl[0] = x[0] - fXbeg;
+ xl[1] = x[1] - fYbeg;
+ xl[2] = x[2] - fZbeg;
-// hix=xl[0]*fXdeli;
-// ratx=hix-int(hix);
-// ix=int(hix);
+ hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
+ ratx=hix-int(hix);
+ ix=int(hix);
-// hiy=xl[1]*fYdeli;
-// raty=hiy-int(hiy);
-// iy=int(hiy);
+ hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
+ raty=hiy-int(hiy);
+ iy=int(hiy);
-// hiz=xl[2]*fZdeli;
-// ratz=hiz-int(hiz);
-// iz=int(hiz);
-
-// ratx1=kone-ratx;
-// raty1=kone-raty;
-// ratz1=kone-ratz;
+ 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;
+}
-// 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(Float_t *x, Float_t *b)
+void AliFieldMap::Field(double *x, double *b) const
{
//
// Use simple interpolation to obtain field at point x
- // Faster version which uses direct access to the TVector
- // instead of Bx(...), By(...), Bz(...) interfaces.
- // The interpolation is done first in Z (4 points),
- // then in Y (2 points), and finally in X direction (one point).
-
- Double_t h=(x[2]-fZbeg)*fZdeli;
- Int_t i=int(h); // ix
- register Int_t index=i;
- Double_t ratz=h-i;
- Double_t ratz1=1-ratz;
-
- h=(TMath::Abs(x[1])-fYbeg)*fYdeli;
- i=int(h); // iy
- index+=(i*fZn);
- Double_t raty=h-i;
- Double_t raty1=1-raty;
-
- h=(TMath::Abs(x[0])-fXbeg)*fXdeli;
- i=int(h); // iz
- index+=(i*fZn*fYn);
- index*=3; //index now points to the beginning of the field data (ix,iy,iz)
- Double_t ratx=h-i;
- Double_t ratx1=1-ratx;
-
- // The way we retrive data from fB should be optimized
- // with respect to the CPU cache
-
- Float_t bx00=(*fB)[index++]*ratz1;
- Float_t by00=(*fB)[index++]*ratz1;
- Float_t bz00=(*fB)[index++]*ratz1;
- bx00+=(*fB)[index++]*ratz;
- by00+=(*fB)[index++]*ratz;
- bz00+=(*fB)[index]*ratz;
-
- index+=(3*fZn-5);
-
- Float_t bx01=(*fB)[index++]*ratz1;
- Float_t by01=(*fB)[index++]*ratz1;
- Float_t bz01=(*fB)[index++]*ratz1;
- bx01+=(*fB)[index++]*ratz;
- by01+=(*fB)[index++]*ratz;
- bz01+=(*fB)[index]*ratz;
-
- index+=(3*fYn*fZn-5);
-
- Float_t bx11=(*fB)[index++]*ratz1;
- Float_t by11=(*fB)[index++]*ratz1;
- Float_t bz11=(*fB)[index++]*ratz1;
- bx11+=(*fB)[index++]*ratz;
- by11+=(*fB)[index++]*ratz;
- bz11+=(*fB)[index]*ratz;
-
- index-=(3*fZn+5);
-
- Float_t bx10=(*fB)[index++]*ratz1;
- Float_t by10=(*fB)[index++]*ratz1;
- Float_t bz10=(*fB)[index++]*ratz1;
- bx10+=(*fB)[index++]*ratz;
- by10+=(*fB)[index++]*ratz;
- bz10+=(*fB)[index]*ratz;
-
- b[0]= (bx00*raty1+bx01*raty)*ratx1 +(bx10*raty1+bx11*raty)*ratx;
- b[1]= (by00*raty1+by01*raty)*ratx1 +(by10*raty1+by11*raty)*ratx;
- b[2]= (bz00*raty1+bz01*raty)*ratx1 +(bz10*raty1+bz11*raty)*ratx;
+ //
+ Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
+ bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
+ const Double_t 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::Copy(AliFieldMap & /* magf */) const
+void AliFieldMap::Copy(TObject & /* magf */) const
{
//
// Copy *this onto magf -- Not implemented
//
- Fatal("Copy","Not implemented!\n");
+ AliFatal("Not implemented!");
}
//_______________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff