]> 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
Alignment framework (C.Cheshkov). More information is available in http://agenda...
[u/mrichter/AliRoot.git] / STEER / AliRieman.cxx
CommitLineData
cabb917f 1#include "TTree.h"
2#include "TMatrixDSym.h"
3#include "TMatrixD.h"
4#include "AliRieman.h"
5#include "TRandom.h"
6
7ClassImp(AliRieman)
8
9
10
11AliRieman::AliRieman(){
12 //
13 // default constructor
14 //
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;
18
19 fCapacity = 0;
20 fN =0;
21 fX = 0;
22 fY = 0;
23 fZ = 0;
24 fSy = 0;
25 fSz = 0;
26 fCovar = 0;
27 fConv = kFALSE;
28}
29
30
31AliRieman::AliRieman(Int_t capacity){
32 //
33 // default constructor
34 //
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;
38
39 fCapacity = capacity;
40 fN =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);
47 fConv = kFALSE;
48}
49
50AliRieman::AliRieman(const AliRieman &rieman):TObject(){
51 //
52 // copy constructor
53 //
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;
58 fN =rieman.fN;
59 fCovar = new TMatrixDSym(*(rieman.fCovar));
60 fConv = rieman.fConv;
61 fX = new Float_t[fN];
62 fY = new Float_t[fN];
63 fZ = new Float_t[fN];
64 fSy = new Float_t[fN];
65 fSz = new Float_t[fN];
66 for (Int_t i=0;i<rieman.fN;i++){
67 fX[i] = rieman.fX[i];
68 fY[i] = rieman.fY[i];
69 fZ[i] = rieman.fZ[i];
70 fSy[i] = rieman.fSy[i];
71 fSz[i] = rieman.fSz[i];
72 }
73}
74
75AliRieman::~AliRieman()
76{
77 delete[]fX;
78 delete[]fY;
79 delete[]fZ;
80 delete[]fSy;
81 delete[]fSz;
82 delete fCovar;
83}
84
85void AliRieman::Reset()
86{
87 fN=0;
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;
91 fConv =kFALSE;
92}
93
94
95void AliRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
96{
97 //
98 // Rieman update
99 //
100 //------------------------------------------------------
101 // XY direction
102 //
103 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
104 //
105 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
106 //
107 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
108 //
109 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
110 //
111 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
112 //
113 // final linear equation : a + u*b +t*c - v =0;
114 //
115 // Minimization :
116 //
117 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
118 //
119 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
120 //
121 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
122 //
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;
125 //
126 // XY part
127 //
128 Double_t t = x*x+y*y;
129 if (t<2) return;
130 t = 1./t;
131 Double_t u = 2.*x*t;
132 Double_t v = 2.*y*t;
133 Double_t error = 2.*sy*t;
134 error *=error;
135 Double_t weight = 1./error;
136 fSumXY[0] +=weight;
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;
140 //
141 // XZ part
142 //
143 weight = 1./sz;
144 fSumXZ[0] +=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;
147 //
148 fN++;
149}
150
151
152
153
154
155void AliRieman::UpdatePol(){
156 //
157 // Rieman update
158 //
159 //
160 if (fN==0) return;
161 for (Int_t i=0;i<6;i++)fParams[i]=0;
162 Int_t conv=0;
163 //
164 // XY part
165 //
166 TMatrixDSym smatrix(3);
167 TMatrixD sums(1,3);
168 //
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;
172
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];
176 smatrix.Invert();
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);
181 }
182 TMatrixD res = sums*smatrix;
183 fParams[0] = res(0,0);
184 fParams[1] = res(0,1);
185 fParams[2] = res(0,2);
186 conv++;
187 }
188 //
189 // XZ part
190 //
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];
195 smatrixz.Invert();
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);
207 }
208 conv++;
209 }
210
211 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
212 //
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
215 //
216 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
217 //
218 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
219 // final linear equation : a + u*b +t*c - v =0;
220 //
221 // fParam[0] = 1/y0
222 // fParam[1] = -x0/y0
223 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
224 if (conv>1) fConv =kTRUE;
225 else
226 fConv=kFALSE;
227}
228
229void AliRieman::Update(){
230 //
231 // Rieman update
232 //
233 //
234 if (fN==0) return;
235 for (Int_t i=0;i<6;i++)fParams[i]=0;
236 Int_t conv=0;
237 //
238 // XY part
239 //
240 TMatrixDSym smatrix(3);
241 TMatrixD sums(1,3);
242 //
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;
246
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];
250 smatrix.Invert();
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);
255 }
256 TMatrixD res = sums*smatrix;
257 fParams[0] = res(0,0);
258 fParams[1] = res(0,1);
259 fParams[2] = res(0,2);
260 conv++;
261 }
262 if (conv==0){
263 fConv = kFALSE; //non converged
264 return;
265 }
266 //
267 // XZ part
268 //
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);
271 Double_t sumXZ[9];
272
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];
278 weight *=weight;
279 Double_t dphi = (phi-phi0)/Rm1;
280 sumXZ[0] +=weight;
281 sumXZ[1] +=weight*dphi;
282 sumXZ[2] +=weight*dphi*dphi;
283 sumXZ[3] +=weight*fZ[i];
284 sumXZ[4] +=weight*dphi*fZ[i];
285
286 }
287
288 TMatrixDSym smatrixz(2);
289 TMatrixD sumsz(1,2);
290 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
291 smatrixz.Invert();
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);
301 }
302 conv++;
303 }
304
305 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
306 //
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
309 //
310 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
311 //
312 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
313 // final linear equation : a + u*b +t*c - v =0;
314 //
315 // fParam[0] = 1/y0
316 // fParam[1] = -x0/y0
317 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
318 if (conv>1) fConv =kTRUE;
319 else
320 fConv=kFALSE;
321}
322
323
324
325Double_t AliRieman::GetYat(Double_t x){
326 if (!fConv) return 0.;
327 Double_t res2 = (x*fParams[0]+fParams[1]);
328 res2*=res2;
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];
333 return res2;
334}
335
336Double_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.;
344 return res;
345}
346
347
348
349Double_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;
358 return res;
359}
360
361Double_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));
367 return res;
368}
369
370
371Double_t AliRieman::GetC(){
372 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
373}
374
375