2 #include "TTreeStream.h"
3 #include "AliTPCExBExact.h"
5 ClassImp(AliTPCExBExact)
7 const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
8 const Double_t AliTPCExBExact::fgkDriftField=40.e3;
10 AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
11 Double_t driftVelocity,
12 Int_t nx,Int_t ny,Int_t nz,Int_t n)
13 : fkMap(0),fkField(bField),fkN(n),
14 fkNX(nx),fkNY(ny),fkNZ(nz),
15 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
18 // The constructor. One has to supply a magnetic field and an (initial)
19 // drift velocity. Since some kind of lookuptable is created the
20 // number of its meshpoints can be supplied.
21 // n sets the number of integration steps to be used when integrating
22 // over the full drift length.
24 fDriftVelocity=driftVelocity;
28 AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
29 Double_t driftVelocity,Int_t n)
30 : fkMap(bFieldMap),fkField(0),fkN(n) {
32 // The constructor. One has to supply a field map and an (initial)
34 // n sets the number of integration steps to be used when integrating
35 // over the full drift length.
37 fDriftVelocity=driftVelocity;
39 fkXMin=bFieldMap->Xmin()
40 -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
42 fkXMax=bFieldMap->Xmax()
43 -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
45 fkYMin=bFieldMap->Ymin()
46 -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
48 fkYMax=bFieldMap->Ymax()
49 -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
51 fkZMax=bFieldMap->Zmax()
52 -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
54 fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!
56 fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
57 fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
58 fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
63 AliTPCExBExact::~AliTPCExBExact() {
65 // destruct the poor object.
70 void AliTPCExBExact::Correct(const Double_t *position,Double_t *corrected) {
71 Double_t r=TMath::Sqrt(position[0]*position[0]+position[1]*position[1]);
72 if (TMath::Abs(position[2])>250.||r<90.||250.<r) {
73 for (Int_t i=0;i<3;++i) corrected[i]=position[i];
76 Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
77 Int_t xi1=static_cast<Int_t>(x);
78 xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
83 Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
84 Int_t yi1=static_cast<Int_t>(y);
85 yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
90 Double_t z=position[2]/fkZMax*(fkNZ-1);
99 Int_t zi1=static_cast<Int_t>(z);
100 zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
103 Double_t dz1=(zi2-z);
105 for (int i=0;i<3;++i)
107 =fLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
108 +fLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
109 +fLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
110 +fLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
111 +fLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
112 +fLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
113 +fLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
114 +fLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
115 // corrected[2]=position[2];
119 void AliTPCExBExact::TestThisBeautifulObject(const char* fileName) {
121 // well, as the name sais...
123 TTreeSRedirector ts(fileName);
125 for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
126 for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
127 for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
131 CalculateDistortion(x,dnl);
132 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
133 Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
135 Double_t phi=TMath::ATan2(x[0],x[1]);
136 Double_t phid=TMath::ATan2(d[0],d[1]);
137 Double_t dphi=phi-phid;
138 if (dphi<0.) dphi+=TMath::TwoPi();
139 if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
140 Double_t drphi=r*dphi;
141 Double_t dx=x[0]-d[0];
142 Double_t dy=x[1]-d[1];
143 Double_t dz=x[2]-d[2];
144 Double_t dnlx=x[0]-dnl[0];
145 Double_t dnly=x[1]-dnl[1];
146 Double_t dnlz=x[2]-dnl[2];
165 void AliTPCExBExact::CreateLookupTable() {
167 // Helper function to fill the lookup table.
169 fLook=new Double_t[fkNX*fkNY*fkNZ*2*3];
171 for (int i=0;i<fkNX;++i) {
172 x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
173 for (int j=0;j<fkNY;++j) {
174 x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
175 for (int k=0;k<fkNZ;++k) {
176 x[2]=1.*fkZMax/(fkNZ-1)*k;
177 x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
178 CalculateDistortion(x,&fLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
180 CalculateDistortion(x,&fLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
186 void AliTPCExBExact::GetE(Double_t *E,const Double_t *x) const {
188 // Helper function returning the E field in SI units (V/m).
192 E[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
195 void AliTPCExBExact::GetB(Double_t *B,const Double_t *x) const {
197 // Helper function returning the B field in SI units (T).
200 // the beautiful m to cm (and the ugly "const_cast") and Double_t
201 // to Float_t read the NRs introduction!:
202 for (int i=0;i<3;++i) xm[i]=x[i]*100.;
207 fkField->Field(xm,Bf);
208 for (int i=0;i<3;++i) B[i]=Bf[i]/10.;
211 void AliTPCExBExact::Motion(const Double_t *x,Double_t,
212 Double_t *dxdt) const {
214 // The differential equation of motion of the electrons.
216 const Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
217 const Double_t tau2=tau*tau;
222 Double_t wx=fgkEM*B[0];
223 Double_t wy=fgkEM*B[1];
224 Double_t wz=fgkEM*B[2];
225 Double_t ex=fgkEM*E[0];
226 Double_t ey=fgkEM*E[1];
227 Double_t ez=fgkEM*E[2];
228 Double_t w2=(wx*wx+wy*wy+wz*wz);
229 dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
230 dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
231 dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
232 Double_t fac=tau/(1.+w2*tau2);
238 void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
239 Double_t *dist) const {
241 // Helper function that calculates one distortion by integration
242 // (only used to fill the lookup table).
244 const Double_t h=0.01*250./fDriftVelocity/fkN;
248 for (int i=0;i<3;++i)
249 xo[i]=xt[i]=x0[i]*0.01;
250 while (TMath::Abs(xt[2])<250.*0.01) {
251 for (int i=0;i<3;++i)
257 Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
258 dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
259 dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
260 // dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
261 dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
262 dist[2]+=(x0[2]<0.?-1:1.)*250.;
269 // reverse the distortion, i.e. get the correction
270 dist[0]=x0[0]-(dist[0]-x0[0]);
271 dist[1]=x0[1]-(dist[1]-x0[1]);
274 void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
276 // An elementary integration step.
277 // (simple Euler Method)
281 for (int i=0;i<3;++i)
284 /* suggestions about how to write it this way are welcome!
285 void DGLStep(void (*f)(const Double_t *x,Double_t t,Double_t *dxdt),
286 Double_t *x,Double_t t,Double_t h,Int_t n) const;