]> 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
Support for automatic pedestal estimation per pad (Marian and Matevz).
[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
53AliRieman::AliRieman(){
54 //
55 // default constructor
56 //
57 for (Int_t i=0;i<6;i++) fParams[i] = 0;
58 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
59 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
0100c5fc 60 for (Int_t i=0;i<6;i++) {
61 fSumPolY[i]=0;
62 fSumPolZ[i]=0;
63 }
8849da04 64 fSumZZ = 0;
cabb917f 65 fCapacity = 0;
66 fN =0;
67 fX = 0;
68 fY = 0;
69 fZ = 0;
70 fSy = 0;
71 fSz = 0;
049179df 72 fChi2 = 0;
73 fChi2Y = 0;
74 fChi2Z = 0;
cabb917f 75 fCovar = 0;
0100c5fc 76 fCovarPolY = 0;
77 fCovarPolZ = 0;
cabb917f 78 fConv = kFALSE;
79}
80
81
82AliRieman::AliRieman(Int_t capacity){
83 //
84 // default constructor
85 //
86 for (Int_t i=0;i<6;i++) fParams[i] = 0;
87 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
88 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
0100c5fc 89 for (Int_t i=0;i<6;i++) {
90 fSumPolY[i]=0;
91 fSumPolZ[i]=0;
92 }
8849da04 93 fSumZZ =0;
cabb917f 94 fCapacity = capacity;
95 fN =0;
049179df 96 fX = new Double_t[fCapacity];
97 fY = new Double_t[fCapacity];
98 fZ = new Double_t[fCapacity];
99 fSy = new Double_t[fCapacity];
100 fSz = new Double_t[fCapacity];
cabb917f 101 fCovar = new TMatrixDSym(6);
0100c5fc 102 fCovarPolY = new TMatrixDSym(3);
103 fCovarPolZ = new TMatrixDSym(2);
049179df 104 fChi2 = 0;
105 fChi2Y = 0;
106 fChi2Z = 0;
cabb917f 107 fConv = kFALSE;
108}
109
049179df 110AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
cabb917f 111 //
112 // copy constructor
113 //
114 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
115 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
116 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
0100c5fc 117 for (Int_t i=0;i<6;i++) {
118 fSumPolY[i]=rieman.fSumPolY[i];
119 fSumPolZ[i]=rieman.fSumPolZ[i];
120 }
8849da04 121 fSumZZ = rieman.fSumZZ;
cabb917f 122 fCapacity = rieman.fN;
123 fN =rieman.fN;
124 fCovar = new TMatrixDSym(*(rieman.fCovar));
0100c5fc 125 fCovarPolY = new TMatrixDSym(*(rieman.fCovarPolY));
126 fCovarPolZ = new TMatrixDSym(*(rieman.fCovarPolZ));
cabb917f 127 fConv = rieman.fConv;
049179df 128 fX = new Double_t[fN];
129 fY = new Double_t[fN];
130 fZ = new Double_t[fN];
131 fSy = new Double_t[fN];
132 fSz = new Double_t[fN];
cabb917f 133 for (Int_t i=0;i<rieman.fN;i++){
134 fX[i] = rieman.fX[i];
135 fY[i] = rieman.fY[i];
136 fZ[i] = rieman.fZ[i];
137 fSy[i] = rieman.fSy[i];
138 fSz[i] = rieman.fSz[i];
139 }
140}
141
142AliRieman::~AliRieman()
143{
049179df 144 //
145 // Destructor
146 //
cabb917f 147 delete[]fX;
148 delete[]fY;
149 delete[]fZ;
150 delete[]fSy;
151 delete[]fSz;
152 delete fCovar;
0100c5fc 153 delete fCovarPolY;
154 delete fCovarPolZ;
cabb917f 155}
156
157void AliRieman::Reset()
158{
049179df 159 //
160 // Reset all the data members
161 //
cabb917f 162 fN=0;
163 for (Int_t i=0;i<6;i++) fParams[i] = 0;
164 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
0100c5fc 165 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
166 for (Int_t i=0;i<6;i++) {
167 fSumPolY[i]=0;
168 fSumPolZ[i]=0;
169 }
8849da04 170 fSumZZ =0;
cabb917f 171 fConv =kFALSE;
172}
173
174
049179df 175void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
cabb917f 176{
177 //
178 // Rieman update
179 //
180 //------------------------------------------------------
181 // XY direction
182 //
183 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
184 //
185 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
186 //
049179df 187 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
cabb917f 188 //
189 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
190 //
191 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
192 //
193 // final linear equation : a + u*b +t*c - v =0;
194 //
195 // Minimization :
196 //
197 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
198 //
199 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
200 //
201 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
202 //
203 if (fN==fCapacity-1) return; // out of capacity
204 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
205 //
206 // XY part
207 //
208 Double_t t = x*x+y*y;
209 if (t<2) return;
210 t = 1./t;
211 Double_t u = 2.*x*t;
212 Double_t v = 2.*y*t;
213 Double_t error = 2.*sy*t;
214 error *=error;
215 Double_t weight = 1./error;
216 fSumXY[0] +=weight;
217 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
218 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
219 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
220 //
221 // XZ part
222 //
8849da04 223 weight = 1./(sz*sz);
cabb917f 224 fSumXZ[0] +=weight;
8849da04 225 Double_t r = TMath::Sqrt(x*x+y*y);
226 fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
227 // Now the auxulary sums
228 fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
229 fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
230 fSumZZ += z*z*weight;
cabb917f 231 //
0100c5fc 232 // sum accumulation for rough error estimates of the track extrapolation error
233 //
234 Double_t maty = 1./(sy*sy);
235 Double_t matz = 1./(sz*sz);
236 for (Int_t i=0;i<5; i++){
237 fSumPolY[i] += maty;
238 fSumPolZ[i] += matz;
239 maty *= x;
240 matz *= x;
241 }
cabb917f 242 fN++;
243}
244
245
cabb917f 246void AliRieman::UpdatePol(){
247 //
248 // Rieman update
249 //
250 //
251 if (fN==0) return;
252 for (Int_t i=0;i<6;i++)fParams[i]=0;
253 Int_t conv=0;
254 //
255 // XY part
256 //
257 TMatrixDSym smatrix(3);
258 TMatrixD sums(1,3);
259 //
260 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
261 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
262 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
263
264 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
265 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
0100c5fc 266 smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
cabb917f 267 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
268 smatrix.Invert();
269 if (smatrix.IsValid()){
270 for (Int_t i=0;i<3;i++)
271 for (Int_t j=0;j<=i;j++){
272 (*fCovar)(i,j)=smatrix(i,j);
273 }
274 TMatrixD res = sums*smatrix;
275 fParams[0] = res(0,0);
276 fParams[1] = res(0,1);
277 fParams[2] = res(0,2);
278 conv++;
279 }
280 //
281 // XZ part
282 //
0100c5fc 283 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
284
285// fSumXZ[1] += fSumXZ[5]*rm1*rm1;
286// fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
287// fSumXZ[4] += fSumXZ[7]*rm1*rm1;
288 Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1;
289 Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
290 Double_t sum4 = fSumXZ[4] + fSumXZ[7]*rm1*rm1;
8849da04 291
292 TMatrixDSym smatrixz(2);
0100c5fc 293 // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
294 //smatrixz(1,0) = fSumXZ[1];
295 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2;
296 smatrixz(1,0) = sum1;
cabb917f 297 smatrixz.Invert();
8849da04 298 TMatrixD sumsxz(1,2);
cabb917f 299 if (smatrixz.IsValid()){
8849da04 300 sumsxz(0,0) = fSumXZ[3];
0100c5fc 301 // sumsxz(0,1) = fSumXZ[4];
302 sumsxz(0,1) = sum4;
8849da04 303 TMatrixD res = sumsxz*smatrixz;
cabb917f 304 fParams[3] = res(0,0);
305 fParams[4] = res(0,1);
8849da04 306 fParams[5] = 0;
307 for (Int_t i=0;i<2;i++)
cabb917f 308 for (Int_t j=0;j<=i;j++){
8849da04 309 (*fCovar)(i+3,j+3)=smatrixz(i,j);
310 }
cabb917f 311 conv++;
312 }
0100c5fc 313 UpdateCovariancePol();
cabb917f 314 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
315 //
316 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
317 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
318 //
319 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
320 //
321 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
322 // final linear equation : a + u*b +t*c - v =0;
323 //
324 // fParam[0] = 1/y0
325 // fParam[1] = -x0/y0
326 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
327 if (conv>1) fConv =kTRUE;
328 else
329 fConv=kFALSE;
330}
331
332void AliRieman::Update(){
333 //
334 // Rieman update
335 //
336 //
337 if (fN==0) return;
338 for (Int_t i=0;i<6;i++)fParams[i]=0;
339 Int_t conv=0;
340 //
341 // XY part
342 //
343 TMatrixDSym smatrix(3);
344 TMatrixD sums(1,3);
345 //
346 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
347 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
348 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
349
350 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
351 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
049179df 352 //
353 smatrix(1,0) = fSumXY[1];
354 smatrix(2,0) = fSumXY[3];
355 smatrix(2,1) = fSumXY[7];
356
cabb917f 357 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
049179df 358 //TDecompChol decomp(smatrix);
359 //decomp.SetTol(1.0e-14);
360 //smatrix =
361 //decomp.Invert();
cabb917f 362 smatrix.Invert();
363 if (smatrix.IsValid()){
364 for (Int_t i=0;i<3;i++)
365 for (Int_t j=0;j<3;j++){
366 (*fCovar)(i,j)=smatrix(i,j);
367 }
368 TMatrixD res = sums*smatrix;
369 fParams[0] = res(0,0);
370 fParams[1] = res(0,1);
371 fParams[2] = res(0,2);
372 conv++;
373 }
374 if (conv==0){
375 fConv = kFALSE; //non converged
376 return;
377 }
049179df 378 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
379 fConv = kFALSE; //non converged
380 return;
381 }
cabb917f 382 //
383 // XZ part
384 //
385 Double_t x0 = -fParams[1]/fParams[0];
049179df 386 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
cabb917f 387 Double_t sumXZ[9];
388
389 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
390 for (Int_t i=0;i<fN;i++){
049179df 391 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
392 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
cabb917f 393 Double_t weight = 1/fSz[i];
394 weight *=weight;
049179df 395 Double_t dphi = (phi-phi0)/rm1;
cabb917f 396 sumXZ[0] +=weight;
397 sumXZ[1] +=weight*dphi;
398 sumXZ[2] +=weight*dphi*dphi;
399 sumXZ[3] +=weight*fZ[i];
400 sumXZ[4] +=weight*dphi*fZ[i];
401
402 }
403
404 TMatrixDSym smatrixz(2);
405 TMatrixD sumsz(1,2);
406 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
049179df 407 smatrixz(1,0) = sumXZ[1];
cabb917f 408 smatrixz.Invert();
409 if (smatrixz.IsValid()){
410 sumsz(0,0) = sumXZ[3];
411 sumsz(0,1) = sumXZ[4];
412 TMatrixD res = sumsz*smatrixz;
413 fParams[3] = res(0,0);
414 fParams[4] = res(0,1);
415 for (Int_t i=0;i<2;i++)
416 for (Int_t j=0;j<2;j++){
417 (*fCovar)(i+3,j+3)=smatrixz(i,j);
418 }
419 conv++;
420 }
0100c5fc 421 UpdateCovariancePol();
cabb917f 422 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
423 //
424 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
049179df 425 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
cabb917f 426 //
427 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
428 //
429 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
430 // final linear equation : a + u*b +t*c - v =0;
431 //
432 // fParam[0] = 1/y0
433 // fParam[1] = -x0/y0
434 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
435 if (conv>1) fConv =kTRUE;
436 else
437 fConv=kFALSE;
049179df 438 fChi2Y = CalcChi2Y();
439 fChi2Z = CalcChi2Z();
440 fChi2 = fChi2Y +fChi2Z;
cabb917f 441}
442
0100c5fc 443void AliRieman::UpdateCovariancePol(){
444 //
445 // covariance matrices for rough extrapolation error estimate
446 // covariance matrix - get error estimates at point x = 0
447 // ! NOTE - error estimates is very rough and it is valid only if dl<<R
448 // when dl is the distance between first and last point
449 //
450 //
451 (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2];
452 (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3];
453 (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4];
454 fCovarPolY->Invert();
455 //
456
457 (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1];
458 (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
459 fCovarPolZ->Invert();
460 //
461}
462
463Double_t AliRieman::GetErrY(Double_t x) const{
464 //
465 // P0' = P0 + P1 * x + P2 * x^2
466 // P1' = P1 + P2 * x
467 // P2' = + P2
468 TMatrixD trans(3,3);
469 trans(0,0) = 1;
470 trans(0,1) = x;
471 trans(0,2) = x*x;
472 trans(1,1) = 1;
473 trans(1,2) = x;
474 trans(2,2) = 1;
475 //
476 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
477 covarX*=trans.T();
478 return covarX(0,0);
479}
480
481Double_t AliRieman::GetErrZ(Double_t x) const{
482 //
483 // assumption error of curvature determination neggligible
484 //
485 // P0' = P0 + P1 * x
486 // P1' = P1
487 TMatrixD trans(2,2);
488 trans(0,0) = 1;
489 trans(0,1) = x;
490 trans(1,1) = 1;
491 //
492 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
493 covarX*=trans.T();
494 return covarX(0,0);
495}
496
497Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
498 //
499 // Get external parameters
500 // + approximative covariance
501 // 0 1 2 3 4 // y - local y
502 // 5 6 7 8 // z - local z
503 // 9 10 11 // snp - local sin fi
504 // 12 13 // tgl - deep angle
505 // 14 // C - curvature
506 Double_t sign = (GetC()>0) ? 1.:-1.;
507
508 params[0] = GetYat(xref);
509 params[1] = GetZat(xref);
510 params[2] = TMath::Sin(TMath::ATan(GetDYat(xref)));
511 params[3] = sign*fParams[4];
512 params[4] = GetC();
513 //
514 // covariance matrix in y
515 //
516 TMatrixD transY(3,3);
517 transY(0,0) = 1;
518 transY(0,1) = xref;
519 transY(0,2) = xref*xref;
520 transY(1,1) = 1;
521 transY(1,2) = xref;
522 transY(2,2) = 1;
523 TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
524 covarY*=transY.T();
525 //
526 TMatrixD transZ(2,2);
527 transZ(0,0) = 1;
528 transZ(0,1) = xref;
529 transZ(1,1) = 1;
530 TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
531 covarZ*=transZ.T();
532 //
533 // C ~ 2*P2 - in rotated frame
534 //
535 covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2);
536 covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0;
537 covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2);
538 covar[12] = covarZ(1,1); covar[13]= 0;
539 covar[14] = 2.*covarY(2,2);
540 //
541 return fConv;
542}
543
544
cabb917f 545
546
049179df 547Double_t AliRieman::GetYat(Double_t x) const {
548 //
549 // get y at given x position
550 //
cabb917f 551 if (!fConv) return 0.;
552 Double_t res2 = (x*fParams[0]+fParams[1]);
553 res2*=res2;
554 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
555 if (res2<0) return 0;
556 res2 = TMath::Sqrt(res2);
557 res2 = (1-res2)/fParams[0];
558 return res2;
559}
560
049179df 561Double_t AliRieman::GetDYat(Double_t x) const {
562 //
563 // get dy/dx at given x position
564 //
cabb917f 565 if (!fConv) return 0.;
566 Double_t x0 = -fParams[1]/fParams[0];
567 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
049179df 568 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
569 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
570 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
cabb917f 571 if (fParams[0]<0) res*=-1.;
572 return res;
573}
574
575
576
049179df 577Double_t AliRieman::GetZat(Double_t x) const {
578 //
579 // get z at given x position
580 //
cabb917f 581 if (!fConv) return 0.;
582 Double_t x0 = -fParams[1]/fParams[0];
583 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
049179df 584 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
585 Double_t phi = TMath::ASin((x-x0)*rm1);
586 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
cabb917f 587 Double_t dphi = (phi-phi0);
049179df 588 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
cabb917f 589 return res;
590}
591
049179df 592Double_t AliRieman::GetDZat(Double_t x) const {
593 //
594 // get dz/dx at given x postion
595 //
cabb917f 596 if (!fConv) return 0.;
597 Double_t x0 = -fParams[1]/fParams[0];
598 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
049179df 599 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
600 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
cabb917f 601 return res;
602}
603
604
8849da04 605//_____________________________________________________________________________
606Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
607{
608 //
609 // Returns position given radius
610 //
611 if (!fConv) return kFALSE;
612 Double_t res2 = (r*fParams[0]+fParams[1]);
613 res2*=res2;
614 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
615 if (res2<0) return kFALSE;
616 res2 = TMath::Sqrt(res2);
617 res2 = (1-res2)/fParams[0];
618
619 if (!fConv) return kFALSE;
620 Double_t x0 = -fParams[1]/fParams[0];
621 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
622 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
623 Double_t phi = TMath::ASin((r-x0)*rm1);
624 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
625 Double_t dphi = (phi-phi0);
626 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
627
628 Double_t sin = TMath::Sin(alpha);
629 Double_t cos = TMath::Cos(alpha);
630 xyz[0] = r*cos - res2*sin;
631 xyz[1] = res2*cos + r*sin;
632 xyz[2] = res;
633
634 return kTRUE;
635}
636
637
049179df 638Double_t AliRieman::GetC() const{
639 //
640 // get curvature
641 //
cabb917f 642 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
643}
644
049179df 645Double_t AliRieman::CalcChi2Y() const{
646 //
647 // calculate chi2 for Y
648 //
649 Double_t sumchi2 = 0;
650 for (Int_t i=0;i<fN;i++){
651 Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
652 sumchi2+=chi2*chi2;
653 }
654 return sumchi2;
655}
656
657
658Double_t AliRieman::CalcChi2Z() const{
659 //
660 // calculate chi2 for Z
661 //
662 Double_t sumchi2 = 0;
663 for (Int_t i=0;i<fN;i++){
664 Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
665 sumchi2+=chi2*chi2;
666 }
667 return sumchi2;
668}
669
670Double_t AliRieman::CalcChi2() const{
671 //
672 // sum chi2 in both coord - supposing Y and Z independent
673 //
674 return CalcChi2Y()+CalcChi2Z();
675}
676
677AliRieman * AliRieman::MakeResiduals() const{
678 //
679 // create residual structure - ONLY for debug purposes
680 //
681 AliRieman *rieman = new AliRieman(fN);
682 for (Int_t i=0;i<fN;i++){
683 rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
684 }
685 return rieman;
686}
687
cabb917f 688