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