]>
Commit | Line | Data |
---|---|---|
cabb917f | 1 | #include "TTree.h" |
2 | #include "TMatrixDSym.h" | |
3 | #include "TMatrixD.h" | |
4 | #include "AliRieman.h" | |
5 | #include "TRandom.h" | |
6 | ||
7 | ClassImp(AliRieman) | |
8 | ||
9 | ||
10 | ||
11 | AliRieman::AliRieman(){ | |
12 | // | |
13 | // default constructor | |
14 | // | |
15 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
16 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
17 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
18 | ||
19 | fCapacity = 0; | |
20 | fN =0; | |
21 | fX = 0; | |
22 | fY = 0; | |
23 | fZ = 0; | |
24 | fSy = 0; | |
25 | fSz = 0; | |
26 | fCovar = 0; | |
27 | fConv = kFALSE; | |
28 | } | |
29 | ||
30 | ||
31 | AliRieman::AliRieman(Int_t capacity){ | |
32 | // | |
33 | // default constructor | |
34 | // | |
35 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
36 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
37 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
38 | ||
39 | fCapacity = capacity; | |
40 | fN =0; | |
41 | fX = new Float_t[fCapacity]; | |
42 | fY = new Float_t[fCapacity]; | |
43 | fZ = new Float_t[fCapacity]; | |
44 | fSy = new Float_t[fCapacity]; | |
45 | fSz = new Float_t[fCapacity]; | |
46 | fCovar = new TMatrixDSym(6); | |
47 | fConv = kFALSE; | |
48 | } | |
49 | ||
50 | AliRieman::AliRieman(const AliRieman &rieman):TObject(){ | |
51 | // | |
52 | // copy constructor | |
53 | // | |
54 | for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i]; | |
55 | for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; | |
56 | for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; | |
57 | fCapacity = rieman.fN; | |
58 | fN =rieman.fN; | |
59 | fCovar = new TMatrixDSym(*(rieman.fCovar)); | |
60 | fConv = rieman.fConv; | |
61 | fX = new Float_t[fN]; | |
62 | fY = new Float_t[fN]; | |
63 | fZ = new Float_t[fN]; | |
64 | fSy = new Float_t[fN]; | |
65 | fSz = new Float_t[fN]; | |
66 | for (Int_t i=0;i<rieman.fN;i++){ | |
67 | fX[i] = rieman.fX[i]; | |
68 | fY[i] = rieman.fY[i]; | |
69 | fZ[i] = rieman.fZ[i]; | |
70 | fSy[i] = rieman.fSy[i]; | |
71 | fSz[i] = rieman.fSz[i]; | |
72 | } | |
73 | } | |
74 | ||
75 | AliRieman::~AliRieman() | |
76 | { | |
77 | delete[]fX; | |
78 | delete[]fY; | |
79 | delete[]fZ; | |
80 | delete[]fSy; | |
81 | delete[]fSz; | |
82 | delete fCovar; | |
83 | } | |
84 | ||
85 | void AliRieman::Reset() | |
86 | { | |
87 | fN=0; | |
88 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
89 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
90 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
91 | fConv =kFALSE; | |
92 | } | |
93 | ||
94 | ||
95 | void AliRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz) | |
96 | { | |
97 | // | |
98 | // Rieman update | |
99 | // | |
100 | //------------------------------------------------------ | |
101 | // XY direction | |
102 | // | |
103 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
104 | // | |
105 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
106 | // | |
107 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
108 | // | |
109 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
110 | // | |
111 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
112 | // | |
113 | // final linear equation : a + u*b +t*c - v =0; | |
114 | // | |
115 | // Minimization : | |
116 | // | |
117 | // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; | |
118 | // | |
119 | // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) | |
120 | // | |
121 | // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t | |
122 | // | |
123 | if (fN==fCapacity-1) return; // out of capacity | |
124 | fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz; | |
125 | // | |
126 | // XY part | |
127 | // | |
128 | Double_t t = x*x+y*y; | |
129 | if (t<2) return; | |
130 | t = 1./t; | |
131 | Double_t u = 2.*x*t; | |
132 | Double_t v = 2.*y*t; | |
133 | Double_t error = 2.*sy*t; | |
134 | error *=error; | |
135 | Double_t weight = 1./error; | |
136 | fSumXY[0] +=weight; | |
137 | fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight; | |
138 | fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight; | |
139 | fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight; | |
140 | // | |
141 | // XZ part | |
142 | // | |
143 | weight = 1./sz; | |
144 | fSumXZ[0] +=weight; | |
145 | fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x; | |
146 | fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z; | |
147 | // | |
148 | fN++; | |
149 | } | |
150 | ||
151 | ||
152 | ||
153 | ||
154 | ||
155 | void AliRieman::UpdatePol(){ | |
156 | // | |
157 | // Rieman update | |
158 | // | |
159 | // | |
160 | if (fN==0) return; | |
161 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
162 | Int_t conv=0; | |
163 | // | |
164 | // XY part | |
165 | // | |
166 | TMatrixDSym smatrix(3); | |
167 | TMatrixD sums(1,3); | |
168 | // | |
169 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
170 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
171 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
172 | ||
173 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
174 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
175 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; | |
176 | smatrix.Invert(); | |
177 | if (smatrix.IsValid()){ | |
178 | for (Int_t i=0;i<3;i++) | |
179 | for (Int_t j=0;j<=i;j++){ | |
180 | (*fCovar)(i,j)=smatrix(i,j); | |
181 | } | |
182 | TMatrixD res = sums*smatrix; | |
183 | fParams[0] = res(0,0); | |
184 | fParams[1] = res(0,1); | |
185 | fParams[2] = res(0,2); | |
186 | conv++; | |
187 | } | |
188 | // | |
189 | // XZ part | |
190 | // | |
191 | TMatrixDSym smatrixz(3); | |
192 | smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2]; | |
193 | smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3]; | |
194 | smatrixz(2,2) = fSumXZ[4]; | |
195 | smatrixz.Invert(); | |
196 | if (smatrixz.IsValid()){ | |
197 | sums(0,0) = fSumXZ[5]; | |
198 | sums(0,1) = fSumXZ[6]; | |
199 | sums(0,2) = fSumXZ[7]; | |
200 | TMatrixD res = sums*smatrixz; | |
201 | fParams[3] = res(0,0); | |
202 | fParams[4] = res(0,1); | |
203 | fParams[5] = res(0,2); | |
204 | for (Int_t i=0;i<3;i++) | |
205 | for (Int_t j=0;j<=i;j++){ | |
206 | (*fCovar)(i+2,j+2)=smatrixz(i,j); | |
207 | } | |
208 | conv++; | |
209 | } | |
210 | ||
211 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
212 | // | |
213 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
214 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
215 | // | |
216 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
217 | // | |
218 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
219 | // final linear equation : a + u*b +t*c - v =0; | |
220 | // | |
221 | // fParam[0] = 1/y0 | |
222 | // fParam[1] = -x0/y0 | |
223 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
224 | if (conv>1) fConv =kTRUE; | |
225 | else | |
226 | fConv=kFALSE; | |
227 | } | |
228 | ||
229 | void AliRieman::Update(){ | |
230 | // | |
231 | // Rieman update | |
232 | // | |
233 | // | |
234 | if (fN==0) return; | |
235 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
236 | Int_t conv=0; | |
237 | // | |
238 | // XY part | |
239 | // | |
240 | TMatrixDSym smatrix(3); | |
241 | TMatrixD sums(1,3); | |
242 | // | |
243 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
244 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
245 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
246 | ||
247 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
248 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
249 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; | |
250 | smatrix.Invert(); | |
251 | if (smatrix.IsValid()){ | |
252 | for (Int_t i=0;i<3;i++) | |
253 | for (Int_t j=0;j<3;j++){ | |
254 | (*fCovar)(i,j)=smatrix(i,j); | |
255 | } | |
256 | TMatrixD res = sums*smatrix; | |
257 | fParams[0] = res(0,0); | |
258 | fParams[1] = res(0,1); | |
259 | fParams[2] = res(0,2); | |
260 | conv++; | |
261 | } | |
262 | if (conv==0){ | |
263 | fConv = kFALSE; //non converged | |
264 | return; | |
265 | } | |
266 | // | |
267 | // XZ part | |
268 | // | |
269 | Double_t x0 = -fParams[1]/fParams[0]; | |
270 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
271 | Double_t sumXZ[9]; | |
272 | ||
273 | for (Int_t i=0;i<9;i++) sumXZ[i]=0; | |
274 | for (Int_t i=0;i<fN;i++){ | |
275 | Double_t phi = TMath::ASin((fX[i]-x0)*Rm1); | |
276 | Double_t phi0 = TMath::ASin((0.-x0)*Rm1); | |
277 | Double_t weight = 1/fSz[i]; | |
278 | weight *=weight; | |
279 | Double_t dphi = (phi-phi0)/Rm1; | |
280 | sumXZ[0] +=weight; | |
281 | sumXZ[1] +=weight*dphi; | |
282 | sumXZ[2] +=weight*dphi*dphi; | |
283 | sumXZ[3] +=weight*fZ[i]; | |
284 | sumXZ[4] +=weight*dphi*fZ[i]; | |
285 | ||
286 | } | |
287 | ||
288 | TMatrixDSym smatrixz(2); | |
289 | TMatrixD sumsz(1,2); | |
290 | smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2]; | |
291 | smatrixz.Invert(); | |
292 | if (smatrixz.IsValid()){ | |
293 | sumsz(0,0) = sumXZ[3]; | |
294 | sumsz(0,1) = sumXZ[4]; | |
295 | TMatrixD res = sumsz*smatrixz; | |
296 | fParams[3] = res(0,0); | |
297 | fParams[4] = res(0,1); | |
298 | for (Int_t i=0;i<2;i++) | |
299 | for (Int_t j=0;j<2;j++){ | |
300 | (*fCovar)(i+3,j+3)=smatrixz(i,j); | |
301 | } | |
302 | conv++; | |
303 | } | |
304 | ||
305 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
306 | // | |
307 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
308 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
309 | // | |
310 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
311 | // | |
312 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
313 | // final linear equation : a + u*b +t*c - v =0; | |
314 | // | |
315 | // fParam[0] = 1/y0 | |
316 | // fParam[1] = -x0/y0 | |
317 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
318 | if (conv>1) fConv =kTRUE; | |
319 | else | |
320 | fConv=kFALSE; | |
321 | } | |
322 | ||
323 | ||
324 | ||
325 | Double_t AliRieman::GetYat(Double_t x){ | |
326 | if (!fConv) return 0.; | |
327 | Double_t res2 = (x*fParams[0]+fParams[1]); | |
328 | res2*=res2; | |
329 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
330 | if (res2<0) return 0; | |
331 | res2 = TMath::Sqrt(res2); | |
332 | res2 = (1-res2)/fParams[0]; | |
333 | return res2; | |
334 | } | |
335 | ||
336 | Double_t AliRieman::GetDYat(Double_t x){ | |
337 | if (!fConv) return 0.; | |
338 | Double_t x0 = -fParams[1]/fParams[0]; | |
339 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0; | |
340 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
341 | if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0; | |
342 | Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0)); | |
343 | if (fParams[0]<0) res*=-1.; | |
344 | return res; | |
345 | } | |
346 | ||
347 | ||
348 | ||
349 | Double_t AliRieman::GetZat(Double_t x){ | |
350 | if (!fConv) return 0.; | |
351 | Double_t x0 = -fParams[1]/fParams[0]; | |
352 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
353 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
354 | Double_t phi = TMath::ASin((x-x0)*Rm1); | |
355 | Double_t phi0 = TMath::ASin((0.-x0)*Rm1); | |
356 | Double_t dphi = (phi-phi0); | |
357 | Double_t res = fParams[3]+fParams[4]*dphi/Rm1; | |
358 | return res; | |
359 | } | |
360 | ||
361 | Double_t AliRieman::GetDZat(Double_t x){ | |
362 | if (!fConv) return 0.; | |
363 | Double_t x0 = -fParams[1]/fParams[0]; | |
364 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
365 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
366 | Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1)); | |
367 | return res; | |
368 | } | |
369 | ||
370 | ||
371 | Double_t AliRieman::GetC(){ | |
372 | return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
373 | } | |
374 | ||
375 |