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

#include "TMath.h"
#include "TTreeStream.h"
#include "AliTPCExBExact.h"

ClassImp(AliTPCExBExact)

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

AliTPCExBExact::AliTPCExBExact()
  : fDriftVelocity(0),
    fkMap(0),fkField(0),fkN(0),
    fkNX(0),fkNY(0),fkNZ(0),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNLook(0),fkLook(0) {
  //
  // purely for I/O
  //
}

AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
			       Double_t driftVelocity,
			       Int_t nx,Int_t ny,Int_t nz,Int_t n)
  : fDriftVelocity(driftVelocity),
    fkMap(0),fkField(bField),fkN(n),
    fkNX(nx),fkNY(ny),fkNZ(nz),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNLook(0),fkLook(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.
  // n sets the number of integration steps to be used when integrating
  // over the full drift length.
  //
  CreateLookupTable();
}

AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
			       Double_t driftVelocity,Int_t n) 
  : fDriftVelocity(driftVelocity),
    fkMap(bFieldMap),fkField(0),fkN(n), 
    fkNX(0),fkNY(0),fkNZ(0),
    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
    fkZMin(-250.),fkZMax(250.),
    fkNLook(0),fkLook(0) {
  //
  // The constructor. One has to supply a field map and an (initial)
  // drift velocity.
  // n sets the number of integration steps to be used when integrating
  // over the full drift length.
  //

  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();
  fkZMax=bFieldMap->Zmax()
    -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
    *bFieldMap->DelZ();
  fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!

  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);

  CreateLookupTable();
}

AliTPCExBExact::~AliTPCExBExact() {
  //
  // destruct the poor object.
  //
  delete[] fkLook;
}

void AliTPCExBExact::Correct(const Double_t *position, Double_t *corrected) {
  //
  // correct for the distortion
  //
  Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
  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=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
  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=position[2]/fkZMax*(fkNZ-1);
  Int_t side;
  if (z>0) {
    side=1;
  }
  else {
    z=-z;
    side=0;
  }
  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_t dz1=(zi2-z);
  
  for (int i=0;i<3;++i)
    corrected[i]
      =fkLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
      +fkLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
      +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
      +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
      +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
      +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
      +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
      +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
  //    corrected[2]=position[2];
}

void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
					     const char* fileName) {
  fkMap=bFieldMap;
  fkField=0;
  TestThisBeautifulObjectGeneric(fileName);
}

void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
					     const char* fileName) {
  fkField=bField;
  fkMap=0;
  TestThisBeautifulObjectGeneric(fileName);
}

void AliTPCExBExact::TestThisBeautifulObjectGeneric(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];
	Double_t dnl[3];
	Correct(x,d);
	CalculateDistortion(x,dnl);
	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];
	Double_t dnlx=x[0]-dnl[0];
	Double_t dnly=x[1]-dnl[1];
	Double_t dnlz=x[2]-dnl[2];
	ts<<"positions"
	  <<"x0="<<x[0]
	  <<"x1="<<x[1]
	  <<"x2="<<x[2]
	  <<"dx="<<dx
	  <<"dy="<<dy
	  <<"dz="<<dz
	  <<"dnlx="<<dnlx
	  <<"dnly="<<dnly
	  <<"dnlz="<<dnlz
	  <<"r="<<r
	  <<"phi="<<phi
	  <<"dr="<<dr
	  <<"drphi="<<drphi
	  <<"\n";
      }
}

void AliTPCExBExact::CreateLookupTable() {
  //
  // Helper function to fill the lookup table.
  //
  fkNLook=fkNX*fkNY*fkNZ*2*3;
  fkLook=new Double_t[fkNLook];
  Double_t x[3];
  for (int i=0;i<fkNX;++i) {
    x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
    for (int j=0;j<fkNY;++j) {
      x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
      for (int k=0;k<fkNZ;++k) {
	x[2]=1.*fkZMax/(fkNZ-1)*k;
	x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
	CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
	x[2]=-x[2];
	CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
      }
    }
  }
}

void AliTPCExBExact::GetE(Double_t *E,const Double_t *x) const {
  //
  // Helper function returning the E field in SI units (V/m).
  //
  E[0]=0.;
  E[1]=0.;
  E[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
}

void AliTPCExBExact::GetB(Double_t *B,const Double_t *x) const {
  //
  // Helper function returning the B field in SI units (T).
  //
  Float_t xm[3];
  // the beautiful m to cm (and the ugly "const_cast") and Double_t 
  // to Float_t read the NRs introduction!:
  for (int i=0;i<3;++i) xm[i]=x[i]*100.;
  Float_t Bf[3];
  if (fkMap!=0)
    fkMap->Field(xm,Bf);
  else
    fkField->Field(xm,Bf);
  for (int i=0;i<3;++i) B[i]=Bf[i]/10.;
}

void AliTPCExBExact::Motion(const Double_t *x,Double_t,
			    Double_t *dxdt) const {
  //
  // The differential equation of motion of the electrons.
  //
  const Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
  const Double_t tau2=tau*tau;
  Double_t E[3];
  Double_t B[3];
  GetE(E,x);
  GetB(B,x);
  Double_t wx=fgkEM*B[0];
  Double_t wy=fgkEM*B[1];
  Double_t wz=fgkEM*B[2];
  Double_t ex=fgkEM*E[0];
  Double_t ey=fgkEM*E[1];
  Double_t ez=fgkEM*E[2];
  Double_t w2=(wx*wx+wy*wy+wz*wz);
  dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
  dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
  dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
  Double_t fac=tau/(1.+w2*tau2);
  dxdt[0]*=fac;
  dxdt[1]*=fac;
  dxdt[2]*=fac;
}

void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
					 Double_t *dist) const {
  //
  // Helper function that calculates one distortion by integration
  // (only used to fill the lookup table).
  //
  const Double_t h=0.01*250./fDriftVelocity/fkN;
  Double_t t=0.;
  Double_t xt[3];
  Double_t xo[3];
  for (int i=0;i<3;++i)
    xo[i]=xt[i]=x0[i]*0.01;
  while (TMath::Abs(xt[2])<250.*0.01) {
    for (int i=0;i<3;++i)
      xo[i]=xt[i];
    DGLStep(xt,t,h);
    t+=h;
  }
  if (t!=0.) {
    Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
    dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
    dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
    //    dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
    dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
    dist[2]+=(x0[2]<0.?-1:1.)*250.;
  }
  else {
    dist[0]=x0[0];
    dist[1]=x0[1];
    dist[2]=x0[2];
  }
  // reverse the distortion, i.e. get the correction
  dist[0]=x0[0]-(dist[0]-x0[0]);
  dist[1]=x0[1]-(dist[1]-x0[1]);
}

void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
  //
  // An elementary integration step.
  // (simple Euler Method)
  //
  Double_t dxdt[3];
  Motion(x,t,dxdt);
  for (int i=0;i<3;++i)
    x[i]+=h*dxdt[i];

  /* suggestions about how to write it this way are welcome!
     void DGLStep(void (*f)(const Double_t *x,Double_t t,Double_t *dxdt),
                   Double_t *x,Double_t t,Double_t h,Int_t n) const;
  */

}
Commit	Line	Data
faf93237	1	#include "TMath.h"
	2	#include "TTreeStream.h"
	3	#include "AliTPCExBExact.h"
	4
	5	ClassImp(AliTPCExBExact)
	6
	7	const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
	8	const Double_t AliTPCExBExact::fgkDriftField=40.e3;
	9
481f877b	10	AliTPCExBExact::AliTPCExBExact()
	11	: fDriftVelocity(0),
	12	fkMap(0),fkField(0),fkN(0),
	13	fkNX(0),fkNY(0),fkNZ(0),
	14	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	15	fkZMin(-250.),fkZMax(250.),
	16	fkNLook(0),fkLook(0) {
	17	//
	18	// purely for I/O
	19	//
	20	}
	21
faf93237	22	AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
	23	Double_t driftVelocity,
	24	Int_t nx,Int_t ny,Int_t nz,Int_t n)
481f877b	25	: fDriftVelocity(driftVelocity),
481f877b	26	fkMap(0),fkField(bField),fkN(n),
faf93237	27	fkNX(nx),fkNY(ny),fkNZ(nz),
faf93237	28	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b	29	fkZMin(-250.),fkZMax(250.),
481f877b	30	fkNLook(0),fkLook(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	// n sets the number of integration steps to be used when integrating
	36	// over the full drift length.
	37	//
faf93237	38	CreateLookupTable();
	39	}
	40
	41	AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
	42	Double_t driftVelocity,Int_t n)
481f877b	43	: fDriftVelocity(driftVelocity),
	44	fkMap(bFieldMap),fkField(0),fkN(n),
	45	fkNX(0),fkNY(0),fkNZ(0),
	46	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	47	fkZMin(-250.),fkZMax(250.),
	48	fkNLook(0),fkLook(0) {
faf93237	49	//
	50	// The constructor. One has to supply a field map and an (initial)
	51	// drift velocity.
	52	// n sets the number of integration steps to be used when integrating
	53	// over the full drift length.
	54	//
faf93237	55
	56	fkXMin=bFieldMap->Xmin()
	57	-TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
	58	*bFieldMap->DelX();
	59	fkXMax=bFieldMap->Xmax()
	60	-TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
	61	*bFieldMap->DelX();
	62	fkYMin=bFieldMap->Ymin()
	63	-TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
	64	*bFieldMap->DelY();
	65	fkYMax=bFieldMap->Ymax()
	66	-TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
	67	*bFieldMap->DelY();
	68	fkZMax=bFieldMap->Zmax()
	69	-TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
	70	*bFieldMap->DelZ();
	71	fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!
	72
	73	fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
	74	fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
	75	fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
	76
	77	CreateLookupTable();
	78	}
	79
	80	AliTPCExBExact::~AliTPCExBExact() {
	81	//
	82	// destruct the poor object.
	83	//
481f877b	84	delete[] fkLook;
faf93237	85	}
faf93237	86
481f877b	87	void AliTPCExBExact::Correct(const Double_t position, Double_t corrected) {
	88	//
	89	// correct for the distortion
	90	//
	91	Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
	92	Int_t xi1=static_cast<Int_t>(x);
	93	xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
	94	Int_t xi2=xi1+1;
	95	Double_t dx=(x-xi1);
	96	Double_t dx1=(xi2-x);
	97
	98	Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
	99	Int_t yi1=static_cast<Int_t>(y);
	100	yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
	101	Int_t yi2=yi1+1;
	102	Double_t dy=(y-yi1);
	103	Double_t dy1=(yi2-y);
	104
	105	Double_t z=position[2]/fkZMax*(fkNZ-1);
	106	Int_t side;
	107	if (z>0) {
	108	side=1;
faf93237	109	}
faf93237	110	else {
481f877b	111	z=-z;
481f877b	112	side=0;
faf93237	113	}
481f877b	114	Int_t zi1=static_cast<Int_t>(z);
	115	zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
	116	Int_t zi2=zi1+1;
	117	Double_t dz=(z-zi1);
	118	Double_t dz1=(zi2-z);
	119
	120	for (int i=0;i<3;++i)
	121	corrected[i]
	122	=fkLook[(((xi1fkNY+yi1)fkNZ+zi1)2+side)3+i]dx1dy1*dz1
	123	+fkLook[(((xi1fkNY+yi1)fkNZ+zi2)2+side)3+i]dx1dy1*dz
	124	+fkLook[(((xi1fkNY+yi2)fkNZ+zi1)2+side)3+i]dx1dy *dz1
	125	+fkLook[(((xi1fkNY+yi2)fkNZ+zi2)2+side)3+i]dx1dy *dz
	126	+fkLook[(((xi2fkNY+yi2)fkNZ+zi1)2+side)3+i]dx dy *dz1
	127	+fkLook[(((xi2fkNY+yi2)fkNZ+zi2)2+side)3+i]dx dy *dz
	128	+fkLook[(((xi2fkNY+yi1)fkNZ+zi1)2+side)3+i]dx dy1*dz1
	129	+fkLook[(((xi2fkNY+yi1)fkNZ+zi2)2+side)3+i]dx dy1*dz ;
	130	// corrected[2]=position[2];
	131	}
	132
	133	void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
	134	const char* fileName) {
	135	fkMap=bFieldMap;
	136	fkField=0;
	137	TestThisBeautifulObjectGeneric(fileName);
	138	}
	139
	140	void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
	141	const char* fileName) {
	142	fkField=bField;
	143	fkMap=0;
	144	TestThisBeautifulObjectGeneric(fileName);
faf93237	145	}
faf93237	146
481f877b	147	void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) {
faf93237	148	//
	149	// well, as the name sais...
	150	//
	151	TTreeSRedirector ts(fileName);
	152	Double_t x[3];
	153	for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
	154	for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
	155	for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
	156	Double_t d[3];
	157	Double_t dnl[3];
	158	Correct(x,d);
	159	CalculateDistortion(x,dnl);
	160	Double_t r=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
	161	Double_t rd=TMath::Sqrt(d[0]d[0]+d[1]d[1]);
	162	Double_t dr=r-rd;
	163	Double_t phi=TMath::ATan2(x[0],x[1]);
	164	Double_t phid=TMath::ATan2(d[0],d[1]);
	165	Double_t dphi=phi-phid;
	166	if (dphi<0.) dphi+=TMath::TwoPi();
	167	if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
	168	Double_t drphi=r*dphi;
	169	Double_t dx=x[0]-d[0];
	170	Double_t dy=x[1]-d[1];
	171	Double_t dz=x[2]-d[2];
	172	Double_t dnlx=x[0]-dnl[0];
	173	Double_t dnly=x[1]-dnl[1];
	174	Double_t dnlz=x[2]-dnl[2];
	175	ts<<"positions"
	176	<<"x0="<<x[0]
	177	<<"x1="<<x[1]
	178	<<"x2="<<x[2]
	179	<<"dx="<<dx
	180	<<"dy="<<dy
	181	<<"dz="<<dz
	182	<<"dnlx="<<dnlx
	183	<<"dnly="<<dnly
	184	<<"dnlz="<<dnlz
	185	<<"r="<<r
	186	<<"phi="<<phi
	187	<<"dr="<<dr
	188	<<"drphi="<<drphi
	189	<<"\n";
	190	}
	191	}
	192
	193	void AliTPCExBExact::CreateLookupTable() {
	194	//
	195	// Helper function to fill the lookup table.
	196	//
481f877b	197	fkNLook=fkNXfkNYfkNZ23;
481f877b	198	fkLook=new Double_t[fkNLook];
faf93237	199	Double_t x[3];
	200	for (int i=0;i<fkNX;++i) {
	201	x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
	202	for (int j=0;j<fkNY;++j) {
	203	x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
	204	for (int k=0;k<fkNZ;++k) {
	205	x[2]=1.fkZMax/(fkNZ-1)k;
	206	x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
481f877b	207	CalculateDistortion(x,&fkLook[(((ifkNY+j)fkNZ+k)2+1)3]);
faf93237	208	x[2]=-x[2];
481f877b	209	CalculateDistortion(x,&fkLook[(((ifkNY+j)fkNZ+k)2+0)3]);
faf93237	210	}
	211	}
	212	}
	213	}
	214
	215	void AliTPCExBExact::GetE(Double_t E,const Double_t x) const {
	216	//
	217	// Helper function returning the E field in SI units (V/m).
	218	//
	219	E[0]=0.;
	220	E[1]=0.;
	221	E[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
	222	}
	223
	224	void AliTPCExBExact::GetB(Double_t B,const Double_t x) const {
	225	//
	226	// Helper function returning the B field in SI units (T).
	227	//
	228	Float_t xm[3];
	229	// the beautiful m to cm (and the ugly "const_cast") and Double_t
	230	// to Float_t read the NRs introduction!:
	231	for (int i=0;i<3;++i) xm[i]=x[i]*100.;
	232	Float_t Bf[3];
	233	if (fkMap!=0)
	234	fkMap->Field(xm,Bf);
	235	else
	236	fkField->Field(xm,Bf);
	237	for (int i=0;i<3;++i) B[i]=Bf[i]/10.;
	238	}
	239
	240	void AliTPCExBExact::Motion(const Double_t *x,Double_t,
	241	Double_t *dxdt) const {
	242	//
	243	// The differential equation of motion of the electrons.
	244	//
	245	const Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
	246	const Double_t tau2=tau*tau;
	247	Double_t E[3];
	248	Double_t B[3];
	249	GetE(E,x);
	250	GetB(B,x);
	251	Double_t wx=fgkEM*B[0];
	252	Double_t wy=fgkEM*B[1];
	253	Double_t wz=fgkEM*B[2];
	254	Double_t ex=fgkEM*E[0];
	255	Double_t ey=fgkEM*E[1];
	256	Double_t ez=fgkEM*E[2];
	257	Double_t w2=(wxwx+wywy+wz*wz);
	258	dxdt[0]=(1.+wxwxtau2)ex+(wztau+wxwytau2)ey+(-wytau+wxwztau2)*ez;
	259	dxdt[1]=(-wztau+wxwytau2)ex+(1.+wywytau2)ey+(wxtau+wywztau2)*ez;
	260	dxdt[2]=(wytau+wxwztau2)ex+(-wxtau+wywztau2)ey+(1.+wzwztau2)*ez;
	261	Double_t fac=tau/(1.+w2*tau2);
	262	dxdt[0]*=fac;
	263	dxdt[1]*=fac;
	264	dxdt[2]*=fac;
	265	}
	266
	267	void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
	268	Double_t *dist) const {
	269	//
	270	// Helper function that calculates one distortion by integration
	271	// (only used to fill the lookup table).
	272	//
	273	const Double_t h=0.01*250./fDriftVelocity/fkN;
274	Double_t t=0.;
275	Double_t xt[3];
276	Double_t xo[3];
277	for (int i=0;i<3;++i)
278	xo[i]=xt[i]=x0[i]*0.01;
279	while (TMath::Abs(xt[2])<250.*0.01) {
280	for (int i=0;i<3;++i)
281	xo[i]=xt[i];
282	DGLStep(xt,t,h);
283	t+=h;
284	}
285	if (t!=0.) {
286	Double_t p=((xt[2]<0.?-1.:1.)250.0.01-xo[2])/(xt[2]-xo[2]);
287	dist[0]=(xo[0]+p(xt[0]-xo[0]))100.;
288	dist[1]=(xo[1]+p(xt[1]-xo[1]))100.;
289	// dist[2]=(xo[2]+p(xt[2]-xo[2]))100.;
290	dist[2]=(x0[2]>0.?-1:1.)(t-h+ph)fDriftVelocity100.;
291	dist[2]+=(x0[2]<0.?-1:1.)*250.;
292	}
293	else {
294	dist[0]=x0[0];
295	dist[1]=x0[1];
296	dist[2]=x0[2];
297	}
298	// reverse the distortion, i.e. get the correction
299	dist[0]=x0[0]-(dist[0]-x0[0]);
300	dist[1]=x0[1]-(dist[1]-x0[1]);
301	}
302
303	void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
304	//
305	// An elementary integration step.
306	// (simple Euler Method)
307	//
308	Double_t dxdt[3];
309	Motion(x,t,dxdt);
310	for (int i=0;i<3;++i)
311	x[i]+=h*dxdt[i];
312
313	/* suggestions about how to write it this way are welcome!
314	void DGLStep(void (f)(const Double_t x,Double_t t,Double_t *dxdt),
315	Double_t *x,Double_t t,Double_t h,Int_t n) const;
316	*/
317
318	}