]>
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 | ||
53 | AliRieman::AliRieman(){ | |
54 | // | |
55 | // default constructor | |
56 | // | |
57 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
58 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
59 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
8849da04 | 60 | fSumZZ = 0; |
cabb917f | 61 | fCapacity = 0; |
62 | fN =0; | |
63 | fX = 0; | |
64 | fY = 0; | |
65 | fZ = 0; | |
66 | fSy = 0; | |
67 | fSz = 0; | |
049179df | 68 | fChi2 = 0; |
69 | fChi2Y = 0; | |
70 | fChi2Z = 0; | |
cabb917f | 71 | fCovar = 0; |
72 | fConv = kFALSE; | |
73 | } | |
74 | ||
75 | ||
76 | AliRieman::AliRieman(Int_t capacity){ | |
77 | // | |
78 | // default constructor | |
79 | // | |
80 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
81 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
82 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
8849da04 | 83 | fSumZZ =0; |
cabb917f | 84 | fCapacity = capacity; |
85 | fN =0; | |
049179df | 86 | fX = new Double_t[fCapacity]; |
87 | fY = new Double_t[fCapacity]; | |
88 | fZ = new Double_t[fCapacity]; | |
89 | fSy = new Double_t[fCapacity]; | |
90 | fSz = new Double_t[fCapacity]; | |
cabb917f | 91 | fCovar = new TMatrixDSym(6); |
049179df | 92 | fChi2 = 0; |
93 | fChi2Y = 0; | |
94 | fChi2Z = 0; | |
cabb917f | 95 | fConv = kFALSE; |
96 | } | |
97 | ||
049179df | 98 | AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){ |
cabb917f | 99 | // |
100 | // copy constructor | |
101 | // | |
102 | for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i]; | |
103 | for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; | |
104 | for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; | |
8849da04 | 105 | fSumZZ = rieman.fSumZZ; |
cabb917f | 106 | fCapacity = rieman.fN; |
107 | fN =rieman.fN; | |
108 | fCovar = new TMatrixDSym(*(rieman.fCovar)); | |
109 | fConv = rieman.fConv; | |
049179df | 110 | fX = new Double_t[fN]; |
111 | fY = new Double_t[fN]; | |
112 | fZ = new Double_t[fN]; | |
113 | fSy = new Double_t[fN]; | |
114 | fSz = new Double_t[fN]; | |
cabb917f | 115 | for (Int_t i=0;i<rieman.fN;i++){ |
116 | fX[i] = rieman.fX[i]; | |
117 | fY[i] = rieman.fY[i]; | |
118 | fZ[i] = rieman.fZ[i]; | |
119 | fSy[i] = rieman.fSy[i]; | |
120 | fSz[i] = rieman.fSz[i]; | |
121 | } | |
122 | } | |
123 | ||
124 | AliRieman::~AliRieman() | |
125 | { | |
049179df | 126 | // |
127 | // Destructor | |
128 | // | |
cabb917f | 129 | delete[]fX; |
130 | delete[]fY; | |
131 | delete[]fZ; | |
132 | delete[]fSy; | |
133 | delete[]fSz; | |
134 | delete fCovar; | |
135 | } | |
136 | ||
137 | void AliRieman::Reset() | |
138 | { | |
049179df | 139 | // |
140 | // Reset all the data members | |
141 | // | |
cabb917f | 142 | fN=0; |
143 | for (Int_t i=0;i<6;i++) fParams[i] = 0; | |
144 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
145 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; | |
8849da04 | 146 | fSumZZ =0; |
cabb917f | 147 | fConv =kFALSE; |
148 | } | |
149 | ||
150 | ||
049179df | 151 | void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz) |
cabb917f | 152 | { |
153 | // | |
154 | // Rieman update | |
155 | // | |
156 | //------------------------------------------------------ | |
157 | // XY direction | |
158 | // | |
159 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
160 | // | |
161 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
162 | // | |
049179df | 163 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 |
cabb917f | 164 | // |
165 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
166 | // | |
167 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
168 | // | |
169 | // final linear equation : a + u*b +t*c - v =0; | |
170 | // | |
171 | // Minimization : | |
172 | // | |
173 | // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; | |
174 | // | |
175 | // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) | |
176 | // | |
177 | // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t | |
178 | // | |
179 | if (fN==fCapacity-1) return; // out of capacity | |
180 | fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz; | |
181 | // | |
182 | // XY part | |
183 | // | |
184 | Double_t t = x*x+y*y; | |
185 | if (t<2) return; | |
186 | t = 1./t; | |
187 | Double_t u = 2.*x*t; | |
188 | Double_t v = 2.*y*t; | |
189 | Double_t error = 2.*sy*t; | |
190 | error *=error; | |
191 | Double_t weight = 1./error; | |
192 | fSumXY[0] +=weight; | |
193 | fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight; | |
194 | fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight; | |
195 | fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight; | |
196 | // | |
197 | // XZ part | |
198 | // | |
8849da04 | 199 | weight = 1./(sz*sz); |
cabb917f | 200 | fSumXZ[0] +=weight; |
8849da04 | 201 | Double_t r = TMath::Sqrt(x*x+y*y); |
202 | fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z; | |
203 | // Now the auxulary sums | |
204 | fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24; | |
205 | fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24); | |
206 | fSumZZ += z*z*weight; | |
cabb917f | 207 | // |
208 | fN++; | |
209 | } | |
210 | ||
211 | ||
212 | ||
213 | ||
214 | ||
215 | void AliRieman::UpdatePol(){ | |
216 | // | |
217 | // Rieman update | |
218 | // | |
219 | // | |
220 | if (fN==0) return; | |
221 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
222 | Int_t conv=0; | |
223 | // | |
224 | // XY part | |
225 | // | |
226 | TMatrixDSym smatrix(3); | |
227 | TMatrixD sums(1,3); | |
228 | // | |
229 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
230 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
231 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
232 | ||
233 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
234 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
235 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; | |
236 | smatrix.Invert(); | |
237 | if (smatrix.IsValid()){ | |
238 | for (Int_t i=0;i<3;i++) | |
239 | for (Int_t j=0;j<=i;j++){ | |
240 | (*fCovar)(i,j)=smatrix(i,j); | |
241 | } | |
242 | TMatrixD res = sums*smatrix; | |
243 | fParams[0] = res(0,0); | |
244 | fParams[1] = res(0,1); | |
245 | fParams[2] = res(0,2); | |
246 | conv++; | |
247 | } | |
248 | // | |
249 | // XZ part | |
250 | // | |
8849da04 | 251 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
252 | fSumXZ[1] += fSumXZ[5]*Rm1*Rm1; | |
253 | fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1; | |
254 | fSumXZ[4] += fSumXZ[7]*Rm1*Rm1; | |
255 | ||
256 | TMatrixDSym smatrixz(2); | |
257 | smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2]; | |
258 | smatrixz(1,0) = fSumXZ[1]; | |
cabb917f | 259 | smatrixz.Invert(); |
8849da04 | 260 | TMatrixD sumsxz(1,2); |
cabb917f | 261 | if (smatrixz.IsValid()){ |
8849da04 | 262 | sumsxz(0,0) = fSumXZ[3]; |
263 | sumsxz(0,1) = fSumXZ[4]; | |
264 | TMatrixD res = sumsxz*smatrixz; | |
cabb917f | 265 | fParams[3] = res(0,0); |
266 | fParams[4] = res(0,1); | |
8849da04 | 267 | fParams[5] = 0; |
268 | for (Int_t i=0;i<2;i++) | |
cabb917f | 269 | for (Int_t j=0;j<=i;j++){ |
8849da04 | 270 | (*fCovar)(i+3,j+3)=smatrixz(i,j); |
271 | } | |
cabb917f | 272 | conv++; |
273 | } | |
cabb917f | 274 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> |
275 | // | |
276 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
277 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
278 | // | |
279 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
280 | // | |
281 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
282 | // final linear equation : a + u*b +t*c - v =0; | |
283 | // | |
284 | // fParam[0] = 1/y0 | |
285 | // fParam[1] = -x0/y0 | |
286 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
287 | if (conv>1) fConv =kTRUE; | |
288 | else | |
289 | fConv=kFALSE; | |
290 | } | |
291 | ||
292 | void AliRieman::Update(){ | |
293 | // | |
294 | // Rieman update | |
295 | // | |
296 | // | |
297 | if (fN==0) return; | |
298 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
299 | Int_t conv=0; | |
300 | // | |
301 | // XY part | |
302 | // | |
303 | TMatrixDSym smatrix(3); | |
304 | TMatrixD sums(1,3); | |
305 | // | |
306 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
307 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
308 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
309 | ||
310 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
311 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
049179df | 312 | // |
313 | smatrix(1,0) = fSumXY[1]; | |
314 | smatrix(2,0) = fSumXY[3]; | |
315 | smatrix(2,1) = fSumXY[7]; | |
316 | ||
cabb917f | 317 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; |
049179df | 318 | //TDecompChol decomp(smatrix); |
319 | //decomp.SetTol(1.0e-14); | |
320 | //smatrix = | |
321 | //decomp.Invert(); | |
cabb917f | 322 | smatrix.Invert(); |
323 | if (smatrix.IsValid()){ | |
324 | for (Int_t i=0;i<3;i++) | |
325 | for (Int_t j=0;j<3;j++){ | |
326 | (*fCovar)(i,j)=smatrix(i,j); | |
327 | } | |
328 | TMatrixD res = sums*smatrix; | |
329 | fParams[0] = res(0,0); | |
330 | fParams[1] = res(0,1); | |
331 | fParams[2] = res(0,2); | |
332 | conv++; | |
333 | } | |
334 | if (conv==0){ | |
335 | fConv = kFALSE; //non converged | |
336 | return; | |
337 | } | |
049179df | 338 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){ |
339 | fConv = kFALSE; //non converged | |
340 | return; | |
341 | } | |
cabb917f | 342 | // |
343 | // XZ part | |
344 | // | |
345 | Double_t x0 = -fParams[1]/fParams[0]; | |
049179df | 346 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.); |
cabb917f | 347 | Double_t sumXZ[9]; |
348 | ||
349 | for (Int_t i=0;i<9;i++) sumXZ[i]=0; | |
350 | for (Int_t i=0;i<fN;i++){ | |
049179df | 351 | Double_t phi = TMath::ASin((fX[i]-x0)*rm1); |
352 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
cabb917f | 353 | Double_t weight = 1/fSz[i]; |
354 | weight *=weight; | |
049179df | 355 | Double_t dphi = (phi-phi0)/rm1; |
cabb917f | 356 | sumXZ[0] +=weight; |
357 | sumXZ[1] +=weight*dphi; | |
358 | sumXZ[2] +=weight*dphi*dphi; | |
359 | sumXZ[3] +=weight*fZ[i]; | |
360 | sumXZ[4] +=weight*dphi*fZ[i]; | |
361 | ||
362 | } | |
363 | ||
364 | TMatrixDSym smatrixz(2); | |
365 | TMatrixD sumsz(1,2); | |
366 | smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2]; | |
049179df | 367 | smatrixz(1,0) = sumXZ[1]; |
cabb917f | 368 | smatrixz.Invert(); |
369 | if (smatrixz.IsValid()){ | |
370 | sumsz(0,0) = sumXZ[3]; | |
371 | sumsz(0,1) = sumXZ[4]; | |
372 | TMatrixD res = sumsz*smatrixz; | |
373 | fParams[3] = res(0,0); | |
374 | fParams[4] = res(0,1); | |
375 | for (Int_t i=0;i<2;i++) | |
376 | for (Int_t j=0;j<2;j++){ | |
377 | (*fCovar)(i+3,j+3)=smatrixz(i,j); | |
378 | } | |
379 | conv++; | |
380 | } | |
381 | ||
382 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
383 | // | |
384 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
049179df | 385 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 |
cabb917f | 386 | // |
387 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
388 | // | |
389 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
390 | // final linear equation : a + u*b +t*c - v =0; | |
391 | // | |
392 | // fParam[0] = 1/y0 | |
393 | // fParam[1] = -x0/y0 | |
394 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
395 | if (conv>1) fConv =kTRUE; | |
396 | else | |
397 | fConv=kFALSE; | |
049179df | 398 | fChi2Y = CalcChi2Y(); |
399 | fChi2Z = CalcChi2Z(); | |
400 | fChi2 = fChi2Y +fChi2Z; | |
cabb917f | 401 | } |
402 | ||
403 | ||
404 | ||
049179df | 405 | Double_t AliRieman::GetYat(Double_t x) const { |
406 | // | |
407 | // get y at given x position | |
408 | // | |
cabb917f | 409 | if (!fConv) return 0.; |
410 | Double_t res2 = (x*fParams[0]+fParams[1]); | |
411 | res2*=res2; | |
412 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
413 | if (res2<0) return 0; | |
414 | res2 = TMath::Sqrt(res2); | |
415 | res2 = (1-res2)/fParams[0]; | |
416 | return res2; | |
417 | } | |
418 | ||
049179df | 419 | Double_t AliRieman::GetDYat(Double_t x) const { |
420 | // | |
421 | // get dy/dx at given x position | |
422 | // | |
cabb917f | 423 | if (!fConv) return 0.; |
424 | Double_t x0 = -fParams[1]/fParams[0]; | |
425 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0; | |
049179df | 426 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
427 | if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0; | |
428 | Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0)); | |
cabb917f | 429 | if (fParams[0]<0) res*=-1.; |
430 | return res; | |
431 | } | |
432 | ||
433 | ||
434 | ||
049179df | 435 | Double_t AliRieman::GetZat(Double_t x) const { |
436 | // | |
437 | // get z at given x position | |
438 | // | |
cabb917f | 439 | if (!fConv) return 0.; |
440 | Double_t x0 = -fParams[1]/fParams[0]; | |
441 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
049179df | 442 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
443 | Double_t phi = TMath::ASin((x-x0)*rm1); | |
444 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
cabb917f | 445 | Double_t dphi = (phi-phi0); |
049179df | 446 | Double_t res = fParams[3]+fParams[4]*dphi/rm1; |
cabb917f | 447 | return res; |
448 | } | |
449 | ||
049179df | 450 | Double_t AliRieman::GetDZat(Double_t x) const { |
451 | // | |
452 | // get dz/dx at given x postion | |
453 | // | |
cabb917f | 454 | if (!fConv) return 0.; |
455 | Double_t x0 = -fParams[1]/fParams[0]; | |
456 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
049179df | 457 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
458 | Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1)); | |
cabb917f | 459 | return res; |
460 | } | |
461 | ||
462 | ||
8849da04 | 463 | //_____________________________________________________________________________ |
464 | Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const | |
465 | { | |
466 | // | |
467 | // Returns position given radius | |
468 | // | |
469 | if (!fConv) return kFALSE; | |
470 | Double_t res2 = (r*fParams[0]+fParams[1]); | |
471 | res2*=res2; | |
472 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
473 | if (res2<0) return kFALSE; | |
474 | res2 = TMath::Sqrt(res2); | |
475 | res2 = (1-res2)/fParams[0]; | |
476 | ||
477 | if (!fConv) return kFALSE; | |
478 | Double_t x0 = -fParams[1]/fParams[0]; | |
479 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
480 | Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
481 | Double_t phi = TMath::ASin((r-x0)*rm1); | |
482 | Double_t phi0 = TMath::ASin((0.-x0)*rm1); | |
483 | Double_t dphi = (phi-phi0); | |
484 | Double_t res = fParams[3]+fParams[4]*dphi/rm1; | |
485 | ||
486 | Double_t sin = TMath::Sin(alpha); | |
487 | Double_t cos = TMath::Cos(alpha); | |
488 | xyz[0] = r*cos - res2*sin; | |
489 | xyz[1] = res2*cos + r*sin; | |
490 | xyz[2] = res; | |
491 | ||
492 | return kTRUE; | |
493 | } | |
494 | ||
495 | ||
049179df | 496 | Double_t AliRieman::GetC() const{ |
497 | // | |
498 | // get curvature | |
499 | // | |
cabb917f | 500 | return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); |
501 | } | |
502 | ||
049179df | 503 | Double_t AliRieman::CalcChi2Y() const{ |
504 | // | |
505 | // calculate chi2 for Y | |
506 | // | |
507 | Double_t sumchi2 = 0; | |
508 | for (Int_t i=0;i<fN;i++){ | |
509 | Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i]; | |
510 | sumchi2+=chi2*chi2; | |
511 | } | |
512 | return sumchi2; | |
513 | } | |
514 | ||
515 | ||
516 | Double_t AliRieman::CalcChi2Z() const{ | |
517 | // | |
518 | // calculate chi2 for Z | |
519 | // | |
520 | Double_t sumchi2 = 0; | |
521 | for (Int_t i=0;i<fN;i++){ | |
522 | Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i]; | |
523 | sumchi2+=chi2*chi2; | |
524 | } | |
525 | return sumchi2; | |
526 | } | |
527 | ||
528 | Double_t AliRieman::CalcChi2() const{ | |
529 | // | |
530 | // sum chi2 in both coord - supposing Y and Z independent | |
531 | // | |
532 | return CalcChi2Y()+CalcChi2Z(); | |
533 | } | |
534 | ||
535 | AliRieman * AliRieman::MakeResiduals() const{ | |
536 | // | |
537 | // create residual structure - ONLY for debug purposes | |
538 | // | |
539 | AliRieman *rieman = new AliRieman(fN); | |
540 | for (Int_t i=0;i<fN;i++){ | |
541 | rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]); | |
542 | } | |
543 | return rieman; | |
544 | } | |
545 | ||
cabb917f | 546 |