]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliRieman.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes required by Effective C++
[u/mrichter/AliRoot.git] / STEER / AliRieman.cxx
CommitLineData
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
49ClassImp(AliRieman)
50
51
52
75e3794b 53AliRieman::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 84AliRieman::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 114AliRieman::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
152AliRieman::~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
167void 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 185void 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 256void 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
342void 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 453void 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
473Double_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
491Double_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
507Bool_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 557Double_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 571Double_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 587Double_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 602Double_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//_____________________________________________________________________________
616Bool_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 648Double_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 655Double_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
668Double_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
680Double_t AliRieman::CalcChi2() const{
681 //
682 // sum chi2 in both coord - supposing Y and Z independent
683 //
684 return CalcChi2Y()+CalcChi2Z();
685}
686
687AliRieman * 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