bug fix

[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx

diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx
index f7d159f3c014f77d0002436b003acf6d22bf414b..fc6a763a4cbd1c43241cb9b6f6b3c2c439714d09 100644 (file)
--- a/TRD/AliTRDseedV1.cxx
+++ b/TRD/AliTRDseedV1.cxx
@@ -83,7 +83,6 @@ AliTRDseedV1::AliTRDseedV1(Int_t det)
   ,fZ(0.)
   ,fS2Y(0.)
   ,fS2Z(0.)
-  ,fC(0.)
   ,fChi2(0.)
 {
   //
@@ -101,6 +100,8 @@ AliTRDseedV1::AliTRDseedV1(Int_t det)
   fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
   fLabels[2]=0;  // number of different labels for tracklet
   memset(fRefCov, 0, 7*sizeof(Double_t));
+  // stand alone curvature
+  fC[0] = 0.; fC[1] = 0.; 
   // covariance matrix [diagonal]
   // default sy = 200um and sz = 2.3 cm 
   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3; 
@@ -130,7 +131,6 @@ AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
   ,fZ(0.)
   ,fS2Y(0.)
   ,fS2Z(0.)
-  ,fC(0.)
   ,fChi2(0.)
 {
   //
@@ -208,7 +208,6 @@ void AliTRDseedV1::Copy(TObject &ref) const
   target.fZ             = fZ;
   target.fS2Y           = fS2Y;
   target.fS2Z           = fS2Z;
-  target.fC             = fC;
   target.fChi2          = fChi2;
   
   memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
@@ -222,6 +221,7 @@ void AliTRDseedV1::Copy(TObject &ref) const
   memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t)); 
   memcpy(target.fLabels, fLabels, 3*sizeof(Int_t)); 
   memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t)); 
+  target.fC[0] = fC[0]; target.fC[1] = fC[1];
   memcpy(target.fCov, fCov, 3*sizeof(Double_t)); 
   
   TObject::Copy(ref);
@@ -271,7 +271,8 @@ void AliTRDseedV1::Reset(Option_t *opt)
   fPt=0.;
   fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
   fS2Y=0.; fS2Z=0.;
-  fC=0.; fChi2 = 0.;
+  fC[0]=0.; fC[1]=0.; 
+  fChi2 = 0.;
 
   memset(fPad, 0, 3*sizeof(Float_t));
   fYref[0] = 0.; fYref[1] = 0.; 
@@ -692,7 +693,7 @@ void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
 }
 
 //____________________________________________________________
-Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
+Int_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
 {
 // Helper function to calculate the square root of the covariance matrix. 
 // The input matrix is stored in the vector c and the result in the vector d. 
@@ -711,6 +712,7 @@ Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
 // Author A.Bercuci <A.Bercuci@gsi.de>
 // Date   Mar 19 2009
 
+  const Double_t kZero(1.e-20);
   Double_t l[2], // eigenvalues
            v[3]; // eigenvectors
   // the secular equation and its solution :
@@ -719,26 +721,33 @@ Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
   // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
   Double_t tr = c[0]+c[2],           // trace
           det = c[0]*c[2]-c[1]*c[1]; // determinant
-  if(TMath::Abs(det)<1.e-20) return -1.;
+  if(TMath::Abs(det)<kZero) return 1;
   Double_t dd = TMath::Sqrt(tr*tr - 4*det);
-  l[0] = .5*(tr + dd);
-  l[1] = .5*(tr - dd);
-  if(l[0]<0. || l[1]<0.) return -1.;
-
+  l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
+  l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
+  if(l[0]<kZero || l[1]<kZero) return 2;
   // the sym V matrix
   // | v00   v10|
   // | v10   v11|
-  Double_t tmp = (l[0]-c[0])/c[1];
-  v[0] = TMath::Sqrt(1./(tmp*tmp+1));
-  v[1] = tmp*v[0];
-  v[2] = v[1]*c[1]/(l[1]-c[2]);
+  Double_t den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
+  if(den<kZero){ // almost diagonal
+    v[0] = TMath::Sign(0., c[1]);
+    v[1] = TMath::Sign(1., (l[0]-c[0]));
+    v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
+  } else {
+    Double_t tmp = 1./TMath::Sqrt(den);
+    v[0] = c[1]* tmp;
+    v[1] = (l[0]-c[0])*tmp;
+    if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
+    else v[2] = v[0]*(l[0]-c[0])/(l[1]-c[2]);
+  }
   // the VD^{1/2}V is: 
   l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
   d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
   d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
   d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
 
-  return 1.;
+  return 0;
 }
 
 //____________________________________________________________
@@ -1337,6 +1346,7 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
 
     // Recalculate cluster error based on tracking information
     c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
+    c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
     sy[n]  = TMath::Sqrt(c->GetSigmaY2());
 
     yc[n]  = recoParam->UseGAUS() ? 
@@ -1352,7 +1362,10 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
   }
 
   // to few clusters
-  if (n < kClmin) return kFALSE; 
+  if (n < kClmin){ 
+    SetErrorMsg(kFitFailed);
+    return kFALSE; 
+  }
 
   // fit XY
   if(!fitterY.Eval()){
@@ -1404,6 +1417,9 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
       qc[n]   = TMath::Abs(c->GetQ());
       xc[n]   = fX0 - c->GetX();
       zc[n]   = c->GetZ();
+      // Recalculate cluster error based on tracking information
+      c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?(z0 - xc[n]*dzdx):-1., dydx);
+      c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
       fitterZ.AddPoint(&xc[n], -qc[n], 1.);
       n--;m++;
     }
@@ -1706,6 +1722,7 @@ void AliTRDseedV1::Print(Option_t *o) const
   AliInfo(Form("Fit | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | ----- |", x, GetY(), TMath::Sqrt(cov[0]), GetZ(), TMath::Sqrt(cov[2]), fYfit[1]));
   AliInfo(Form("Ref | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | %5.2f |", x, fYref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[0]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]))
   AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
+  if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
   AliInfo(Form("dEdx [a.u.]    = %f / %f / %f / %f / %f/ %f / %f / %f", fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7]));
   AliInfo(Form("PID            = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));