[u/mrichter/AliRoot.git] / STEER / AliFieldMap.cxx

/**************************************************************************
 * 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 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 "AliLog.h"
#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;
  AliInfo(Form("Reading Magnetic Field Map %s from file %s",
	       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 *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
  //
    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(TObject & /* magf */) const
{
  //
  // Copy *this onto magf -- Not implemented
  //
  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;
    }
}
Commit	Line	Data
84737f5e	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
acd84897	16	/* $Id$ */
84737f5e	17
116cbefd	18	//-----------------------------------------------------------------------
90dd3681	19	//
af7ba10c	20	// Class to handle the field
	21	// I/O and interpolation
	22	// of the field map to return the value
	23	// of the magnetic field in an arbitrary position
84737f5e	24	// Author: Andreas Morsch <andreas.morsch@cern.ch>
90dd3681	25	//
116cbefd	26	//-----------------------------------------------------------------------
84737f5e	27
b09247a2	28	#include <cstdlib>
a1e17193	29	#include <TClass.h>
116cbefd	30	#include <TSystem.h>
116cbefd	31
21bf7095	32	#include "AliLog.h"
84737f5e	33	#include "AliFieldMap.h"
84737f5e	34
	35	ClassImp(AliFieldMap)
	36
e2afb3b6	37	//_______________________________________________________________________
	38	AliFieldMap::AliFieldMap():
	39	fXbeg(0),
	40	fYbeg(0),
	41	fZbeg(0),
	42	fXend(0),
	43	fYend(0),
	44	fZend(0),
	45	fXdel(0),
	46	fYdel(0),
	47	fZdel(0),
	48	fXdeli(0),
	49	fYdeli(0),
	50	fZdeli(0),
	51	fXn(0),
	52	fYn(0),
	53	fZn(0),
	54	fWriteEnable(0),
	55	fB(0)
84737f5e	56	{
	57	//
	58	// Standard constructor
	59	//
7b6cddfa	60	SetWriteEnable();
84737f5e	61	}
84737f5e	62
e2afb3b6	63	//_______________________________________________________________________
	64	AliFieldMap::AliFieldMap(const char name, const char title):
	65	TNamed(name,title),
	66	fXbeg(0),
	67	fYbeg(0),
	68	fZbeg(0),
	69	fXend(0),
	70	fYend(0),
	71	fZend(0),
	72	fXdel(0),
	73	fYdel(0),
	74	fZdel(0),
	75	fXdeli(0),
	76	fYdeli(0),
	77	fZdeli(0),
	78	fXn(0),
	79	fYn(0),
	80	fZn(0),
	81	fWriteEnable(0),
	82	fB(0)
84737f5e	83	{
	84	//
	85	// Standard constructor
	86	//
84737f5e	87	ReadField();
7b6cddfa	88	SetWriteEnable();
84737f5e	89	}
84737f5e	90
e2afb3b6	91	//_______________________________________________________________________
84737f5e	92	AliFieldMap::~AliFieldMap()
84737f5e	93	{
e2afb3b6	94	//
	95	// Destructor
	96	//
84737f5e	97	delete fB;
	98	}
	99
e2afb3b6	100	//_______________________________________________________________________
	101	AliFieldMap::AliFieldMap(const AliFieldMap &map):
	102	TNamed(map),
	103	fXbeg(0),
	104	fYbeg(0),
	105	fZbeg(0),
	106	fXend(0),
	107	fYend(0),
	108	fZend(0),
	109	fXdel(0),
	110	fYdel(0),
	111	fZdel(0),
	112	fXdeli(0),
	113	fYdeli(0),
	114	fZdeli(0),
	115	fXn(0),
	116	fYn(0),
	117	fZn(0),
	118	fWriteEnable(0),
	119	fB(0)
84737f5e	120	{
	121	//
	122	// Copy constructor
	123	//
	124	map.Copy(*this);
	125	}
	126
e2afb3b6	127	//_______________________________________________________________________
84737f5e	128	void AliFieldMap::ReadField()
	129	{
	130	//
	131	// Method to read the magnetic field map from file
	132	//
	133	FILE* magfile;
e2afb3b6	134	// FILE* endf = fopen("end.table", "r");
e2afb3b6	135	// FILE* out = fopen("out", "w");
84737f5e	136
	137	Int_t ix, iy, iz, ipx, ipy, ipz;
	138	Float_t bx, by, bz;
	139	char *fname = 0;
21bf7095	140	AliInfo(Form("Reading Magnetic Field Map %s from file %s",
21bf7095	141	fName.Data(),fTitle.Data()));
84737f5e	142
	143	fname = gSystem->ExpandPathName(fTitle.Data());
	144	magfile = fopen(fname,"r");
	145	delete [] fname;
	146
	147	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	148	&fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
	149
	150	fXdeli = 1./fXdel;
	151	fYdeli = 1./fYdel;
	152	fZdeli = 1./fZdel;
	153	fXend = fXbeg + (fXn-1)*fXdel;
	154	fYend = fYbeg + (fYn-1)*fYdel;
	155	fZend = fZbeg + (fZn-1)*fZdel;
	156
	157	Int_t nDim = fXnfYnfZn;
	158
e2afb3b6	159	// Float_t x,y,z,b;
84737f5e	160
	161	fB = new TVector(3*nDim);
	162	if (magfile) {
	163	for (ix = 0; ix < fXn; ix++) {
	164	ipx=ix3(fZn*fYn);
	165	for (iy = 0; iy < fYn; iy++) {
	166	ipy=ipx+iy3fZn;
	167	for (iz = 0; iz < fZn; iz++) {
	168	ipz=ipy+iz*3;
	169
7b6cddfa	170	if (iz == -1) {
	171	// fscanf(endf,"%f %f %f", &bx,&by,&bz);
	172	} else if (iz > -1) {
84737f5e	173	fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa	174	} else {
84737f5e	175	continue;
7b6cddfa	176	}
7b6cddfa	177
84737f5e	178	// fscanf(magfile,"%f %f %f %f %f %f %f ",
	179	// &x, &y, &z, &bx,&by,&bz, &b);
	180	// fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
	181
	182	(fB)(ipz+2) = 10.bz;
	183	(fB)(ipz+1) = 10.by;
	184	(fB)(ipz ) = 10.bx;
	185	} //iz
	186	} // iy
	187	} // ix
	188	/*
	189	//
	190	// this part for interpolation between z = 700 and 720 cm to get
	191	// z = 710 cm
	192	//
	193	for (ix = 0; ix < fXn; ix++) {
	194	ipx=ix3(fZn*fYn);
	195	for (iy = 0; iy < fYn; iy++) {
	196	ipy=ipx+iy3fZn;
	197	Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
	198	Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
	199	Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;
	200	ipz=ipy+3;
	201	(*fB)(ipz+2) = bzz;
	202	(*fB)(ipz+1) = byy;
	203	(*fB)(ipz ) = bxx;
	204	} // iy
	205	} // ix
	206	*/
	207	} else {
	208	printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
	209	exit(1);
	210	} // if mafile
	211	}
	212
5d8718b8	213	//_______________________________________________________________________
ff66b122	214	void AliFieldMap::Field(float x, float b) const
	215	{
	216	//
	217	// Use simple interpolation to obtain field at point x
	218	//
	219	float ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	220	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	221	const float kone=1;
	222	Int_t ix, iy, iz;
	223	b[0]=b[1]=b[2]=0;
	224	//
	225
	226	xl[0] = x[0] - fXbeg;
	227	xl[1] = x[1] - fYbeg;
	228	xl[2] = x[2] - fZbeg;
	229
	230	hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
	231	ratx=hix-int(hix);
	232	ix=int(hix);
	233
	234	hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
	235	raty=hiy-int(hiy);
	236	iy=int(hiy);
	237
	238	hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
	239	ratz=hiz-int(hiz);
	240	iz=int(hiz);
	241
	242	ratx1=kone-ratx;
	243	raty1=kone-raty;
	244	ratz1=kone-ratz;
	245
	246	if (!fB) return;
	247
	248	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	249	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	250	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	251	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	252	bhz = blyhz raty1+bhyhz raty;
	253	blz = blylz raty1+bhylz raty;
	254	b[0] = blz ratz1+bhz ratz;
	255	//
	256	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	257	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	258	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	259	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	260	bhz = blyhz raty1+bhyhz raty;
	261	blz = blylz raty1+bhylz raty;
	262	b[1] = blz ratz1+bhz ratz;
	263	//
	264	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	265	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	266	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	267	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	268	bhz = blyhz raty1+bhyhz raty;
	269	blz = blylz raty1+bhylz raty;
	270	b[2] = blz ratz1+bhz ratz;
	271	}
	272
	273	//_______________________________________________________________________
	274	void AliFieldMap::Field(double x, double b) const
5d8718b8	275	{
	276	//
	277	// Use simple interpolation to obtain field at point x
90dd3681	278	//
	279	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	280	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	281	const Double_t kone=1;
	282	Int_t ix, iy, iz;
	283	b[0]=b[1]=b[2]=0;
	284	//
	285
415bda4f	286	xl[0] = x[0] - fXbeg;
	287	xl[1] = x[1] - fYbeg;
	288	xl[2] = x[2] - fZbeg;
90dd3681	289
cb3dd940	290	hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
90dd3681	291	ratx=hix-int(hix);
	292	ix=int(hix);
	293
cb3dd940	294	hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
90dd3681	295	raty=hiy-int(hiy);
	296	iy=int(hiy);
	297
cb3dd940	298	hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
90dd3681	299	ratz=hiz-int(hiz);
	300	iz=int(hiz);
	301
	302	ratx1=kone-ratx;
	303	raty1=kone-raty;
	304	ratz1=kone-ratz;
	305
39b41d10	306	if (!fB) return;
39b41d10	307
90dd3681	308	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	309	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	310	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	311	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	312	bhz = blyhz raty1+bhyhz raty;
	313	blz = blylz raty1+bhylz raty;
	314	b[0] = blz ratz1+bhz ratz;
	315	//
	316	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	317	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	318	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	319	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	320	bhz = blyhz raty1+bhyhz raty;
	321	blz = blylz raty1+bhylz raty;
	322	b[1] = blz ratz1+bhz ratz;
	323	//
	324	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	325	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	326	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	327	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	328	bhz = blyhz raty1+bhyhz raty;
	329	blz = blylz raty1+bhylz raty;
	330	b[2] = blz ratz1+bhz ratz;
84737f5e	331	}
84737f5e	332
e2afb3b6	333	//_______________________________________________________________________
6c4904c2	334	void AliFieldMap::Copy(TObject & /* magf */) const
84737f5e	335	{
	336	//
	337	// Copy *this onto magf -- Not implemented
	338	//
21bf7095	339	AliFatal("Not implemented!");
84737f5e	340	}
84737f5e	341
e2afb3b6	342	//_______________________________________________________________________
84737f5e	343	AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
	344	{
	345	magf.Copy(*this);
	346	return *this;
	347	}
7b6cddfa	348
e2afb3b6	349	//_______________________________________________________________________
7b6cddfa	350	void AliFieldMap::Streamer(TBuffer &R__b)
	351	{
	352	// Stream an object of class AliFieldMap.
f9aab227	353	TVector* save = 0;
b4cfe01e	354
	355	if (R__b.IsReading()) {
	356	AliFieldMap::Class()->ReadBuffer(R__b, this);
	357	} else {
	358	if (!fWriteEnable) {
	359	save = fB;
	360	fB = 0;
	361	}
	362	AliFieldMap::Class()->WriteBuffer(R__b, this);
	363	if (!fWriteEnable) fB = save;
	364	}
7b6cddfa	365	}