Classes for space correction
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Jun 2007 14:23:07 +0000 (14:23 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Jun 2007 14:23:07 +0000 (14:23 +0000)
virtual class AliCorrector
+ExB correctors
(Magnus, Marian)

TPC/AliCorrector.cxx [new file with mode: 0644]
TPC/AliCorrector.h [new file with mode: 0644]
TPC/AliTPCExB.cxx [new file with mode: 0644]
TPC/AliTPCExB.h [new file with mode: 0644]
TPC/AliTPCExBExact.cxx [new file with mode: 0644]
TPC/AliTPCExBExact.h [new file with mode: 0644]
TPC/AliTPCExBFirst.cxx [new file with mode: 0644]
TPC/AliTPCExBFirst.h [new file with mode: 0644]
TPC/TPCbaseLinkDef.h
TPC/libTPCbase.pkg

diff --git a/TPC/AliCorrector.cxx b/TPC/AliCorrector.cxx
new file mode 100644 (file)
index 0000000..867511e
--- /dev/null
@@ -0,0 +1,4 @@
+#include "AliCorrector.h"
+
+ClassImp(AliCorrector)
+
diff --git a/TPC/AliCorrector.h b/TPC/AliCorrector.h
new file mode 100644 (file)
index 0000000..c5c42ec
--- /dev/null
@@ -0,0 +1,18 @@
+#ifndef ALI_CORRECTOR_H
+#define ALI_CORRECTOR_H
+
+#include "TObject.h"
+
+class AliCorrector:public TObject {
+public:
+  virtual ~AliCorrector() {};
+  virtual void Correct(const Double_t *position,Double_t *corrected)=0;
+  virtual void CorrectInverse(const Double_t *position,Double_t *corrected) {
+    Correct(position,corrected);
+    for (Int_t i=0;i<3;++i)
+      corrected[i]=position[i]-(corrected[i]-position[i]);
+  }
+  ClassDef(AliCorrector,1)
+};
+
+#endif
diff --git a/TPC/AliTPCExB.cxx b/TPC/AliTPCExB.cxx
new file mode 100644 (file)
index 0000000..cff7236
--- /dev/null
@@ -0,0 +1,4 @@
+#include "AliTPCExB.h"
+
+ClassImp(AliTPCExB)
+
diff --git a/TPC/AliTPCExB.h b/TPC/AliTPCExB.h
new file mode 100644 (file)
index 0000000..e8fd263
--- /dev/null
@@ -0,0 +1,18 @@
+#ifndef ALITPC_EXB
+#define ALITPC_EXB
+
+#include "AliCorrector.h"
+
+class AliTPCExB:public AliCorrector {
+public:
+  virtual ~AliTPCExB() {};
+  void SetDriftVelocity(Double_t driftVelocity) {
+    fDriftVelocity=driftVelocity;
+  };
+protected:
+  Double_t fDriftVelocity; // The electron drift velocity.
+  
+  ClassDef(AliTPCExB,1)
+};
+
+#endif
diff --git a/TPC/AliTPCExBExact.cxx b/TPC/AliTPCExBExact.cxx
new file mode 100644 (file)
index 0000000..053ecb2
--- /dev/null
@@ -0,0 +1,289 @@
+#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(const AliMagF *bField,
+                              Double_t driftVelocity,
+                              Int_t nx,Int_t ny,Int_t nz,Int_t n)
+  : fkMap(0),fkField(bField),fkN(n),
+    fkNX(nx),fkNY(ny),fkNZ(nz),
+    fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(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.
+  // n sets the number of integration steps to be used when integrating
+  // over the full drift length.
+  //
+  fDriftVelocity=driftVelocity;
+  CreateLookupTable();
+}
+
+AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
+                              Double_t driftVelocity,Int_t n) 
+  : fkMap(bFieldMap),fkField(0),fkN(n) {
+  //
+  // 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.
+  //
+  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();
+  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[] fLook;
+}
+
+void AliTPCExBExact::Correct(const Double_t *position,Double_t *corrected) {
+  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 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]
+       =fLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
+       +fLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
+       +fLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
+       +fLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
+       +fLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
+       +fLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
+       +fLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
+       +fLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
+    //    corrected[2]=position[2];
+  }
+}
+
+void AliTPCExBExact::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];
+       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.
+  //
+  fLook=new Double_t[fkNX*fkNY*fkNZ*2*3];
+  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,&fLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
+       x[2]=-x[2];
+       CalculateDistortion(x,&fLook[(((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;
+  */
+
+}
diff --git a/TPC/AliTPCExBExact.h b/TPC/AliTPCExBExact.h
new file mode 100644 (file)
index 0000000..c9a86f3
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef ALITPC_EXB_EXACT
+#define ALITPC_EXB_EXACT
+
+#include "AliTPCExB.h"
+#include "AliFieldMap.h"
+#include "AliMagF.h"
+
+class AliTPCExBExact:public AliTPCExB  {
+public:
+  AliTPCExBExact(const AliFieldMap *bFieldMap,Double_t driftVelocity,
+                Int_t n=100);
+  AliTPCExBExact(const AliMagF *bField,Double_t driftVelocity,Int_t n=100,
+                Int_t nx=100,Int_t ny=100,Int_t nz=100);
+  virtual ~AliTPCExBExact();
+  virtual void Correct(const Double_t *position,Double_t *corrected);
+  void TestThisBeautifulObject(const char* fileName);
+private:
+  AliTPCExBExact& operator=(const AliTPCExBExact&); // don't assign me
+  AliTPCExBExact(const AliTPCExBExact&); // don't copy me
+  void CreateLookupTable();
+  void GetE(Double_t *E,const Double_t *x) const;
+  void GetB(Double_t *B,const Double_t *x) const;
+  void Motion(const Double_t *x,Double_t t,Double_t *dxdt) const;
+  void CalculateDistortion(const Double_t *x,Double_t *dist) const;
+  void DGLStep(Double_t *x,Double_t t,Double_t h) const;
+  const AliFieldMap *fkMap; // the magnetic field map as supplied by the user
+  const AliMagF *fkField;   // the magnetic field as supplied by the user
+  Int_t fkN;       // max number of integration steps
+  Int_t fkNX;      // field mesh points in x direction
+  Int_t fkNY;      // field mesh points in y direction
+  Int_t fkNZ;      // field mesh points in z direction
+  Double_t fkXMin; // the first grid point in x direction
+  Double_t fkXMax; // the last grid point in x direction
+  Double_t fkYMin; // the first grid point in y direction
+  Double_t fkYMax; // the last grid point in y direction
+  Double_t fkZMin; // the first grid point in z direction
+  Double_t fkZMax; // the last grid point in z direction
+  Double_t *fLook; // the great lookup table
+  static const Double_t fgkEM; // elementary charge over electron mass (C/kg)
+  static const Double_t fgkDriftField; // the TPC drift field (V/m) (modulus)
+
+  ClassDef(AliTPCExBExact,1)
+};
+
+#endif
diff --git a/TPC/AliTPCExBFirst.cxx b/TPC/AliTPCExBFirst.cxx
new file mode 100644 (file)
index 0000000..e793246
--- /dev/null
@@ -0,0 +1,266 @@
+#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;
+  }
+}
diff --git a/TPC/AliTPCExBFirst.h b/TPC/AliTPCExBFirst.h
new file mode 100644 (file)
index 0000000..c003ea0
--- /dev/null
@@ -0,0 +1,40 @@
+#ifndef ALITPC_EXB_FIRST_H
+#define ALITPC_EXB_FIRST_H
+
+#include "AliTPCExB.h"
+#include "AliFieldMap.h"
+#include "AliMagF.h"
+
+class AliTPCExBFirst:public AliTPCExB {
+public:
+  AliTPCExBFirst(const AliFieldMap *bFieldMap,Double_t driftVelocity);
+  AliTPCExBFirst(const AliMagF *bField,Double_t driftVelocity,
+                Int_t nx=100,Int_t ny=100,Int_t nz=100);
+  virtual ~AliTPCExBFirst();
+  virtual void Correct(const Double_t *position,Double_t *corrected);
+  void TestThisBeautifulObject(const char* fileName);
+private:
+  AliTPCExBFirst& operator=(const AliTPCExBFirst&); // don't assign me
+  AliTPCExBFirst(const AliTPCExBFirst&); // don't copy me
+  void ConstructCommon(const AliFieldMap *bFieldMap,const AliMagF *bField);
+  void GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
+                    Double_t *Bx,Double_t *By) const;
+  Int_t fkNX;        // field mesh points in x direction
+  Int_t fkNY;        // field mesh points in y direction
+  Int_t fkNZ;        // field mesh points in z direction
+  Double_t fkXMin;   // the first grid point in x direction
+  Double_t fkXMax;   // the last grid point in x direction
+  Double_t fkYMin;   // the first grid point in y direction
+  Double_t fkYMax;   // the last grid point in y direction
+  Double_t fkZMin;   // the first grid point in z direction
+  Double_t fkZMax;   // the last grid point in z direction
+  Double_t *fMeanBx; // the mean field in x direction upto a certain z value
+  Double_t *fMeanBy; // the mean field in y direction upto a certain z value
+  Double_t fMeanBz;  // the mean field in z direction inside the TPC volume
+  static const Double_t fgkEM; // elementary charge over electron mass (C/kg)
+  static const Double_t fgkDriftField; // the TPC drift field (V/m) (modulus)
+
+  ClassDef(AliTPCExBFirst,1)
+};
+
+#endif
index 56abca3..35fffa2 100644 (file)
 #pragma link C++ class AliTPCCalibViewerGUI+;
 
 
+#pragma link C++ class AliCorrector+;
+#pragma link C++ class AliTPCExBExact+;
+#pragma link C++ class AliTPCExBFirst+ ;
+#pragma link C++ class AliTPCExB+;
+
+
+
 #endif
 
index c0fd513..6e6d8da 100644 (file)
@@ -11,7 +11,8 @@ SRCS:=  AliSegmentID.cxx  AliSegmentArray.cxx AliDigits.cxx AliH2F.cxx \
        AliTPCLaserTracks.cxx AliTPCSensorTemp.cxx  AliTPCSensorTempArray.cxx \
        AliTPCCalibPedestal.cxx AliTPCCalibSignal.cxx   AliTPCCalibCE.cxx \
         AliTPCPreprocessor.cxx  \
-        AliTPCCalibViewer.cxx  AliTPCCalibViewerGUI.cxx
+        AliTPCCalibViewer.cxx  AliTPCCalibViewerGUI.cxx \
+        AliCorrector.cxx  AliTPCExB.cxx  AliTPCExBExact.cxx AliTPCExBFirst.cxx