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