coverity fix
[u/mrichter/AliRoot.git] / TPC / AliTPCExBExact.cxx
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
24 #include "TMath.h"
25 #include "TTreeStream.h"
26 #include "AliMagF.h"
27 #include "AliTPCExBExact.h"
28
29 ClassImp(AliTPCExBExact)
30
31 const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
32 const Double_t AliTPCExBExact::fgkDriftField=-40.e3;
33
34 AliTPCExBExact::AliTPCExBExact()
35   : fDriftVelocity(0),
36     //fkMap(0),
37     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
47 AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
48                                Double_t driftVelocity,
49                                Int_t nx,Int_t ny,Int_t nz,Int_t n)
50   : fDriftVelocity(driftVelocity),
51     //fkMap(0),
52     fkField(bField),fkN(n),
53     fkNX(nx),fkNY(ny),fkNZ(nz),
54     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
55     fkZMin(-250.),fkZMax(250.),
56     fkNLook(0),fkLook(0) {
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   //
64   CreateLookupTable();
65 }
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 AliTPCExBExact::~AliTPCExBExact() {
109   //
110   // destruct the poor object.
111   //
112   delete[] fkLook;
113 }
114
115 void 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;
137   }
138   else {
139     z=-z;
140     side=0;
141   }
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
161 /*
162 void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
163                                              const char* fileName) {
164   //
165   // Have a look at the common part "TestThisBeautifulObjectGeneric".
166   //
167   fkMap=bFieldMap;
168   fkField=0;
169   TestThisBeautifulObjectGeneric(fileName);
170 }
171 */
172
173 void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
174                                              const char* fileName) {
175   //
176   // Have a look at the common part "TestThisBeautifulObjectGeneric".
177   //
178   fkField=bField;
179   //fkMap=0;
180   TestThisBeautifulObjectGeneric(fileName);
181 }
182
183 void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) {
184   //
185   // Well, as the name sais... it tests the object.
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];
211         Double_t b[3];
212         GetB(b,x);
213         ts<<"positions"
214           <<"bx="<<b[0]
215           <<"by="<<b[1]
216           <<"bz="<<b[2]
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
234 void AliTPCExBExact::CreateLookupTable() {
235   //
236   // Helper function to fill the lookup table.
237   //
238   fkNLook=fkNX*fkNY*fkNZ*2*3;
239   fkLook=new Double_t[fkNLook];
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
248         CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
249         x[2]=-x[2];
250         CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
251       }
252     }
253   }
254 }
255
256 void AliTPCExBExact::GetE(Double_t *e,const Double_t *x) const {
257   //
258   // Helper function returning the E field in SI units (V/m).
259   //
260   e[0]=0.;
261   e[1]=0.;
262   e[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
263 }
264
265 void AliTPCExBExact::GetB(Double_t *b,const Double_t *x) const {
266   //
267   // Helper function returning the B field in SI units (T).
268   //
269   Double_t xm[3];
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.;
273   Double_t bf[3];
274   //if (fkMap!=0)
275   //  fkMap->Field(xm,bf);
276   //else
277   ((AliMagF*)fkField)->Field(xm,bf);
278   for (int i=0;i<3;++i) b[i]=bf[i]/10.;
279 }
280
281 void AliTPCExBExact::Motion(const Double_t *x,Double_t,
282                             Double_t *dxdt) const {
283   //
284   // The differential equation of motion of the electrons.
285   //
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];
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
308 void 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   //
314   Double_t h=0.01*250./fDriftVelocity/fkN;
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
344 void 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 }