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