From 4dcdc747c237bc1bed58af4e4d78eb7c00a804be Mon Sep 17 00:00:00 2001
From: cvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Thu, 5 Jul 2007 11:23:31 +0000
Subject: [PATCH] New method to calculate correctly the track extrapolation on
 a reference frame. Or, in other words, the method finds the point of the
 closest approach between two space-points

---
 STEER/AliTrackPointArray.cxx | 59 ++++++++++++++++++++++++++++++++++++
 STEER/AliTrackPointArray.h   |  1 +
 2 files changed, 60 insertions(+)

diff --git a/STEER/AliTrackPointArray.cxx b/STEER/AliTrackPointArray.cxx
index 92b2b15db0f..3b843d7f29f 100644
--- a/STEER/AliTrackPointArray.cxx
+++ b/STEER/AliTrackPointArray.cxx
@@ -23,6 +23,7 @@
 //////////////////////////////////////////////////////////////////////////////
 
 #include <TMath.h>
+#include <TMatrixD.h>
 #include <TMatrixDSym.h>
 
 #include "AliTrackPointArray.h"
@@ -299,6 +300,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
 {
diff --git a/STEER/AliTrackPointArray.h b/STEER/AliTrackPointArray.h
index b7c6aa48145..6e604f35229 100644
--- a/STEER/AliTrackPointArray.h
+++ b/STEER/AliTrackPointArray.h
@@ -38,6 +38,7 @@ class AliTrackPoint : public TObject {
   UShort_t GetVolumeID() const { return fVolumeID; }
 
   Float_t  GetResidual(const AliTrackPoint &p, Bool_t weighted = kFALSE) const;
+  Bool_t   GetPCA(const AliTrackPoint &p, AliTrackPoint &out) const;
 
   Float_t  GetAngle() const;
   AliTrackPoint& Rotate(Float_t alpha) const;
-- 
2.43.0