]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliTrackFitterRieman.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 / AliTrackFitterRieman.cxx
CommitLineData
98937d93 1#include "TMatrixDSym.h"
2#include "TMatrixD.h"
3#include "AliTrackFitterRieman.h"
4
5ClassImp(AliTrackFitterRieman)
6
7AliTrackFitterRieman::AliTrackFitterRieman():
8 AliTrackFitter()
9{
10 //
11 // default constructor
12 //
13 fAlpha = 0.;
14 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
15 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
16 fConv = kFALSE;
17}
18
19
20AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
21 AliTrackFitter(array,owner)
22{
23 //
24 // Constructor
25 //
26 fAlpha = 0.;
27 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
28 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
29 fConv = kFALSE;
30}
31
32AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
33 AliTrackFitter(rieman)
34{
35 //
36 // copy constructor
37 //
38 fAlpha = rieman.fAlpha;
39 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
40 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
41 fConv = rieman.fConv;
42}
43
44//_____________________________________________________________________________
45AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
46{
47 // assignment operator
48 //
49 if(this==&rieman) return *this;
50 ((AliTrackFitter *)this)->operator=(rieman);
51
52 fAlpha = rieman.fAlpha;
53 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
54 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
55 fConv = rieman.fConv;
56
57 return *this;
58}
59
60AliTrackFitterRieman::~AliTrackFitterRieman()
61{
62 // destructor
63 //
64}
65
66void AliTrackFitterRieman::Reset()
67{
68 // Reset the track parameters and
69 // rieman sums
70 AliTrackFitter::Reset();
71 fAlpha = 0.;
72 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
73 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
74 fConv =kFALSE;
75}
76
77Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
78 AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
79 AliAlignObj::ELayerID layerRangeMin,
80 AliAlignObj::ELayerID layerRangeMax)
81{
82 // Fit the track points. The method takes as an input
83 // the id (volid) of the volume to be skipped from fitting.
84 // The following two parameters are used to define the
85 // range of volumes to be used in the fitting
86 // As a result two AliTrackPointArray's obects are filled.
87 // The first one contains the space points with
88 // volume id = volid. The second array of points represents
89 // the track extrapolation corresponding to the space points
90 // in the first array. The two arrays can be used to find
91 // the residuals in the volid and consequently construct a
92 // chi2 function to be minimized during the alignment
93 // procedures. For the moment the track extrapolation is taken
94 // as follows: in XY plane - at the CDA between track circle
95 // and the space point; in Z - the track extrapolation on the Z
96 // plane defined by the space point.
97
98 pVolId = pTrack = 0x0;
99 fConv = kFALSE;
100
101 Int_t npoints = fPoints->GetNPoints();
102 if (npoints < 3) return kFALSE;
103
104 AliTrackPoint p;
105 fPoints->GetPoint(p,0);
106 fAlpha = TMath::ATan2(p.GetY(),p.GetX());
107 Double_t sin = TMath::Sin(fAlpha);
108 Double_t cos = TMath::Cos(fAlpha);
109
110 Int_t npVolId = 0;
111 Int_t npused = 0;
112 Int_t *pindex = new Int_t[npoints];
113 for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
114 {
115 fPoints->GetPoint(p,ipoint);
116 UShort_t iVolId = p.GetVolumeID();
117 if (iVolId == volId) {
118 pindex[npVolId] = ipoint;
119 npVolId++;
120 }
121 if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
122 iVolId >= AliAlignObj::LayerToVolUID(layerRangeMax,0)) continue;
123 Float_t x = p.GetX()*cos + p.GetY()*sin;
124 Float_t y = p.GetY()*cos - p.GetX()*sin;
125 AddPoint(x,y,p.GetZ(),1,1);
126 npused++;
127 }
128
129 if (npused < 3) {
130 delete [] pindex;
131 return kFALSE;
132 }
133
134 Update();
135
136 if (!fConv) {
137 delete [] pindex;
138 return kFALSE;
139 }
140
141 if ((fParams[0] == 0) ||
142 ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
143 delete [] pindex;
144 return kFALSE;
145 }
146
147 pVolId = new AliTrackPointArray(npVolId);
148 pTrack = new AliTrackPointArray(npVolId);
149 AliTrackPoint p2;
150 for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
151 {
152 Int_t index = pindex[ipoint];
153 fPoints->GetPoint(p,index);
154 if (GetPCA(p,p2)) {
155 pVolId->AddPoint(ipoint,&p);
156 pTrack->AddPoint(ipoint,&p2);
157 }
158 }
159
160 delete [] pindex;
161
162 return kTRUE;
163}
164
165void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
166{
167 //
168 // Rieman update
169 //
170 //------------------------------------------------------
171 // XY direction
172 //
173 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
174 //
175 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
176 //
177 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
178 //
179 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
180 //
181 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
182 //
183 // final linear equation : a + u*b +t*c - v =0;
184 //
185 // Minimization :
186 //
187 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
188 //
189 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
190 //
191 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
192 //
193 //
194 // XY part
195 //
196 Double_t t = x*x+y*y;
197 if (t<2) return;
198 t = 1./t;
199 Double_t u = 2.*x*t;
200 Double_t v = 2.*y*t;
201 Double_t error = 2.*sy*t;
202 error *=error;
203 Double_t weight = 1./error;
204 fSumXY[0] +=weight;
205 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
206 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
207 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
208 //
209 // XZ part
210 //
211 weight = 1./sz;
212 fSumXZ[0] +=weight;
213 fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
214 fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
215}
216
217void AliTrackFitterRieman::Update(){
218 //
219 // Rieman update
220 //
221 //
222 for (Int_t i=0;i<6;i++)fParams[i]=0;
223 Int_t conv=0;
224 //
225 // XY part
226 //
227 TMatrixDSym smatrix(3);
228 TMatrixD sums(1,3);
229 //
230 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
231 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
232 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
233
234 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
235 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
236 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
237 smatrix.Invert();
238 if (smatrix.IsValid()){
239 for (Int_t i=0;i<3;i++)
240 for (Int_t j=0;j<=i;j++){
241 (*fCov)(i,j)=smatrix(i,j);
242 }
243 TMatrixD res = sums*smatrix;
244 fParams[0] = res(0,0);
245 fParams[1] = res(0,1);
246 fParams[2] = res(0,2);
247 conv++;
248 }
249 //
250 // XZ part
251 //
252 TMatrixDSym smatrixz(3);
253 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
254 smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
255 smatrixz(2,2) = fSumXZ[4];
256 smatrixz.Invert();
257 if (smatrixz.IsValid()){
258 sums(0,0) = fSumXZ[5];
259 sums(0,1) = fSumXZ[6];
260 sums(0,2) = fSumXZ[7];
261 TMatrixD res = sums*smatrixz;
262 fParams[3] = res(0,0);
263 fParams[4] = res(0,1);
264 fParams[5] = res(0,2);
265 for (Int_t i=0;i<3;i++)
266 for (Int_t j=0;j<=i;j++){
267 (*fCov)(i+2,j+2)=smatrixz(i,j);
268 }
269 conv++;
270 }
271
272 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
273 //
274 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
275 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
276 //
277 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
278 //
279 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
280 // final linear equation : a + u*b +t*c - v =0;
281 //
282 // fParam[0] = 1/y0
283 // fParam[1] = -x0/y0
284 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
285 if (conv>1) fConv =kTRUE;
286 else
287 fConv=kFALSE;
288}
289
290Double_t AliTrackFitterRieman::GetYat(Double_t x){
291 if (!fConv) return 0.;
292 Double_t res2 = (x*fParams[0]+fParams[1]);
293 res2*=res2;
294 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
295 if (res2<0) return 0;
296 res2 = TMath::Sqrt(res2);
297 res2 = (1-res2)/fParams[0];
298 return res2;
299}
300
301Double_t AliTrackFitterRieman::GetDYat(Double_t x){
302 if (!fConv) return 0.;
303 Double_t x0 = -fParams[1]/fParams[0];
304 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
305 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
306 if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
307 Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
308 if (fParams[0]<0) res*=-1.;
309 return res;
310}
311
312
313
314Double_t AliTrackFitterRieman::GetZat(Double_t x){
315 if (!fConv) return 0.;
316 Double_t x0 = -fParams[1]/fParams[0];
317 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
318 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
319 Double_t phi = TMath::ASin((x-x0)*Rm1);
320 Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
321 Double_t dphi = (phi-phi0);
322 Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
323 return res;
324}
325
326Double_t AliTrackFitterRieman::GetDZat(Double_t x){
327 if (!fConv) return 0.;
328 Double_t x0 = -fParams[1]/fParams[0];
329 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
330 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
331 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
332 return res;
333}
334
335
336Double_t AliTrackFitterRieman::GetC(){
337 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
338}
339
340Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){
341 if (!fConv) return kFALSE;
342 Double_t res2 = (r*fParams[0]+fParams[1]);
343 res2*=res2;
344 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
345 if (res2<0) return kFALSE;
346 res2 = TMath::Sqrt(res2);
347 res2 = (1-res2)/fParams[0];
348
349 if (!fConv) return kFALSE;
350 Double_t x0 = -fParams[1]/fParams[0];
351 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
352 Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
353 Double_t phi = TMath::ASin((r-x0)*Rm1);
354 Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
355 Double_t dphi = (phi-phi0);
356 Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
357
358 Double_t sin = TMath::Sin(fAlpha);
359 Double_t cos = TMath::Cos(fAlpha);
360 xyz[0] = r*cos - res2*sin;
361 xyz[1] = res2*cos + r*sin;
362 xyz[2] = res;
363
364 return kTRUE;
365}
366
367Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
368{
369 // Get the closest to a given spacepoint track trajectory point
370 // Look for details in the description of the Fit() method
371
372 if (!fConv) return kFALSE;
373
374 // First X and Y coordinates
375 Double_t sin = TMath::Sin(fAlpha);
376 Double_t cos = TMath::Cos(fAlpha);
377 // fParam[0] = 1/y0
378 // fParam[1] = -x0/y0
379 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
380 if (fParams[0] == 0) return kFALSE;
381 Double_t x0 = -fParams[1]/fParams[0]*cos - 1./fParams[0]*sin;
382 Double_t y0 = 1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
383 if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
384 Double_t R = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
385 fParams[0];
386
387 Double_t x = p.GetX();
388 Double_t y = p.GetY();
389 Double_t dR = TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
390 Double_t xprime = TMath::Abs(R)*(x-x0)/dR + x0;
391 Double_t yprime = TMath::Abs(R)*(y-y0)/dR + y0;
392
393 // Now Z coordinate
394 Double_t phi = TMath::ASin((x-x0)/R);
395 Double_t phi0 = TMath::ASin((0.-x0)/R);
396 Double_t dphi = (phi-phi0);
397 Double_t zprime = fParams[3]+fParams[4]*dphi*R;
398
399 p2.SetXYZ(xprime,yprime,zprime);
400
401 return kTRUE;
402}