- Three classes by MinJung Kweon AliHFEpriVtx, AliHFEsecVtx and AliHFEmcQA for primar...
[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 <TMath.h>
46 #include <TMatrixD.h>
47
48 #include "AliRieman.h"
49
50 ClassImp(AliRieman)
51
52
53
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 {
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;
78   for (Int_t i=0;i<6;i++) {
79     fSumPolY[i]=0;
80     fSumPolZ[i]=0;
81   }
82 }
83
84
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 {
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;
109   for (Int_t i=0;i<6;i++) {
110     fSumPolY[i]=0;
111     fSumPolZ[i]=0;
112   }
113 }
114
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 {
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];
140   for (Int_t i=0;i<6;i++) {
141     fSumPolY[i]=rieman.fSumPolY[i];
142     fSumPolZ[i]=rieman.fSumPolZ[i];
143   }
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 {
155   //
156   // Destructor
157   //
158   delete[]fX;
159   delete[]fY;
160   delete[]fZ;
161   delete[]fSy;
162   delete[]fSz;
163   delete fCovar;
164   delete fCovarPolY;
165   delete fCovarPolZ;
166 }
167
168 void AliRieman::Reset()
169 {
170   //
171   // Reset all the data members
172   //
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;
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   }
181   fSumZZ =0;
182   fConv =kFALSE;
183 }
184
185
186 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
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   //
198   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
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   //
234   weight = 1./(sz*sz);
235   fSumXZ[0] +=weight;
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;
242   //
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   }
253   fN++;
254 }
255
256
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];
277   smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
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   //
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;
302   
303   TMatrixDSym     smatrixz(2);
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;
308   smatrixz.Invert();
309   TMatrixD        sumsxz(1,2);
310   if (smatrixz.IsValid()){
311     sumsxz(0,0)    = fSumXZ[3];
312     //    sumsxz(0,1)    = fSumXZ[4];
313     sumsxz(0,1)    = sum4;
314     TMatrixD res = sumsxz*smatrixz;
315     fParams[3] = res(0,0);
316     fParams[4] = res(0,1);
317     fParams[5] = 0;
318     for (Int_t i=0;i<2;i++)
319       for (Int_t j=0;j<=i;j++){
320         (*fCovar)(i+3,j+3)=smatrixz(i,j);
321       }
322     conv++;
323   }
324   UpdateCovariancePol();
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];
363   //
364   smatrix(1,0) = fSumXY[1];
365   smatrix(2,0) = fSumXY[3];
366   smatrix(2,1) = fSumXY[7];
367
368   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
369   //TDecompChol decomp(smatrix);
370   //decomp.SetTol(1.0e-14);
371   //smatrix = 
372   //decomp.Invert();
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   }
389   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
390     fConv = kFALSE;   //non converged
391     return;
392   }
393   //
394   // XZ part
395   //
396   Double_t x0 = -fParams[1]/fParams[0];
397   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.); 
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++){
402     Double_t phi  = TMath::ASin((fX[i]-x0)*rm1);
403     Double_t phi0 = TMath::ASin((0.-x0)*rm1);
404     Double_t weight = 1/fSz[i];
405     weight *=weight;
406     Double_t dphi = (phi-phi0)/rm1;
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];
418   smatrixz(1,0) = sumXZ[1];
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   }
432   UpdateCovariancePol();
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;  ==>
436   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
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;
449   fChi2Y = CalcChi2Y();
450   fChi2Z = CalcChi2Z();
451   fChi2  = fChi2Y +fChi2Z;
452 }
453
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
556
557
558 Double_t AliRieman::GetYat(Double_t x) const {
559   //
560   // get y at given x position
561   //
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
572 Double_t AliRieman::GetDYat(Double_t x) const {
573   //
574   // get dy/dx at given x position
575   //
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;
579   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
580   Double_t arg = (1./rm1-(x-x0))*(1./rm1+(x-x0));
581   if ( arg <= 0) return 0;
582   Double_t res = (x-x0)/TMath::Sqrt(arg);
583   if (fParams[0]<0) res*=-1.;
584   return res;
585 }
586
587
588
589 Double_t AliRieman::GetZat(Double_t x) const {
590   //
591   // get z at given x position
592   //
593   if (!fConv) return 0.;
594   Double_t x0 = -fParams[1]/fParams[0];
595   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
596   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
597   Double_t phi  = TMath::ASin((x-x0)*rm1);
598   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
599   Double_t dphi = (phi-phi0);
600   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
601   return res;
602 }
603
604 Double_t AliRieman::GetDZat(Double_t x) const {
605   //
606   // get dz/dx at given x postion
607   //
608   if (!fConv) return 0.;
609   Double_t x0 = -fParams[1]/fParams[0]; 
610   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
611   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
612   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
613   return res;
614 }
615
616
617 //_____________________________________________________________________________
618 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
619 {
620   //
621   // Returns position given radius
622   //
623   if (!fConv) return kFALSE;
624   Double_t res2 = (r*fParams[0]+fParams[1]);
625   res2*=res2;
626   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
627   if (res2<0) return kFALSE;
628   res2 = TMath::Sqrt(res2);
629   res2 = (1-res2)/fParams[0];
630
631   if (!fConv) return kFALSE;
632   Double_t x0 = -fParams[1]/fParams[0];
633   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
634   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
635   Double_t phi  = TMath::ASin((r-x0)*rm1);
636   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
637   Double_t dphi = (phi-phi0);
638   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
639
640   Double_t sin = TMath::Sin(alpha);
641   Double_t cos = TMath::Cos(alpha);
642   xyz[0] = r*cos - res2*sin;
643   xyz[1] = res2*cos + r*sin;
644   xyz[2] = res;
645
646   return kTRUE;
647 }
648
649
650 Double_t AliRieman::GetC() const{
651   //
652   // get curvature
653   //
654   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
655 }
656
657 Double_t AliRieman::CalcChi2Y() const{ 
658   //
659   // calculate chi2 for Y
660   //
661   Double_t sumchi2 = 0;
662   for (Int_t i=0;i<fN;i++){
663     Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
664     sumchi2+=chi2*chi2;    
665   }  
666   return sumchi2;
667 }
668
669
670 Double_t AliRieman::CalcChi2Z() const{
671   //
672   // calculate chi2 for Z
673   //
674   Double_t sumchi2 = 0;
675   for (Int_t i=0;i<fN;i++){
676     Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
677     sumchi2+=chi2*chi2;    
678   }  
679   return sumchi2;
680 }
681
682 Double_t AliRieman::CalcChi2() const{
683   //
684   // sum chi2 in both coord - supposing Y and Z independent
685   //
686   return CalcChi2Y()+CalcChi2Z();
687 }
688
689 AliRieman * AliRieman::MakeResiduals() const{
690   //
691   // create residual structure - ONLY for debug purposes
692   //
693   AliRieman *rieman = new AliRieman(fN);  
694   for (Int_t i=0;i<fN;i++){
695     rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
696   }
697   return rieman;
698 }
699
700