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