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