Replaced AliInfo with AliDebug

[u/mrichter/AliRoot.git] / STEER / AliTrackFitterKalman.cxx

diff --git a/STEER/AliTrackFitterKalman.cxx b/STEER/AliTrackFitterKalman.cxx
index 5251c96a22153d40653f1324ab85d6a0d1f4c4b1..71523b3d76c1566f240c886d09b69e34402a020e 100755 (executable)
--- a/STEER/AliTrackFitterKalman.cxx
+++ b/STEER/AliTrackFitterKalman.cxx
@@ -32,7 +32,7 @@
 #include <TMatrixD.h>
 #include <TMatrixDSym.h>
 
-//#include "AliLog.h"
+#include "AliLog.h"
 #include "AliTrackFitterKalman.h"
 
 ClassImp(AliTrackFitterKalman)
@@ -42,8 +42,7 @@ const Double_t AliTrackFitterKalman::fgkMaxChi2=33.;
 AliTrackFitterKalman::
 AliTrackFitterKalman(AliTrackPointArray *array, Bool_t owner):
   AliTrackFitter(array,owner),
-  fMaxChi2(fgkMaxChi2),
-  fSeed(kFALSE)
+  fMaxChi2(fgkMaxChi2)
 {
   //
   // Constructor
@@ -55,7 +54,7 @@ Bool_t AliTrackFitterKalman::Begin(Int_t first, Int_t last) {
   // Make a seed out of the track points with the indices "first" and "last". 
   // This is the "default" seed.
   //
-  if (fSeed) return kTRUE;
+  fPoints->Sort();
   AliTrackPoint fp,lp;
   fPoints->GetPoint(fp,first); fPoints->GetPoint(lp,last);
   return MakeSeed(&fp,&lp);
@@ -69,17 +68,19 @@ MakeSeed(const AliTrackPoint *p0, const AliTrackPoint *p1) {
   //
   Double_t x0=p0->GetX(), y0=p0->GetY(), z0=p0->GetZ();
   Double_t dx=p1->GetX()-x0, dy=p1->GetY()-y0, dz=p1->GetZ()-z0;
+  if (dy==0.) { 
+     AliError("Seeds perpendicular to Y axis are not allowed !"); 
+     return kFALSE; 
+  }
   Double_t vx=dx/dy;
   Double_t vz=dz/dy;
   Double_t par[5]={x0,y0,z0,vx,vz};
 
   const Float_t *cv0=p0->GetCov();
   const Float_t *cv1=p1->GetCov();
-  Double_t rdx2=(cv0[0]+cv1[0])/dx/dx;
   Double_t rdy2=(cv0[3]+cv1[3])/dy/dy;
-  Double_t rdz2=(cv0[5]+cv1[5])/dz/dz;
-  Double_t svx2=vx*vx*(rdx2+rdy2);     
-  Double_t svz2=vz*vz*(rdz2+rdy2);     
+  Double_t svx2=(cv0[0]+cv1[0])/dy/dy + vx*vx*rdy2;     
+  Double_t svz2=(cv0[5]+cv1[5])/dy/dy + vz*vz*rdy2;     
   Double_t cov[15]={
      cv0[0],
      cv0[1],cv0[3],
@@ -179,6 +180,7 @@ Bool_t AliTrackFitterKalman::Update(const AliTrackPoint *p, Double_t chi2) {
   v(0,0)=cov[0]+c(0,0); v(0,1)=cov[1]+c(0,1); v(0,2)=cov[2]+c(0,2);
   v(1,0)=cov[1]+c(1,0); v(1,1)=cov[3]+c(1,1); v(1,2)=cov[4]+c(1,2);
   v(2,0)=cov[2]+c(2,0); v(2,1)=cov[4]+c(2,1); v(2,2)=cov[5]+c(2,2);
+  if (TMath::Abs(v.Determinant()) < 1.e-33) return kFALSE;
   v.Invert();
 
   TMatrixD ch(5,3);
@@ -240,13 +242,13 @@ Bool_t AliTrackFitterKalman::GetPCA(const AliTrackPoint &p, AliTrackPoint &i) co
   TMatrixDSym tC(3);
   {
   Double_t sig2=s*s/12.;
-  TMatrixDSym &cv=*fCov;
-  tC(0,0) = cv(0,0) + 2*s*cv(3,0) + s*s*cv(3,3) + vx*vx*sig2;
-  tC(1,0) = cv(1,0) + s*cv(1,3) + vx*sig2; 
-  tC(1,1) = cv(1,1) + sig2;
-  tC(2,0) = cv(2,0) + s*(cv(4,0) + cv(2,3)) + s*s*cv(4,3) + vx*vz*sig2;
-  tC(2,1) = cv(2,1) + s*cv(4,1) + vz*sig2;
-  tC(2,2) = cv(2,2) + 2*s*cv(4,2) + s*s*cv(4,4) + vz*vz*sig2;
+
+  tC(0,0) = vx*vx*sig2;
+  tC(1,0) = vx*sig2; 
+  tC(1,1) = sig2;
+  tC(2,0) = vx*vz*sig2;
+  tC(2,1) = vz*sig2;
+  tC(2,2) = vz*vz*sig2;
 
   tC(0,1) = tC(1,0); tC(0,2) = tC(2,0);
                      tC(1,2) = tC(2,1);
@@ -267,6 +269,7 @@ Bool_t AliTrackFitterKalman::GetPCA(const AliTrackPoint &p, AliTrackPoint &i) co
 
   TMatrixDSym tmW(tC);
   tmW+=mC;
+  if (TMath::Abs(tmW.Determinant()) < 1.e-33) return kFALSE;
   tmW.Invert();
 
   TMatrixD mW(tC,TMatrixD::kMult,tmW);
@@ -278,6 +281,14 @@ Bool_t AliTrackFitterKalman::GetPCA(const AliTrackPoint &p, AliTrackPoint &i) co
 
   TMatrixD iC(tC,TMatrixD::kMult,tmW);
   iC*=mC;
+  Double_t sig2=1.;  // Releasing the matrix by 1 cm along the track direction 
+  iC(0,0) += vx*vx*sig2;
+  iC(0,1) += vx*sig2; 
+  iC(1,1) += sig2;
+  iC(0,2) += vx*vz*sig2;
+  iC(1,2) += vz*sig2;
+  iC(2,2) += vz*vz*sig2;
+
 
   Float_t cov[6]={
     iC(0,0), iC(0,1), iC(0,2),
@@ -316,5 +327,4 @@ AliTrackFitterKalman::SetSeed(const Double_t par[5], const Double_t cov[15]) {
   cv(0,3)=cv(3,0); cv(1,3)=cv(3,1); cv(2,3)=cv(3,2);
   cv(0,4)=cv(4,0); cv(1,4)=cv(4,1); cv(2,4)=cv(4,2); cv(3,4)=cv(4,3);
 
-  fSeed=kTRUE;
 }