1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //////////////////////////////////////////////////////////////////////////////
27 #include <TMatrixDSym.h>
28 #include <TGeoMatrix.h>
29 #include <TMatrixDSymEigen.h>
31 #include "AliTrackPointArray.h"
33 ClassImp(AliTrackPointArray)
35 //______________________________________________________________________________
36 AliTrackPointArray::AliTrackPointArray() :
54 //______________________________________________________________________________
55 AliTrackPointArray::AliTrackPointArray(Int_t npoints):
59 fX(new Float_t[npoints]),
60 fY(new Float_t[npoints]),
61 fZ(new Float_t[npoints]),
62 fCharge(new Float_t[npoints]),
63 fDriftTime(new Float_t[npoints]),
64 fChargeRatio(new Float_t[npoints]),
65 fClusterType(new Int_t[npoints]),
66 fIsExtra(new Bool_t[npoints]),
68 fCov(new Float_t[fSize]),
69 fVolumeID(new UShort_t[npoints])
73 for (Int_t ip=0; ip<npoints;ip++){
83 for (Int_t icov=0;icov<6; icov++)
88 //______________________________________________________________________________
89 AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
91 fSorted(array.fSorted),
92 fNPoints(array.fNPoints),
93 fX(new Float_t[fNPoints]),
94 fY(new Float_t[fNPoints]),
95 fZ(new Float_t[fNPoints]),
96 fCharge(new Float_t[fNPoints]),
97 fDriftTime(new Float_t[fNPoints]),
98 fChargeRatio(new Float_t[fNPoints]),
99 fClusterType(new Int_t[fNPoints]),
100 fIsExtra(new Bool_t[fNPoints]),
102 fCov(new Float_t[fSize]),
103 fVolumeID(new UShort_t[fNPoints])
107 memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
108 memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
109 memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
111 memcpy(fCharge,array.fCharge,fNPoints*sizeof(Float_t));
113 memset(fCharge, 0, fNPoints*sizeof(Float_t));
115 if (array.fDriftTime) {
116 memcpy(fDriftTime,array.fDriftTime,fNPoints*sizeof(Float_t));
118 memset(fDriftTime, 0, fNPoints*sizeof(Float_t));
120 if (array.fChargeRatio) {
121 memcpy(fChargeRatio,array.fChargeRatio,fNPoints*sizeof(Float_t));
123 memset(fChargeRatio, 0, fNPoints*sizeof(Float_t));
125 if (array.fClusterType) {
126 memcpy(fClusterType,array.fClusterType,fNPoints*sizeof(Int_t));
128 memset(fClusterType, 0, fNPoints*sizeof(Int_t));
130 if (array.fIsExtra) {
131 memcpy(fIsExtra,array.fIsExtra,fNPoints*sizeof(Bool_t));
133 memset(fIsExtra, 0, fNPoints*sizeof(Bool_t));
135 memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
136 memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
139 //_____________________________________________________________________________
140 AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& array)
142 // assignment operator
144 if(this==&array) return *this;
145 ((TObject *)this)->operator=(array);
147 fSorted = array.fSorted;
148 fNPoints = array.fNPoints;
151 fX = new Float_t[fNPoints];
153 fY = new Float_t[fNPoints];
155 fZ = new Float_t[fNPoints];
157 fCharge = new Float_t[fNPoints];
158 delete [] fDriftTime;
159 fDriftTime = new Float_t[fNPoints];
160 delete [] fChargeRatio;
161 fChargeRatio = new Float_t[fNPoints];
162 delete [] fClusterType;
163 fClusterType = new Int_t[fNPoints];
165 fIsExtra = new Bool_t[fNPoints];
167 fVolumeID = new UShort_t[fNPoints];
169 fCov = new Float_t[fSize];
170 memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
171 memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
172 memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
173 memcpy(fCharge,array.fCharge,fNPoints*sizeof(Float_t));
174 memcpy(fDriftTime,array.fDriftTime,fNPoints*sizeof(Float_t));
175 memcpy(fChargeRatio,array.fChargeRatio,fNPoints*sizeof(Float_t));
176 memcpy(fClusterType,array.fClusterType,fNPoints*sizeof(Int_t));
177 memcpy(fIsExtra,array.fIsExtra,fNPoints*sizeof(Bool_t));
178 memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
179 memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
184 //______________________________________________________________________________
185 AliTrackPointArray::~AliTrackPointArray()
193 delete [] fDriftTime;
194 delete [] fChargeRatio;
195 delete [] fClusterType;
202 //______________________________________________________________________________
203 Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
205 // Add a point to the array at position i
207 if (i >= fNPoints) return kFALSE;
211 fCharge[i] = p->GetCharge();
212 fDriftTime[i] = p->GetDriftTime();
213 fChargeRatio[i]=p->GetChargeRatio();
214 fClusterType[i]=p->GetClusterType();
215 fIsExtra[i] = p->IsExtra();
216 fVolumeID[i] = p->GetVolumeID();
217 memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
222 //______________________________________________________________________________
223 Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
225 // Get the point at position i
227 if (i >= fNPoints) return kFALSE;
228 p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
229 p.SetVolumeID(fVolumeID[i]);
230 p.SetCharge(fCharge ? fCharge[i] : 0);
231 p.SetDriftTime(fDriftTime ? fDriftTime[i] : 0);
232 p.SetChargeRatio(fChargeRatio ? fChargeRatio[i] : 0);
233 p.SetClusterType(fClusterType ? fClusterType[i] : 0);
234 p.SetExtra(fIsExtra ? fIsExtra[i] : kFALSE);
238 //______________________________________________________________________________
239 Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
241 // This method checks if the array
242 // has at least one hit in the detector
243 // volume defined by volid
244 Bool_t check = kFALSE;
245 for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
246 if (fVolumeID[ipoint] == volid) check = kTRUE;
251 //______________________________________________________________________________
252 void AliTrackPointArray::Sort(Bool_t down)
254 // Sort the array by the values of Y-coordinate of the track points.
255 // The order is given by "down".
256 // Optimized more for maintenance rather than for speed.
260 Int_t *index=new Int_t[fNPoints];
261 AliTrackPointArray a(*this);
262 TMath::Sort(fNPoints,a.GetY(),index,down);
265 for (Int_t i = 0; i < fNPoints; i++) {
266 a.GetPoint(p,index[i]);
274 //_____________________________________________________________________________
275 void AliTrackPointArray::Print(Option_t *) const
277 // Print the space-point coordinates and modules info
278 for (int i=0;i<fNPoints;i++) {
279 printf("#%3d VID %5d XYZ:%+9.3f/%+9.3f/%+9.3f |q: %+7.2f |DT: %+8.1f| ChR: %+.2e |Cl: %d %s\n",
280 i,fVolumeID[i],fX[i],fY[i],fZ[i],fCharge[i],fDriftTime[i],
281 fChargeRatio[i],fClusterType[i],fIsExtra[i] ? "|E":"");
286 ClassImp(AliTrackPoint)
288 //______________________________________________________________________________
289 AliTrackPoint::AliTrackPoint() :
301 // Default constructor
303 memset(fCov,0,6*sizeof(Float_t));
307 //______________________________________________________________________________
308 AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid, Float_t charge, Float_t drifttime,Float_t chargeratio, Int_t clutyp) :
324 SetDriftTime(drifttime);
325 SetChargeRatio(chargeratio);
326 SetClusterType(clutyp);
330 //______________________________________________________________________________
331 AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid, Float_t charge, Float_t drifttime,Float_t chargeratio, Int_t clutyp) :
345 SetXYZ(xyz[0],xyz[1],xyz[2],cov);
347 SetDriftTime(drifttime);
348 SetChargeRatio(chargeratio);
350 SetClusterType(clutyp);
353 //______________________________________________________________________________
354 AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
368 SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
369 SetCharge(p.fCharge);
370 SetDriftTime(p.fDriftTime);
371 SetChargeRatio(p.fChargeRatio);
372 SetClusterType(p.fClusterType);
373 SetExtra(p.fIsExtra);
374 SetVolumeID(p.fVolumeID);
377 //_____________________________________________________________________________
378 AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
380 // assignment operator
382 if(this==&p) return *this;
383 ((TObject *)this)->operator=(p);
385 SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
386 SetCharge(p.fCharge);
387 SetDriftTime(p.fDriftTime);
388 SetChargeRatio(p.fChargeRatio);
389 SetClusterType(p.fClusterType);
390 SetExtra(p.fIsExtra);
391 SetVolumeID(p.fVolumeID);
396 //______________________________________________________________________________
397 void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
399 // Set XYZ coordinates and their cov matrix
405 memcpy(fCov,cov,6*sizeof(Float_t));
408 //______________________________________________________________________________
409 void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
411 // Set XYZ coordinates and their cov matrix
413 SetXYZ(xyz[0],xyz[1],xyz[2],cov);
416 //______________________________________________________________________________
417 void AliTrackPoint::SetCov(const Float_t *cov)
419 // Set XYZ cov matrix
422 memcpy(fCov,cov,6*sizeof(Float_t));
425 //______________________________________________________________________________
426 void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
432 memcpy(cov,fCov,6*sizeof(Float_t));
435 //______________________________________________________________________________
436 Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
438 // This method calculates the track to space-point residuals. The track
439 // interpolation is also stored as AliTrackPoint. Using the option
440 // 'weighted' one can calculate the residual either with or without
441 // taking into account the covariance matrix of the space-point and
442 // track interpolation. The second case the residual becomes a pull.
447 Float_t xyz[3],xyzp[3];
450 res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
451 (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
452 (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
455 Float_t xyz[3],xyzp[3];
456 Float_t cov[6],covp[6];
459 mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
460 mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
461 mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
463 TMatrixDSym mcovp(3);
464 mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
465 mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
466 mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
467 TMatrixDSym msum = mcov + mcovp;
469 // mcov.Print(); mcovp.Print(); msum.Print();
470 if (msum.IsValid()) {
471 for (Int_t i = 0; i < 3; i++)
472 for (Int_t j = 0; j < 3; j++)
473 res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
480 //_____________________________________________________________________________
481 Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
484 // Get the intersection point between this point and
485 // the point "p" belongs to.
486 // The result is stored as a point 'out'
487 // return kFALSE in case of failure.
497 const Float_t *cv=GetCov();
498 tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
499 tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
500 tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
510 const Float_t *cv=p.GetCov();
511 mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
512 mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
513 mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
519 if (!tmW.IsValid()) return kFALSE;
521 TMatrixD mW(tC,TMatrixD::kMult,tmW);
522 TMatrixD tW(mC,TMatrixD::kMult,tmW);
524 TMatrixD mi(mW,TMatrixD::kMult,m);
525 TMatrixD ti(tW,TMatrixD::kMult,t);
528 TMatrixD iC(tC,TMatrixD::kMult,tmW);
531 out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
532 UShort_t id=p.GetVolumeID();
538 //______________________________________________________________________________
539 Float_t AliTrackPoint::GetAngle() const
541 // The method uses the covariance matrix of
542 // the space-point in order to extract the
543 // orientation of the detector plane.
544 // The rotation in XY plane only is calculated.
546 Float_t phi= TMath::ATan2(TMath::Sqrt(fCov[0]),TMath::Sqrt(fCov[3]));
548 phi = TMath::Pi() - phi;
549 if ((fY-fX) < 0) phi += TMath::Pi();
552 if ((fX+fY) < 0) phi += TMath::Pi();
559 //______________________________________________________________________________
560 Bool_t AliTrackPoint::GetRotMatrix(TGeoRotation& rot) const
562 // Returns the orientation of the
563 // sensitive layer (using cluster
564 // covariance matrix).
565 // Assumes that cluster has errors only in the layer's plane.
566 // Return value is kTRUE in case of success.
570 const Float_t *cov=GetCov();
571 mcov(0,0)=cov[0]; mcov(0,1)=cov[1]; mcov(0,2)=cov[2];
572 mcov(1,0)=cov[1]; mcov(1,1)=cov[3]; mcov(1,2)=cov[4];
573 mcov(2,0)=cov[2]; mcov(2,1)=cov[4]; mcov(2,2)=cov[5];
576 TMatrixDSymEigen eigen(mcov);
577 TMatrixD eigenMatrix = eigen.GetEigenVectors();
579 rot.SetMatrix(eigenMatrix.GetMatrixArray());
585 //_____________________________________________________________________________
586 AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
588 // Transform the space-point coordinates
589 // and covariance matrix from global to
590 // local (detector plane) coordinate system
591 // XY plane rotation only
593 static AliTrackPoint p;
596 Float_t xyz[3],cov[6];
599 Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
601 Float_t newxyz[3],newcov[6];
602 newxyz[0] = cos*xyz[0] + sin*xyz[1];
603 newxyz[1] = cos*xyz[1] - sin*xyz[0];
606 newcov[0] = cov[0]*cos*cos+
609 newcov[1] = cov[1]*(cos*cos-sin*sin)+
610 (cov[3]-cov[0])*sin*cos;
611 newcov[2] = cov[2]*cos+
613 newcov[3] = cov[0]*sin*sin-
616 newcov[4] = cov[4]*cos-
620 p.SetXYZ(newxyz,newcov);
621 p.SetVolumeID(GetVolumeID());
626 //_____________________________________________________________________________
627 AliTrackPoint& AliTrackPoint::MasterToLocal() const
629 // Transform the space-point coordinates
630 // and the covariance matrix from the
631 // (master) to the local (tracking)
634 Float_t alpha = GetAngle();
635 return Rotate(alpha);
638 //_____________________________________________________________________________
639 void AliTrackPoint::Print(Option_t *) const
641 // Print the space-point coordinates and
644 printf("VolumeID=%d\n", GetVolumeID());
645 printf("X = %12.6f Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
646 printf("Y = %12.6f Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
647 printf("Z = %12.6f Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
648 printf("Charge = %f\n", fCharge);
649 printf("Drift Time = %f\n", fDriftTime);
650 printf("Charge Ratio = %f\n", fChargeRatio);
651 printf("Cluster Type = %d\n", fClusterType);
652 if(fIsExtra) printf("This is an extra point\n");
657 //________________________________
658 void AliTrackPoint::SetAlignCovMatrix(const TMatrixDSym& alignparmtrx){
659 // Add the uncertainty on the cluster position due to alignment
660 // (using the 6x6 AliAlignObj Cov. Matrix alignparmtrx) to the already
661 // present Cov. Matrix
669 cov(1,0)=cov(0,1)=fCov[1];
670 cov(2,0)=cov(0,2)=fCov[2];
672 cov(2,1)=cov(1,2)=fCov[4];
675 jacob(0,0) = 1; jacob(1,0) = 0; jacob(2,0) = 0;
676 jacob(0,1) = 0; jacob(1,1) = 1; jacob(2,1) = 0;
677 jacob(0,2) = 0; jacob(1,2) = 0; jacob(2,2) = 1;
678 jacob(0,3) = 0; jacob(1,3) =-fZ; jacob(2,3) = fY;
679 jacob(0,4) = fZ; jacob(1,4) = 0; jacob(2,4) =-fX;
680 jacob(0,5) = -fY; jacob(1,5) = fX; jacob(2,5) = 0;
682 TMatrixD jacobT=jacob.T();jacob.T();
684 coval=jacob*alignparmtrx*jacobT+cov;
687 newcov[0]=coval(0,0);
688 newcov[1]=coval(1,0);
689 newcov[2]=coval(2,0);
690 newcov[3]=coval(1,1);
691 newcov[4]=coval(2,1);
692 newcov[5]=coval(2,2);
694 SetXYZ(fX,fY,fZ,newcov);