Record changes.
[u/mrichter/AliRoot.git] / TPC / AliTPCExBExact.cxx
CommitLineData
faf93237 1#include "TMath.h"
2#include "TTreeStream.h"
3#include "AliTPCExBExact.h"
4
5ClassImp(AliTPCExBExact)
6
7const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
8const Double_t AliTPCExBExact::fgkDriftField=40.e3;
9
481f877b 10AliTPCExBExact::AliTPCExBExact()
11 : fDriftVelocity(0),
12 fkMap(0),fkField(0),fkN(0),
13 fkNX(0),fkNY(0),fkNZ(0),
14 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
15 fkZMin(-250.),fkZMax(250.),
16 fkNLook(0),fkLook(0) {
17 //
18 // purely for I/O
19 //
20}
21
faf93237 22AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
23 Double_t driftVelocity,
24 Int_t nx,Int_t ny,Int_t nz,Int_t n)
481f877b 25 : fDriftVelocity(driftVelocity),
26 fkMap(0),fkField(bField),fkN(n),
faf93237 27 fkNX(nx),fkNY(ny),fkNZ(nz),
28 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b 29 fkZMin(-250.),fkZMax(250.),
30 fkNLook(0),fkLook(0) {
faf93237 31 //
32 // The constructor. One has to supply a magnetic field and an (initial)
33 // drift velocity. Since some kind of lookuptable is created the
34 // number of its meshpoints can be supplied.
35 // n sets the number of integration steps to be used when integrating
36 // over the full drift length.
37 //
faf93237 38 CreateLookupTable();
39}
40
41AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
42 Double_t driftVelocity,Int_t n)
481f877b 43 : fDriftVelocity(driftVelocity),
44 fkMap(bFieldMap),fkField(0),fkN(n),
45 fkNX(0),fkNY(0),fkNZ(0),
46 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
47 fkZMin(-250.),fkZMax(250.),
48 fkNLook(0),fkLook(0) {
faf93237 49 //
50 // The constructor. One has to supply a field map and an (initial)
51 // drift velocity.
52 // n sets the number of integration steps to be used when integrating
53 // over the full drift length.
54 //
faf93237 55
56 fkXMin=bFieldMap->Xmin()
57 -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
58 *bFieldMap->DelX();
59 fkXMax=bFieldMap->Xmax()
60 -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
61 *bFieldMap->DelX();
62 fkYMin=bFieldMap->Ymin()
63 -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
64 *bFieldMap->DelY();
65 fkYMax=bFieldMap->Ymax()
66 -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
67 *bFieldMap->DelY();
68 fkZMax=bFieldMap->Zmax()
69 -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
70 *bFieldMap->DelZ();
71 fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!
72
73 fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
74 fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
75 fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
76
77 CreateLookupTable();
78}
79
80AliTPCExBExact::~AliTPCExBExact() {
81 //
82 // destruct the poor object.
83 //
481f877b 84 delete[] fkLook;
faf93237 85}
86
481f877b 87void AliTPCExBExact::Correct(const Double_t *position, Double_t *corrected) {
88 //
89 // correct for the distortion
90 //
91 Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
92 Int_t xi1=static_cast<Int_t>(x);
93 xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
94 Int_t xi2=xi1+1;
95 Double_t dx=(x-xi1);
96 Double_t dx1=(xi2-x);
97
98 Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
99 Int_t yi1=static_cast<Int_t>(y);
100 yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
101 Int_t yi2=yi1+1;
102 Double_t dy=(y-yi1);
103 Double_t dy1=(yi2-y);
104
105 Double_t z=position[2]/fkZMax*(fkNZ-1);
106 Int_t side;
107 if (z>0) {
108 side=1;
faf93237 109 }
110 else {
481f877b 111 z=-z;
112 side=0;
faf93237 113 }
481f877b 114 Int_t zi1=static_cast<Int_t>(z);
115 zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
116 Int_t zi2=zi1+1;
117 Double_t dz=(z-zi1);
118 Double_t dz1=(zi2-z);
119
120 for (int i=0;i<3;++i)
121 corrected[i]
122 =fkLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
123 +fkLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
124 +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
125 +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
126 +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
127 +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
128 +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
129 +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
130 // corrected[2]=position[2];
131}
132
133void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
134 const char* fileName) {
135 fkMap=bFieldMap;
136 fkField=0;
137 TestThisBeautifulObjectGeneric(fileName);
138}
139
140void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
141 const char* fileName) {
142 fkField=bField;
143 fkMap=0;
144 TestThisBeautifulObjectGeneric(fileName);
faf93237 145}
146
481f877b 147void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) {
faf93237 148 //
149 // well, as the name sais...
150 //
151 TTreeSRedirector ts(fileName);
152 Double_t x[3];
153 for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
154 for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
155 for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
156 Double_t d[3];
157 Double_t dnl[3];
158 Correct(x,d);
159 CalculateDistortion(x,dnl);
160 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
161 Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
162 Double_t dr=r-rd;
163 Double_t phi=TMath::ATan2(x[0],x[1]);
164 Double_t phid=TMath::ATan2(d[0],d[1]);
165 Double_t dphi=phi-phid;
166 if (dphi<0.) dphi+=TMath::TwoPi();
167 if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
168 Double_t drphi=r*dphi;
169 Double_t dx=x[0]-d[0];
170 Double_t dy=x[1]-d[1];
171 Double_t dz=x[2]-d[2];
172 Double_t dnlx=x[0]-dnl[0];
173 Double_t dnly=x[1]-dnl[1];
174 Double_t dnlz=x[2]-dnl[2];
175 ts<<"positions"
176 <<"x0="<<x[0]
177 <<"x1="<<x[1]
178 <<"x2="<<x[2]
179 <<"dx="<<dx
180 <<"dy="<<dy
181 <<"dz="<<dz
182 <<"dnlx="<<dnlx
183 <<"dnly="<<dnly
184 <<"dnlz="<<dnlz
185 <<"r="<<r
186 <<"phi="<<phi
187 <<"dr="<<dr
188 <<"drphi="<<drphi
189 <<"\n";
190 }
191}
192
193void AliTPCExBExact::CreateLookupTable() {
194 //
195 // Helper function to fill the lookup table.
196 //
481f877b 197 fkNLook=fkNX*fkNY*fkNZ*2*3;
198 fkLook=new Double_t[fkNLook];
faf93237 199 Double_t x[3];
200 for (int i=0;i<fkNX;++i) {
201 x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
202 for (int j=0;j<fkNY;++j) {
203 x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
204 for (int k=0;k<fkNZ;++k) {
205 x[2]=1.*fkZMax/(fkNZ-1)*k;
206 x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
481f877b 207 CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
faf93237 208 x[2]=-x[2];
481f877b 209 CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
faf93237 210 }
211 }
212 }
213}
214
215void AliTPCExBExact::GetE(Double_t *E,const Double_t *x) const {
216 //
217 // Helper function returning the E field in SI units (V/m).
218 //
219 E[0]=0.;
220 E[1]=0.;
221 E[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
222}
223
224void AliTPCExBExact::GetB(Double_t *B,const Double_t *x) const {
225 //
226 // Helper function returning the B field in SI units (T).
227 //
228 Float_t xm[3];
229 // the beautiful m to cm (and the ugly "const_cast") and Double_t
230 // to Float_t read the NRs introduction!:
231 for (int i=0;i<3;++i) xm[i]=x[i]*100.;
232 Float_t Bf[3];
233 if (fkMap!=0)
234 fkMap->Field(xm,Bf);
235 else
236 fkField->Field(xm,Bf);
237 for (int i=0;i<3;++i) B[i]=Bf[i]/10.;
238}
239
240void AliTPCExBExact::Motion(const Double_t *x,Double_t,
241 Double_t *dxdt) const {
242 //
243 // The differential equation of motion of the electrons.
244 //
245 const Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
246 const Double_t tau2=tau*tau;
247 Double_t E[3];
248 Double_t B[3];
249 GetE(E,x);
250 GetB(B,x);
251 Double_t wx=fgkEM*B[0];
252 Double_t wy=fgkEM*B[1];
253 Double_t wz=fgkEM*B[2];
254 Double_t ex=fgkEM*E[0];
255 Double_t ey=fgkEM*E[1];
256 Double_t ez=fgkEM*E[2];
257 Double_t w2=(wx*wx+wy*wy+wz*wz);
258 dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
259 dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
260 dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
261 Double_t fac=tau/(1.+w2*tau2);
262 dxdt[0]*=fac;
263 dxdt[1]*=fac;
264 dxdt[2]*=fac;
265}
266
267void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
268 Double_t *dist) const {
269 //
270 // Helper function that calculates one distortion by integration
271 // (only used to fill the lookup table).
272 //
273 const Double_t h=0.01*250./fDriftVelocity/fkN;
274 Double_t t=0.;
275 Double_t xt[3];
276 Double_t xo[3];
277 for (int i=0;i<3;++i)
278 xo[i]=xt[i]=x0[i]*0.01;
279 while (TMath::Abs(xt[2])<250.*0.01) {
280 for (int i=0;i<3;++i)
281 xo[i]=xt[i];
282 DGLStep(xt,t,h);
283 t+=h;
284 }
285 if (t!=0.) {
286 Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
287 dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
288 dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
289 // dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
290 dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
291 dist[2]+=(x0[2]<0.?-1:1.)*250.;
292 }
293 else {
294 dist[0]=x0[0];
295 dist[1]=x0[1];
296 dist[2]=x0[2];
297 }
298 // reverse the distortion, i.e. get the correction
299 dist[0]=x0[0]-(dist[0]-x0[0]);
300 dist[1]=x0[1]-(dist[1]-x0[1]);
301}
302
303void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
304 //
305 // An elementary integration step.
306 // (simple Euler Method)
307 //
308 Double_t dxdt[3];
309 Motion(x,t,dxdt);
310 for (int i=0;i<3;++i)
311 x[i]+=h*dxdt[i];
312
313 /* suggestions about how to write it this way are welcome!
314 void DGLStep(void (*f)(const Double_t *x,Double_t t,Double_t *dxdt),
315 Double_t *x,Double_t t,Double_t h,Int_t n) const;
316 */
317
318}