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