Line data Source code
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 : /// \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.
23 :
24 : #include "TMath.h"
25 : #include "TTreeStream.h"
26 : #include "AliMagF.h"
27 : #include "AliTPCExBExact.h"
28 :
29 : /// \cond CLASSIMP
30 24 : ClassImp(AliTPCExBExact)
31 : /// \endcond
32 :
33 : const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
34 : const Double_t AliTPCExBExact::fgkDriftField=-40.e3;
35 :
36 0 : AliTPCExBExact::AliTPCExBExact()
37 0 : : fDriftVelocity(0),
38 : //fkMap(0),
39 0 : fkField(0),fkN(0),
40 0 : fkNX(0),fkNY(0),fkNZ(0),
41 0 : fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
42 0 : fkZMin(-250.),fkZMax(250.),
43 0 : fkNLook(0),fkLook(0) {
44 : /// purely for I/O
45 :
46 0 : }
47 :
48 0 : AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
49 : Double_t driftVelocity,
50 : Int_t nx,Int_t ny,Int_t nz,Int_t n)
51 0 : : fDriftVelocity(driftVelocity),
52 : //fkMap(0),
53 0 : fkField(bField),fkN(n),
54 0 : fkNX(nx),fkNY(ny),fkNZ(nz),
55 0 : fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
56 0 : fkZMin(-250.),fkZMax(250.),
57 0 : fkNLook(0),fkLook(0) {
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 :
64 0 : CreateLookupTable();
65 0 : }
66 :
67 : /*
68 : AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
69 : Double_t driftVelocity,Int_t n)
70 : : fDriftVelocity(driftVelocity),
71 : fkMap(bFieldMap),fkField(0),fkN(n),
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) {
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 : //
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 : }
106 : */
107 :
108 0 : AliTPCExBExact::~AliTPCExBExact() {
109 : /// destruct the poor object.
110 :
111 0 : delete[] fkLook;
112 0 : }
113 :
114 : void AliTPCExBExact::Correct(const Double_t *position, Double_t *corrected) {
115 : /// correct for the distortion
116 :
117 0 : Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
118 0 : Int_t xi1=static_cast<Int_t>(x);
119 0 : xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
120 0 : Int_t xi2=xi1+1;
121 0 : Double_t dx=(x-xi1);
122 0 : Double_t dx1=(xi2-x);
123 :
124 0 : Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
125 0 : Int_t yi1=static_cast<Int_t>(y);
126 0 : yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
127 0 : Int_t yi2=yi1+1;
128 0 : Double_t dy=(y-yi1);
129 0 : Double_t dy1=(yi2-y);
130 :
131 0 : Double_t z=position[2]/fkZMax*(fkNZ-1);
132 : Int_t side;
133 0 : if (z>0) {
134 : side=1;
135 0 : }
136 : else {
137 0 : z=-z;
138 : side=0;
139 : }
140 0 : Int_t zi1=static_cast<Int_t>(z);
141 0 : zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
142 0 : Int_t zi2=zi1+1;
143 0 : Double_t dz=(z-zi1);
144 0 : Double_t dz1=(zi2-z);
145 :
146 0 : for (int i=0;i<3;++i)
147 0 : corrected[i]
148 0 : =fkLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
149 0 : +fkLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
150 0 : +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
151 0 : +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
152 0 : +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
153 0 : +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
154 0 : +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
155 0 : +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
156 : // corrected[2]=position[2];
157 0 : }
158 :
159 : /*
160 : void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
161 : const char* fileName) {
162 : //
163 : // Have a look at the common part "TestThisBeautifulObjectGeneric".
164 : //
165 : fkMap=bFieldMap;
166 : fkField=0;
167 : TestThisBeautifulObjectGeneric(fileName);
168 : }
169 : */
170 :
171 : void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
172 : const char* fileName) {
173 : /// Have a look at the common part "TestThisBeautifulObjectGeneric".
174 :
175 0 : fkField=bField;
176 : //fkMap=0;
177 0 : TestThisBeautifulObjectGeneric(fileName);
178 0 : }
179 :
180 : void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) {
181 : /// Well, as the name sais... it tests the object.
182 :
183 0 : TTreeSRedirector ts(fileName);
184 0 : Double_t x[3];
185 0 : for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
186 0 : for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
187 0 : for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
188 0 : Double_t d[3];
189 0 : Double_t dnl[3];
190 0 : Correct(x,d);
191 0 : CalculateDistortion(x,dnl);
192 0 : Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
193 0 : Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
194 0 : Double_t dr=r-rd;
195 0 : Double_t phi=TMath::ATan2(x[0],x[1]);
196 0 : Double_t phid=TMath::ATan2(d[0],d[1]);
197 0 : Double_t dphi=phi-phid;
198 0 : if (dphi<0.) dphi+=TMath::TwoPi();
199 0 : if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
200 0 : Double_t drphi=r*dphi;
201 0 : Double_t dx=x[0]-d[0];
202 0 : Double_t dy=x[1]-d[1];
203 0 : Double_t dz=x[2]-d[2];
204 0 : Double_t dnlx=x[0]-dnl[0];
205 0 : Double_t dnly=x[1]-dnl[1];
206 0 : Double_t dnlz=x[2]-dnl[2];
207 0 : Double_t b[3];
208 0 : GetB(b,x);
209 0 : ts<<"positions"
210 0 : <<"bx="<<b[0]
211 0 : <<"by="<<b[1]
212 0 : <<"bz="<<b[2]
213 0 : <<"x0="<<x[0]
214 0 : <<"x1="<<x[1]
215 0 : <<"x2="<<x[2]
216 0 : <<"dx="<<dx
217 0 : <<"dy="<<dy
218 0 : <<"dz="<<dz
219 0 : <<"dnlx="<<dnlx
220 0 : <<"dnly="<<dnly
221 0 : <<"dnlz="<<dnlz
222 0 : <<"r="<<r
223 0 : <<"phi="<<phi
224 0 : <<"dr="<<dr
225 0 : <<"drphi="<<drphi
226 0 : <<"\n";
227 0 : }
228 0 : }
229 :
230 : void AliTPCExBExact::CreateLookupTable() {
231 : /// Helper function to fill the lookup table.
232 :
233 0 : fkNLook=fkNX*fkNY*fkNZ*2*3;
234 0 : fkLook=new Double_t[fkNLook];
235 0 : Double_t x[3];
236 0 : for (int i=0;i<fkNX;++i) {
237 0 : x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
238 0 : for (int j=0;j<fkNY;++j) {
239 0 : x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
240 0 : for (int k=0;k<fkNZ;++k) {
241 0 : x[2]=1.*fkZMax/(fkNZ-1)*k;
242 0 : x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
243 0 : CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
244 0 : x[2]=-x[2];
245 0 : CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
246 : }
247 : }
248 : }
249 0 : }
250 :
251 : void AliTPCExBExact::GetE(Double_t *e,const Double_t *x) const {
252 : /// Helper function returning the E field in SI units (V/m).
253 :
254 0 : e[0]=0.;
255 0 : e[1]=0.;
256 0 : e[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
257 0 : }
258 :
259 : void AliTPCExBExact::GetB(Double_t *b,const Double_t *x) const {
260 : /// Helper function returning the B field in SI units (T).
261 :
262 0 : Double_t xm[3];
263 : // the beautiful m to cm (and the ugly "const_cast") and Double_t
264 : // to Float_t read the NRs introduction!:
265 0 : for (int i=0;i<3;++i) xm[i]=x[i]*100.;
266 0 : Double_t bf[3];
267 : //if (fkMap!=0)
268 : // fkMap->Field(xm,bf);
269 : //else
270 0 : ((AliMagF*)fkField)->Field(xm,bf);
271 0 : for (int i=0;i<3;++i) b[i]=bf[i]/10.;
272 0 : }
273 :
274 : void AliTPCExBExact::Motion(const Double_t *x,Double_t,
275 : Double_t *dxdt) const {
276 : /// The differential equation of motion of the electrons.
277 :
278 0 : Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
279 0 : Double_t tau2=tau*tau;
280 0 : Double_t e[3];
281 0 : Double_t b[3];
282 0 : GetE(e,x);
283 0 : GetB(b,x);
284 0 : Double_t wx=fgkEM*b[0];
285 0 : Double_t wy=fgkEM*b[1];
286 0 : Double_t wz=fgkEM*b[2];
287 0 : Double_t ex=fgkEM*e[0];
288 0 : Double_t ey=fgkEM*e[1];
289 0 : Double_t ez=fgkEM*e[2];
290 0 : Double_t w2=(wx*wx+wy*wy+wz*wz);
291 0 : dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
292 0 : dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
293 0 : dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
294 0 : Double_t fac=tau/(1.+w2*tau2);
295 0 : dxdt[0]*=fac;
296 0 : dxdt[1]*=fac;
297 0 : dxdt[2]*=fac;
298 0 : }
299 :
300 : void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
301 : Double_t *dist) const {
302 : /// Helper function that calculates one distortion by integration
303 : /// (only used to fill the lookup table).
304 :
305 0 : Double_t h=0.01*250./fDriftVelocity/fkN;
306 : Double_t t=0.;
307 0 : Double_t xt[3];
308 0 : Double_t xo[3];
309 0 : for (int i=0;i<3;++i)
310 0 : xo[i]=xt[i]=x0[i]*0.01;
311 0 : while (TMath::Abs(xt[2])<250.*0.01) {
312 0 : for (int i=0;i<3;++i)
313 0 : xo[i]=xt[i];
314 0 : DGLStep(xt,t,h);
315 0 : t+=h;
316 : }
317 0 : if (t!=0.) {
318 0 : Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
319 0 : dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
320 0 : dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
321 : // dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
322 0 : dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
323 0 : dist[2]+=(x0[2]<0.?-1:1.)*250.;
324 0 : }
325 : else {
326 0 : dist[0]=x0[0];
327 0 : dist[1]=x0[1];
328 0 : dist[2]=x0[2];
329 : }
330 : // reverse the distortion, i.e. get the correction
331 0 : dist[0]=x0[0]-(dist[0]-x0[0]);
332 0 : dist[1]=x0[1]-(dist[1]-x0[1]);
333 0 : }
334 :
335 : void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
336 : /// An elementary integration step.
337 : /// (simple Euler Method)
338 :
339 0 : Double_t dxdt[3];
340 0 : Motion(x,t,dxdt);
341 0 : for (int i=0;i<3;++i)
342 0 : 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 0 : }
|
