]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliTrackPointArray.cxx
protect QA agains non TH objects
[u/mrichter/AliRoot.git] / STEER / AliTrackPointArray.cxx
index 92b2b15db0ff3628f9bc1a09c3807ec98ef92356..4256745386b449792bef4c5f6c29851b6100ced5 100644 (file)
@@ -23,6 +23,7 @@
 //////////////////////////////////////////////////////////////////////////////
 
 #include <TMath.h>
+#include <TMatrixD.h>
 #include <TMatrixDSym.h>
 
 #include "AliTrackPointArray.h"
@@ -32,6 +33,7 @@ ClassImp(AliTrackPointArray)
 //______________________________________________________________________________
 AliTrackPointArray::AliTrackPointArray() :
   TObject(),
+  fSorted(kFALSE),
   fNPoints(0),
   fX(0),
   fY(0),
@@ -45,6 +47,7 @@ AliTrackPointArray::AliTrackPointArray() :
 //______________________________________________________________________________
 AliTrackPointArray::AliTrackPointArray(Int_t npoints):
   TObject(),
+  fSorted(kFALSE),
   fNPoints(npoints),
   fX(new Float_t[npoints]),
   fY(new Float_t[npoints]),
@@ -55,11 +58,20 @@ AliTrackPointArray::AliTrackPointArray(Int_t npoints):
 {
   // Constructor
   //
+  for (Int_t ip=0; ip<npoints;ip++){
+    fX[ip]=0;
+    fY[ip]=0;
+    fZ[ip]=0;
+    fVolumeID[ip]=0;
+    for (Int_t icov=0;icov<6; icov++)
+      fCov[6*ip+icov]=0;
+  }
 }
 
 //______________________________________________________________________________
 AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
   TObject(array),
+  fSorted(array.fSorted),
   fNPoints(array.fNPoints),
   fX(new Float_t[fNPoints]),
   fY(new Float_t[fNPoints]),
@@ -85,12 +97,18 @@ AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& arr
   if(this==&array) return *this;
   ((TObject *)this)->operator=(array);
 
+  fSorted = array.fSorted;
   fNPoints = array.fNPoints;
   fSize = array.fSize;
+  delete [] fX;
   fX = new Float_t[fNPoints];
+  delete [] fY;
   fY = new Float_t[fNPoints];
+  delete [] fZ;
   fZ = new Float_t[fNPoints];
+  delete [] fVolumeID;
   fVolumeID = new UShort_t[fNPoints];
+  delete [] fCov;
   fCov = new Float_t[fSize];
   memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
   memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
@@ -128,6 +146,7 @@ Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
   return kTRUE;
 }
 
+
 //______________________________________________________________________________
 Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
 {
@@ -152,6 +171,29 @@ Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
   return check;
 }
 
+//______________________________________________________________________________
+void AliTrackPointArray::Sort(Bool_t down)
+{
+  // Sort the array by the values of Y-coordinate of the track points.
+  // The order is given by "down".
+  // Optimized more for maintenance rather than for speed.
+  if (fSorted) return;
+
+  Int_t *index=new Int_t[fNPoints];
+  AliTrackPointArray a(*this);
+  TMath::Sort(fNPoints,a.GetY(),index,down);
+  AliTrackPoint p;
+  for (Int_t i = 0; i < fNPoints; i++) {
+    a.GetPoint(p,index[i]);
+    AddPoint(i,&p);
+  }
+
+  delete[] index;
+  fSorted=kTRUE;
+}
+
 ClassImp(AliTrackPoint)
 
 //______________________________________________________________________________
@@ -299,6 +341,64 @@ Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) cons
   return res;
 }
 
+//_____________________________________________________________________________
+Bool_t AliTrackPoint::GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const
+{
+  //
+  // Get the intersection point between this point and
+  // the point "p" belongs to.
+  // The result is stored as a point 'out'
+  // return kFALSE in case of failure.
+  out.SetXYZ(0,0,0);
+
+  TMatrixD t(3,1);
+  t(0,0)=GetX();
+  t(1,0)=GetY();
+  t(2,0)=GetZ();
+  TMatrixDSym tC(3);
+  {
+  const Float_t *cv=GetCov();
+  tC(0,0)=cv[0]; tC(0,1)=cv[1]; tC(0,2)=cv[2];
+  tC(1,0)=cv[1]; tC(1,1)=cv[3]; tC(1,2)=cv[4];
+  tC(2,0)=cv[2]; tC(2,1)=cv[4]; tC(2,2)=cv[5];
+  }
+
+  TMatrixD m(3,1);
+  m(0,0)=p.GetX();
+  m(1,0)=p.GetY();
+  m(2,0)=p.GetZ();
+  TMatrixDSym mC(3);
+  {
+  const Float_t *cv=p.GetCov();
+  mC(0,0)=cv[0]; mC(0,1)=cv[1]; mC(0,2)=cv[2];
+  mC(1,0)=cv[1]; mC(1,1)=cv[3]; mC(1,2)=cv[4];
+  mC(2,0)=cv[2]; mC(2,1)=cv[4]; mC(2,2)=cv[5];
+  }
+
+  TMatrixDSym tmW(tC);
+  tmW+=mC;
+  tmW.Invert();
+  if (!tmW.IsValid()) return kFALSE; 
+
+  TMatrixD mW(tC,TMatrixD::kMult,tmW);
+  TMatrixD tW(mC,TMatrixD::kMult,tmW);
+
+  TMatrixD mi(mW,TMatrixD::kMult,m);
+  TMatrixD ti(tW,TMatrixD::kMult,t);
+  ti+=mi;
+
+  TMatrixD iC(tC,TMatrixD::kMult,tmW);
+  iC*=mC;
+
+  out.SetXYZ(ti(0,0),ti(1,0),ti(2,0));
+  UShort_t id=p.GetVolumeID();
+  out.SetVolumeID(id);
+
+  return kTRUE;
+}
+
 //______________________________________________________________________________
 Float_t AliTrackPoint::GetAngle() const
 {
@@ -385,3 +485,46 @@ void AliTrackPoint::Print(Option_t *) const
   printf("Z = %12.6f    Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
 
 }
+
+
+//________________________________
+void AliTrackPoint::SetAlignCovMatrix(const TMatrixDSym alignparmtrx){
+  // Add the uncertainty on the cluster position due to alignment
+  // (using the 6x6 AliAlignObj Cov. Matrix alignparmtrx) to the already
+  // present Cov. Matrix 
+
+  TMatrixDSym cov(3);
+  TMatrixD coval(3,3);
+  TMatrixD jacob(3,6);
+  Float_t newcov[6];
+
+  cov(0,0)=fCov[0];
+  cov(1,0)=cov(0,1)=fCov[1];
+  cov(2,0)=cov(0,2)=fCov[2];
+  cov(1,1)=fCov[3];
+  cov(2,1)=cov(1,2)=fCov[4];
+  cov(2,2)=fCov[5];
+  
+  jacob(0,0) = 1;      jacob(1,0) = 0;       jacob(2,0) = 0;
+  jacob(0,1) = 0;      jacob(1,1) = 1;       jacob(2,1) = 0;
+  jacob(0,2) = 0;      jacob(1,2) = 0;       jacob(2,2) = 1;
+  jacob(0,3) = 0;      jacob(1,3) =-fZ;      jacob(2,3) = fY;
+  jacob(0,4) = fZ;     jacob(1,4) = 0;       jacob(2,4) =-fX;
+  jacob(0,5) = -fY;    jacob(1,5) = fX;      jacob(2,5) = 0;
+  
+  TMatrixD jacobT=jacob.T();jacob.T();
+  
+  coval=jacob*alignparmtrx*jacobT+cov;
+
+
+  newcov[0]=coval(0,0);
+  newcov[1]=coval(1,0);
+  newcov[2]=coval(2,0);
+  newcov[3]=coval(1,1);
+  newcov[4]=coval(2,1);
+  newcov[5]=coval(2,2);
+  SetXYZ(fX,fY,fZ,newcov);
+
+}
+