/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $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 // Author: Andreas Morsch //----------------------------------------------------------------------- #include #include #include "AliFieldMap.h" ClassImp(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 // SetWriteEnable(); } //_______________________________________________________________________ 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 // ReadField(); SetWriteEnable(); } //_______________________________________________________________________ AliFieldMap::~AliFieldMap() { // // Destructor // delete fB; } //_______________________________________________________________________ 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"); 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()); fname = gSystem->ExpandPathName(fTitle.Data()); magfile = fopen(fname,"r"); delete [] fname; fscanf(magfile,"%d %d %d %f %f %f %f %f %f", &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel); fXdeli = 1./fXdel; fYdeli = 1./fYdel; fZdeli = 1./fZdel; fXend = fXbeg + (fXn-1)*fXdel; fYend = fYbeg + (fYn-1)*fYdel; fZend = fZbeg + (fZn-1)*fZdel; Int_t nDim = fXn*fYn*fZn; // Float_t x,y,z,b; fB = new TVector(3*nDim); if (magfile) { for (ix = 0; ix < fXn; ix++) { ipx=ix*3*(fZn*fYn); for (iy = 0; iy < fYn; iy++) { ipy=ipx+iy*3*fZn; 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) { fscanf(magfile," %f %f %f", &bx, &by, &bz); } 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); (*fB)(ipz+2) = 10.*bz; (*fB)(ipz+1) = 10.*by; (*fB)(ipz ) = 10.*bx; } //iz } // iy } // ix /* // // this part for interpolation between z = 700 and 720 cm to get // z = 710 cm // for (ix = 0; ix < fXn; ix++) { ipx=ix*3*(fZn*fYn); for (iy = 0; iy < fYn; iy++) { ipy=ipx+iy*3*fZn; Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.; Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.; Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.; ipz=ipy+3; (*fB)(ipz+2) = bzz; (*fB)(ipz+1) = byy; (*fB)(ipz ) = bxx; } // iy } // ix */ } else { printf("%s: File %s not found !\n",ClassName(),fTitle.Data()); exit(1); } // if mafile } //_______________________________________________________________________ // 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; // // // xl[0]=TMath::Abs(x[0])-fXbeg; // xl[1]=TMath::Abs(x[1])-fYbeg; // xl[2]=x[2]-fZbeg; // 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); // 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; // } //_______________________________________________________________________ void AliFieldMap::Field(Float_t *x, Float_t *b) { // // 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; } //_______________________________________________________________________ void AliFieldMap::Copy(AliFieldMap & /* magf */) const { // // Copy *this onto magf -- Not implemented // Fatal("Copy","Not implemented!\n"); } //_______________________________________________________________________ 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; } }