]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDseedV1.cxx
bug fix
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx
index 441cbee35b7ed0e561a2c9100b85964a9077aec4..fc6a763a4cbd1c43241cb9b6f6b3c2c439714d09 100644 (file)
@@ -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.; 
@@ -531,6 +532,20 @@ Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
   return p;
 }
 
+//____________________________________________________________________
+Float_t AliTRDseedV1::GetOccupancyTB() const
+{
+// Returns procentage of TB occupied by clusters
+
+  Int_t n(0);
+  AliTRDcluster *c(NULL);
+  for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
+    if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
+    n++;
+  }
+
+  return Float_t(n)/AliTRDtrackerV1::GetNTimeBins();
+}
 
 //____________________________________________________________________
 Float_t* AliTRDseedV1::GetProbability(Bool_t force)
@@ -678,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. 
@@ -697,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 :
@@ -705,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;
 }
 
 //____________________________________________________________
@@ -867,24 +890,26 @@ Bool_t    AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t
 // Author : Alexandru Bercuci <A.Bercuci@gsi.de>
 // Debug  : level >3
 
-  if(!fkReconstructor->GetRecoParam() ){
+  const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
+
+  if(!recoParam){
     AliError("Tracklets can not be used without a valid RecoParam.");
     return kFALSE;
   }
   // Initialize reco params for this tracklet
   // 1. first time bin in the drift region
   Int_t t0 = 14;
-  Int_t kClmin = Int_t(fkReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
+  Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
 
-  Double_t sysCov[5]; fkReconstructor->GetRecoParam()->GetSysCovMatrix(sysCov);        
+  Double_t sysCov[5]; recoParam->GetSysCovMatrix(sysCov);      
   Double_t s2yTrk= fRefCov[0], 
            s2yCl = 0., 
            s2zCl = GetPadLength()*GetPadLength()/12., 
            syRef = TMath::Sqrt(s2yTrk),
            t2    = GetTilt()*GetTilt();
   //define roads
-  Double_t kroady = 1., //fkReconstructor->GetRecoParam() ->GetRoad1y();
-           kroadz = GetPadLength() * fkReconstructor->GetRecoParam()->GetRoadzMultiplicator() + 1.;
+  Double_t kroady = 1., //recoParam->GetRoad1y();
+           kroadz = GetPadLength() * recoParam->GetRoadzMultiplicator() + 1.;
   // define probing cluster (the perfect cluster) and default calibration
   Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
   AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
@@ -997,7 +1022,7 @@ Bool_t     AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t
     }
     lr = ir; if(nr>=3) break;
   }
-  if(fkReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
+  if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
     TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
     UChar_t stat(0);
     if(IsKink()) SETBIT(stat, 1);
@@ -1293,6 +1318,7 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
 
   Int_t n = 0;
   AliTRDcluster *c=NULL, **jc = &fClusters[0];
+  const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
     xc[ic]  = -1.;
     yc[ic]  = 999.;
@@ -1320,9 +1346,10 @@ 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]   = fkReconstructor->GetRecoParam()->UseGAUS() ? 
+    yc[n]  = recoParam->UseGAUS() ? 
       c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
     zc[n]   = c->GetZ();
     //optional tilt correction
@@ -1335,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()){
@@ -1387,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++;
     }
@@ -1689,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]));