Change THnSparseD to THnF in truncated mean related calibration code to reduce the...

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

diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx
index 4bfe661664105deea2d8ce47b37ee19a0ef7b934..7f166c88a568d057275f8afcd0cb0d6630e8286a 100644 (file)
--- a/TRD/AliTRDtrackerV1.cxx
+++ b/TRD/AliTRDtrackerV1.cxx
@@ -77,7 +77,7 @@ TLinearFitter* AliTRDtrackerV1::fgTiltedRieman = NULL;
 TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = NULL;
 
 //____________________________________________________________________
-AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) 
+AliTRDtrackerV1::AliTRDtrackerV1(const AliTRDReconstructor *rec) 
   :AliTracker()
   ,fkReconstructor(NULL)
   ,fkRecoParam(NULL)
@@ -351,12 +351,12 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
     if (expectedClr<0){      
       seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop);
       continue;
-    }
-
-    if(expectedClr){
+    } else {
       nFound++;  
-      // computes PID for track
+      // compute PID
       track.CookPID();
+      //compute MC label
+      track.CookLabel(1. - AliTRDReconstructor::GetLabelFraction());
       // update calibration references using this track
       if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track);
       // save calibration object
@@ -370,12 +370,10 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
       track.UpdateESDtrack(seed);
     }
 
-    if ((TMath::Abs(track.GetC(track.GetBz()) - p4) / TMath::Abs(p4) < 0.2) ||(track.Pt() > 0.8)) {
-
-      // Make backup for back propagation
+    // Make backup for back propagation
+    if ((TMath::Abs(track.GetC(track.GetBz()) - p4) / TMath::Abs(p4) < 0.2) || (track.Pt() > 0.8)) {
       Int_t foundClr = track.GetNumberOfClusters();
       if (foundClr >= foundMin) {
-        track.CookLabel(1. - AliTRDReconstructor::GetLabelFraction());
         //if(track.GetBackupTrack()) UseClusters(track.GetBackupTrack());
 
         // Sign only gold tracks
@@ -735,7 +733,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     // rough estimate of the entry point
     if (!t.GetProlongation(fR[ily], y, z)){
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kProlongation);
+      t.SetErrStat(AliTRDtrackV1::kProlongation);
       AliDebug(4, Form("Failed Rough Prolongation to ly[%d] x[%7.2f] y[%7.2f] z[%7.2f]", ily, fR[ily], y, z));
       break;
     }
@@ -756,23 +754,23 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       // propagate to the default radial position
       if(fR[ily] > (AliTRDReconstructor::GetMaxStep() + t.GetX()) && !PropagateToX(t, fR[ily], AliTRDReconstructor::GetMaxStep())){
         n=-1; 
-        t.SetStatus(AliTRDtrackV1::kPropagation);
+        t.SetErrStat(AliTRDtrackV1::kPropagation);
         AliDebug(4, "Failed Propagation [Missing Geometry]");
         break;
       }
       if(!AdjustSector(&t)){
         n=-1; 
-        t.SetStatus(AliTRDtrackV1::kAdjustSector);
+        t.SetErrStat(AliTRDtrackV1::kAdjustSector);
         AliDebug(4, "Failed Adjust Sector [Missing Geometry]");
         break;
       }
       if(TMath::Abs(t.GetSnp()) > AliTRDReconstructor::GetMaxSnp()){
         n=-1; 
-        t.SetStatus(AliTRDtrackV1::kSnp);
+        t.SetErrStat(AliTRDtrackV1::kSnp);
         AliDebug(4, "Failed Max Snp [Missing Geometry]");
         break;
       }
-      t.SetStatus(AliTRDtrackV1::kGeometry, ily);
+      t.SetErrStat(AliTRDtrackV1::kGeometry, ily);
       continue;
     }
 
@@ -786,19 +784,19 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     x = glb[0] - AliTRDReconstructor::GetMaxStep();
     if(x > (AliTRDReconstructor::GetMaxStep() + t.GetX()) && !PropagateToX(t, x, AliTRDReconstructor::GetMaxStep())){
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kPropagation);
+      t.SetErrStat(AliTRDtrackV1::kPropagation);
       AliDebug(4, Form("Failed Initial Propagation to x[%7.2f]", x));
       break;
     }
     if(!AdjustSector(&t)){
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kAdjustSector);
+      t.SetErrStat(AliTRDtrackV1::kAdjustSector);
       AliDebug(4, "Failed Adjust Sector Start");
       break;
     }
     if(TMath::Abs(t.GetSnp()) > AliTRDReconstructor::GetMaxSnp()) {
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kSnp);
+      t.SetErrStat(AliTRDtrackV1::kSnp);
       AliDebug(4, Form("Failed Max Snp[%f] MaxSnp[%f]", t.GetSnp(), AliTRDReconstructor::GetMaxSnp()));
       break;
     }
@@ -814,7 +812,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     if(doRecalculate){
       det = AliTRDgeometry::GetDetector(ily, stk, sm);
       if(!(matrix = fGeom->GetClusterMatrix(det))){ 
-        t.SetStatus(AliTRDtrackV1::kGeometry, ily);
+        t.SetErrStat(AliTRDtrackV1::kGeometry, ily);
         AliDebug(4, Form("Failed Geometry Matrix ly[%d]", ily));
         continue;
       }
@@ -825,12 +823,12 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     // check if track is well inside fiducial volume 
     if (!t.GetProlongation(x+AliTRDReconstructor::GetMaxStep(), y, z)) {
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kProlongation);
+      t.SetErrStat(AliTRDtrackV1::kProlongation);
       AliDebug(4, Form("Failed Prolongation to x[%7.2f] y[%7.2f] z[%7.2f]", x+AliTRDReconstructor::GetMaxStep(), y, z));
       break;
     }
     if(fGeom->IsOnBoundary(det, y, z, .5)){ 
-      t.SetStatus(AliTRDtrackV1::kBoundary, ily);
+      t.SetErrStat(AliTRDtrackV1::kBoundary, ily);
       AliDebug(4, "Failed Track on Boundary");
       continue;
     }
@@ -840,24 +838,24 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       AliDebug(3, Form("Building tracklet det[%d]", det));
       // check data in supermodule
       if(!fTrSec[sm].GetNChambers()){ 
-        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        t.SetErrStat(AliTRDtrackV1::kNoClusters, ily);
         AliDebug(4, "Failed NoClusters");
         continue;
       }
       if(fTrSec[sm].GetX(ily) < 1.){ 
-        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        t.SetErrStat(AliTRDtrackV1::kNoClusters, ily);
         AliDebug(4, "Failed NoX");
         continue;
       }
       
       // check data in chamber
       if(!(chamber = fTrSec[sm].GetChamber(stk, ily))){ 
-        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        t.SetErrStat(AliTRDtrackV1::kNoClusters, ily);
         AliDebug(4, "Failed No Detector");
         continue;
       }
       if(chamber->GetNClusters() < fgNTimeBins*fkRecoParam ->GetFindableClusters()){ 
-        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        t.SetErrStat(AliTRDtrackV1::kNoClusters, ily);
         AliDebug(4, "Failed Not Enough Clusters in Detector");
         continue;
       }      
@@ -871,12 +869,16 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       ptrTracklet->SetX0(glb[0]+driftLength);
       if(!ptrTracklet->Init(&t)){
         n=-1; 
-        t.SetStatus(AliTRDtrackV1::kTrackletInit);
+        t.SetErrStat(AliTRDtrackV1::kTrackletInit);
         AliDebug(4, "Failed Tracklet Init");
         break;
       }
-      if(!ptrTracklet->AttachClusters(chamber, kTRUE, t.Charge()>0?kTRUE:kFALSE, fEventInFile)){
-        t.SetStatus(AliTRDtrackV1::kNoAttach, ily);
+      // Select attachment base on track to B field sign not only track charge which is buggy
+      // mark kFALSE same sign tracks and kTRUE opposite sign tracks
+      // A.Bercuci 3.11.2011
+      Float_t prod(t.GetBz()*t.Charge());
+      if(!ptrTracklet->AttachClusters(chamber, kTRUE, prod<0.?kTRUE:kFALSE, fEventInFile)){
+        t.SetErrStat(AliTRDtrackV1::kNoAttach, ily);
         if(debugLevel>3){
           AliTRDseedV1 trackletCp(*ptrTracklet);
           UChar_t status(t.GetStatusTRD(ily));
@@ -890,7 +892,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       }
       AliDebug(3, Form("Number of Clusters in Tracklet: %d", ptrTracklet->GetN()));
       if(ptrTracklet->GetN() < fgNTimeBins*fkRecoParam->GetFindableClusters()){
-        t.SetStatus(AliTRDtrackV1::kNoClustersTracklet, ily);
+        t.SetErrStat(AliTRDtrackV1::kNoClustersTracklet, ily);
         if(debugLevel>3){
           AliTRDseedV1 trackletCp(*ptrTracklet);
           UChar_t status(t.GetStatusTRD(ily));
@@ -911,26 +913,26 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     // 0 : no correction
     // 2 : pseudo tilt correction
     if(!ptrTracklet->FitRobust(t.Charge()>0?kTRUE:kFALSE)){
-      t.SetStatus(AliTRDtrackV1::kNoFit, ily);
+      t.SetErrStat(AliTRDtrackV1::kNoFit, ily);
       AliDebug(4, "Failed Tracklet Fit");
       continue;
     } 
     x = ptrTracklet->GetX(); //GetX0();
     if(x > (AliTRDReconstructor::GetMaxStep() + t.GetX()) && !PropagateToX(t, x, AliTRDReconstructor::GetMaxStep())) {
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kPropagation);
+      t.SetErrStat(AliTRDtrackV1::kPropagation);
       AliDebug(4, Form("Failed Propagation to Tracklet x[%7.2f]", x));
       break;
     }
     if(!AdjustSector(&t)) {
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kAdjustSector);
+      t.SetErrStat(AliTRDtrackV1::kAdjustSector);
       AliDebug(4, "Failed Adjust Sector");
       break;
     }
     if(TMath::Abs(t.GetSnp()) > AliTRDReconstructor::GetMaxSnp()) {
       n=-1; 
-      t.SetStatus(AliTRDtrackV1::kSnp);
+      t.SetErrStat(AliTRDtrackV1::kSnp);
       AliDebug(4, Form("Failed Max Snp[%f] MaxSnp[%f]", t.GetSnp(), AliTRDReconstructor::GetMaxSnp()));
       break;
     }
@@ -939,7 +941,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     Double_t chi2 = ((AliExternalTrackParam)t).GetPredictedChi2(p, cov);
     // update Kalman with the TRD measurement
     if(chi2>1e+10){ // TODO
-      t.SetStatus(AliTRDtrackV1::kChi2, ily);
+      t.SetErrStat(AliTRDtrackV1::kChi2, ily);
       if(debugLevel > 2){
         UChar_t status(t.GetStatusTRD());
         AliTRDseedV1  trackletCp(*ptrTracklet);
@@ -962,7 +964,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     if(kUseTRD){
       if(!((AliExternalTrackParam&)t).Update(p, cov)) {
         n=-1; 
-        t.SetStatus(AliTRDtrackV1::kUpdate);
+        t.SetErrStat(AliTRDtrackV1::kUpdate);
         if(debugLevel > 2){
           UChar_t status(t.GetStatusTRD());
           AliTRDseedV1  trackletCp(*ptrTracklet);
@@ -984,7 +986,10 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
   
 
     // register tracklet with the tracker and track
-    ptrTracklet->Update(&t);
+    // Save inside the tracklet the track parameters BEFORE track update.
+    // Commented out their overwriting AFTER track update
+    // A.Bercuci 3.11.2011
+    //ptrTracklet->Update(&t); 
     ptrTracklet = SetTracklet(ptrTracklet);
     Int_t index(fTracklets->GetEntriesFast()-1);
     t.SetTracklet(ptrTracklet, index);
@@ -1156,7 +1161,10 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub
   Double_t b = fitter->GetParameter(1);
   Double_t curvature = a/TMath::Sqrt(b*b + 1);
 
-  Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
+  Float_t chi2track = 0.0;
+  if (nPoints > 0) {
+    chi2track = fitter->GetChisquare()/Double_t(nPoints);
+  }
   for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++)
     tracklets[ip].SetC(curvature, 1);
 
@@ -1804,8 +1812,8 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
   //
 
   // Current track X-position
-  Double_t xpos = t.GetX(),
-           mass = t.GetMass();
+  Double_t xpos = t.GetX()/*,
+           mass = t.GetMass()*/;
 
   // Direction: inward or outward
   Double_t dir  = (xpos < xToGo) ? 1.0 : -1.0;
@@ -1843,7 +1851,7 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
     // Propagate the track to the X-position after the next step
     if (!t.PropagateTo(x, param[1], param[0]*param[4])) return 0;
 
-    // Correct for mean material budget
+/*    // Correct for mean material budget
     Double_t dEdx(0.),
              bg(t.GetP()/mass);
     if(AliLog::GetDebugLevel("TRD", "AliTRDtrackerV1")>=3){
@@ -1877,7 +1885,7 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
     }
     if(AliLog::GetDebugLevel("TRD", "AliTRDtrackerV1")>=2) printf("D-AliTRDtrackerV1::PropagateTo(): dEdx(bg=%e, m=%e)= %e[GeV/cm]\n", bg, mass, dEdx);
     if (!t.CorrectForMeanMaterialdEdx(param[1], dir*param[0]*param[4], mass, dEdx)) return 0;
-
+*/
     // Rotate the track if necessary
     if(!AdjustSector(&t)) return 0;
 
@@ -2102,7 +2110,7 @@ Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *const track)
 
 
 //____________________________________________________________________
-AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *const track, Int_t p, Int_t &idx)
+AliTRDseedV1* AliTRDtrackerV1::GetTracklet(const AliTRDtrackV1 *const track, Int_t p, Int_t &idx)
 {
   // Find tracklet for TRD track <track>
   // Parameters
@@ -3151,7 +3159,9 @@ Bool_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRD
     lQuality[jLayer] = bseed[jLayer].GetQuality(kTRUE);
     quality    += lQuality[jLayer];
   }
-  quality /= rLayers;
+  if (rLayers > 0) {
+    quality /= rLayers;
+  }
   AliDebug(2, Form("Start N[%d] Q[%f] chi2[%f]", rLayers, quality, chi2));
 
   for (Int_t iter = 0; iter < 4; iter++) {
@@ -4078,7 +4088,7 @@ Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::Eval(){
 }
 
 //_____________________________________________________________________________
-void AliTRDtrackerV1::AliTRDtrackFitterRieman::UpdateFitters(AliTRDseedV1 * const tracklet){
+void AliTRDtrackerV1::AliTRDtrackFitterRieman::UpdateFitters(const AliTRDseedV1 * const tracklet){
   //
   // Does the transformations and updates the fitters
   // The following transformation is applied