]>
Commit | Line | Data |
---|---|---|
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 | ||
22 | #include "TMath.h" | |
23 | #include "AliFieldMap.h" | |
24 | #include "AliMagF.h" | |
25 | #include "TTreeStream.h" | |
26 | #include "AliTPCExBFirst.h" | |
27 | ||
28 | ClassImp(AliTPCExBFirst) | |
29 | ||
30 | const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31; | |
31 | const Double_t AliTPCExBFirst::fgkDriftField=40.e3; | |
32 | ||
33 | AliTPCExBFirst::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 | ||
45 | AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField, | |
46 | Double_t driftVelocity, | |
47 | Int_t nx,Int_t ny,Int_t nz) | |
48 | : fDriftVelocity(driftVelocity), | |
49 | fkNX(nx),fkNY(ny),fkNZ(nz), | |
50 | fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.), | |
51 | fkZMin(-250.),fkZMax(250.), | |
52 | fkNMean(0), | |
53 | fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) { | |
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 | // | |
59 | ConstructCommon(0,bField); | |
60 | } | |
61 | ||
62 | AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap, | |
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.) { | |
70 | // | |
71 | // The constructor. One has to supply a field map and an (initial) | |
72 | // drift velocity. | |
73 | // | |
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 | ||
101 | AliTPCExBFirst::~AliTPCExBFirst() { | |
102 | // | |
103 | // destruct the poor object. | |
104 | // | |
105 | delete[] fkMeanBx; | |
106 | delete[] fkMeanBy; | |
107 | } | |
108 | ||
109 | void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) { | |
110 | // | |
111 | // correct for the distortion | |
112 | // | |
113 | Double_t bx,by; | |
114 | GetMeanFields(position[0],position[1],position[2],&bx,&by); | |
115 | if (position[2]>0.) { | |
116 | Double_t bxe,bye; | |
117 | GetMeanFields(position[0],position[1],250.,&bxe,&bye); | |
118 | if (position[2]!=250.) { | |
119 | bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]); | |
120 | by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]); | |
121 | } | |
122 | else { | |
123 | bx=bxe; | |
124 | by=bye; | |
125 | } | |
126 | } | |
127 | ||
128 | Double_t mu=fDriftVelocity/fgkDriftField; | |
129 | Double_t wt=mu*fkMeanBz; | |
130 | ||
131 | corrected[0]=mu*(wt*bx-by)/(1.+wt*wt); | |
132 | corrected[1]=mu*(wt*by+bx)/(1.+wt*wt); | |
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]; | |
146 | } | |
147 | ||
148 | void 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 | ||
186 | void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap, | |
187 | const AliMagF *bField) { | |
188 | // | |
189 | // THIS IS PRIVATE! (a helper for the constructor) | |
190 | // | |
191 | fkNMean=fkNX*fkNY*fkNZ; | |
192 | fkMeanBx=new Double_t[fkNMean]; | |
193 | fkMeanBy=new Double_t[fkNMean]; | |
194 | ||
195 | Double_t x[3]; | |
196 | Double_t nBz=0; | |
197 | fkMeanBz=0.; | |
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); | |
202 | Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]); | |
203 | Double_t bx=0.,by=0.; | |
204 | for (int k=0;k<fkNZ;++k) { | |
205 | x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1); | |
206 | Float_t b[3]; | |
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) | |
212 | bFieldMap->Field(xt,b); | |
213 | else | |
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.; | |
221 | ++nBz; | |
222 | } | |
223 | } | |
224 | } | |
225 | } | |
226 | fkMeanBz/=nBz; | |
227 | } | |
228 | ||
229 | void 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 | ||
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; | |
262 | Int_t zi0=zi1-1; | |
263 | double snmx,snmy; | |
264 | if (zi0>=0) { | |
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; | |
273 | } | |
274 | else | |
275 | snmx=snmy=0.; | |
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; | |
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 | } |