Fixes for coverity.
[u/mrichter/AliRoot.git] / STEER / STEER / AliTrackFitterStraight.cxx
CommitLineData
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
24ClassImp(AliTrackFitterStraight)
25
26AliTrackFitterStraight::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
42AliTrackFitterStraight::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
57AliTrackFitterStraight::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//_____________________________________________________________________________
73AliTrackFitterStraight &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
91AliTrackFitterStraight::~AliTrackFitterStraight()
92{
93 // destructor
94 //
95}
96
97void 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
111Bool_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
207void 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 231Bool_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
294Double_t AliTrackFitterStraight::GetYat(Double_t x) const {
295 if (!fConv) return 0.;
296 return (fParams[0]+x*fParams[1]);
297}
298
299Double_t AliTrackFitterStraight::GetDYat(Double_t x) const {
300 if (!fConv) return 0.;
301 return fParams[1]+0.*x;
302}
303
304
305
306Double_t AliTrackFitterStraight::GetZat(Double_t x) const {
307 if (!fConv) return 0.;
308 return (fParams[2]+x*fParams[3]);
309}
310
311Double_t AliTrackFitterStraight::GetDZat(Double_t x) const {
312 if (!fConv) return 0.;
313 return fParams[3]+0.*x;
314}
315
316Bool_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
330Bool_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}