]> 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
Releasing the covariance matrix of intersection point along the track direction
[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
3010c308 44#include <TMatrixDSym.h>
45#include <TMath.h>
46#include <TMatrixD.h>
47
cabb917f 48#include "AliRieman.h"
cabb917f 49
50ClassImp(AliRieman)
51
52
53
75e3794b 54AliRieman::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 85AliRieman::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 115AliRieman::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
153AliRieman::~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
168void 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 186void 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 257void 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
343void 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 454void 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
474Double_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
492Double_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
508Bool_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 558Double_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 572Double_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 588Double_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 603Double_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//_____________________________________________________________________________
617Bool_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 649Double_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 656Double_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
669Double_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
681Double_t AliRieman::CalcChi2() const{
682 //
683 // sum chi2 in both coord - supposing Y and Z independent
684 //
685 return CalcChi2Y()+CalcChi2Z();
686}
687
688AliRieman * 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