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