]>
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" |
f7a1cc68 | 39 | //#include "AliFieldMap.h" |
cf585711 | 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; | |
a187c771 | 47 | const Double_t AliTPCExBFirst::fgkDriftField=-40.e3; |
faf93237 | 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 | // | |
f7a1cc68 | 76 | // ConstructCommon(0,bField); |
77 | ConstructCommon(bField); | |
2abfc1e6 | 78 | SetInstance(this); |
faf93237 | 79 | } |
80 | ||
f7a1cc68 | 81 | /* |
faf93237 | 82 | AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap, |
481f877b | 83 | Double_t driftVelocity) |
84 | : fDriftVelocity(driftVelocity), | |
85 | fkNX(0),fkNY(0),fkNZ(0), | |
86 | fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.), | |
87 | fkZMin(-250.),fkZMax(250.), | |
88 | fkNMean(0), | |
89 | fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) { | |
faf93237 | 90 | // |
91 | // The constructor. One has to supply a field map and an (initial) | |
92 | // drift velocity. | |
93 | // | |
2abfc1e6 | 94 | SetInstance(this); |
faf93237 | 95 | fkXMin=bFieldMap->Xmin() |
96 | -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX()) | |
97 | *bFieldMap->DelX(); | |
98 | fkXMax=bFieldMap->Xmax() | |
99 | -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX()) | |
100 | *bFieldMap->DelX(); | |
101 | fkYMin=bFieldMap->Ymin() | |
102 | -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY()) | |
103 | *bFieldMap->DelY(); | |
104 | fkYMax=bFieldMap->Ymax() | |
105 | -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY()) | |
106 | *bFieldMap->DelY(); | |
107 | fkZMin=bFieldMap->Zmin() | |
108 | -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ()) | |
109 | *bFieldMap->DelZ(); | |
110 | fkZMax=bFieldMap->Zmax() | |
111 | -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ()) | |
112 | *bFieldMap->DelZ(); | |
113 | ||
114 | fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1); | |
115 | fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1); | |
116 | fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1); | |
117 | ||
118 | ConstructCommon(bFieldMap,0); | |
119 | } | |
f7a1cc68 | 120 | */ |
faf93237 | 121 | |
122 | AliTPCExBFirst::~AliTPCExBFirst() { | |
123 | // | |
124 | // destruct the poor object. | |
125 | // | |
481f877b | 126 | delete[] fkMeanBx; |
127 | delete[] fkMeanBy; | |
faf93237 | 128 | } |
129 | ||
130 | void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) { | |
131 | // | |
132 | // correct for the distortion | |
133 | // | |
cf585711 | 134 | Double_t bx,by; |
135 | GetMeanFields(position[0],position[1],position[2],&bx,&by); | |
481f877b | 136 | if (position[2]>0.) { |
cf585711 | 137 | Double_t bxe,bye; |
138 | GetMeanFields(position[0],position[1],250.,&bxe,&bye); | |
481f877b | 139 | if (position[2]!=250.) { |
cf585711 | 140 | bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]); |
141 | by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]); | |
faf93237 | 142 | } |
143 | else { | |
cf585711 | 144 | bx=bxe; |
145 | by=bye; | |
faf93237 | 146 | } |
faf93237 | 147 | } |
481f877b | 148 | |
149 | Double_t mu=fDriftVelocity/fgkDriftField; | |
150 | Double_t wt=mu*fkMeanBz; | |
151 | ||
cf585711 | 152 | corrected[0]=mu*(wt*bx-by)/(1.+wt*wt); |
153 | corrected[1]=mu*(wt*by+bx)/(1.+wt*wt); | |
481f877b | 154 | |
155 | if (position[2]>0.) { | |
156 | corrected[0]*=(250.-position[2]); | |
157 | corrected[1]*=(250.-position[2]); | |
158 | } | |
159 | else { | |
160 | corrected[0]*=(-250.-position[2]); | |
161 | corrected[1]*=(-250.-position[2]); | |
162 | } | |
163 | ||
164 | corrected[0]=position[0]-corrected[0]; | |
165 | corrected[1]=position[1]-corrected[1]; | |
166 | corrected[2]=position[2]; | |
faf93237 | 167 | } |
168 | ||
169 | void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) { | |
170 | // | |
171 | // well, as the name sais... | |
172 | // | |
173 | TTreeSRedirector ts(fileName); | |
174 | Double_t x[3]; | |
175 | for (x[0]=-250.;x[0]<=250.;x[0]+=10.) | |
176 | for (x[1]=-250.;x[1]<=250.;x[1]+=10.) | |
177 | for (x[2]=-250.;x[2]<=250.;x[2]+=10.) { | |
178 | Double_t d[3]; | |
179 | Correct(x,d); | |
180 | Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]); | |
181 | Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]); | |
182 | Double_t dr=r-rd; | |
183 | Double_t phi=TMath::ATan2(x[0],x[1]); | |
184 | Double_t phid=TMath::ATan2(d[0],d[1]); | |
185 | Double_t dphi=phi-phid; | |
186 | if (dphi<0.) dphi+=TMath::TwoPi(); | |
187 | if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi; | |
188 | Double_t drphi=r*dphi; | |
189 | Double_t dx=x[0]-d[0]; | |
190 | Double_t dy=x[1]-d[1]; | |
191 | Double_t dz=x[2]-d[2]; | |
192 | ts<<"positions" | |
193 | <<"x0="<<x[0] | |
194 | <<"x1="<<x[1] | |
195 | <<"x2="<<x[2] | |
196 | <<"dx="<<dx | |
197 | <<"dy="<<dy | |
198 | <<"dz="<<dz | |
199 | <<"r="<<r | |
200 | <<"phi="<<phi | |
201 | <<"dr="<<dr | |
202 | <<"drphi="<<drphi | |
203 | <<"\n"; | |
204 | } | |
205 | } | |
206 | ||
f7a1cc68 | 207 | |
208 | void AliTPCExBFirst::ConstructCommon(//const AliFieldMap *bFieldMap, | |
faf93237 | 209 | const AliMagF *bField) { |
210 | // | |
211 | // THIS IS PRIVATE! (a helper for the constructor) | |
212 | // | |
481f877b | 213 | fkNMean=fkNX*fkNY*fkNZ; |
cf585711 | 214 | fkMeanBx=new Double_t[fkNMean]; |
215 | fkMeanBy=new Double_t[fkNMean]; | |
faf93237 | 216 | |
217 | Double_t x[3]; | |
218 | Double_t nBz=0; | |
481f877b | 219 | fkMeanBz=0.; |
faf93237 | 220 | for (int i=0;i<fkNX;++i) { |
221 | x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1); | |
222 | for (int j=0;j<fkNY;++j) { | |
223 | x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1); | |
cf585711 | 224 | Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]); |
225 | Double_t bx=0.,by=0.; | |
faf93237 | 226 | for (int k=0;k<fkNZ;++k) { |
227 | x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1); | |
f7a1cc68 | 228 | Double_t b[3]; |
faf93237 | 229 | // the x is not const in the Field function... |
f7a1cc68 | 230 | Double_t xt[3]; |
faf93237 | 231 | for (int l=0;l<3;++l) xt[l]=x[l]; |
232 | // that happens due to the lack of a sophisticated class design: | |
f7a1cc68 | 233 | // if (bFieldMap!=0) |
234 | // bFieldMap->Field(xt,b); | |
235 | // else | |
236 | ((AliMagF*)bField)->Field(xt,b); | |
cf585711 | 237 | bx+=b[0]/10.; |
238 | by+=b[1]/10.; | |
239 | fkMeanBx[(k*fkNY+j)*fkNX+i]=bx; | |
240 | fkMeanBy[(k*fkNY+j)*fkNX+i]=by; | |
241 | if (90.<=r&&r<=250.) { | |
242 | fkMeanBz+=b[2]/10.; | |
faf93237 | 243 | ++nBz; |
244 | } | |
245 | } | |
246 | } | |
247 | } | |
481f877b | 248 | fkMeanBz/=nBz; |
faf93237 | 249 | } |
250 | ||
f7a1cc68 | 251 | |
faf93237 | 252 | void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz, |
253 | Double_t *Bx,Double_t *By) const { | |
254 | // | |
255 | // THIS IS PRIVATE! (calculates the mean field utilising a lookup table) | |
256 | // | |
257 | Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin); | |
258 | Int_t xi1=static_cast<Int_t>(x); | |
259 | xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0); | |
260 | Int_t xi2=xi1+1; | |
261 | Double_t dx=(x-xi1); | |
262 | Double_t dx1=(xi2-x); | |
263 | ||
264 | Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin); | |
265 | Int_t yi1=static_cast<Int_t>(y); | |
266 | yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0); | |
267 | Int_t yi2=yi1+1; | |
268 | Double_t dy=(y-yi1); | |
269 | Double_t dy1=(yi2-y); | |
270 | ||
271 | Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin); | |
272 | Int_t zi1=static_cast<Int_t>(z); | |
273 | zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0); | |
274 | Int_t zi2=zi1+1; | |
275 | Double_t dz=(z-zi1); | |
276 | ||
481f877b | 277 | double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1 |
278 | +fkMeanBx[yi2*fkNX+xi1]*dx1*dy | |
82e0477c | 279 | +fkMeanBx[yi1*fkNX+xi2]*dx *dy1 |
280 | +fkMeanBx[yi2*fkNX+xi2]*dx *dy; | |
481f877b | 281 | double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1 |
282 | +fkMeanBy[yi2*fkNX+xi1]*dx1*dy | |
82e0477c | 283 | +fkMeanBy[yi1*fkNX+xi2]*dx *dy1 |
284 | +fkMeanBy[yi2*fkNX+xi2]*dx *dy; | |
faf93237 | 285 | Int_t zi0=zi1-1; |
286 | double snmx,snmy; | |
287 | if (zi0>=0) { | |
481f877b | 288 | snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
289 | +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 290 | +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
291 | +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
481f877b | 292 | snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
293 | +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 294 | +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
295 | +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
faf93237 | 296 | } |
297 | else | |
298 | snmx=snmy=0.; | |
481f877b | 299 | double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
300 | +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 301 | +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
302 | +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
481f877b | 303 | double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
304 | +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 305 | +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
306 | +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
481f877b | 307 | double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
308 | +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 309 | +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
310 | +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
481f877b | 311 | double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1 |
312 | +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy | |
82e0477c | 313 | +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy1 |
314 | +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy; | |
315 | ||
316 | ||
317 | ||
faf93237 | 318 | *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx); |
319 | *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy); | |
320 | //TODO: make this nice | |
321 | if (TMath::Abs(z)>0.001) { | |
322 | *Bx/=z; | |
323 | *By/=z; | |
324 | } | |
325 | } |