[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(const AliMagF *bField,
			       Double_t driftVelocity,
			       Int_t nx,Int_t ny,Int_t nz)
  : fkNX(nx),fkNY(ny),fkNZ(nz),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.) {
  //
  // 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.
  //
  fDriftVelocity=driftVelocity;
  ConstructCommon(0,bField);
}

AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
			       Double_t driftVelocity) {
  //
  // The constructor. One has to supply a field map and an (initial)
  // drift velocity.
  //
  fDriftVelocity=driftVelocity;

  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[] fMeanBx;
  delete[] fMeanBy;
}

void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
  //
  // correct for the distortion
  //
  Double_t r=TMath::Sqrt(position[0]*position[0]+position[1]*position[1]);
  if (TMath::Abs(position[2])>250.||r<90.||250.<r) {
    for (Int_t i=0;i<3;++i) corrected[i]=position[i];
  }
  else {
    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*fMeanBz;

    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)
  //
  fMeanBx=new Double_t[fkNX*fkNY*fkNZ];
  fMeanBy=new Double_t[fkNX*fkNY*fkNZ];

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