]> 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
Efficient C++ initialization of data members
[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;
60
61 fCapacity = 0;
62 fN =0;
63 fX = 0;
64 fY = 0;
65 fZ = 0;
66 fSy = 0;
67 fSz = 0;
049179df 68 fChi2 = 0;
69 fChi2Y = 0;
70 fChi2Z = 0;
cabb917f 71 fCovar = 0;
72 fConv = kFALSE;
73}
74
75
76AliRieman::AliRieman(Int_t capacity){
77 //
78 // default constructor
79 //
80 for (Int_t i=0;i<6;i++) fParams[i] = 0;
81 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
82 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
83
84 fCapacity = capacity;
85 fN =0;
049179df 86 fX = new Double_t[fCapacity];
87 fY = new Double_t[fCapacity];
88 fZ = new Double_t[fCapacity];
89 fSy = new Double_t[fCapacity];
90 fSz = new Double_t[fCapacity];
cabb917f 91 fCovar = new TMatrixDSym(6);
049179df 92 fChi2 = 0;
93 fChi2Y = 0;
94 fChi2Z = 0;
cabb917f 95 fConv = kFALSE;
96}
97
049179df 98AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
cabb917f 99 //
100 // copy constructor
101 //
102 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
103 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
104 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
105 fCapacity = rieman.fN;
106 fN =rieman.fN;
107 fCovar = new TMatrixDSym(*(rieman.fCovar));
108 fConv = rieman.fConv;
049179df 109 fX = new Double_t[fN];
110 fY = new Double_t[fN];
111 fZ = new Double_t[fN];
112 fSy = new Double_t[fN];
113 fSz = new Double_t[fN];
cabb917f 114 for (Int_t i=0;i<rieman.fN;i++){
115 fX[i] = rieman.fX[i];
116 fY[i] = rieman.fY[i];
117 fZ[i] = rieman.fZ[i];
118 fSy[i] = rieman.fSy[i];
119 fSz[i] = rieman.fSz[i];
120 }
121}
122
123AliRieman::~AliRieman()
124{
049179df 125 //
126 // Destructor
127 //
cabb917f 128 delete[]fX;
129 delete[]fY;
130 delete[]fZ;
131 delete[]fSy;
132 delete[]fSz;
133 delete fCovar;
134}
135
136void AliRieman::Reset()
137{
049179df 138 //
139 // Reset all the data members
140 //
cabb917f 141 fN=0;
142 for (Int_t i=0;i<6;i++) fParams[i] = 0;
143 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
144 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
145 fConv =kFALSE;
146}
147
148
049179df 149void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
cabb917f 150{
151 //
152 // Rieman update
153 //
154 //------------------------------------------------------
155 // XY direction
156 //
157 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
158 //
159 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
160 //
049179df 161 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
cabb917f 162 //
163 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
164 //
165 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
166 //
167 // final linear equation : a + u*b +t*c - v =0;
168 //
169 // Minimization :
170 //
171 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
172 //
173 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
174 //
175 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
176 //
177 if (fN==fCapacity-1) return; // out of capacity
178 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
179 //
180 // XY part
181 //
182 Double_t t = x*x+y*y;
183 if (t<2) return;
184 t = 1./t;
185 Double_t u = 2.*x*t;
186 Double_t v = 2.*y*t;
187 Double_t error = 2.*sy*t;
188 error *=error;
189 Double_t weight = 1./error;
190 fSumXY[0] +=weight;
191 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
192 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
193 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
194 //
195 // XZ part
196 //
197 weight = 1./sz;
198 fSumXZ[0] +=weight;
199 fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
200 fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
201 //
202 fN++;
203}
204
205
206
207
208
209void AliRieman::UpdatePol(){
210 //
211 // Rieman update
212 //
213 //
214 if (fN==0) return;
215 for (Int_t i=0;i<6;i++)fParams[i]=0;
216 Int_t conv=0;
217 //
218 // XY part
219 //
220 TMatrixDSym smatrix(3);
221 TMatrixD sums(1,3);
222 //
223 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
224 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
225 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
226
227 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
228 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
229 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
230 smatrix.Invert();
231 if (smatrix.IsValid()){
232 for (Int_t i=0;i<3;i++)
233 for (Int_t j=0;j<=i;j++){
234 (*fCovar)(i,j)=smatrix(i,j);
235 }
236 TMatrixD res = sums*smatrix;
237 fParams[0] = res(0,0);
238 fParams[1] = res(0,1);
239 fParams[2] = res(0,2);
240 conv++;
241 }
242 //
243 // XZ part
244 //
245 TMatrixDSym smatrixz(3);
246 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
247 smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
248 smatrixz(2,2) = fSumXZ[4];
249 smatrixz.Invert();
250 if (smatrixz.IsValid()){
251 sums(0,0) = fSumXZ[5];
252 sums(0,1) = fSumXZ[6];
253 sums(0,2) = fSumXZ[7];
254 TMatrixD res = sums*smatrixz;
255 fParams[3] = res(0,0);
256 fParams[4] = res(0,1);
257 fParams[5] = res(0,2);
258 for (Int_t i=0;i<3;i++)
259 for (Int_t j=0;j<=i;j++){
260 (*fCovar)(i+2,j+2)=smatrixz(i,j);
261 }
262 conv++;
263 }
264
265 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
266 //
267 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
268 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
269 //
270 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
271 //
272 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
273 // final linear equation : a + u*b +t*c - v =0;
274 //
275 // fParam[0] = 1/y0
276 // fParam[1] = -x0/y0
277 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
278 if (conv>1) fConv =kTRUE;
279 else
280 fConv=kFALSE;
281}
282
283void AliRieman::Update(){
284 //
285 // Rieman update
286 //
287 //
288 if (fN==0) return;
289 for (Int_t i=0;i<6;i++)fParams[i]=0;
290 Int_t conv=0;
291 //
292 // XY part
293 //
294 TMatrixDSym smatrix(3);
295 TMatrixD sums(1,3);
296 //
297 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
298 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
299 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
300
301 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
302 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
049179df 303 //
304 smatrix(1,0) = fSumXY[1];
305 smatrix(2,0) = fSumXY[3];
306 smatrix(2,1) = fSumXY[7];
307
cabb917f 308 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
049179df 309 //TDecompChol decomp(smatrix);
310 //decomp.SetTol(1.0e-14);
311 //smatrix =
312 //decomp.Invert();
cabb917f 313 smatrix.Invert();
314 if (smatrix.IsValid()){
315 for (Int_t i=0;i<3;i++)
316 for (Int_t j=0;j<3;j++){
317 (*fCovar)(i,j)=smatrix(i,j);
318 }
319 TMatrixD res = sums*smatrix;
320 fParams[0] = res(0,0);
321 fParams[1] = res(0,1);
322 fParams[2] = res(0,2);
323 conv++;
324 }
325 if (conv==0){
326 fConv = kFALSE; //non converged
327 return;
328 }
049179df 329 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
330 fConv = kFALSE; //non converged
331 return;
332 }
cabb917f 333 //
334 // XZ part
335 //
336 Double_t x0 = -fParams[1]/fParams[0];
049179df 337 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
cabb917f 338 Double_t sumXZ[9];
339
340 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
341 for (Int_t i=0;i<fN;i++){
049179df 342 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
343 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
cabb917f 344 Double_t weight = 1/fSz[i];
345 weight *=weight;
049179df 346 Double_t dphi = (phi-phi0)/rm1;
cabb917f 347 sumXZ[0] +=weight;
348 sumXZ[1] +=weight*dphi;
349 sumXZ[2] +=weight*dphi*dphi;
350 sumXZ[3] +=weight*fZ[i];
351 sumXZ[4] +=weight*dphi*fZ[i];
352
353 }
354
355 TMatrixDSym smatrixz(2);
356 TMatrixD sumsz(1,2);
357 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
049179df 358 smatrixz(1,0) = sumXZ[1];
cabb917f 359 smatrixz.Invert();
360 if (smatrixz.IsValid()){
361 sumsz(0,0) = sumXZ[3];
362 sumsz(0,1) = sumXZ[4];
363 TMatrixD res = sumsz*smatrixz;
364 fParams[3] = res(0,0);
365 fParams[4] = res(0,1);
366 for (Int_t i=0;i<2;i++)
367 for (Int_t j=0;j<2;j++){
368 (*fCovar)(i+3,j+3)=smatrixz(i,j);
369 }
370 conv++;
371 }
372
373 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
374 //
375 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
049179df 376 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
cabb917f 377 //
378 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
379 //
380 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
381 // final linear equation : a + u*b +t*c - v =0;
382 //
383 // fParam[0] = 1/y0
384 // fParam[1] = -x0/y0
385 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
386 if (conv>1) fConv =kTRUE;
387 else
388 fConv=kFALSE;
049179df 389 fChi2Y = CalcChi2Y();
390 fChi2Z = CalcChi2Z();
391 fChi2 = fChi2Y +fChi2Z;
cabb917f 392}
393
394
395
049179df 396Double_t AliRieman::GetYat(Double_t x) const {
397 //
398 // get y at given x position
399 //
cabb917f 400 if (!fConv) return 0.;
401 Double_t res2 = (x*fParams[0]+fParams[1]);
402 res2*=res2;
403 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
404 if (res2<0) return 0;
405 res2 = TMath::Sqrt(res2);
406 res2 = (1-res2)/fParams[0];
407 return res2;
408}
409
049179df 410Double_t AliRieman::GetDYat(Double_t x) const {
411 //
412 // get dy/dx at given x position
413 //
cabb917f 414 if (!fConv) return 0.;
415 Double_t x0 = -fParams[1]/fParams[0];
416 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
049179df 417 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
418 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
419 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
cabb917f 420 if (fParams[0]<0) res*=-1.;
421 return res;
422}
423
424
425
049179df 426Double_t AliRieman::GetZat(Double_t x) const {
427 //
428 // get z at given x position
429 //
cabb917f 430 if (!fConv) return 0.;
431 Double_t x0 = -fParams[1]/fParams[0];
432 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
049179df 433 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
434 Double_t phi = TMath::ASin((x-x0)*rm1);
435 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
cabb917f 436 Double_t dphi = (phi-phi0);
049179df 437 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
cabb917f 438 return res;
439}
440
049179df 441Double_t AliRieman::GetDZat(Double_t x) const {
442 //
443 // get dz/dx at given x postion
444 //
cabb917f 445 if (!fConv) return 0.;
446 Double_t x0 = -fParams[1]/fParams[0];
447 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
049179df 448 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
449 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
cabb917f 450 return res;
451}
452
453
049179df 454Double_t AliRieman::GetC() const{
455 //
456 // get curvature
457 //
cabb917f 458 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
459}
460
049179df 461Double_t AliRieman::CalcChi2Y() const{
462 //
463 // calculate chi2 for Y
464 //
465 Double_t sumchi2 = 0;
466 for (Int_t i=0;i<fN;i++){
467 Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
468 sumchi2+=chi2*chi2;
469 }
470 return sumchi2;
471}
472
473
474Double_t AliRieman::CalcChi2Z() const{
475 //
476 // calculate chi2 for Z
477 //
478 Double_t sumchi2 = 0;
479 for (Int_t i=0;i<fN;i++){
480 Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
481 sumchi2+=chi2*chi2;
482 }
483 return sumchi2;
484}
485
486Double_t AliRieman::CalcChi2() const{
487 //
488 // sum chi2 in both coord - supposing Y and Z independent
489 //
490 return CalcChi2Y()+CalcChi2Z();
491}
492
493AliRieman * AliRieman::MakeResiduals() const{
494 //
495 // create residual structure - ONLY for debug purposes
496 //
497 AliRieman *rieman = new AliRieman(fN);
498 for (Int_t i=0;i<fN;i++){
499 rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
500 }
501 return rieman;
502}
503
cabb917f 504