2 #include "TMatrixDSym.h"
11 AliRieman::AliRieman(){
13 // default constructor
15 for (Int_t i=0;i<6;i++) fParams[i] = 0;
16 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
17 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
31 AliRieman::AliRieman(Int_t capacity){
33 // default constructor
35 for (Int_t i=0;i<6;i++) fParams[i] = 0;
36 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
37 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
41 fX = new Float_t[fCapacity];
42 fY = new Float_t[fCapacity];
43 fZ = new Float_t[fCapacity];
44 fSy = new Float_t[fCapacity];
45 fSz = new Float_t[fCapacity];
46 fCovar = new TMatrixDSym(6);
50 AliRieman::AliRieman(const AliRieman &rieman):TObject(){
54 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
55 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
56 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
57 fCapacity = rieman.fN;
59 fCovar = new TMatrixDSym(*(rieman.fCovar));
64 fSy = new Float_t[fN];
65 fSz = new Float_t[fN];
66 for (Int_t i=0;i<rieman.fN;i++){
70 fSy[i] = rieman.fSy[i];
71 fSz[i] = rieman.fSz[i];
75 AliRieman::~AliRieman()
85 void AliRieman::Reset()
88 for (Int_t i=0;i<6;i++) fParams[i] = 0;
89 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
90 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
95 void AliRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
100 //------------------------------------------------------
103 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
105 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
107 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
109 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
111 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
113 // final linear equation : a + u*b +t*c - v =0;
117 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
119 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
121 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
123 if (fN==fCapacity-1) return; // out of capacity
124 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
128 Double_t t = x*x+y*y;
133 Double_t error = 2.*sy*t;
135 Double_t weight = 1./error;
137 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
138 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
139 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
145 fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
146 fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
155 void AliRieman::UpdatePol(){
161 for (Int_t i=0;i<6;i++)fParams[i]=0;
166 TMatrixDSym smatrix(3);
169 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
170 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
171 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
173 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
174 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
175 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
177 if (smatrix.IsValid()){
178 for (Int_t i=0;i<3;i++)
179 for (Int_t j=0;j<=i;j++){
180 (*fCovar)(i,j)=smatrix(i,j);
182 TMatrixD res = sums*smatrix;
183 fParams[0] = res(0,0);
184 fParams[1] = res(0,1);
185 fParams[2] = res(0,2);
191 TMatrixDSym smatrixz(3);
192 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
193 smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
194 smatrixz(2,2) = fSumXZ[4];
196 if (smatrixz.IsValid()){
197 sums(0,0) = fSumXZ[5];
198 sums(0,1) = fSumXZ[6];
199 sums(0,2) = fSumXZ[7];
200 TMatrixD res = sums*smatrixz;
201 fParams[3] = res(0,0);
202 fParams[4] = res(0,1);
203 fParams[5] = res(0,2);
204 for (Int_t i=0;i<3;i++)
205 for (Int_t j=0;j<=i;j++){
206 (*fCovar)(i+2,j+2)=smatrixz(i,j);
211 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
213 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
214 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
216 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
218 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
219 // final linear equation : a + u*b +t*c - v =0;
222 // fParam[1] = -x0/y0
223 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
224 if (conv>1) fConv =kTRUE;
229 void AliRieman::Update(){
235 for (Int_t i=0;i<6;i++)fParams[i]=0;
240 TMatrixDSym smatrix(3);
243 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
244 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
245 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
247 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
248 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
249 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
251 if (smatrix.IsValid()){
252 for (Int_t i=0;i<3;i++)
253 for (Int_t j=0;j<3;j++){
254 (*fCovar)(i,j)=smatrix(i,j);
256 TMatrixD res = sums*smatrix;
257 fParams[0] = res(0,0);
258 fParams[1] = res(0,1);
259 fParams[2] = res(0,2);
263 fConv = kFALSE; //non converged
269 Double_t x0 = -fParams[1]/fParams[0];
270 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
273 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
274 for (Int_t i=0;i<fN;i++){
275 Double_t phi = TMath::ASin((fX[i]-x0)*Rm1);
276 Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
277 Double_t weight = 1/fSz[i];
279 Double_t dphi = (phi-phi0)/Rm1;
281 sumXZ[1] +=weight*dphi;
282 sumXZ[2] +=weight*dphi*dphi;
283 sumXZ[3] +=weight*fZ[i];
284 sumXZ[4] +=weight*dphi*fZ[i];
288 TMatrixDSym smatrixz(2);
290 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
292 if (smatrixz.IsValid()){
293 sumsz(0,0) = sumXZ[3];
294 sumsz(0,1) = sumXZ[4];
295 TMatrixD res = sumsz*smatrixz;
296 fParams[3] = res(0,0);
297 fParams[4] = res(0,1);
298 for (Int_t i=0;i<2;i++)
299 for (Int_t j=0;j<2;j++){
300 (*fCovar)(i+3,j+3)=smatrixz(i,j);
305 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
307 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
308 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
310 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
312 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
313 // final linear equation : a + u*b +t*c - v =0;
316 // fParam[1] = -x0/y0
317 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
318 if (conv>1) fConv =kTRUE;
325 Double_t AliRieman::GetYat(Double_t x){
326 if (!fConv) return 0.;
327 Double_t res2 = (x*fParams[0]+fParams[1]);
329 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
330 if (res2<0) return 0;
331 res2 = TMath::Sqrt(res2);
332 res2 = (1-res2)/fParams[0];
336 Double_t AliRieman::GetDYat(Double_t x){
337 if (!fConv) return 0.;
338 Double_t x0 = -fParams[1]/fParams[0];
339 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
340 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
341 if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
342 Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
343 if (fParams[0]<0) res*=-1.;
349 Double_t AliRieman::GetZat(Double_t x){
350 if (!fConv) return 0.;
351 Double_t x0 = -fParams[1]/fParams[0];
352 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
353 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
354 Double_t phi = TMath::ASin((x-x0)*Rm1);
355 Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
356 Double_t dphi = (phi-phi0);
357 Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
361 Double_t AliRieman::GetDZat(Double_t x){
362 if (!fConv) return 0.;
363 Double_t x0 = -fParams[1]/fParams[0];
364 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
365 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
366 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
371 Double_t AliRieman::GetC(){
372 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);