* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-*/
+/* $Id$ */
+//-----------------------------------------------------------------------
//
+// 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"
-#include "TSystem.h"
ClassImp(AliFieldMap)
-//________________________________________
-AliFieldMap::AliFieldMap()
+//_______________________________________________________________________
+AliFieldMap::AliFieldMap():
+ fXbeg(0),
+ fYbeg(0),
+ fZbeg(0),
+ fXend(0),
+ fYend(0),
+ fZend(0),
+ fXdel(0),
+ fYdel(0),
+ fZdel(0),
+ fXdeli(0),
+ fYdeli(0),
+ fZdeli(0),
+ fXn(0),
+ fYn(0),
+ fZn(0),
+ fWriteEnable(0),
+ fB(0)
{
//
// Standard constructor
//
- fB = 0;
+ SetWriteEnable();
}
-AliFieldMap::AliFieldMap(const char *name, const char *title)
- : TNamed(name,title)
+//_______________________________________________________________________
+AliFieldMap::AliFieldMap(const char *name, const char *title):
+ TNamed(name,title),
+ fXbeg(0),
+ fYbeg(0),
+ fZbeg(0),
+ fXend(0),
+ fYend(0),
+ fZend(0),
+ fXdel(0),
+ fYdel(0),
+ fZdel(0),
+ fXdeli(0),
+ fYdeli(0),
+ fZdeli(0),
+ fXn(0),
+ fYn(0),
+ fZn(0),
+ fWriteEnable(0),
+ fB(0)
{
//
// Standard constructor
//
- fB = 0;
ReadField();
+ SetWriteEnable();
}
+//_______________________________________________________________________
AliFieldMap::~AliFieldMap()
{
-//
-// Destructor
-//
+ //
+ // Destructor
+ //
delete fB;
}
-//________________________________________
-AliFieldMap::AliFieldMap(const AliFieldMap &map)
+//_______________________________________________________________________
+AliFieldMap::AliFieldMap(const AliFieldMap &map):
+ TNamed(map),
+ fXbeg(0),
+ fYbeg(0),
+ fZbeg(0),
+ fXend(0),
+ fYend(0),
+ fZend(0),
+ fXdel(0),
+ fYdel(0),
+ fZdel(0),
+ fXdeli(0),
+ fYdeli(0),
+ fZdeli(0),
+ fXn(0),
+ fYn(0),
+ fZn(0),
+ fWriteEnable(0),
+ fB(0)
{
//
// Copy constructor
map.Copy(*this);
}
-//________________________________________
+//_______________________________________________________________________
void AliFieldMap::ReadField()
{
//
// Method to read the magnetic field map from file
//
FILE* magfile;
- FILE* endf = fopen("end.table", "r");
- FILE* out = fopen("out", "w");
+ // FILE* endf = fopen("end.table", "r");
+ // FILE* out = fopen("out", "w");
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");
Int_t nDim = fXn*fYn*fZn;
-// Float_t x,y,z,b;
+ // Float_t x,y,z,b;
fB = new TVector(3*nDim);
if (magfile) {
for (iz = 0; iz < fZn; iz++) {
ipz=ipy+iz*3;
- if (iz == -1)
- fscanf(endf,"%f %f %f", &bx,&by,&bz);
- else if (iz > -1)
+ if (iz == -1) {
+// fscanf(endf,"%f %f %f", &bx,&by,&bz);
+ } else if (iz > -1) {
fscanf(magfile," %f %f %f", &bx, &by, &bz);
- else
+ } else {
continue;
-
+ }
+
// fscanf(magfile,"%f %f %f %f %f %f %f ",
// &x, &y, &z, &bx,&by,&bz, &b);
// fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
} // if mafile
}
-void AliFieldMap::Field(Float_t *x, Float_t *b)
+//_______________________________________________________________________
+void AliFieldMap::Field(float *x, float *b) const
{
-//
-// Use simple interpolation to obtain field at point x
-//
+ //
+ // 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
+ //
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]=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;
+ hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
ratx=hix-int(hix);
ix=int(hix);
- hiy=xl[1]*fYdeli;
+ hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
raty=hiy-int(hiy);
iy=int(hiy);
- hiz=xl[2]*fZdeli;
+ hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
ratz=hiz-int(hiz);
iz=int(hiz);
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;
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!");
}
-//________________________________________
+//_______________________________________________________________________
AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
{
magf.Copy(*this);
return *this;
}
+
+//_______________________________________________________________________
+void AliFieldMap::Streamer(TBuffer &R__b)
+{
+ // Stream an object of class AliFieldMap.
+ TVector* save = 0;
+
+ if (R__b.IsReading()) {
+ AliFieldMap::Class()->ReadBuffer(R__b, this);
+ } else {
+ if (!fWriteEnable) {
+ save = fB;
+ fB = 0;
+ }
+ AliFieldMap::Class()->WriteBuffer(R__b, this);
+ if (!fWriteEnable) fB = save;
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff