From faf9323738d83434592a89e8b4439f5b9130db9c Mon Sep 17 00:00:00 2001 From: marian Date: Wed, 6 Jun 2007 14:23:07 +0000 Subject: [PATCH] Classes for space correction virtual class AliCorrector +ExB correctors (Magnus, Marian) --- TPC/AliCorrector.cxx | 4 + TPC/AliCorrector.h | 18 +++ TPC/AliTPCExB.cxx | 4 + TPC/AliTPCExB.h | 18 +++ TPC/AliTPCExBExact.cxx | 289 +++++++++++++++++++++++++++++++++++++++++ TPC/AliTPCExBExact.h | 45 +++++++ TPC/AliTPCExBFirst.cxx | 266 +++++++++++++++++++++++++++++++++++++ TPC/AliTPCExBFirst.h | 40 ++++++ TPC/TPCbaseLinkDef.h | 7 + TPC/libTPCbase.pkg | 3 +- 10 files changed, 693 insertions(+), 1 deletion(-) create mode 100644 TPC/AliCorrector.cxx create mode 100644 TPC/AliCorrector.h create mode 100644 TPC/AliTPCExB.cxx create mode 100644 TPC/AliTPCExB.h create mode 100644 TPC/AliTPCExBExact.cxx create mode 100644 TPC/AliTPCExBExact.h create mode 100644 TPC/AliTPCExBFirst.cxx create mode 100644 TPC/AliTPCExBFirst.h diff --git a/TPC/AliCorrector.cxx b/TPC/AliCorrector.cxx new file mode 100644 index 00000000000..867511e51cc --- /dev/null +++ b/TPC/AliCorrector.cxx @@ -0,0 +1,4 @@ +#include "AliCorrector.h" + +ClassImp(AliCorrector) + diff --git a/TPC/AliCorrector.h b/TPC/AliCorrector.h new file mode 100644 index 00000000000..c5c42ec18b7 --- /dev/null +++ b/TPC/AliCorrector.h @@ -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 index 00000000000..cff72361408 --- /dev/null +++ b/TPC/AliTPCExB.cxx @@ -0,0 +1,4 @@ +#include "AliTPCExB.h" + +ClassImp(AliTPCExB) + diff --git a/TPC/AliTPCExB.h b/TPC/AliTPCExB.h new file mode 100644 index 00000000000..e8fd263046e --- /dev/null +++ b/TPC/AliTPCExB.h @@ -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 index 00000000000..053ecb234ff --- /dev/null +++ b/TPC/AliTPCExBExact.cxx @@ -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((fkXMax-fkXMin)/bFieldMap->DelX()+1.1); + fkNY=static_cast((fkYMax-fkYMin)/bFieldMap->DelY()+1.1); + fkNZ=static_cast((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.(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(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(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="<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 index 00000000000..c9a86f3a76e --- /dev/null +++ b/TPC/AliTPCExBExact.h @@ -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 index 00000000000..e7932464d8a --- /dev/null +++ b/TPC/AliTPCExBFirst.cxx @@ -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((fkXMax-fkXMin)/bFieldMap->DelX()+1.1); + fkNY=static_cast((fkYMax-fkYMin)/bFieldMap->DelY()+1.1); + fkNZ=static_cast((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.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="<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(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(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(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 index 00000000000..c003ea0f675 --- /dev/null +++ b/TPC/AliTPCExBFirst.h @@ -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 diff --git a/TPC/TPCbaseLinkDef.h b/TPC/TPCbaseLinkDef.h index 56abca3f8b7..35fffa20c11 100644 --- a/TPC/TPCbaseLinkDef.h +++ b/TPC/TPCbaseLinkDef.h @@ -51,5 +51,12 @@ #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 diff --git a/TPC/libTPCbase.pkg b/TPC/libTPCbase.pkg index c0fd5136ebc..6e6d8da2033 100644 --- a/TPC/libTPCbase.pkg +++ b/TPC/libTPCbase.pkg @@ -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 -- 2.43.0