]>
Commit | Line | Data |
---|---|---|
1 | #include "TMatrixDSym.h" | |
2 | #include "TMatrixD.h" | |
3 | #include "AliTrackFitterStraight.h" | |
4 | ||
5 | ClassImp(AliTrackFitterStraight) | |
6 | ||
7 | AliTrackFitterStraight::AliTrackFitterStraight(): | |
8 | AliTrackFitter() | |
9 | { | |
10 | // | |
11 | // default constructor | |
12 | // | |
13 | fAlpha = 0.; | |
14 | for (Int_t i=0;i<5;i++) fSumXY[i] = 0; | |
15 | fSumYY = 0; | |
16 | for (Int_t i=0;i<5;i++) fSumXZ[i] = 0; | |
17 | fSumZZ = 0; | |
18 | fNUsed = 0; | |
19 | fConv = kFALSE; | |
20 | } | |
21 | ||
22 | ||
23 | AliTrackFitterStraight::AliTrackFitterStraight(AliTrackPointArray *array, Bool_t owner): | |
24 | AliTrackFitter(array,owner) | |
25 | { | |
26 | // | |
27 | // Constructor | |
28 | // | |
29 | fAlpha = 0.; | |
30 | for (Int_t i=0;i<5;i++) fSumXY[i] = 0; | |
31 | fSumYY = 0; | |
32 | for (Int_t i=0;i<5;i++) fSumXZ[i] = 0; | |
33 | fSumZZ = 0; | |
34 | fNUsed = 0; | |
35 | fConv = kFALSE; | |
36 | } | |
37 | ||
38 | AliTrackFitterStraight::AliTrackFitterStraight(const AliTrackFitterStraight &fitter): | |
39 | AliTrackFitter(fitter) | |
40 | { | |
41 | // | |
42 | // copy constructor | |
43 | // | |
44 | fAlpha = fitter.fAlpha; | |
45 | for (Int_t i=0;i<5;i++) fSumXY[i] = fitter.fSumXY[i]; | |
46 | fSumYY = fitter.fSumYY; | |
47 | for (Int_t i=0;i<5;i++) fSumXZ[i] = fitter.fSumXZ[i]; | |
48 | fSumZZ = fitter.fSumZZ; | |
49 | fNUsed = fitter.fNUsed; | |
50 | fConv = fitter.fConv; | |
51 | } | |
52 | ||
53 | //_____________________________________________________________________________ | |
54 | AliTrackFitterStraight &AliTrackFitterStraight::operator =(const AliTrackFitterStraight& fitter) | |
55 | { | |
56 | // assignment operator | |
57 | // | |
58 | if(this==&fitter) return *this; | |
59 | ((AliTrackFitter *)this)->operator=(fitter); | |
60 | ||
61 | fAlpha = fitter.fAlpha; | |
62 | for (Int_t i=0;i<5;i++) fSumXY[i] = fitter.fSumXY[i]; | |
63 | fSumYY = fitter.fSumYY; | |
64 | for (Int_t i=0;i<5;i++) fSumXZ[i] = fitter.fSumXZ[i]; | |
65 | fSumZZ = fitter.fSumZZ; | |
66 | fNUsed = fitter.fNUsed; | |
67 | fConv = fitter.fConv; | |
68 | ||
69 | return *this; | |
70 | } | |
71 | ||
72 | AliTrackFitterStraight::~AliTrackFitterStraight() | |
73 | { | |
74 | // destructor | |
75 | // | |
76 | } | |
77 | ||
78 | void AliTrackFitterStraight::Reset() | |
79 | { | |
80 | // Reset the track parameters and | |
81 | // sums | |
82 | AliTrackFitter::Reset(); | |
83 | fAlpha = 0.; | |
84 | for (Int_t i=0;i<5;i++) fSumXY[i] = 0; | |
85 | fSumYY = 0; | |
86 | for (Int_t i=0;i<5;i++) fSumXZ[i] = 0; | |
87 | fSumZZ = 0; | |
88 | fNUsed = 0; | |
89 | fConv =kFALSE; | |
90 | } | |
91 | ||
92 | Bool_t AliTrackFitterStraight::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, | |
93 | AliAlignObj::ELayerID layerRangeMin, | |
94 | AliAlignObj::ELayerID layerRangeMax) | |
95 | { | |
96 | // Fit the track points. The method takes as an input | |
97 | // the set of id's (volids) of the volumes in which | |
98 | // one wants to calculate the residuals. | |
99 | // The following parameters are used to define the | |
100 | // range of volumes to be used in the fitting | |
101 | // As a result two AliTrackPointArray's obects are filled. | |
102 | // The first one contains the space points with | |
103 | // volume id's from volids list. The second array of points represents | |
104 | // the track extrapolations corresponding to the space points | |
105 | // in the first array. The two arrays can be used to find | |
106 | // the residuals in the volids and consequently construct a | |
107 | // chi2 function to be minimized during the alignment | |
108 | // procedures. For the moment the track extrapolation is taken | |
109 | // at the space-point reference plane. The reference plane is | |
110 | // found using the covariance matrix of the point | |
111 | // (assuming sigma(x)=0 at the reference coordinate system. | |
112 | ||
113 | Reset(); | |
114 | ||
115 | Int_t npoints = fPoints->GetNPoints(); | |
116 | if (npoints < 2) return kFALSE; | |
117 | ||
118 | Bool_t isAlphaCalc = kFALSE; | |
119 | AliTrackPoint p,plocal; | |
120 | ||
121 | Int_t npVolId = 0; | |
122 | fNUsed = 0; | |
123 | Int_t *pindex = new Int_t[npoints]; | |
124 | for (Int_t ipoint = 0; ipoint < npoints; ipoint++) | |
125 | { | |
126 | fPoints->GetPoint(p,ipoint); | |
127 | UShort_t iVolId = p.GetVolumeID(); | |
128 | if (FindVolId(volIds,iVolId)) { | |
129 | pindex[npVolId] = ipoint; | |
130 | npVolId++; | |
131 | } | |
132 | if (volIdsFit != 0x0) { | |
133 | if (!FindVolId(volIdsFit,iVolId)) continue; | |
134 | } | |
135 | else { | |
136 | if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) || | |
137 | iVolId > AliAlignObj::LayerToVolUID(layerRangeMax, | |
138 | AliAlignObj::LayerSize(layerRangeMax))) continue; | |
139 | } | |
140 | if (!isAlphaCalc) { | |
141 | fAlpha = p.GetAngle(); | |
142 | isAlphaCalc = kTRUE; | |
143 | } | |
144 | plocal = p.Rotate(fAlpha); | |
145 | AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(), | |
146 | TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5])); | |
147 | fNUsed++; | |
148 | } | |
149 | ||
150 | if (fNUsed < 2) { | |
151 | delete [] pindex; | |
152 | return kFALSE; | |
153 | } | |
154 | ||
155 | Update(); | |
156 | ||
157 | if (!fConv) { | |
158 | delete [] pindex; | |
159 | return kFALSE; | |
160 | } | |
161 | ||
162 | if (fNUsed < fMinNPoints ) { | |
163 | delete [] pindex; | |
164 | return kFALSE; | |
165 | } | |
166 | ||
167 | fPVolId = new AliTrackPointArray(npVolId); | |
168 | fPTrack = new AliTrackPointArray(npVolId); | |
169 | AliTrackPoint p2; | |
170 | for (Int_t ipoint = 0; ipoint < npVolId; ipoint++) | |
171 | { | |
172 | Int_t index = pindex[ipoint]; | |
173 | fPoints->GetPoint(p,index); | |
174 | if (GetPCA(p,p2)) { | |
175 | Float_t xyz[3],xyz2[3]; | |
176 | p.GetXYZ(xyz); p2.GetXYZ(xyz2); | |
177 | // printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]); | |
178 | fPVolId->AddPoint(ipoint,&p); | |
179 | fPTrack->AddPoint(ipoint,&p2); | |
180 | } | |
181 | } | |
182 | ||
183 | delete [] pindex; | |
184 | ||
185 | return kTRUE; | |
186 | } | |
187 | ||
188 | void AliTrackFitterStraight::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz) | |
189 | { | |
190 | // Straight track fitter | |
191 | // The method add a point to the sums | |
192 | // used to extract track parameters | |
193 | ||
194 | // | |
195 | // XY part | |
196 | // | |
197 | Double_t weight = 1./(sy*sy); | |
198 | fSumXY[0] +=weight; | |
199 | fSumXY[1] +=x*weight; fSumXY[2] +=x*x*weight; | |
200 | fSumXY[3] +=y*weight; fSumXY[4] +=x*y*weight; | |
201 | fSumYY += y*y*weight; | |
202 | // | |
203 | // XZ part | |
204 | // | |
205 | weight = 1./(sz*sz); | |
206 | fSumXZ[0] +=weight; | |
207 | fSumXZ[1] +=x*weight; fSumXZ[2] +=x*x*weight; | |
208 | fSumXZ[3] +=z*weight; fSumXZ[4] +=x*z*weight; | |
209 | fSumZZ += z*z*weight; | |
210 | } | |
211 | ||
212 | void AliTrackFitterStraight::Update(){ | |
213 | // | |
214 | // Track fitter update | |
215 | // | |
216 | // | |
217 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
218 | fChi2 = 0; | |
219 | fNdf = 0; | |
220 | Int_t conv=0; | |
221 | // | |
222 | // XY part | |
223 | // | |
224 | TMatrixDSym smatrix(2); | |
225 | TMatrixD sums(1,2); | |
226 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[2]; | |
227 | smatrix(0,1) = fSumXY[1]; smatrix(1,0)=fSumXY[1]; | |
228 | sums(0,0) = fSumXY[3]; sums(0,1) =fSumXY[4]; | |
229 | smatrix.Invert(); | |
230 | if (smatrix.IsValid()){ | |
231 | for (Int_t i=0;i<2;i++) | |
232 | for (Int_t j=0;j<=i;j++){ | |
233 | (*fCov)(i,j)=smatrix(i,j); | |
234 | } | |
235 | TMatrixD res = sums*smatrix; | |
236 | fParams[0] = res(0,0); | |
237 | fParams[1] = res(0,1); | |
238 | TMatrixD tmp = res*sums.T(); | |
239 | fChi2 += fSumYY - tmp(0,0); | |
240 | fNdf += fNUsed - 2; | |
241 | conv++; | |
242 | } | |
243 | // | |
244 | // XZ part | |
245 | // | |
246 | TMatrixDSym smatrixz(2); | |
247 | TMatrixD sumsxz(1,2); | |
248 | smatrixz(0,0) = fSumXZ[0]; smatrixz(1,1) = fSumXZ[2]; | |
249 | smatrixz(0,1) = fSumXZ[1]; smatrixz(1,0) = fSumXZ[1]; | |
250 | sumsxz(0,0) = fSumXZ[3]; sumsxz(0,1) = fSumXZ[4]; | |
251 | smatrixz.Invert(); | |
252 | if (smatrixz.IsValid()){ | |
253 | TMatrixD res = sumsxz*smatrixz; | |
254 | fParams[2] = res(0,0); | |
255 | fParams[3] = res(0,1); | |
256 | fParams[4] = fParams[5] = 0; | |
257 | for (Int_t i=0;i<2;i++) | |
258 | for (Int_t j=0;j<=i;j++){ | |
259 | (*fCov)(i+2,j+2)=smatrixz(i,j); | |
260 | } | |
261 | TMatrixD tmp = res*sumsxz.T(); | |
262 | fChi2 += fSumZZ - tmp(0,0); | |
263 | fNdf += fNUsed - 2; | |
264 | conv++; | |
265 | } | |
266 | ||
267 | if (conv>1) | |
268 | fConv =kTRUE; | |
269 | else | |
270 | fConv=kFALSE; | |
271 | } | |
272 | ||
273 | Double_t AliTrackFitterStraight::GetYat(Double_t x) const { | |
274 | if (!fConv) return 0.; | |
275 | return (fParams[0]+x*fParams[1]); | |
276 | } | |
277 | ||
278 | Double_t AliTrackFitterStraight::GetDYat(Double_t x) const { | |
279 | if (!fConv) return 0.; | |
280 | return fParams[1]+0.*x; | |
281 | } | |
282 | ||
283 | ||
284 | ||
285 | Double_t AliTrackFitterStraight::GetZat(Double_t x) const { | |
286 | if (!fConv) return 0.; | |
287 | return (fParams[2]+x*fParams[3]); | |
288 | } | |
289 | ||
290 | Double_t AliTrackFitterStraight::GetDZat(Double_t x) const { | |
291 | if (!fConv) return 0.; | |
292 | return fParams[3]+0.*x; | |
293 | } | |
294 | ||
295 | Bool_t AliTrackFitterStraight::GetXYZat(Double_t r, Float_t *xyz) const { | |
296 | if (!fConv) return kFALSE; | |
297 | Double_t y = (fParams[0]+r*fParams[1]); | |
298 | Double_t z = (fParams[2]+r*fParams[3]); | |
299 | ||
300 | Double_t sin = TMath::Sin(fAlpha); | |
301 | Double_t cos = TMath::Cos(fAlpha); | |
302 | xyz[0] = r*cos - y*sin; | |
303 | xyz[1] = y*cos + r*sin; | |
304 | xyz[2] = z; | |
305 | ||
306 | return kTRUE; | |
307 | } | |
308 | ||
309 | Bool_t AliTrackFitterStraight::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const | |
310 | { | |
311 | // Get the closest to a given spacepoint track trajectory point | |
312 | // Look for details in the description of the Fit() method | |
313 | ||
314 | if (!fConv) return kFALSE; | |
315 | ||
316 | // First X and Y coordinates | |
317 | Double_t sin = TMath::Sin(fAlpha); | |
318 | Double_t cos = TMath::Cos(fAlpha); | |
319 | // Track parameters in the global coordinate system | |
320 | Double_t x0 = -fParams[0]*sin; | |
321 | Double_t y0 = fParams[0]*cos; | |
322 | if ((cos - fParams[1]*sin) == 0) return kFALSE; | |
323 | Double_t dydx = (fParams[1]*cos + sin)/(cos - fParams[1]*sin); | |
324 | ||
325 | // Define space-point refence plane | |
326 | Double_t alphap = p.GetAngle(); | |
327 | Double_t sinp = TMath::Sin(alphap); | |
328 | Double_t cosp = TMath::Cos(alphap); | |
329 | Double_t x = p.GetX()*cosp + p.GetY()*sinp; | |
330 | // Double_t y = p.GetY()*cosp - p.GetX()*sinp; | |
331 | Double_t x0p= x0*cosp + y0*sinp; | |
332 | Double_t y0p= y0*cosp - x0*sinp; | |
333 | if ((cos + dydx*sin) == 0) return kFALSE; | |
334 | Double_t dydxp = (dydx*cos - sin)/(cos + dydx*sin); | |
335 | Double_t yprime = y0p + dydxp*(x-x0p); | |
336 | ||
337 | // Back to the global coordinate system | |
338 | Double_t xsecond = x*cosp - yprime*sinp; | |
339 | Double_t ysecond = yprime*cosp + x*sinp; | |
340 | ||
341 | // Now Z coordinate and track angles | |
342 | Double_t x2 = xsecond*cos + ysecond*sin; | |
343 | Double_t zsecond = GetZat(x2); | |
344 | Double_t dydx2 = fParams[1]; | |
345 | Double_t dzdx = fParams[3]; | |
346 | ||
347 | // Fill the cov matrix of the track extrapolation point | |
348 | Double_t cov[6] = {0,0,0,0,0,0}; | |
349 | Double_t sigmax = 100*100.; | |
350 | cov[0] = sigmax; cov[1] = sigmax*dydx2; cov[2] = sigmax*dzdx; | |
351 | cov[3] = sigmax*dydx2*dydx2; cov[4] = sigmax*dydx2*dzdx; | |
352 | cov[5] = sigmax*dzdx*dzdx; | |
353 | ||
354 | Float_t newcov[6]; | |
355 | newcov[0] = cov[0]*cos*cos- | |
356 | 2*cov[1]*sin*cos+ | |
357 | cov[3]*sin*sin; | |
358 | newcov[1] = cov[1]*(cos*cos-sin*sin)- | |
359 | (cov[3]-cov[0])*sin*cos; | |
360 | newcov[2] = cov[2]*cos- | |
361 | cov[4]*sin; | |
362 | newcov[3] = cov[0]*sin*sin+ | |
363 | 2*cov[1]*sin*cos+ | |
364 | cov[3]*cos*cos; | |
365 | newcov[4] = cov[4]*cos+ | |
366 | cov[2]*sin; | |
367 | newcov[5] = cov[5]; | |
368 | ||
369 | p2.SetXYZ(xsecond,ysecond,zsecond,newcov); | |
370 | ||
371 | return kTRUE; | |
372 | } |