]>
Commit | Line | Data |
---|---|---|
049179df | 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 | /* $Id$ */ | |
17 | // Class for global helix fit of a track | |
18 | // Author: M.Ivanov | |
19 | // The method uses the following idea: | |
20 | //------------------------------------------------------ | |
21 | // XY direction | |
22 | // | |
23 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
24 | // | |
25 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
26 | // | |
27 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
28 | // | |
29 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
30 | // | |
31 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
32 | // | |
33 | // final linear equation : a + u*b +t*c - v =0; | |
34 | // | |
35 | // Minimization : | |
36 | // | |
37 | // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; | |
38 | // | |
39 | // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) | |
40 | // | |
41 | // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t | |
42 | ||
43 | ||
cabb917f | 44 | #include "TMatrixDSym.h" |
049179df | 45 | //#include "TDecompChol.h" |
cabb917f | 46 | #include "TMatrixD.h" |
47 | #include "AliRieman.h" | |
cabb917f | 48 | |
49 | ClassImp(AliRieman) | |
50 | ||
51 | ||
52 | ||
75e3794b | 53 | AliRieman::AliRieman() : |
54 | TObject(), | |
55 | fCapacity(0), | |
56 | fN(0), | |
57 | fX(0), | |
58 | fY(0), | |
59 | fZ(0), | |
60 | fSy(0), | |
61 | fSz(0), | |
62 | fCovar(0), | |
63 | fCovarPolY(0), | |
64 | fCovarPolZ(0), | |
65 | fSumZZ(0), | |
66 | fChi2(0), | |
67 | fChi2Y(0), | |
68 | fChi2Z(0), | |
69 | fConv(kFALSE) | |
70 | { | |
cabb917f | 71 | // |
72 | // default constructor | |
73 | // | |
74 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
75 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
76 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
0100c5fc | 77 | for (Int_t i=0;i<6;i++) { |
78 | fSumPolY[i]=0; | |
79 | fSumPolZ[i]=0; | |
80 | } | |
cabb917f | 81 | } |
82 | ||
83 | ||
75e3794b | 84 | AliRieman::AliRieman(Int_t capacity) : |
85 | TObject(), | |
86 | fCapacity(capacity), | |
87 | fN(0), | |
88 | fX(new Double_t[fCapacity]), | |
89 | fY(new Double_t[fCapacity]), | |
90 | fZ(new Double_t[fCapacity]), | |
91 | fSy(new Double_t[fCapacity]), | |
92 | fSz(new Double_t[fCapacity]), | |
93 | fCovar(new TMatrixDSym(6)), | |
94 | fCovarPolY(new TMatrixDSym(3)), | |
95 | fCovarPolZ(new TMatrixDSym(2)), | |
96 | fSumZZ(0), | |
97 | fChi2(0), | |
98 | fChi2Y(0), | |
99 | fChi2Z(0), | |
100 | fConv(kFALSE) | |
101 | { | |
cabb917f | 102 | // |
103 | // default constructor | |
104 | // | |
105 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
106 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
107 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
0100c5fc | 108 | for (Int_t i=0;i<6;i++) { |
109 | fSumPolY[i]=0; | |
110 | fSumPolZ[i]=0; | |
111 | } | |
cabb917f | 112 | } |
113 | ||
75e3794b | 114 | AliRieman::AliRieman(const AliRieman &rieman): |
115 | TObject(rieman), | |
116 | fCapacity(rieman.fN), | |
117 | fN(rieman.fN), | |
118 | fX(new Double_t[fN]), | |
119 | fY(new Double_t[fN]), | |
120 | fZ(new Double_t[fN]), | |
121 | fSy(new Double_t[fN]), | |
122 | fSz(new Double_t[fN]), | |
123 | fCovar(new TMatrixDSym(*(rieman.fCovar))), | |
124 | fCovarPolY(new TMatrixDSym(*(rieman.fCovarPolY))), | |
125 | fCovarPolZ(new TMatrixDSym(*(rieman.fCovarPolZ))), | |
126 | fSumZZ(rieman.fSumZZ), | |
127 | fChi2(rieman.fChi2), | |
128 | fChi2Y(rieman.fChi2Y), | |
129 | fChi2Z(rieman.fChi2Z), | |
130 | fConv(rieman.fConv) | |
131 | ||
132 | { | |
cabb917f | 133 | // |
134 | // copy constructor | |
135 | // | |
136 | for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i]; | |
137 | for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; | |
138 | for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; | |
0100c5fc | 139 | for (Int_t i=0;i<6;i++) { |
140 | fSumPolY[i]=rieman.fSumPolY[i]; | |
141 | fSumPolZ[i]=rieman.fSumPolZ[i]; | |
142 | } | |
cabb917f | 143 | for (Int_t i=0;i<rieman.fN;i++){ |
144 | fX[i] = rieman.fX[i]; | |
145 | fY[i] = rieman.fY[i]; | |
146 | fZ[i] = rieman.fZ[i]; | |
147 | fSy[i] = rieman.fSy[i]; | |
148 | fSz[i] = rieman.fSz[i]; | |
149 | } | |
150 | } | |
151 | ||
152 | AliRieman::~AliRieman() | |
153 | { | |
049179df | 154 | // |
155 | // Destructor | |
156 | // | |
cabb917f | 157 | delete[]fX; |
158 | delete[]fY; | |
159 | delete[]fZ; | |
160 | delete[]fSy; | |
161 | delete[]fSz; | |
162 | delete fCovar; | |
0100c5fc | 163 | delete fCovarPolY; |
164 | delete fCovarPolZ; | |
cabb917f | 165 | } |
166 | ||
167 | void AliRieman::Reset() | |
168 | { | |
049179df | 169 | // |
170 | // Reset all the data members | |
171 | // | |
cabb917f | 172 | fN=0; |
173 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
174 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
0100c5fc | 175 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; |
176 | for (Int_t i=0;i<6;i++) { | |
177 | fSumPolY[i]=0; | |
178 | fSumPolZ[i]=0; | |
179 | } | |
8849da04 | 180 | fSumZZ =0; |
cabb917f | 181 | fConv =kFALSE; |
182 | } | |
183 | ||
184 | ||
049179df | 185 | void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz) |
cabb917f | 186 | { |
187 | // | |
188 | // Rieman update | |
189 | // | |
190 | //------------------------------------------------------ | |
191 | // XY direction | |
192 | // | |
193 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
194 | // | |
195 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
196 | // | |
049179df | 197 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 |
cabb917f | 198 | // |
199 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
200 | // | |
201 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
202 | // | |
203 | // final linear equation : a + u*b +t*c - v =0; | |
204 | // | |
205 | // Minimization : | |
206 | // | |
207 | // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; | |
208 | // | |
209 | // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) | |
210 | // | |
211 | // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t | |
212 | // | |
213 | if (fN==fCapacity-1) return; // out of capacity | |
214 | fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz; | |
215 | // | |
216 | // XY part | |
217 | // | |
218 | Double_t t = x*x+y*y; | |
219 | if (t<2) return; | |
220 | t = 1./t; | |
221 | Double_t u = 2.*x*t; | |
222 | Double_t v = 2.*y*t; | |
223 | Double_t error = 2.*sy*t; | |
224 | error *=error; | |
225 | Double_t weight = 1./error; | |
226 | fSumXY[0] +=weight; | |
227 | fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight; | |
228 | fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight; | |
229 | fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight; | |
230 | // | |
231 | // XZ part | |
232 | // | |
8849da04 | 233 | weight = 1./(sz*sz); |
cabb917f | 234 | fSumXZ[0] +=weight; |
8849da04 | 235 | Double_t r = TMath::Sqrt(x*x+y*y); |
236 | fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z; | |
237 | // Now the auxulary sums | |
238 | fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24; | |
239 | fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24); | |
240 | fSumZZ += z*z*weight; | |
cabb917f | 241 | // |
0100c5fc | 242 | // sum accumulation for rough error estimates of the track extrapolation error |
243 | // | |
244 | Double_t maty = 1./(sy*sy); | |
245 | Double_t matz = 1./(sz*sz); | |
246 | for (Int_t i=0;i<5; i++){ | |
247 | fSumPolY[i] += maty; | |
248 | fSumPolZ[i] += matz; | |
249 | maty *= x; | |
250 | matz *= x; | |
251 | } | |
cabb917f | 252 | fN++; |
253 | } | |
254 | ||
255 | ||
cabb917f | 256 | void AliRieman::UpdatePol(){ |
257 | // | |
258 | // Rieman update | |
259 | // | |
260 | // | |
261 | if (fN==0) return; | |
262 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
263 | Int_t conv=0; | |
264 | // | |
265 | // XY part | |
266 | // | |
267 | TMatrixDSym smatrix(3); | |
268 | TMatrixD sums(1,3); | |
269 | // | |
270 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
271 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
272 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
273 | ||
274 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
275 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
0100c5fc | 276 | smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7]; |
cabb917f | 277 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; |
278 | smatrix.Invert(); | |
279 | if (smatrix.IsValid()){ | |
280 | for (Int_t i=0;i<3;i++) | |
281 | for (Int_t j=0;j<=i;j++){ | |
282 | (*fCovar)(i,j)=smatrix(i,j); | |
283 | } | |
284 | TMatrixD res = sums*smatrix; | |
285 | fParams[0] = res(0,0); | |
286 | fParams[1] = res(0,1); | |
287 | fParams[2] = res(0,2); | |
288 | conv++; | |
289 | } | |
290 | // | |
291 | // XZ part | |
292 | // | |
0100c5fc | 293 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
294 | ||
295 | // fSumXZ[1] += fSumXZ[5]*rm1*rm1; | |
296 | // fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1; | |
297 | // fSumXZ[4] += fSumXZ[7]*rm1*rm1; | |
298 | Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1; | |
299 | Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1; | |
300 | Double_t sum4 = fSumXZ[4] + fSumXZ[7]*rm1*rm1; | |
8849da04 | 301 | |
302 | TMatrixDSym smatrixz(2); | |
0100c5fc | 303 | // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2]; |
304 | //smatrixz(1,0) = fSumXZ[1]; | |
305 | smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2; | |
306 | smatrixz(1,0) = sum1; | |
cabb917f | 307 | smatrixz.Invert(); |
8849da04 | 308 | TMatrixD sumsxz(1,2); |
cabb917f | 309 | if (smatrixz.IsValid()){ |
8849da04 | 310 | sumsxz(0,0) = fSumXZ[3]; |
0100c5fc | 311 | // sumsxz(0,1) = fSumXZ[4]; |
312 | sumsxz(0,1) = sum4; | |
8849da04 | 313 | TMatrixD res = sumsxz*smatrixz; |
cabb917f | 314 | fParams[3] = res(0,0); |
315 | fParams[4] = res(0,1); | |
8849da04 | 316 | fParams[5] = 0; |
317 | for (Int_t i=0;i<2;i++) | |
cabb917f | 318 | for (Int_t j=0;j<=i;j++){ |
8849da04 | 319 | (*fCovar)(i+3,j+3)=smatrixz(i,j); |
320 | } | |
cabb917f | 321 | conv++; |
322 | } | |
0100c5fc | 323 | UpdateCovariancePol(); |
cabb917f | 324 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> |
325 | // | |
326 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
327 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
328 | // | |
329 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
330 | // | |
331 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
332 | // final linear equation : a + u*b +t*c - v =0; | |
333 | // | |
334 | // fParam[0] = 1/y0 | |
335 | // fParam[1] = -x0/y0 | |
336 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
337 | if (conv>1) fConv =kTRUE; | |
338 | else | |
339 | fConv=kFALSE; | |
340 | } | |
341 | ||
342 | void AliRieman::Update(){ | |
343 | // | |
344 | // Rieman update | |
345 | // | |
346 | // | |
347 | if (fN==0) return; | |
348 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
349 | Int_t conv=0; | |
350 | // | |
351 | // XY part | |
352 | // | |
353 | TMatrixDSym smatrix(3); | |
354 | TMatrixD sums(1,3); | |
355 | // | |
356 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
357 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
358 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
359 | ||
360 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
361 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
049179df | 362 | // |
363 | smatrix(1,0) = fSumXY[1]; | |
364 | smatrix(2,0) = fSumXY[3]; | |
365 | smatrix(2,1) = fSumXY[7]; | |
366 | ||
cabb917f | 367 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; |
049179df | 368 | //TDecompChol decomp(smatrix); |
369 | //decomp.SetTol(1.0e-14); | |
370 | //smatrix = | |
371 | //decomp.Invert(); | |
cabb917f | 372 | smatrix.Invert(); |
373 | if (smatrix.IsValid()){ | |
374 | for (Int_t i=0;i<3;i++) | |
375 | for (Int_t j=0;j<3;j++){ | |
376 | (*fCovar)(i,j)=smatrix(i,j); | |
377 | } | |
378 | TMatrixD res = sums*smatrix; | |
379 | fParams[0] = res(0,0); | |
380 | fParams[1] = res(0,1); | |
381 | fParams[2] = res(0,2); | |
382 | conv++; | |
383 | } | |
384 | if (conv==0){ | |
385 | fConv = kFALSE; //non converged | |
386 | return; | |
387 | } | |
049179df | 388 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){ |
389 | fConv = kFALSE; //non converged | |
390 | return; | |
391 | } | |
cabb917f | 392 | // |
393 | // XZ part | |
394 | // | |
395 | Double_t x0 = -fParams[1]/fParams[0]; | |
049179df | 396 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.); |
cabb917f | 397 | Double_t sumXZ[9]; |
398 | ||
399 | for (Int_t i=0;i<9;i++) sumXZ[i]=0; | |
400 | for (Int_t i=0;i<fN;i++){ | |
049179df | 401 | Double_t phi = TMath::ASin((fX[i]-x0)*rm1); |
402 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
cabb917f | 403 | Double_t weight = 1/fSz[i]; |
404 | weight *=weight; | |
049179df | 405 | Double_t dphi = (phi-phi0)/rm1; |
cabb917f | 406 | sumXZ[0] +=weight; |
407 | sumXZ[1] +=weight*dphi; | |
408 | sumXZ[2] +=weight*dphi*dphi; | |
409 | sumXZ[3] +=weight*fZ[i]; | |
410 | sumXZ[4] +=weight*dphi*fZ[i]; | |
411 | ||
412 | } | |
413 | ||
414 | TMatrixDSym smatrixz(2); | |
415 | TMatrixD sumsz(1,2); | |
416 | smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2]; | |
049179df | 417 | smatrixz(1,0) = sumXZ[1]; |
cabb917f | 418 | smatrixz.Invert(); |
419 | if (smatrixz.IsValid()){ | |
420 | sumsz(0,0) = sumXZ[3]; | |
421 | sumsz(0,1) = sumXZ[4]; | |
422 | TMatrixD res = sumsz*smatrixz; | |
423 | fParams[3] = res(0,0); | |
424 | fParams[4] = res(0,1); | |
425 | for (Int_t i=0;i<2;i++) | |
426 | for (Int_t j=0;j<2;j++){ | |
427 | (*fCovar)(i+3,j+3)=smatrixz(i,j); | |
428 | } | |
429 | conv++; | |
430 | } | |
0100c5fc | 431 | UpdateCovariancePol(); |
cabb917f | 432 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> |
433 | // | |
434 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
049179df | 435 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 |
cabb917f | 436 | // |
437 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
438 | // | |
439 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
440 | // final linear equation : a + u*b +t*c - v =0; | |
441 | // | |
442 | // fParam[0] = 1/y0 | |
443 | // fParam[1] = -x0/y0 | |
444 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
445 | if (conv>1) fConv =kTRUE; | |
446 | else | |
447 | fConv=kFALSE; | |
049179df | 448 | fChi2Y = CalcChi2Y(); |
449 | fChi2Z = CalcChi2Z(); | |
450 | fChi2 = fChi2Y +fChi2Z; | |
cabb917f | 451 | } |
452 | ||
0100c5fc | 453 | void AliRieman::UpdateCovariancePol(){ |
454 | // | |
455 | // covariance matrices for rough extrapolation error estimate | |
456 | // covariance matrix - get error estimates at point x = 0 | |
457 | // ! NOTE - error estimates is very rough and it is valid only if dl<<R | |
458 | // when dl is the distance between first and last point | |
459 | // | |
460 | // | |
461 | (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2]; | |
462 | (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3]; | |
463 | (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4]; | |
464 | fCovarPolY->Invert(); | |
465 | // | |
466 | ||
467 | (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1]; | |
468 | (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2]; | |
469 | fCovarPolZ->Invert(); | |
470 | // | |
471 | } | |
472 | ||
473 | Double_t AliRieman::GetErrY(Double_t x) const{ | |
474 | // | |
475 | // P0' = P0 + P1 * x + P2 * x^2 | |
476 | // P1' = P1 + P2 * x | |
477 | // P2' = + P2 | |
478 | TMatrixD trans(3,3); | |
479 | trans(0,0) = 1; | |
480 | trans(0,1) = x; | |
481 | trans(0,2) = x*x; | |
482 | trans(1,1) = 1; | |
483 | trans(1,2) = x; | |
484 | trans(2,2) = 1; | |
485 | // | |
486 | TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY); | |
487 | covarX*=trans.T(); | |
488 | return covarX(0,0); | |
489 | } | |
490 | ||
491 | Double_t AliRieman::GetErrZ(Double_t x) const{ | |
492 | // | |
493 | // assumption error of curvature determination neggligible | |
494 | // | |
495 | // P0' = P0 + P1 * x | |
496 | // P1' = P1 | |
497 | TMatrixD trans(2,2); | |
498 | trans(0,0) = 1; | |
499 | trans(0,1) = x; | |
500 | trans(1,1) = 1; | |
501 | // | |
502 | TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ); | |
503 | covarX*=trans.T(); | |
504 | return covarX(0,0); | |
505 | } | |
506 | ||
507 | Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){ | |
508 | // | |
509 | // Get external parameters | |
510 | // + approximative covariance | |
511 | // 0 1 2 3 4 // y - local y | |
512 | // 5 6 7 8 // z - local z | |
513 | // 9 10 11 // snp - local sin fi | |
514 | // 12 13 // tgl - deep angle | |
515 | // 14 // C - curvature | |
516 | Double_t sign = (GetC()>0) ? 1.:-1.; | |
517 | ||
518 | params[0] = GetYat(xref); | |
519 | params[1] = GetZat(xref); | |
520 | params[2] = TMath::Sin(TMath::ATan(GetDYat(xref))); | |
521 | params[3] = sign*fParams[4]; | |
522 | params[4] = GetC(); | |
523 | // | |
524 | // covariance matrix in y | |
525 | // | |
526 | TMatrixD transY(3,3); | |
527 | transY(0,0) = 1; | |
528 | transY(0,1) = xref; | |
529 | transY(0,2) = xref*xref; | |
530 | transY(1,1) = 1; | |
531 | transY(1,2) = xref; | |
532 | transY(2,2) = 1; | |
533 | TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY); | |
534 | covarY*=transY.T(); | |
535 | // | |
536 | TMatrixD transZ(2,2); | |
537 | transZ(0,0) = 1; | |
538 | transZ(0,1) = xref; | |
539 | transZ(1,1) = 1; | |
540 | TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ); | |
541 | covarZ*=transZ.T(); | |
542 | // | |
543 | // C ~ 2*P2 - in rotated frame | |
544 | // | |
545 | covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2); | |
546 | covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0; | |
547 | covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2); | |
548 | covar[12] = covarZ(1,1); covar[13]= 0; | |
549 | covar[14] = 2.*covarY(2,2); | |
550 | // | |
551 | return fConv; | |
552 | } | |
553 | ||
554 | ||
cabb917f | 555 | |
556 | ||
049179df | 557 | Double_t AliRieman::GetYat(Double_t x) const { |
558 | // | |
559 | // get y at given x position | |
560 | // | |
cabb917f | 561 | if (!fConv) return 0.; |
562 | Double_t res2 = (x*fParams[0]+fParams[1]); | |
563 | res2*=res2; | |
564 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
565 | if (res2<0) return 0; | |
566 | res2 = TMath::Sqrt(res2); | |
567 | res2 = (1-res2)/fParams[0]; | |
568 | return res2; | |
569 | } | |
570 | ||
049179df | 571 | Double_t AliRieman::GetDYat(Double_t x) const { |
572 | // | |
573 | // get dy/dx at given x position | |
574 | // | |
cabb917f | 575 | if (!fConv) return 0.; |
576 | Double_t x0 = -fParams[1]/fParams[0]; | |
577 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0; | |
049179df | 578 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
579 | if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0; | |
580 | Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0)); | |
cabb917f | 581 | if (fParams[0]<0) res*=-1.; |
582 | return res; | |
583 | } | |
584 | ||
585 | ||
586 | ||
049179df | 587 | Double_t AliRieman::GetZat(Double_t x) const { |
588 | // | |
589 | // get z at given x position | |
590 | // | |
cabb917f | 591 | if (!fConv) return 0.; |
592 | Double_t x0 = -fParams[1]/fParams[0]; | |
593 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
049179df | 594 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
595 | Double_t phi = TMath::ASin((x-x0)*rm1); | |
596 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
cabb917f | 597 | Double_t dphi = (phi-phi0); |
049179df | 598 | Double_t res = fParams[3]+fParams[4]*dphi/rm1; |
cabb917f | 599 | return res; |
600 | } | |
601 | ||
049179df | 602 | Double_t AliRieman::GetDZat(Double_t x) const { |
603 | // | |
604 | // get dz/dx at given x postion | |
605 | // | |
cabb917f | 606 | if (!fConv) return 0.; |
607 | Double_t x0 = -fParams[1]/fParams[0]; | |
608 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
049179df | 609 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
610 | Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1)); | |
cabb917f | 611 | return res; |
612 | } | |
613 | ||
614 | ||
8849da04 | 615 | //_____________________________________________________________________________ |
616 | Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const | |
617 | { | |
618 | // | |
619 | // Returns position given radius | |
620 | // | |
621 | if (!fConv) return kFALSE; | |
622 | Double_t res2 = (r*fParams[0]+fParams[1]); | |
623 | res2*=res2; | |
624 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
625 | if (res2<0) return kFALSE; | |
626 | res2 = TMath::Sqrt(res2); | |
627 | res2 = (1-res2)/fParams[0]; | |
628 | ||
629 | if (!fConv) return kFALSE; | |
630 | Double_t x0 = -fParams[1]/fParams[0]; | |
631 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
632 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
633 | Double_t phi = TMath::ASin((r-x0)*rm1); | |
634 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
635 | Double_t dphi = (phi-phi0); | |
636 | Double_t res = fParams[3]+fParams[4]*dphi/rm1; | |
637 | ||
638 | Double_t sin = TMath::Sin(alpha); | |
639 | Double_t cos = TMath::Cos(alpha); | |
640 | xyz[0] = r*cos - res2*sin; | |
641 | xyz[1] = res2*cos + r*sin; | |
642 | xyz[2] = res; | |
643 | ||
644 | return kTRUE; | |
645 | } | |
646 | ||
647 | ||
049179df | 648 | Double_t AliRieman::GetC() const{ |
649 | // | |
650 | // get curvature | |
651 | // | |
cabb917f | 652 | return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
653 | } | |
654 | ||
049179df | 655 | Double_t AliRieman::CalcChi2Y() const{ |
656 | // | |
657 | // calculate chi2 for Y | |
658 | // | |
659 | Double_t sumchi2 = 0; | |
660 | for (Int_t i=0;i<fN;i++){ | |
661 | Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i]; | |
662 | sumchi2+=chi2*chi2; | |
663 | } | |
664 | return sumchi2; | |
665 | } | |
666 | ||
667 | ||
668 | Double_t AliRieman::CalcChi2Z() const{ | |
669 | // | |
670 | // calculate chi2 for Z | |
671 | // | |
672 | Double_t sumchi2 = 0; | |
673 | for (Int_t i=0;i<fN;i++){ | |
674 | Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i]; | |
675 | sumchi2+=chi2*chi2; | |
676 | } | |
677 | return sumchi2; | |
678 | } | |
679 | ||
680 | Double_t AliRieman::CalcChi2() const{ | |
681 | // | |
682 | // sum chi2 in both coord - supposing Y and Z independent | |
683 | // | |
684 | return CalcChi2Y()+CalcChi2Z(); | |
685 | } | |
686 | ||
687 | AliRieman * AliRieman::MakeResiduals() const{ | |
688 | // | |
689 | // create residual structure - ONLY for debug purposes | |
690 | // | |
691 | AliRieman *rieman = new AliRieman(fN); | |
692 | for (Int_t i=0;i<fN;i++){ | |
693 | rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]); | |
694 | } | |
695 | return rieman; | |
696 | } | |
697 | ||
cabb917f | 698 |