Completely Updated (Mario Sitta)
[u/mrichter/AliRoot.git] / STEER / AliTrackPointArray.cxx
CommitLineData
98937d93 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//////////////////////////////////////////////////////////////////////////////
17// Class AliTrackPointArray //
18// This class contains the ESD track space-points which are used during //
19// the alignment procedures. Each space-point consist of 3 coordinates //
20// (and their errors) and the index of the sub-detector which contains //
21// the space-point. //
22// cvetan.cheshkov@cern.ch 3/11/2005 //
23//////////////////////////////////////////////////////////////////////////////
24
46ae650f 25#include <TMath.h>
4dcdc747 26#include <TMatrixD.h>
46ae650f 27#include <TMatrixDSym.h>
28
98937d93 29#include "AliTrackPointArray.h"
30
31ClassImp(AliTrackPointArray)
32
33//______________________________________________________________________________
fe12e09c 34AliTrackPointArray::AliTrackPointArray() :
35 TObject(),
bead9796 36 fSorted(kFALSE),
fe12e09c 37 fNPoints(0),
38 fX(0),
39 fY(0),
40 fZ(0),
41 fSize(0),
42 fCov(0),
43 fVolumeID(0)
98937d93 44{
98937d93 45}
46
47//______________________________________________________________________________
48AliTrackPointArray::AliTrackPointArray(Int_t npoints):
fe12e09c 49 TObject(),
bead9796 50 fSorted(kFALSE),
fe12e09c 51 fNPoints(npoints),
52 fX(new Float_t[npoints]),
53 fY(new Float_t[npoints]),
54 fZ(new Float_t[npoints]),
55 fSize(6*npoints),
56 fCov(new Float_t[fSize]),
57 fVolumeID(new UShort_t[npoints])
98937d93 58{
59 // Constructor
60 //
98937d93 61}
62
63//______________________________________________________________________________
64AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
fe12e09c 65 TObject(array),
bead9796 66 fSorted(array.fSorted),
fe12e09c 67 fNPoints(array.fNPoints),
68 fX(new Float_t[fNPoints]),
69 fY(new Float_t[fNPoints]),
70 fZ(new Float_t[fNPoints]),
71 fSize(array.fSize),
72 fCov(new Float_t[fSize]),
73 fVolumeID(new UShort_t[fNPoints])
98937d93 74{
75 // Copy constructor
76 //
98937d93 77 memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
78 memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
79 memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
80 memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
81 memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
82}
83
84//_____________________________________________________________________________
85AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& array)
86{
87 // assignment operator
88 //
89 if(this==&array) return *this;
90 ((TObject *)this)->operator=(array);
91
bead9796 92 fSorted = array.fSorted;
98937d93 93 fNPoints = array.fNPoints;
94 fSize = array.fSize;
6989bff3 95 delete [] fX;
98937d93 96 fX = new Float_t[fNPoints];
6989bff3 97 delete [] fY;
98937d93 98 fY = new Float_t[fNPoints];
6989bff3 99 delete [] fZ;
98937d93 100 fZ = new Float_t[fNPoints];
6989bff3 101 delete [] fVolumeID;
98937d93 102 fVolumeID = new UShort_t[fNPoints];
6989bff3 103 delete [] fCov;
98937d93 104 fCov = new Float_t[fSize];
105 memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
106 memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
107 memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
108 memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
109 memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
110
111 return *this;
112}
113
114//______________________________________________________________________________
115AliTrackPointArray::~AliTrackPointArray()
116{
117 // Destructor
118 //
119 delete [] fX;
120 delete [] fY;
121 delete [] fZ;
122 delete [] fVolumeID;
123 delete [] fCov;
124}
125
126
127//______________________________________________________________________________
128Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
129{
130 // Add a point to the array at position i
131 //
132 if (i >= fNPoints) return kFALSE;
133 fX[i] = p->GetX();
134 fY[i] = p->GetY();
135 fZ[i] = p->GetZ();
136 fVolumeID[i] = p->GetVolumeID();
137 memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
138 return kTRUE;
139}
140
ec9e17f9 141
98937d93 142//______________________________________________________________________________
143Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
144{
145 // Get the point at position i
146 //
147 if (i >= fNPoints) return kFALSE;
148 p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
149 p.SetVolumeID(fVolumeID[i]);
150 return kTRUE;
151}
152
153//______________________________________________________________________________
154Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
155{
156 // This method checks if the array
157 // has at least one hit in the detector
158 // volume defined by volid
159 Bool_t check = kFALSE;
160 for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
161 if (fVolumeID[ipoint] == volid) check = kTRUE;
162
163 return check;
164}
165
bead9796 166//______________________________________________________________________________
167void AliTrackPointArray::Sort(Bool_t down)
168{
169 // Sort the array by the values of Y-coordinate of the track points.
170 // The order is given by "down".
171 // Optimized more for maintenance rather than for speed.
172
173 if (fSorted) return;
174
175 Int_t *index=new Int_t[fNPoints];
176 AliTrackPointArray a(*this);
177 TMath::Sort(fNPoints,a.GetY(),index,down);
178
179 AliTrackPoint p;
180 for (Int_t i = 0; i < fNPoints; i++) {
181 a.GetPoint(p,index[i]);
182 AddPoint(i,&p);
183 }
184
185 delete[] index;
186 fSorted=kTRUE;
187}
188
98937d93 189ClassImp(AliTrackPoint)
190
191//______________________________________________________________________________
fe12e09c 192AliTrackPoint::AliTrackPoint() :
193 TObject(),
194 fX(0),
195 fY(0),
196 fZ(0),
197 fVolumeID(0)
98937d93 198{
199 // Default constructor
200 //
98937d93 201 memset(fCov,0,6*sizeof(Float_t));
202}
203
204
205//______________________________________________________________________________
fe12e09c 206AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid) :
207 TObject(),
208 fX(0),
209 fY(0),
210 fZ(0),
211 fVolumeID(0)
98937d93 212{
213 // Constructor
214 //
215 SetXYZ(x,y,z,cov);
216 SetVolumeID(volid);
217}
218
219//______________________________________________________________________________
fe12e09c 220AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid) :
221 TObject(),
222 fX(0),
223 fY(0),
224 fZ(0),
225 fVolumeID(0)
98937d93 226{
227 // Constructor
228 //
229 SetXYZ(xyz[0],xyz[1],xyz[2],cov);
230 SetVolumeID(volid);
231}
232
233//______________________________________________________________________________
234AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
fe12e09c 235 TObject(p),
236 fX(0),
237 fY(0),
238 fZ(0),
239 fVolumeID(0)
98937d93 240{
241 // Copy constructor
242 //
243 SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
244 SetVolumeID(p.fVolumeID);
245}
246
247//_____________________________________________________________________________
248AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
249{
250 // assignment operator
251 //
252 if(this==&p) return *this;
253 ((TObject *)this)->operator=(p);
254
255 SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
256 SetVolumeID(p.fVolumeID);
257
258 return *this;
259}
260
261//______________________________________________________________________________
262void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
263{
264 // Set XYZ coordinates and their cov matrix
265 //
266 fX = x;
267 fY = y;
268 fZ = z;
269 if (cov)
270 memcpy(fCov,cov,6*sizeof(Float_t));
271}
272
273//______________________________________________________________________________
274void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
275{
276 // Set XYZ coordinates and their cov matrix
277 //
278 SetXYZ(xyz[0],xyz[1],xyz[2],cov);
279}
280
281//______________________________________________________________________________
282void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
283{
284 xyz[0] = fX;
285 xyz[1] = fY;
286 xyz[2] = fZ;
287 if (cov)
288 memcpy(cov,fCov,6*sizeof(Float_t));
289}
46ae650f 290
291//______________________________________________________________________________
292Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
293{
294 // This method calculates the track to space-point residuals. The track
295 // interpolation is also stored as AliTrackPoint. Using the option
296 // 'weighted' one can calculate the residual either with or without
297 // taking into account the covariance matrix of the space-point and
298 // track interpolation. The second case the residual becomes a pull.
299
300 Float_t res = 0;
301
302 if (!weighted) {
303 Float_t xyz[3],xyzp[3];
304 GetXYZ(xyz);
305 p.GetXYZ(xyzp);
306 res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
307 (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
308 (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
309 }
310 else {
311 Float_t xyz[3],xyzp[3];
312 Float_t cov[6],covp[6];
313 GetXYZ(xyz,cov);
314 TMatrixDSym mcov(3);
315 mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
316 mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
317 mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
318 p.GetXYZ(xyzp,covp);
319 TMatrixDSym mcovp(3);
320 mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
321 mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
322 mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
323 TMatrixDSym msum = mcov + mcovp;
324 msum.Invert();
cc345ce3 325 // mcov.Print(); mcovp.Print(); msum.Print();
46ae650f 326 if (msum.IsValid()) {
327 for (Int_t i = 0; i < 3; i++)
328 for (Int_t j = 0; j < 3; j++)
329 res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
330 }
331 }
332
333 return res;
334}
335
4dcdc747 336//_____________________________________________________________________________
337Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
338{
339 //
340 // Get the intersection point between this point and
341 // the point "p" belongs to.
342 // The result is stored as a point 'out'
343 // return kFALSE in case of failure.
344 out.SetXYZ(0,0,0);
345
346 TMatrixD t(3,1);
347 t(0,0)=GetX();
348 t(1,0)=GetY();
349 t(2,0)=GetZ();
350
351 TMatrixDSym tC(3);
352 {
353 const Float_t *cv=GetCov();
354 tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
355 tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
356 tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
357 }
358
359 TMatrixD m(3,1);
360 m(0,0)=p.GetX();
361 m(1,0)=p.GetY();
362 m(2,0)=p.GetZ();
363
364 TMatrixDSym mC(3);
365 {
366 const Float_t *cv=p.GetCov();
367 mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
368 mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
369 mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
370 }
371
372 TMatrixDSym tmW(tC);
373 tmW+=mC;
374 tmW.Invert();
375 if (!tmW.IsValid()) return kFALSE;
376
377 TMatrixD mW(tC,TMatrixD::kMult,tmW);
378 TMatrixD tW(mC,TMatrixD::kMult,tmW);
379
380 TMatrixD mi(mW,TMatrixD::kMult,m);
381 TMatrixD ti(tW,TMatrixD::kMult,t);
382 ti+=mi;
383
384 TMatrixD iC(tC,TMatrixD::kMult,tmW);
385 iC*=mC;
386
387 out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
388 UShort_t id=p.GetVolumeID();
389 out.SetVolumeID(id);
390
391 return kTRUE;
392}
393
46ae650f 394//______________________________________________________________________________
395Float_t AliTrackPoint::GetAngle() const
396{
397 // The method uses the covariance matrix of
398 // the space-point in order to extract the
399 // orientation of the detector plane.
400 // The rotation in XY plane only is calculated.
401
8e52c1a8 402 Float_t phi= TMath::ATan2(TMath::Sqrt(fCov[0]),TMath::Sqrt(fCov[3]));
403 if (fCov[1] > 0) {
404 phi = TMath::Pi() - phi;
405 if ((fY-fX) < 0) phi += TMath::Pi();
406 }
46ae650f 407 else {
8e52c1a8 408 if ((fX+fY) < 0) phi += TMath::Pi();
409 }
410
411 return phi;
412
46ae650f 413}
414
415//_____________________________________________________________________________
416AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
417{
418 // Transform the space-point coordinates
419 // and covariance matrix from global to
420 // local (detector plane) coordinate system
421 // XY plane rotation only
422
423 static AliTrackPoint p;
424 p = *this;
425
426 Float_t xyz[3],cov[6];
427 GetXYZ(xyz,cov);
428
429 Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
430
431 Float_t newxyz[3],newcov[6];
432 newxyz[0] = cos*xyz[0] + sin*xyz[1];
433 newxyz[1] = cos*xyz[1] - sin*xyz[0];
434 newxyz[2] = xyz[2];
435
436 newcov[0] = cov[0]*cos*cos+
437 2*cov[1]*sin*cos+
438 cov[3]*sin*sin;
439 newcov[1] = cov[1]*(cos*cos-sin*sin)+
440 (cov[3]-cov[0])*sin*cos;
441 newcov[2] = cov[2]*cos+
442 cov[4]*sin;
443 newcov[3] = cov[0]*sin*sin-
444 2*cov[1]*sin*cos+
445 cov[3]*cos*cos;
446 newcov[4] = cov[4]*cos-
447 cov[2]*sin;
448 newcov[5] = cov[5];
449
450 p.SetXYZ(newxyz,newcov);
451 p.SetVolumeID(GetVolumeID());
452
453 return p;
454}
455
456//_____________________________________________________________________________
457AliTrackPoint& AliTrackPoint::MasterToLocal() const
458{
459 // Transform the space-point coordinates
460 // and the covariance matrix from the
461 // (master) to the local (tracking)
462 // coordinate system
463
464 Float_t alpha = GetAngle();
465 return Rotate(alpha);
466}
467
468//_____________________________________________________________________________
469void AliTrackPoint::Print(Option_t *) const
470{
471 // Print the space-point coordinates and
472 // covariance matrix
473
474 printf("VolumeID=%d\n", GetVolumeID());
475 printf("X = %12.6f Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
476 printf("Y = %12.6f Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
477 printf("Z = %12.6f Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
478
479}
ec9e17f9 480
481
482//________________________________
483void AliTrackPoint::SetAlignCovMatrix(const TMatrixDSym alignparmtrx){
484 // Add the uncertainty on the cluster position due to alignment
485 // (using the 6x6 AliAlignObj Cov. Matrix alignparmtrx) to the already
486 // present Cov. Matrix
487
488 TMatrixDSym cov(3);
489 TMatrixD coval(3,3);
490 TMatrixD jacob(3,6);
491 Float_t newcov[6];
492
493 cov(0,0)=fCov[0];
494 cov(1,0)=cov(0,1)=fCov[1];
495 cov(2,0)=cov(0,2)=fCov[2];
496 cov(1,1)=fCov[3];
497 cov(2,1)=cov(1,2)=fCov[4];
498 cov(2,2)=fCov[5];
499
500 jacob(0,0) = 1; jacob(1,0) = 0; jacob(2,0) = 0;
501 jacob(0,1) = 0; jacob(1,1) = 1; jacob(2,1) = 0;
502 jacob(0,2) = 0; jacob(1,2) = 0; jacob(2,2) = 1;
503 jacob(0,3) = 0; jacob(1,3) =-fZ; jacob(2,3) = fY;
504 jacob(0,4) = fZ; jacob(1,4) = 0; jacob(2,4) =-fX;
505 jacob(0,5) = fY; jacob(1,5) = fX; jacob(2,5) = 0;
506
507 TMatrixD jacobT=jacob.T();jacob.T();
508
509 coval=jacob*alignparmtrx*jacobT+cov;
510
511
512 newcov[0]=coval(0,0);
513 newcov[1]=coval(1,0);
514 newcov[2]=coval(2,0);
515 newcov[3]=coval(1,1);
516 newcov[4]=coval(2,1);
517 newcov[5]=coval(2,2);
518
519 SetXYZ(fX,fY,fZ,newcov);
520
521}
bead9796 522