[u/mrichter/AliRoot.git] / TPC / AliTPCExBFirst.cxx

#include "TMath.h"
#include "TTreeStream.h"
#include "AliTPCExBFirst.h"

ClassImp(AliTPCExBFirst)

const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
const Double_t AliTPCExBFirst::fgkDriftField=40.e3;

AliTPCExBFirst::AliTPCExBFirst()
  : fDriftVelocity(0),
    fkNX(0),fkNY(0),fkNZ(0),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNMean(0),
    fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
  //
  // purely for I/O
  //
}

AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
			       Double_t driftVelocity,
			       Int_t nx,Int_t ny,Int_t nz)
  : fDriftVelocity(driftVelocity),
    fkNX(nx),fkNY(ny),fkNZ(nz),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNMean(0),
    fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
  //
  // The constructor. One has to supply a magnetic field and an (initial)
  // drift velocity. Since some kind of lookuptable is created the
  // number of its meshpoints can be supplied.
  //
  ConstructCommon(0,bField);
}

AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
			       Double_t driftVelocity) 
  : fDriftVelocity(driftVelocity),
    fkNX(0),fkNY(0),fkNZ(0),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNMean(0),
    fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
  //
  // The constructor. One has to supply a field map and an (initial)
  // drift velocity.
  //

  fkXMin=bFieldMap->Xmin()
    -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
    *bFieldMap->DelX();
  fkXMax=bFieldMap->Xmax()
    -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
    *bFieldMap->DelX();
  fkYMin=bFieldMap->Ymin()
    -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
    *bFieldMap->DelY();
  fkYMax=bFieldMap->Ymax()
    -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
    *bFieldMap->DelY();
  fkZMin=bFieldMap->Zmin()
    -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
    *bFieldMap->DelZ();
  fkZMax=bFieldMap->Zmax()
    -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
    *bFieldMap->DelZ();

  fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
  fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
  fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);

  ConstructCommon(bFieldMap,0);
}

AliTPCExBFirst::~AliTPCExBFirst() { 
  //
  // destruct the poor object.
  //
  delete[] fkMeanBx;
  delete[] fkMeanBy;
}

void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
  //
  // correct for the distortion
  //
  Double_t Bx,By;
  GetMeanFields(position[0],position[1],position[2],&Bx,&By);
  if (position[2]>0.) {
    Double_t Bxe,Bye;
    GetMeanFields(position[0],position[1],250.,&Bxe,&Bye);
    if (position[2]!=250.) {
      Bx=(500.*Bxe-(position[2]+250.)*Bx)/(250.-position[2]);
      By=(500.*Bye-(position[2]+250.)*By)/(250.-position[2]);
    }
    else {
      Bx=Bxe;
      By=Bye;
    }
  }
  
  Double_t mu=fDriftVelocity/fgkDriftField;
  Double_t wt=mu*fkMeanBz;
  
  corrected[0]=mu*(wt*Bx-By)/(1.+wt*wt);
  corrected[1]=mu*(wt*By+Bx)/(1.+wt*wt);
  
  if (position[2]>0.) {
    corrected[0]*=(250.-position[2]);
    corrected[1]*=(250.-position[2]);
  }
  else {
    corrected[0]*=(-250.-position[2]);
    corrected[1]*=(-250.-position[2]);
  }

  corrected[0]=position[0]-corrected[0];
  corrected[1]=position[1]-corrected[1];
  corrected[2]=position[2];
}

void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
  //
  // well, as the name sais...
  //
  TTreeSRedirector ts(fileName);
  Double_t x[3];
  for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
    for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
      for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
	Double_t d[3];
	Correct(x,d);
	Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
	Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
	Double_t dr=r-rd;
	Double_t phi=TMath::ATan2(x[0],x[1]);
	Double_t phid=TMath::ATan2(d[0],d[1]);
	Double_t dphi=phi-phid;
	if (dphi<0.) dphi+=TMath::TwoPi();
	if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
	Double_t drphi=r*dphi;
	Double_t dx=x[0]-d[0];
	Double_t dy=x[1]-d[1];
	Double_t dz=x[2]-d[2];
	ts<<"positions"
	  <<"x0="<<x[0]
	  <<"x1="<<x[1]
	  <<"x2="<<x[2]
	  <<"dx="<<dx
	  <<"dy="<<dy
	  <<"dz="<<dz
	  <<"r="<<r
	  <<"phi="<<phi
	  <<"dr="<<dr
	  <<"drphi="<<drphi
	  <<"\n";
      }
}

void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
				     const AliMagF *bField) {
  //
  // THIS IS PRIVATE! (a helper for the constructor)
  //
  fkNMean=fkNX*fkNY*fkNZ;
  fkMeanBx=new Float_t[fkNMean];
  fkMeanBy=new Float_t[fkNMean];

  Double_t x[3];
  Double_t nBz=0;
  fkMeanBz=0.;
  for (int i=0;i<fkNX;++i) {
    x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
    for (int j=0;j<fkNY;++j) {
      x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
      Double_t R=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
      Double_t Bx=0.,By=0.;
      for (int k=0;k<fkNZ;++k) {
	x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
	Float_t B[3];
	// the x is not const in the Field function...
	Float_t xt[3];
	for (int l=0;l<3;++l) xt[l]=x[l];
	// that happens due to the lack of a sophisticated class design:
	if (bFieldMap!=0)
	  bFieldMap->Field(xt,B);
	else 
	  bField->Field(xt,B);
	Bx+=B[0]/10.;
	By+=B[1]/10.;
	/*old
	fkMeanBx[(k*fkNY+j)*fkNX+i]=Bx/(k+1);
	fkMeanBy[(k*fkNY+j)*fkNX+i]=By/(k+1);
	*/
	fkMeanBx[(k*fkNY+j)*fkNX+i]=Bx;
	fkMeanBy[(k*fkNY+j)*fkNX+i]=By;
	if (90.<=R&&R<=250.) {
	  fkMeanBz+=B[2]/10.;
	  ++nBz;
	}
      }
    }
  }
  fkMeanBz/=nBz;
}

void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
				   Double_t *Bx,Double_t *By) const {
  //
  // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
  //
  Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
  Int_t xi1=static_cast<Int_t>(x);
  xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
  Int_t xi2=xi1+1;
  Double_t dx=(x-xi1);
  Double_t dx1=(xi2-x);

  Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
  Int_t yi1=static_cast<Int_t>(y);
  yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
  Int_t yi2=yi1+1;
  Double_t dy=(y-yi1);
  Double_t dy1=(yi2-y);
  
  Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
  Int_t zi1=static_cast<Int_t>(z);
  zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
  Int_t zi2=zi1+1;
  Double_t dz=(z-zi1);

  double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
	    +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
            +fkMeanBx[yi1*fkNX+xi2]*dx *dy
            +fkMeanBx[yi2*fkNX+xi2]*dx *dy1;
  double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
	    +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
            +fkMeanBy[yi1*fkNX+xi2]*dx *dy
            +fkMeanBy[yi2*fkNX+xi2]*dx *dy1;
  Int_t zi0=zi1-1;
  double snmx,snmy;
  if (zi0>=0) {
    snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
        +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
        +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
        +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
    snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
        +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
        +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
        +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
  }
  else
    snmx=snmy=0.;
  double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 	    +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
            +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
            +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
  double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 	    +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
            +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
            +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
  double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
	     +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
             +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
             +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
  double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
	     +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
             +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
             +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
  *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
  *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
  //TODO: make this nice
  if (TMath::Abs(z)>0.001) {
    *Bx/=z;
    *By/=z;
  }
}
Commit	Line	Data
faf93237	1	#include "TMath.h"
	2	#include "TTreeStream.h"
	3	#include "AliTPCExBFirst.h"
	4
	5	ClassImp(AliTPCExBFirst)
	6
	7	const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
	8	const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
	9
481f877b	10	AliTPCExBFirst::AliTPCExBFirst()
	11	: fDriftVelocity(0),
	12	fkNX(0),fkNY(0),fkNZ(0),
	13	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	14	fkZMin(-250.),fkZMax(250.),
	15	fkNMean(0),
	16	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
	17	//
	18	// purely for I/O
	19	//
	20	}
	21
faf93237	22	AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
	23	Double_t driftVelocity,
	24	Int_t nx,Int_t ny,Int_t nz)
481f877b	25	: fDriftVelocity(driftVelocity),
481f877b	26	fkNX(nx),fkNY(ny),fkNZ(nz),
faf93237	27	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b	28	fkZMin(-250.),fkZMax(250.),
	29	fkNMean(0),
	30	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237	31	//
	32	// The constructor. One has to supply a magnetic field and an (initial)
	33	// drift velocity. Since some kind of lookuptable is created the
	34	// number of its meshpoints can be supplied.
	35	//
faf93237	36	ConstructCommon(0,bField);
	37	}
	38
	39	AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
481f877b	40	Double_t driftVelocity)
	41	: fDriftVelocity(driftVelocity),
	42	fkNX(0),fkNY(0),fkNZ(0),
	43	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	44	fkZMin(-250.),fkZMax(250.),
	45	fkNMean(0),
	46	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237	47	//
	48	// The constructor. One has to supply a field map and an (initial)
	49	// drift velocity.
	50	//
faf93237	51
	52	fkXMin=bFieldMap->Xmin()
	53	-TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
	54	*bFieldMap->DelX();
	55	fkXMax=bFieldMap->Xmax()
	56	-TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
	57	*bFieldMap->DelX();
	58	fkYMin=bFieldMap->Ymin()
	59	-TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
	60	*bFieldMap->DelY();
	61	fkYMax=bFieldMap->Ymax()
	62	-TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
	63	*bFieldMap->DelY();
	64	fkZMin=bFieldMap->Zmin()
	65	-TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
	66	*bFieldMap->DelZ();
	67	fkZMax=bFieldMap->Zmax()
	68	-TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
	69	*bFieldMap->DelZ();
	70
	71	fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
	72	fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
	73	fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
	74
	75	ConstructCommon(bFieldMap,0);
	76	}
	77
	78	AliTPCExBFirst::~AliTPCExBFirst() {
	79	//
	80	// destruct the poor object.
	81	//
481f877b	82	delete[] fkMeanBx;
481f877b	83	delete[] fkMeanBy;
faf93237	84	}
	85
	86	void AliTPCExBFirst::Correct(const Double_t position,Double_t corrected) {
	87	//
	88	// correct for the distortion
	89	//
481f877b	90	Double_t Bx,By;
	91	GetMeanFields(position[0],position[1],position[2],&Bx,&By);
	92	if (position[2]>0.) {
	93	Double_t Bxe,Bye;
	94	GetMeanFields(position[0],position[1],250.,&Bxe,&Bye);
	95	if (position[2]!=250.) {
	96	Bx=(500.Bxe-(position[2]+250.)Bx)/(250.-position[2]);
	97	By=(500.Bye-(position[2]+250.)By)/(250.-position[2]);
faf93237	98	}
faf93237	99	else {
481f877b	100	Bx=Bxe;
481f877b	101	By=Bye;
faf93237	102	}
faf93237	103	}
481f877b	104
	105	Double_t mu=fDriftVelocity/fgkDriftField;
	106	Double_t wt=mu*fkMeanBz;
	107
	108	corrected[0]=mu(wtBx-By)/(1.+wt*wt);
	109	corrected[1]=mu(wtBy+Bx)/(1.+wt*wt);
	110
	111	if (position[2]>0.) {
	112	corrected[0]*=(250.-position[2]);
	113	corrected[1]*=(250.-position[2]);
	114	}
	115	else {
	116	corrected[0]*=(-250.-position[2]);
	117	corrected[1]*=(-250.-position[2]);
	118	}
	119
	120	corrected[0]=position[0]-corrected[0];
	121	corrected[1]=position[1]-corrected[1];
	122	corrected[2]=position[2];
faf93237	123	}
	124
	125	void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
	126	//
	127	// well, as the name sais...
	128	//
	129	TTreeSRedirector ts(fileName);
	130	Double_t x[3];
	131	for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
	132	for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
	133	for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
	134	Double_t d[3];
	135	Correct(x,d);
	136	Double_t r=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
	137	Double_t rd=TMath::Sqrt(d[0]d[0]+d[1]d[1]);
	138	Double_t dr=r-rd;
	139	Double_t phi=TMath::ATan2(x[0],x[1]);
	140	Double_t phid=TMath::ATan2(d[0],d[1]);
	141	Double_t dphi=phi-phid;
	142	if (dphi<0.) dphi+=TMath::TwoPi();
	143	if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
	144	Double_t drphi=r*dphi;
	145	Double_t dx=x[0]-d[0];
	146	Double_t dy=x[1]-d[1];
	147	Double_t dz=x[2]-d[2];
	148	ts<<"positions"
	149	<<"x0="<<x[0]
	150	<<"x1="<<x[1]
	151	<<"x2="<<x[2]
	152	<<"dx="<<dx
	153	<<"dy="<<dy
	154	<<"dz="<<dz
	155	<<"r="<<r
	156	<<"phi="<<phi
	157	<<"dr="<<dr
	158	<<"drphi="<<drphi
	159	<<"\n";
	160	}
	161	}
	162
	163	void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
	164	const AliMagF *bField) {
	165	//
	166	// THIS IS PRIVATE! (a helper for the constructor)
	167	//
481f877b	168	fkNMean=fkNXfkNYfkNZ;
0eab9089	169	fkMeanBx=new Float_t[fkNMean];
0eab9089	170	fkMeanBy=new Float_t[fkNMean];
faf93237	171
	172	Double_t x[3];
	173	Double_t nBz=0;
481f877b	174	fkMeanBz=0.;
faf93237	175	for (int i=0;i<fkNX;++i) {
	176	x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
	177	for (int j=0;j<fkNY;++j) {
	178	x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
	179	Double_t R=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
	180	Double_t Bx=0.,By=0.;
	181	for (int k=0;k<fkNZ;++k) {
	182	x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
	183	Float_t B[3];
	184	// the x is not const in the Field function...
	185	Float_t xt[3];
	186	for (int l=0;l<3;++l) xt[l]=x[l];
	187	// that happens due to the lack of a sophisticated class design:
	188	if (bFieldMap!=0)
	189	bFieldMap->Field(xt,B);
	190	else
	191	bField->Field(xt,B);
	192	Bx+=B[0]/10.;
	193	By+=B[1]/10.;
	194	/*old
481f877b	195	fkMeanBx[(kfkNY+j)fkNX+i]=Bx/(k+1);
481f877b	196	fkMeanBy[(kfkNY+j)fkNX+i]=By/(k+1);
faf93237	197	*/
481f877b	198	fkMeanBx[(kfkNY+j)fkNX+i]=Bx;
481f877b	199	fkMeanBy[(kfkNY+j)fkNX+i]=By;
faf93237	200	if (90.<=R&&R<=250.) {
481f877b	201	fkMeanBz+=B[2]/10.;
faf93237	202	++nBz;
	203	}
	204	}
	205	}
	206	}
481f877b	207	fkMeanBz/=nBz;
faf93237	208	}
	209
	210	void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
	211	Double_t Bx,Double_t By) const {
	212	//
	213	// THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
	214	//
	215	Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
	216	Int_t xi1=static_cast<Int_t>(x);
	217	xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
	218	Int_t xi2=xi1+1;
	219	Double_t dx=(x-xi1);
	220	Double_t dx1=(xi2-x);
	221
	222	Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
	223	Int_t yi1=static_cast<Int_t>(y);
	224	yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
	225	Int_t yi2=yi1+1;
	226	Double_t dy=(y-yi1);
	227	Double_t dy1=(yi2-y);
	228
	229	Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
	230	Int_t zi1=static_cast<Int_t>(z);
	231	zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
	232	Int_t zi2=zi1+1;
	233	Double_t dz=(z-zi1);
	234
481f877b	235	double s0x=fkMeanBx[yi1fkNX+xi1]dx1*dy1
	236	+fkMeanBx[yi2fkNX+xi1]dx1*dy
	237	+fkMeanBx[yi1fkNX+xi2]dx *dy
	238	+fkMeanBx[yi2fkNX+xi2]dx *dy1;
	239	double s0y=fkMeanBy[yi1fkNX+xi1]dx1*dy1
	240	+fkMeanBy[yi2fkNX+xi1]dx1*dy
	241	+fkMeanBy[yi1fkNX+xi2]dx *dy
	242	+fkMeanBy[yi2fkNX+xi2]dx *dy1;
faf93237	243	Int_t zi0=zi1-1;
	244	double snmx,snmy;
	245	if (zi0>=0) {
481f877b	246	snmx=fkMeanBx[(zi0fkNY+yi1)fkNX+xi1]dx1dy1
	247	+fkMeanBx[(zi0fkNY+yi2)fkNX+xi1]dx1dy
	248	+fkMeanBx[(zi0fkNY+yi1)fkNX+xi2]dx dy
	249	+fkMeanBx[(zi0fkNY+yi2)fkNX+xi2]dx dy1;
	250	snmy=fkMeanBy[(zi0fkNY+yi1)fkNX+xi1]dx1dy1
	251	+fkMeanBy[(zi0fkNY+yi2)fkNX+xi1]dx1dy
	252	+fkMeanBy[(zi0fkNY+yi1)fkNX+xi2]dx dy
	253	+fkMeanBy[(zi0fkNY+yi2)fkNX+xi2]dx dy1;
faf93237	254	}
	255	else
	256	snmx=snmy=0.;
481f877b	257	double snx=fkMeanBx[(zi1fkNY+yi1)fkNX+xi1]dx1dy1
	258	+fkMeanBx[(zi1fkNY+yi2)fkNX+xi1]dx1dy
	259	+fkMeanBx[(zi1fkNY+yi1)fkNX+xi2]dx dy
	260	+fkMeanBx[(zi1fkNY+yi2)fkNX+xi2]dx dy1;
	261	double sny=fkMeanBy[(zi1fkNY+yi1)fkNX+xi1]dx1dy1
	262	+fkMeanBy[(zi1fkNY+yi2)fkNX+xi1]dx1dy
	263	+fkMeanBy[(zi1fkNY+yi1)fkNX+xi2]dx dy
	264	+fkMeanBy[(zi1fkNY+yi2)fkNX+xi2]dx dy1;
	265	double snpx=fkMeanBx[(zi2fkNY+yi1)fkNX+xi1]dx1dy1
	266	+fkMeanBx[(zi2fkNY+yi2)fkNX+xi1]dx1dy
	267	+fkMeanBx[(zi2fkNY+yi1)fkNX+xi2]dx dy
	268	+fkMeanBx[(zi2fkNY+yi2)fkNX+xi2]dx dy1;
	269	double snpy=fkMeanBy[(zi2fkNY+yi1)fkNX+xi1]dx1dy1
	270	+fkMeanBy[(zi2fkNY+yi2)fkNX+xi1]dx1dy
	271	+fkMeanBy[(zi2fkNY+yi1)fkNX+xi2]dx dy
	272	+fkMeanBy[(zi2fkNY+yi2)fkNX+xi2]dx dy1;
faf93237	273	Bx=0.5(((snpx-2.snx+snmx)dz+2.(snx-snmx))dz+snx-s0x+snmx);
	274	By=0.5(((snpy-2.sny+snmy)dz+2.(sny-snmy))dz+sny-s0y+snmy);
	275	//TODO: make this nice
	276	if (TMath::Abs(z)>0.001) {
	277	*Bx/=z;
	278	*By/=z;
	279	}
	280	}