X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCExBExact.cxx;h=ae9336e6490f3054379255259c205e4294a0903e;hb=e81dc11276498f9f0adc95ec4f6a0943a552eae5;hp=053ecb234ff38ab23c6f1e5a7e3e3e48b62472ba;hpb=faf9323738d83434592a89e8b4439f5b9130db9c;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCExBExact.cxx b/TPC/AliTPCExBExact.cxx index 053ecb234ff..ae9336e6490 100644 --- a/TPC/AliTPCExBExact.cxx +++ b/TPC/AliTPCExBExact.cxx @@ -1,5 +1,30 @@ +/************************************************************************** + * 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 "AliFieldMap.h" +#include "AliMagF.h" #include "AliTPCExBExact.h" ClassImp(AliTPCExBExact) @@ -7,13 +32,27 @@ 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) - : fkMap(0),fkField(bField),fkN(n), + : fDriftVelocity(driftVelocity), + fkMap(0),fkField(bField),fkN(n), fkNX(nx),fkNY(ny),fkNZ(nz), fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.), - fkZMax(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 @@ -21,20 +60,23 @@ AliTPCExBExact::AliTPCExBExact(const AliMagF *bField, // 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) { + : 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. // - fDriftVelocity=driftVelocity; fkXMin=bFieldMap->Xmin() -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX()) @@ -64,61 +106,78 @@ AliTPCExBExact::~AliTPCExBExact() { // // destruct the poor object. // - delete[] fLook; + delete[] fkLook; } -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 { - Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1); - 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=(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); + 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); - 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] + =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]; +} - 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 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::TestThisBeautifulObject(const char* fileName) { +void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) { // - // well, as the name sais... + // Well, as the name sais... it tests the object. // TTreeSRedirector ts(fileName); Double_t x[3]; @@ -144,7 +203,12 @@ void AliTPCExBExact::TestThisBeautifulObject(const char* fileName) { 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="<Field(xm,Bf); + fkMap->Field(xm,bf); else - fkField->Field(xm,Bf); - for (int i=0;i<3;++i) B[i]=Bf[i]/10.; + 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, @@ -213,18 +278,18 @@ void AliTPCExBExact::Motion(const Double_t *x,Double_t, // // 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 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; @@ -241,7 +306,7 @@ void AliTPCExBExact::CalculateDistortion(const Double_t *x0, // 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 h=0.01*250./fDriftVelocity/fkN; Double_t t=0.; Double_t xt[3]; Double_t xo[3];