Adding AliTPCTracklet to the repository (M.Mager)
[u/mrichter/AliRoot.git] / TPC / AliTPCExBFirst.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 aproximate calculation of the ExB
18// effect. Therefore the drift ODE is Taylor expanded and only the first
19// order part is taken.
20////
21
faf93237 22#include "TMath.h"
cf585711 23#include "AliFieldMap.h"
24#include "AliMagF.h"
faf93237 25#include "TTreeStream.h"
26#include "AliTPCExBFirst.h"
27
28ClassImp(AliTPCExBFirst)
29
30const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
31const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
32
481f877b 33AliTPCExBFirst::AliTPCExBFirst()
34 : fDriftVelocity(0),
35 fkNX(0),fkNY(0),fkNZ(0),
36 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
37 fkZMin(-250.),fkZMax(250.),
38 fkNMean(0),
39 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
40 //
41 // purely for I/O
42 //
43}
44
faf93237 45AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
46 Double_t driftVelocity,
47 Int_t nx,Int_t ny,Int_t nz)
481f877b 48 : fDriftVelocity(driftVelocity),
49 fkNX(nx),fkNY(ny),fkNZ(nz),
faf93237 50 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b 51 fkZMin(-250.),fkZMax(250.),
52 fkNMean(0),
53 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237 54 //
55 // The constructor. One has to supply a magnetic field and an (initial)
56 // drift velocity. Since some kind of lookuptable is created the
57 // number of its meshpoints can be supplied.
58 //
faf93237 59 ConstructCommon(0,bField);
60}
61
62AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
481f877b 63 Double_t driftVelocity)
64 : fDriftVelocity(driftVelocity),
65 fkNX(0),fkNY(0),fkNZ(0),
66 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
67 fkZMin(-250.),fkZMax(250.),
68 fkNMean(0),
69 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237 70 //
71 // The constructor. One has to supply a field map and an (initial)
72 // drift velocity.
73 //
faf93237 74
75 fkXMin=bFieldMap->Xmin()
76 -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
77 *bFieldMap->DelX();
78 fkXMax=bFieldMap->Xmax()
79 -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
80 *bFieldMap->DelX();
81 fkYMin=bFieldMap->Ymin()
82 -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
83 *bFieldMap->DelY();
84 fkYMax=bFieldMap->Ymax()
85 -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
86 *bFieldMap->DelY();
87 fkZMin=bFieldMap->Zmin()
88 -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
89 *bFieldMap->DelZ();
90 fkZMax=bFieldMap->Zmax()
91 -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
92 *bFieldMap->DelZ();
93
94 fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
95 fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
96 fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
97
98 ConstructCommon(bFieldMap,0);
99}
100
101AliTPCExBFirst::~AliTPCExBFirst() {
102 //
103 // destruct the poor object.
104 //
481f877b 105 delete[] fkMeanBx;
106 delete[] fkMeanBy;
faf93237 107}
108
109void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
110 //
111 // correct for the distortion
112 //
cf585711 113 Double_t bx,by;
114 GetMeanFields(position[0],position[1],position[2],&bx,&by);
481f877b 115 if (position[2]>0.) {
cf585711 116 Double_t bxe,bye;
117 GetMeanFields(position[0],position[1],250.,&bxe,&bye);
481f877b 118 if (position[2]!=250.) {
cf585711 119 bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]);
120 by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]);
faf93237 121 }
122 else {
cf585711 123 bx=bxe;
124 by=bye;
faf93237 125 }
faf93237 126 }
481f877b 127
128 Double_t mu=fDriftVelocity/fgkDriftField;
129 Double_t wt=mu*fkMeanBz;
130
cf585711 131 corrected[0]=mu*(wt*bx-by)/(1.+wt*wt);
132 corrected[1]=mu*(wt*by+bx)/(1.+wt*wt);
481f877b 133
134 if (position[2]>0.) {
135 corrected[0]*=(250.-position[2]);
136 corrected[1]*=(250.-position[2]);
137 }
138 else {
139 corrected[0]*=(-250.-position[2]);
140 corrected[1]*=(-250.-position[2]);
141 }
142
143 corrected[0]=position[0]-corrected[0];
144 corrected[1]=position[1]-corrected[1];
145 corrected[2]=position[2];
faf93237 146}
147
148void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
149 //
150 // well, as the name sais...
151 //
152 TTreeSRedirector ts(fileName);
153 Double_t x[3];
154 for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
155 for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
156 for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
157 Double_t d[3];
158 Correct(x,d);
159 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
160 Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
161 Double_t dr=r-rd;
162 Double_t phi=TMath::ATan2(x[0],x[1]);
163 Double_t phid=TMath::ATan2(d[0],d[1]);
164 Double_t dphi=phi-phid;
165 if (dphi<0.) dphi+=TMath::TwoPi();
166 if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
167 Double_t drphi=r*dphi;
168 Double_t dx=x[0]-d[0];
169 Double_t dy=x[1]-d[1];
170 Double_t dz=x[2]-d[2];
171 ts<<"positions"
172 <<"x0="<<x[0]
173 <<"x1="<<x[1]
174 <<"x2="<<x[2]
175 <<"dx="<<dx
176 <<"dy="<<dy
177 <<"dz="<<dz
178 <<"r="<<r
179 <<"phi="<<phi
180 <<"dr="<<dr
181 <<"drphi="<<drphi
182 <<"\n";
183 }
184}
185
186void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
187 const AliMagF *bField) {
188 //
189 // THIS IS PRIVATE! (a helper for the constructor)
190 //
481f877b 191 fkNMean=fkNX*fkNY*fkNZ;
cf585711 192 fkMeanBx=new Double_t[fkNMean];
193 fkMeanBy=new Double_t[fkNMean];
faf93237 194
195 Double_t x[3];
196 Double_t nBz=0;
481f877b 197 fkMeanBz=0.;
faf93237 198 for (int i=0;i<fkNX;++i) {
199 x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
200 for (int j=0;j<fkNY;++j) {
201 x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
cf585711 202 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
203 Double_t bx=0.,by=0.;
faf93237 204 for (int k=0;k<fkNZ;++k) {
205 x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
cf585711 206 Float_t b[3];
faf93237 207 // the x is not const in the Field function...
208 Float_t xt[3];
209 for (int l=0;l<3;++l) xt[l]=x[l];
210 // that happens due to the lack of a sophisticated class design:
211 if (bFieldMap!=0)
cf585711 212 bFieldMap->Field(xt,b);
faf93237 213 else
cf585711 214 bField->Field(xt,b);
215 bx+=b[0]/10.;
216 by+=b[1]/10.;
217 fkMeanBx[(k*fkNY+j)*fkNX+i]=bx;
218 fkMeanBy[(k*fkNY+j)*fkNX+i]=by;
219 if (90.<=r&&r<=250.) {
220 fkMeanBz+=b[2]/10.;
faf93237 221 ++nBz;
222 }
223 }
224 }
225 }
481f877b 226 fkMeanBz/=nBz;
faf93237 227}
228
229void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
230 Double_t *Bx,Double_t *By) const {
231 //
232 // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
233 //
234 Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
235 Int_t xi1=static_cast<Int_t>(x);
236 xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
237 Int_t xi2=xi1+1;
238 Double_t dx=(x-xi1);
239 Double_t dx1=(xi2-x);
240
241 Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
242 Int_t yi1=static_cast<Int_t>(y);
243 yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
244 Int_t yi2=yi1+1;
245 Double_t dy=(y-yi1);
246 Double_t dy1=(yi2-y);
247
248 Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
249 Int_t zi1=static_cast<Int_t>(z);
250 zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
251 Int_t zi2=zi1+1;
252 Double_t dz=(z-zi1);
253
481f877b 254 double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
255 +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
256 +fkMeanBx[yi1*fkNX+xi2]*dx *dy
257 +fkMeanBx[yi2*fkNX+xi2]*dx *dy1;
258 double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
259 +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
260 +fkMeanBy[yi1*fkNX+xi2]*dx *dy
261 +fkMeanBy[yi2*fkNX+xi2]*dx *dy1;
faf93237 262 Int_t zi0=zi1-1;
263 double snmx,snmy;
264 if (zi0>=0) {
481f877b 265 snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
266 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
267 +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
268 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
269 snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
270 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
271 +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
272 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
faf93237 273 }
274 else
275 snmx=snmy=0.;
481f877b 276 double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
277 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
278 +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
279 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
280 double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
281 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
282 +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
283 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
284 double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
285 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
286 +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
287 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
288 double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
289 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
290 +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
291 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
faf93237 292 *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
293 *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
294 //TODO: make this nice
295 if (TMath::Abs(z)>0.001) {
296 *Bx/=z;
297 *By/=z;
298 }
299}