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