[u/mrichter/AliRoot.git] / TPC / AliTPCExBFirst.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 aproximate calculation of the ExB
// effect. Therefore the drift ODE is Taylor expanded and only the first
// order part is taken.
//
// 
// The ExB correction map is stored in the calib DB
// To test current version:
/*
  char *storage = "local://OCDBres"
  Int_t RunNumber=0;
  AliCDBManager::Instance()->SetDefaultStorage(storage);
  AliCDBManager::Instance()->SetRun(RunNumber) 
  AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();

  //  exb->TestExB("exb.root");
  // TFile f("exb.root");
  //positions->Draw("drphi");
*/


#include "TMath.h"
//#include "AliFieldMap.h"
#include "AliMagF.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
  //
  SetInstance(this);
}

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);
  ConstructCommon(bField);
  SetInstance(this);
}

/*
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.
  //
  SetInstance(this);
  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 Double_t[fkNMean];
  fkMeanBy=new Double_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);
	Double_t b[3];
	// the x is not const in the Field function...
	Double_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 
	((AliMagF*)bField)->Field(xt,b);
	bx+=b[0]/10.;
	by+=b[1]/10.;
	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 *dy1
            +fkMeanBx[yi2*fkNX+xi2]*dx *dy;
  double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
	    +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
            +fkMeanBy[yi1*fkNX+xi2]*dx *dy1
            +fkMeanBy[yi2*fkNX+xi2]*dx *dy;
  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 *dy1
        +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy;
    snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
        +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
        +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy1
        +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy;
  }
  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 *dy1
            +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy;
  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 *dy1
            +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy;
  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 *dy1
             +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy;
  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 *dy1
             +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy;


  *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
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
2abfc1e6	16	//
cf585711	17	// This implementation AliTPCExB is using an aproximate calculation of the ExB
	18	// effect. Therefore the drift ODE is Taylor expanded and only the first
	19	// order part is taken.
2abfc1e6	20	//
	21	//
	22	// The ExB correction map is stored in the calib DB
	23	// To test current version:
	24	/*
	25	char *storage = "local://OCDBres"
	26	Int_t RunNumber=0;
	27	AliCDBManager::Instance()->SetDefaultStorage(storage);
	28	AliCDBManager::Instance()->SetRun(RunNumber)
	29	AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();
	30
	31	// exb->TestExB("exb.root");
	32	// TFile f("exb.root");
	33	//positions->Draw("drphi");
	34	*/
	35
	36
cf585711	37
faf93237	38	#include "TMath.h"
f7a1cc68	39	//#include "AliFieldMap.h"
cf585711	40	#include "AliMagF.h"
faf93237	41	#include "TTreeStream.h"
	42	#include "AliTPCExBFirst.h"
	43
	44	ClassImp(AliTPCExBFirst)
	45
	46	const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
a187c771	47	const Double_t AliTPCExBFirst::fgkDriftField=-40.e3;
faf93237	48
481f877b	49	AliTPCExBFirst::AliTPCExBFirst()
	50	: fDriftVelocity(0),
	51	fkNX(0),fkNY(0),fkNZ(0),
	52	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	53	fkZMin(-250.),fkZMax(250.),
	54	fkNMean(0),
	55	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
	56	//
	57	// purely for I/O
	58	//
2abfc1e6	59	SetInstance(this);
481f877b	60	}
481f877b	61
faf93237	62	AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
	63	Double_t driftVelocity,
	64	Int_t nx,Int_t ny,Int_t nz)
481f877b	65	: fDriftVelocity(driftVelocity),
481f877b	66	fkNX(nx),fkNY(ny),fkNZ(nz),
faf93237	67	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b	68	fkZMin(-250.),fkZMax(250.),
	69	fkNMean(0),
	70	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237	71	//
	72	// The constructor. One has to supply a magnetic field and an (initial)
	73	// drift velocity. Since some kind of lookuptable is created the
	74	// number of its meshpoints can be supplied.
	75	//
f7a1cc68	76	// ConstructCommon(0,bField);
f7a1cc68	77	ConstructCommon(bField);
2abfc1e6	78	SetInstance(this);
faf93237	79	}
faf93237	80
f7a1cc68	81	/*
faf93237	82	AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
481f877b	83	Double_t driftVelocity)
	84	: fDriftVelocity(driftVelocity),
	85	fkNX(0),fkNY(0),fkNZ(0),
	86	fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
	87	fkZMin(-250.),fkZMax(250.),
	88	fkNMean(0),
	89	fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237	90	//
	91	// The constructor. One has to supply a field map and an (initial)
	92	// drift velocity.
	93	//
2abfc1e6	94	SetInstance(this);
faf93237	95	fkXMin=bFieldMap->Xmin()
	96	-TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
	97	*bFieldMap->DelX();
	98	fkXMax=bFieldMap->Xmax()
	99	-TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
	100	*bFieldMap->DelX();
	101	fkYMin=bFieldMap->Ymin()
	102	-TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
	103	*bFieldMap->DelY();
	104	fkYMax=bFieldMap->Ymax()
	105	-TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
	106	*bFieldMap->DelY();
	107	fkZMin=bFieldMap->Zmin()
	108	-TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
	109	*bFieldMap->DelZ();
	110	fkZMax=bFieldMap->Zmax()
	111	-TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
	112	*bFieldMap->DelZ();
	113
	114	fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
	115	fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
	116	fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
	117
	118	ConstructCommon(bFieldMap,0);
	119	}
f7a1cc68	120	*/
faf93237	121
	122	AliTPCExBFirst::~AliTPCExBFirst() {
	123	//
	124	// destruct the poor object.
	125	//
481f877b	126	delete[] fkMeanBx;
481f877b	127	delete[] fkMeanBy;
faf93237	128	}
	129
	130	void AliTPCExBFirst::Correct(const Double_t position,Double_t corrected) {
	131	//
	132	// correct for the distortion
	133	//
cf585711	134	Double_t bx,by;
cf585711	135	GetMeanFields(position[0],position[1],position[2],&bx,&by);
481f877b	136	if (position[2]>0.) {
cf585711	137	Double_t bxe,bye;
cf585711	138	GetMeanFields(position[0],position[1],250.,&bxe,&bye);
481f877b	139	if (position[2]!=250.) {
cf585711	140	bx=(500.bxe-(position[2]+250.)bx)/(250.-position[2]);
cf585711	141	by=(500.bye-(position[2]+250.)by)/(250.-position[2]);
faf93237	142	}
faf93237	143	else {
cf585711	144	bx=bxe;
cf585711	145	by=bye;
faf93237	146	}
faf93237	147	}
481f877b	148
	149	Double_t mu=fDriftVelocity/fgkDriftField;
	150	Double_t wt=mu*fkMeanBz;
	151
cf585711	152	corrected[0]=mu(wtbx-by)/(1.+wt*wt);
cf585711	153	corrected[1]=mu(wtby+bx)/(1.+wt*wt);
481f877b	154
	155	if (position[2]>0.) {
	156	corrected[0]*=(250.-position[2]);
	157	corrected[1]*=(250.-position[2]);
	158	}
	159	else {
	160	corrected[0]*=(-250.-position[2]);
	161	corrected[1]*=(-250.-position[2]);
	162	}
	163
	164	corrected[0]=position[0]-corrected[0];
	165	corrected[1]=position[1]-corrected[1];
	166	corrected[2]=position[2];
faf93237	167	}
	168
	169	void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
	170	//
	171	// well, as the name sais...
	172	//
	173	TTreeSRedirector ts(fileName);
	174	Double_t x[3];
	175	for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
	176	for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
	177	for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
	178	Double_t d[3];
	179	Correct(x,d);
	180	Double_t r=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
	181	Double_t rd=TMath::Sqrt(d[0]d[0]+d[1]d[1]);
	182	Double_t dr=r-rd;
	183	Double_t phi=TMath::ATan2(x[0],x[1]);
	184	Double_t phid=TMath::ATan2(d[0],d[1]);
	185	Double_t dphi=phi-phid;
	186	if (dphi<0.) dphi+=TMath::TwoPi();
	187	if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
	188	Double_t drphi=r*dphi;
	189	Double_t dx=x[0]-d[0];
	190	Double_t dy=x[1]-d[1];
	191	Double_t dz=x[2]-d[2];
	192	ts<<"positions"
	193	<<"x0="<<x[0]
	194	<<"x1="<<x[1]
	195	<<"x2="<<x[2]
	196	<<"dx="<<dx
	197	<<"dy="<<dy
	198	<<"dz="<<dz
	199	<<"r="<<r
	200	<<"phi="<<phi
	201	<<"dr="<<dr
	202	<<"drphi="<<drphi
	203	<<"\n";
	204	}
	205	}
	206
f7a1cc68	207
f7a1cc68	208	void AliTPCExBFirst::ConstructCommon(//const AliFieldMap *bFieldMap,
faf93237	209	const AliMagF *bField) {
	210	//
	211	// THIS IS PRIVATE! (a helper for the constructor)
	212	//
481f877b	213	fkNMean=fkNXfkNYfkNZ;
cf585711	214	fkMeanBx=new Double_t[fkNMean];
cf585711	215	fkMeanBy=new Double_t[fkNMean];
faf93237	216
	217	Double_t x[3];
	218	Double_t nBz=0;
481f877b	219	fkMeanBz=0.;
faf93237	220	for (int i=0;i<fkNX;++i) {
	221	x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
	222	for (int j=0;j<fkNY;++j) {
	223	x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
cf585711	224	Double_t r=TMath::Sqrt(x[0]x[0]+x[1]x[1]);
cf585711	225	Double_t bx=0.,by=0.;
faf93237	226	for (int k=0;k<fkNZ;++k) {
faf93237	227	x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
f7a1cc68	228	Double_t b[3];
faf93237	229	// the x is not const in the Field function...
f7a1cc68	230	Double_t xt[3];
faf93237	231	for (int l=0;l<3;++l) xt[l]=x[l];
faf93237	232	// that happens due to the lack of a sophisticated class design:
f7a1cc68	233	// if (bFieldMap!=0)
	234	// bFieldMap->Field(xt,b);
	235	// else
	236	((AliMagF*)bField)->Field(xt,b);
cf585711	237	bx+=b[0]/10.;
	238	by+=b[1]/10.;
	239	fkMeanBx[(kfkNY+j)fkNX+i]=bx;
	240	fkMeanBy[(kfkNY+j)fkNX+i]=by;
	241	if (90.<=r&&r<=250.) {
	242	fkMeanBz+=b[2]/10.;
faf93237	243	++nBz;
	244	}
	245	}
	246	}
	247	}
481f877b	248	fkMeanBz/=nBz;
faf93237	249	}
faf93237	250
f7a1cc68	251
faf93237	252	void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
	253	Double_t Bx,Double_t By) const {
	254	//
	255	// THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
	256	//
	257	Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
	258	Int_t xi1=static_cast<Int_t>(x);
	259	xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
	260	Int_t xi2=xi1+1;
	261	Double_t dx=(x-xi1);
	262	Double_t dx1=(xi2-x);
	263
	264	Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
	265	Int_t yi1=static_cast<Int_t>(y);
	266	yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
	267	Int_t yi2=yi1+1;
	268	Double_t dy=(y-yi1);
	269	Double_t dy1=(yi2-y);
	270
	271	Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
	272	Int_t zi1=static_cast<Int_t>(z);
	273	zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
	274	Int_t zi2=zi1+1;
	275	Double_t dz=(z-zi1);
	276
481f877b	277	double s0x=fkMeanBx[yi1fkNX+xi1]dx1*dy1
481f877b	278	+fkMeanBx[yi2fkNX+xi1]dx1*dy
82e0477c	279	+fkMeanBx[yi1fkNX+xi2]dx *dy1
82e0477c	280	+fkMeanBx[yi2fkNX+xi2]dx *dy;
481f877b	281	double s0y=fkMeanBy[yi1fkNX+xi1]dx1*dy1
481f877b	282	+fkMeanBy[yi2fkNX+xi1]dx1*dy
82e0477c	283	+fkMeanBy[yi1fkNX+xi2]dx *dy1
82e0477c	284	+fkMeanBy[yi2fkNX+xi2]dx *dy;
faf93237	285	Int_t zi0=zi1-1;
	286	double snmx,snmy;
	287	if (zi0>=0) {
481f877b	288	snmx=fkMeanBx[(zi0fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	289	+fkMeanBx[(zi0fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	290	+fkMeanBx[(zi0fkNY+yi1)fkNX+xi2]dx dy1
82e0477c	291	+fkMeanBx[(zi0fkNY+yi2)fkNX+xi2]dx dy;
481f877b	292	snmy=fkMeanBy[(zi0fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	293	+fkMeanBy[(zi0fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	294	+fkMeanBy[(zi0fkNY+yi1)fkNX+xi2]dx dy1
82e0477c	295	+fkMeanBy[(zi0fkNY+yi2)fkNX+xi2]dx dy;
faf93237	296	}
	297	else
	298	snmx=snmy=0.;
481f877b	299	double snx=fkMeanBx[(zi1fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	300	+fkMeanBx[(zi1fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	301	+fkMeanBx[(zi1fkNY+yi1)fkNX+xi2]dx dy1
82e0477c	302	+fkMeanBx[(zi1fkNY+yi2)fkNX+xi2]dx dy;
481f877b	303	double sny=fkMeanBy[(zi1fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	304	+fkMeanBy[(zi1fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	305	+fkMeanBy[(zi1fkNY+yi1)fkNX+xi2]dx dy1
82e0477c	306	+fkMeanBy[(zi1fkNY+yi2)fkNX+xi2]dx dy;
481f877b	307	double snpx=fkMeanBx[(zi2fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	308	+fkMeanBx[(zi2fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	309	+fkMeanBx[(zi2fkNY+yi1)fkNX+xi2]dx dy1
82e0477c	310	+fkMeanBx[(zi2fkNY+yi2)fkNX+xi2]dx dy;
481f877b	311	double snpy=fkMeanBy[(zi2fkNY+yi1)fkNX+xi1]dx1dy1
481f877b	312	+fkMeanBy[(zi2fkNY+yi2)fkNX+xi1]dx1dy
82e0477c	313	+fkMeanBy[(zi2fkNY+yi1)fkNX+xi2]dx dy1
	314	+fkMeanBy[(zi2fkNY+yi2)fkNX+xi2]dx dy;
	315
	316
	317
faf93237	318	Bx=0.5(((snpx-2.snx+snmx)dz+2.(snx-snmx))dz+snx-s0x+snmx);
	319	By=0.5(((snpy-2.sny+snmy)dz+2.(sny-snmy))dz+sny-s0y+snmy);
	320	//TODO: make this nice
	321	if (TMath::Abs(z)>0.001) {
	322	*Bx/=z;
	323	*By/=z;
	324	}
	325	}