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