Changes required by Effective C++
[u/mrichter/AliRoot.git] / STEER / AliRieman.cxx
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
44 #include "TMatrixDSym.h"
45 //#include "TDecompChol.h"
46 #include "TMatrixD.h"
47 #include "AliRieman.h"
48
49 ClassImp(AliRieman)
50
51
52
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 {
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;
77   for (Int_t i=0;i<6;i++) {
78     fSumPolY[i]=0;
79     fSumPolZ[i]=0;
80   }
81 }
82
83
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 {
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;
108   for (Int_t i=0;i<6;i++) {
109     fSumPolY[i]=0;
110     fSumPolZ[i]=0;
111   }
112 }
113
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 {
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];
139   for (Int_t i=0;i<6;i++) {
140     fSumPolY[i]=rieman.fSumPolY[i];
141     fSumPolZ[i]=rieman.fSumPolZ[i];
142   }
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 {
154   //
155   // Destructor
156   //
157   delete[]fX;
158   delete[]fY;
159   delete[]fZ;
160   delete[]fSy;
161   delete[]fSz;
162   delete fCovar;
163   delete fCovarPolY;
164   delete fCovarPolZ;
165 }
166
167 void AliRieman::Reset()
168 {
169   //
170   // Reset all the data members
171   //
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;
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   }
180   fSumZZ =0;
181   fConv =kFALSE;
182 }
183
184
185 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
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   //
197   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
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   //
233   weight = 1./(sz*sz);
234   fSumXZ[0] +=weight;
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;
241   //
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   }
252   fN++;
253 }
254
255
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];
276   smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
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   //
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;
301   
302   TMatrixDSym     smatrixz(2);
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;
307   smatrixz.Invert();
308   TMatrixD        sumsxz(1,2);
309   if (smatrixz.IsValid()){
310     sumsxz(0,0)    = fSumXZ[3];
311     //    sumsxz(0,1)    = fSumXZ[4];
312     sumsxz(0,1)    = sum4;
313     TMatrixD res = sumsxz*smatrixz;
314     fParams[3] = res(0,0);
315     fParams[4] = res(0,1);
316     fParams[5] = 0;
317     for (Int_t i=0;i<2;i++)
318       for (Int_t j=0;j<=i;j++){
319         (*fCovar)(i+3,j+3)=smatrixz(i,j);
320       }
321     conv++;
322   }
323   UpdateCovariancePol();
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];
362   //
363   smatrix(1,0) = fSumXY[1];
364   smatrix(2,0) = fSumXY[3];
365   smatrix(2,1) = fSumXY[7];
366
367   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
368   //TDecompChol decomp(smatrix);
369   //decomp.SetTol(1.0e-14);
370   //smatrix = 
371   //decomp.Invert();
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   }
388   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
389     fConv = kFALSE;   //non converged
390     return;
391   }
392   //
393   // XZ part
394   //
395   Double_t x0 = -fParams[1]/fParams[0];
396   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.); 
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++){
401     Double_t phi  = TMath::ASin((fX[i]-x0)*rm1);
402     Double_t phi0 = TMath::ASin((0.-x0)*rm1);
403     Double_t weight = 1/fSz[i];
404     weight *=weight;
405     Double_t dphi = (phi-phi0)/rm1;
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];
417   smatrixz(1,0) = sumXZ[1];
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   }
431   UpdateCovariancePol();
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;  ==>
435   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
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;
448   fChi2Y = CalcChi2Y();
449   fChi2Z = CalcChi2Z();
450   fChi2  = fChi2Y +fChi2Z;
451 }
452
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
555
556
557 Double_t AliRieman::GetYat(Double_t x) const {
558   //
559   // get y at given x position
560   //
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
571 Double_t AliRieman::GetDYat(Double_t x) const {
572   //
573   // get dy/dx at given x position
574   //
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;
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));
581   if (fParams[0]<0) res*=-1.;
582   return res;
583 }
584
585
586
587 Double_t AliRieman::GetZat(Double_t x) const {
588   //
589   // get z at given x position
590   //
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;
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);
597   Double_t dphi = (phi-phi0);
598   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
599   return res;
600 }
601
602 Double_t AliRieman::GetDZat(Double_t x) const {
603   //
604   // get dz/dx at given x postion
605   //
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;
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));
611   return res;
612 }
613
614
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
648 Double_t AliRieman::GetC() const{
649   //
650   // get curvature
651   //
652   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
653 }
654
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
698