[u/mrichter/AliRoot.git] / TPC / AliTPCExBExact.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.                  *
 **************************************************************************/

////
// This implementation AliTPCExB is using an "exact" calculation of the ExB
// effect. That is, it uses the drift ODE to calculate the distortion
// without any further assumption.
// Due to the numerical integration of the ODE, there are some numerical
// uncertencies involed.
////

#include "TMath.h"
#include "TTreeStream.h"
#include "AliMagF.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) {
  //
  // Have a look at the common part "TestThisBeautifulObjectGeneric".
  //
  fkMap=bFieldMap;
  fkField=0;
  TestThisBeautifulObjectGeneric(fileName);
}
*/

void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
					     const char* fileName) {
  //
  // Have a look at the common part "TestThisBeautifulObjectGeneric".
  //
  fkField=bField;
  //fkMap=0;
  TestThisBeautifulObjectGeneric(fileName);
}

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