Fixes for report #57776 Negative eigenvalues of the track covariance matrix in AliITS...

diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx
index e273da62d88452748acc88b2b4c1607c10f415e6..5bf6da1d3f2b7c9d452f0558d55db0410b48a8d6 100644 (file)
--- a/ITS/AliITStrackV2.cxx
+++ b/ITS/AliITStrackV2.cxx
@@ -389,12 +389,38 @@ Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
   //------------------------------------------------------------------
   //This function improves angular track parameters
   //------------------------------------------------------------------
+  //Store the initail track parameters
+
+  return kTRUE; //PH temporary switched off
+
+  Double_t x = GetX();
+  Double_t alpha = GetAlpha();
+  Double_t par[] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()};
+  Double_t cov[] = {
+    GetSigmaY2(),
+    GetSigmaZY(),
+    GetSigmaZ2(),
+    GetSigmaSnpY(),
+    GetSigmaSnpZ(),
+    GetSigmaSnp2(),
+    GetSigmaTglY(),
+    GetSigmaTglZ(),
+    GetSigmaTglSnp(),
+    GetSigmaTgl2(),
+    GetSigma1PtY(),
+    GetSigma1PtZ(),
+    GetSigma1PtSnp(),
+    GetSigma1PtTgl(),
+    GetSigma1Pt2()
+  }; 
+
+
   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 //Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
   Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
   Double_t zv = xyz[2];                // local frame
 
-  Double_t dy = Par(0) - yv, dz = Par(1) - zv;
+  Double_t dy = par[0] - yv, dz = par[1] - zv;
   Double_t r2=GetX()*GetX() + dy*dy;
   Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
   Double_t beta2=p2/(p2 + GetMass()*GetMass());
@@ -410,28 +436,31 @@ Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
     Double_t sigma2p = theta2*(1.-GetSnp())*(1.+GetSnp())*(1. + GetTgl()*GetTgl());
     Double_t ovSqr2 = 1./TMath::Sqrt(r2);
     Double_t tfact = ovSqr2*(1.-dy*ovSqr2)*(1.+dy*ovSqr2);
-    sigma2p += Cov(0)*tfact*tfact;
+    sigma2p += cov[0]*tfact*tfact;
     sigma2p += ers[1]*ers[1]/r2;
-    sigma2p += 0.25*Cov(14)*cnv*cnv*GetX()*GetX();
-    Double_t eps2p=sigma2p/(Cov(5) + sigma2p);
-    Par(0) += Cov(3)/(Cov(5) + sigma2p)*(parp - GetSnp());
-    Par(2)  = eps2p*GetSnp() + (1 - eps2p)*parp;
-    Cov(5) *= eps2p;
-    Cov(3) *= eps2p;
+    sigma2p += 0.25*cov[14]*cnv*cnv*GetX()*GetX();
+    Double_t eps2p=sigma2p/(cov[5] + sigma2p);
+    par[0] += cov[3]/(cov[5] + sigma2p)*(parp - GetSnp());
+    par[2]  = eps2p*GetSnp() + (1 - eps2p)*parp;
+    cov[5] *= eps2p;
+    cov[3] *= eps2p;
   }
   {
     Double_t parl=0.5*GetC()*dz/TMath::ASin(0.5*GetC()*TMath::Sqrt(r2));
     Double_t sigma2l=theta2;
-    sigma2l += Cov(2)/r2 + Cov(0)*dy*dy*dz*dz/(r2*r2*r2);
+    sigma2l += cov[2]/r2 + cov[0]*dy*dy*dz*dz/(r2*r2*r2);
     sigma2l += ers[2]*ers[2]/r2;
-    Double_t eps2l = sigma2l/(Cov(9) + sigma2l);
-    Par(1) += Cov(7 )/(Cov(9) + sigma2l)*(parl - Par(3));
-    Par(4) += Cov(13)/(Cov(9) + sigma2l)*(parl - Par(3));
-    Par(3)  = eps2l*Par(3) + (1-eps2l)*parl;
-    Cov(9) *= eps2l; 
-    Cov(13)*= eps2l; 
-    Cov(7) *= eps2l; 
+    Double_t eps2l = sigma2l/(cov[9] + sigma2l);
+    par[1] += cov[7 ]/(cov[9] + sigma2l)*(parl - par[3]);
+    par[4] += cov[13]/(cov[9] + sigma2l)*(parl - par[3]);
+    par[3]  = eps2l*par[3] + (1-eps2l)*parl;
+    cov[9] *= eps2l; 
+    cov[13]*= eps2l; 
+    cov[7] *= eps2l; 
   }
+
+  Set(x,alpha,par,cov);
+
   if (!Invariant()) return kFALSE;
 
   return kTRUE;

diff --git a/STEER/AliExternalTrackParam.cxx b/STEER/AliExternalTrackParam.cxx
index a6f0f677ecbcdf606bd7aab59af4283abc51b512..9fe6d955787d86c448cc96600e8061e4843deeca 100644 (file)
--- a/STEER/AliExternalTrackParam.cxx
+++ b/STEER/AliExternalTrackParam.cxx
@@ -25,6 +25,9 @@
 // are implemented.
 // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                            //
 ///////////////////////////////////////////////////////////////////////////////
+#include <cassert>
+
+#include <TVectorD.h>
 #include <TMatrixDSym.h>
 #include <TPolyMarker3D.h>
 #include <TVector3.h>
@@ -62,6 +65,7 @@ AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track)
   //
   for (Int_t i = 0; i < 5; i++) fP[i] = track.fP[i];
   for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i];
+  CheckCovariance();
 }
 
 //_____________________________________________________________________________
@@ -78,6 +82,7 @@ AliExternalTrackParam& AliExternalTrackParam::operator=(const AliExternalTrackPa
 
     for (Int_t i = 0; i < 5; i++) fP[i] = trkPar.fP[i];
     for (Int_t i = 0; i < 15; i++) fC[i] = trkPar.fC[i];
+    CheckCovariance();
   }
 
   return *this;
@@ -96,6 +101,7 @@ AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha,
   //
   for (Int_t i = 0; i < 5; i++)  fP[i] = param[i];
   for (Int_t i = 0; i < 15; i++) fC[i] = covar[i];
+  CheckCovariance();
 }
 
 //_____________________________________________________________________________
@@ -234,6 +240,8 @@ void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
   fC[13] = b1/b3-b2*fC[8]/b3;
   fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
 
+  CheckCovariance();
+
   return;
 }
 
@@ -258,6 +266,7 @@ void AliExternalTrackParam::AddCovariance(const Double_t c[15]) {
     fC[3] +=c[3];  fC[4] +=c[4];  fC[5] +=c[5];
     fC[6] +=c[6];  fC[7] +=c[7];  fC[8] +=c[8];  fC[9] +=c[9];
     fC[10]+=c[10]; fC[11]+=c[11]; fC[12]+=c[12]; fC[13]+=c[13]; fC[14]+=c[14];
+    CheckCovariance();
 }
 
 
@@ -413,6 +422,8 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterial
   fC44 += cC44;
   fP4  *= cP4;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -479,6 +490,8 @@ Bool_t AliExternalTrackParam::CorrectForMaterial
   fC44 += cC44;
   fP4  *= cP4;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -639,6 +652,8 @@ Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
   fC40 *= ca;
   fC42 *= rr;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -711,6 +726,8 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   fC32 += b32;
   fC42 += b42;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -1107,6 +1124,8 @@ Bool_t AliExternalTrackParam::Update(Double_t p[2], Double_t cov[3]) {
 
   fC44-=k40*c04+k41*c14; 
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -1387,7 +1406,6 @@ Double_t b[3], Double_t maxd, Double_t dz[2], Double_t covar[3]) {
   return kTRUE;
 }
 
-
 void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
   //----------------------------------------------------------------
   // This function returns a unit vector along the track direction
@@ -1911,6 +1929,7 @@ Bool_t AliExternalTrackParam::PropagateToBxByBz(Double_t xk, const Double_t b[3]
   fC32 += b32;
   fC42 += b42;
 
+  CheckCovariance();
   
   // Appoximate step length
   Double_t step=dx*TMath::Abs(r2 + f2*(f1+f2)/(r1+r2));
@@ -2076,3 +2095,74 @@ Bool_t AliExternalTrackParam::Translate(Double_t *vTrasl,Double_t *covV){
 
   return kTRUE;
  }
+
+void AliExternalTrackParam::CheckCovariance() {
+
+  // This function forces the diagonal elements of the covariance matrix to be positive.
+  // In case the diagonal element is bigger than the maximal allowed value, it is set to
+  // the limit and the off-diagonal elements that correspond to it are set to zero.
+
+    fC[0] = TMath::Abs(fC[0]);
+    if (fC[0]>kC0max) {
+      fC[0] = kC0max;
+      fC[1] = 0;
+      fC[3] = 0;
+      fC[6] = 0;
+      fC[10] = 0;
+    }
+    fC[2] = TMath::Abs(fC[2]);
+    if (fC[2]>kC2max) {
+      fC[2] = kC2max;
+      fC[1] = 0;
+      fC[4] = 0;
+      fC[7] = 0;
+      fC[11] = 0;
+    }
+    fC[5] = TMath::Abs(fC[5]);
+    if (fC[5]>kC5max) {
+      fC[5] = kC5max;
+      fC[3] = 0;
+      fC[4] = 0;
+      fC[8] = 0;
+      fC[12] = 0;
+    }
+    fC[9] = TMath::Abs(fC[9]);
+    if (fC[9]>kC9max) {
+      fC[9] = kC9max;
+      fC[6] = 0;
+      fC[7] = 0;
+      fC[8] = 0;
+      fC[13] = 0;
+    }
+    fC[14] = TMath::Abs(fC[14]);
+    if (fC[14]>kC14max) {
+      fC[14] = kC14max;
+      fC[10] = 0;
+      fC[11] = 0;
+      fC[12] = 0;
+      fC[13] = 0;
+    }
+    
+    // The part below is used for tests and normally is commented out    
+//     TMatrixDSym m(5);
+//     TVectorD eig(5);
+    
+//     m(0,0)=fC[0];
+//     m(1,0)=fC[1];  m(1,1)=fC[2];
+//     m(2,0)=fC[3];  m(2,1)=fC[4];  m(2,2)=fC[5];
+//     m(3,0)=fC[6];  m(3,1)=fC[7];  m(3,2)=fC[8];  m(3,3)=fC[9];
+//     m(4,0)=fC[10]; m(4,1)=fC[11]; m(4,2)=fC[12]; m(4,3)=fC[13]; m(4,4)=fC[14];
+    
+//     m(0,1)=m(1,0);
+//     m(0,2)=m(2,0); m(1,2)=m(2,1);
+//     m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
+//     m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
+//     m.EigenVectors(eig);
+
+//     //    assert(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0);
+//     if (!(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0)) {
+//       AliWarning("Negative eigenvalues of the covariance matrix!");
+//       this->Print();
+//       eig.Print();
+//     }
+}

diff --git a/STEER/AliExternalTrackParam.h b/STEER/AliExternalTrackParam.h
index 112ffb85baeb97af448228ffcffd415b31b36759..87489e980d51dc299e1da5c683e2b4792173ed34 100644 (file)
--- a/STEER/AliExternalTrackParam.h
+++ b/STEER/AliExternalTrackParam.h
@@ -28,6 +28,12 @@ const Double_t kMostProbablePt=0.35;
 class AliVVertex;
 class TPolyMarker3D; 
 
+const Double_t kC0max=100*100, // SigmaY<=100cm
+               kC2max=100*100, // SigmaZ<=100cm
+               kC5max=1*1,     // SigmaSin<=1
+               kC9max=1*1,     // SigmaTan<=1
+               kC14max=100*100; // Sigma1/Pt<=100 1/GeV
+
 class AliExternalTrackParam: public AliVTrack {
  public:
   AliExternalTrackParam();
@@ -46,6 +52,9 @@ class AliExternalTrackParam: public AliVTrack {
     fX=x; fAlpha=alpha;
     for (Int_t i = 0; i < 5; i++)  fP[i] = param[i];
     for (Int_t i = 0; i < 15; i++) fC[i] = covar[i];
+
+    CheckCovariance();
+
   }
 
   void Set(Double_t xyz[3],Double_t pxpypz[3],Double_t cv[21],Short_t sign);
@@ -204,7 +213,10 @@ class AliExternalTrackParam: public AliVTrack {
 
   virtual Bool_t Translate(Double_t *vTrasl,Double_t *covV);
 
- protected:
+  void CheckCovariance();
+
+/*  protected: */
+ private:
   Double_t &Par(Int_t i) {return fP[i];}
   Double_t &Cov(Int_t i) {return fC[i];}
   
@@ -223,11 +235,26 @@ inline void AliExternalTrackParam::ResetCovariance(Double_t s2) {
   //
   // Reset the covarince matrix to "something big"
   //
-    fC[0]*= s2;
-    fC[1] = 0.;  fC[2]*= s2;
-    fC[3] = 0.;  fC[4] = 0.;  fC[5]*= s2;
-    fC[6] = 0.;  fC[7] = 0.;  fC[8] = 0.;  fC[9]*= s2;
-    fC[10]= 0.;  fC[11]= 0.;  fC[12]= 0.;  fC[13]= 0.;  fC[14]*=s2;
+
+  s2 = TMath::Abs(s2);
+  Double_t fC0=fC[0]*s2,
+           fC2=fC[2]*s2,
+           fC5=fC[5]*s2,
+           fC9=fC[9]*s2,
+           fC14=fC[14]*s2;
+ 
+  if (fC0>kC0max)  fC0 = kC0max;
+  if (fC2>kC2max)  fC2 = kC2max;
+  if (fC5>kC5max)  fC5 = kC5max;
+  if (fC9>kC9max)  fC9 = kC9max;
+  if (fC14>kC14max)  fC14 = kC14max;
+
+
+    fC[0] = fC0;
+    fC[1] = 0.;  fC[2] = fC2;
+    fC[3] = 0.;  fC[4] = 0.;  fC[5] = fC5;
+    fC[6] = 0.;  fC[7] = 0.;  fC[8] = 0.;  fC[9] = fC9;
+    fC[10]= 0.;  fC[11]= 0.;  fC[12]= 0.;  fC[13]= 0.;  fC[14] = fC14;
 }