AliRoot Core  3dc7879 (3dc7879)
AliRieman.cxx
Go to the documentation of this file.
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 
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<5;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<5;i++) {
110  fSumPolY[i]=0;
111  fSumPolZ[i]=0;
112  }
113 }
114 
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<5;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 
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 
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<5;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 
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  }
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 
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  }
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 
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 
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 
Double_t * fX
Definition: AliRieman.h:58
void Update()
Definition: AliRieman.cxx:343
TMatrixDSym * fCovarPolZ
Definition: AliRieman.h:66
Double_t CalcChi2() const
Definition: AliRieman.cxx:682
Bool_t fConv
Definition: AliRieman.h:75
Double_t fChi2Y
Definition: AliRieman.h:73
Double_t fSumXY[9]
Definition: AliRieman.h:67
Double_t GetErrY(Double_t x) const
Definition: AliRieman.cxx:474
Double_t * fZ
Definition: AliRieman.h:60
Double_t * fY
Definition: AliRieman.h:59
void AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
Definition: AliRieman.cxx:186
Double_t fSumPolY[5]
Definition: AliRieman.h:69
Double_t GetC() const
Definition: AliRieman.cxx:650
Double_t fChi2Z
Definition: AliRieman.h:74
Bool_t GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
Definition: AliRieman.cxx:618
bool trans(const AliFMDIndex &x, const AliFMDIndex &y, const AliFMDIndex &z)
Definition: TestIndex.C:94
Double_t * fSz
Definition: AliRieman.h:62
Double_t fSumXZ[9]
Definition: AliRieman.h:68
Double_t chi2
Definition: AnalyzeLaser.C:7
Bool_t GetExternalParameters(Double_t xref, Double_t *params, Double_t *covar)
Definition: AliRieman.cxx:508
Double_t GetYat(Double_t x) const
Definition: AliRieman.cxx:558
Double_t GetErrZ(Double_t x) const
Definition: AliRieman.cxx:492
Int_t fN
Definition: AliRieman.h:57
Double_t GetDYat(Double_t x) const
Definition: AliRieman.cxx:572
TMatrixDSym * fCovarPolY
Definition: AliRieman.h:65
Double_t fChi2
Definition: AliRieman.h:72
Double_t CalcChi2Z() const
Definition: AliRieman.cxx:670
void Reset()
Definition: AliRieman.cxx:168
Double_t GetDZat(Double_t x) const
Definition: AliRieman.cxx:604
Double_t fParams[6]
Definition: AliRieman.h:63
TMatrixDSym * fCovar
Definition: AliRieman.h:64
AliRieman * MakeResiduals() const
Definition: AliRieman.cxx:689
void UpdateCovariancePol()
Definition: AliRieman.cxx:454
Double_t fSumZZ
Definition: AliRieman.h:71
Double_t * fSy
Definition: AliRieman.h:61
Int_t fCapacity
Definition: AliRieman.h:56
void UpdatePol()
Definition: AliRieman.cxx:257
Double_t fSumPolZ[5]
Definition: AliRieman.h:70
Double_t CalcChi2Y() const
Definition: AliRieman.cxx:657
Double_t GetZat(Double_t x) const
Definition: AliRieman.cxx:589
void res(Char_t i)
Definition: Resolution.C:2
for(Int_t itree=0;itree< arrInputTreesDistortionCalib->GetEntriesFast();++itree)