]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDtrackerV1.cxx
update from Prabhat
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerV1.cxx
index 6afdb846085661776859166e4fc839bec6df7c3e..3dbadd98800202a14b7d34539dce5b6b1b7162d2 100644 (file)
@@ -61,6 +61,7 @@ ClassImp(AliTRDtrackerV1)
 ClassImp(AliTRDtrackerV1::AliTRDLeastSquare)
 ClassImp(AliTRDtrackerV1::AliTRDtrackFitterRieman)
 
+AliTRDtrackerV1::ETRDtrackerV1BetheBloch AliTRDtrackerV1::fgBB = AliTRDtrackerV1::kGeant;
 Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
   0.5112, 0.5112, 0.5112, 0.0786, 0.0786,
   0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
@@ -76,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)
@@ -132,7 +133,7 @@ AliTRDtrackerV1::~AliTRDtrackerV1()
   //
   // Destructor
   //
-  
+
   if(fgRieman) delete fgRieman; fgRieman = NULL;
   if(fgTiltedRieman) delete fgTiltedRieman; fgTiltedRieman = NULL;
   if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained; fgTiltedRiemanConstrained = NULL;
@@ -140,7 +141,8 @@ AliTRDtrackerV1::~AliTRDtrackerV1()
   if(fTracksESD){ fTracksESD->Delete(); delete fTracksESD; }
   if(fTracks) {fTracks->Delete(); delete fTracks;}
   if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
-  if(fClusters) {
+  if(IsClustersOwner() && fClusters) {
+    AliInfo(Form("tracker[%p] removing %d own clusters @ %p", (void*)this, fClusters->GetEntries(), (void*)fClusters));
     fClusters->Delete(); delete fClusters;
   }
   if(fGeom) delete fGeom;
@@ -313,8 +315,21 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
   
     // Check the seed status
     ULong_t status = seed->GetStatus();
-    if ((status & AliESDtrack::kTPCout) == 0) continue;
     if ((status & AliESDtrack::kTRDout) != 0) continue;
+    if ((status & AliESDtrack::kTPCout)){
+      AliDebug(3, Form("Prolongate seed[%2d] which is TPC.", iSeed));
+      // set steering parameters for TPC
+      //fkRecoParam->SetTrackParam(kTPC);
+/*    } else {
+      if ((status & AliESDtrack::kITSout)){
+        AliDebug(3, Form("Prolongate seed[%2d] which is ITS.", iSeed));
+        // set steering parameters for ITS
+        //fkRecoParam->SetTrackParam(kITS);
+        // rotate
+        Float_t  globalToTracking = AliTRDgeometry::GetAlpha()*(Int_t(seed->GetAlpha()/AliTRDgeometry::GetAlpha()) + (seed->GetAlpha()>0. ? 0.5 : -0.5));
+        if(!seed->Rotate(globalToTracking)) continue;
+      } else continue;*/
+    } else continue;
 
     // Propagate to the entrance in the TRD mother volume
     track.~AliTRDtrackV1();
@@ -350,12 +365,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
@@ -369,12 +384,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
@@ -611,6 +624,24 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
     Double_t cov[3]; tracklet->GetCovAt(x, cov);
     Double_t p[2] = { tracklet->GetY(), tracklet->GetZ()};
     Double_t chi2 = ((AliExternalTrackParam)t).GetPredictedChi2(p, cov);
+    if(fkReconstructor->IsDebugStreaming()){
+      Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+      TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
+      AliExternalTrackParam param0(t);
+      AliExternalTrackParam param1(t);
+      param1.Update(p, cov);
+      TVectorD vcov(3,cov);
+      TVectorD vpar(3,p);
+      cstreamer << "FollowProlongationInfo"
+               << "EventNumber="       << eventNumber
+               << "iplane="<<iplane
+               << "vcov.="<<&vcov
+               << "vpar.="<<&vpar
+               << "tracklet.="      << tracklet
+               << "param0.="           << &param0
+               << "param1.="           << &param1
+               << "\n";
+    }
     if (chi2 < 1e+10 && ((AliExternalTrackParam&)t).Update(p, cov)){ 
       // Register info to track
       t.SetNumberOfClusters();
@@ -727,14 +758,14 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
 
   // Loop through the TRD layers
   TGeoHMatrix *matrix = NULL;
-  Double_t x, y, z;
+  Double_t x(0.), y(0.), z(0.);
   for (Int_t ily=startLayer, sm=-1, stk=-1, det=-1; ily < AliTRDgeometry::kNlayer; ily++) {
     AliDebug(2, Form("Propagate to x[%d] = %7.2f", ily, fR[ily]));
 
     // 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;
     }
@@ -745,7 +776,6 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     stk = fGeom->GetStack(z, ily);
     det = stk>=0 ? AliTRDgeometry::GetDetector(ily, stk, sm) : -1;
     matrix = det>=0 ? fGeom->GetClusterMatrix(det) : NULL;
-    AliDebug(3, Form("Propagate to det[%3d]", det));
 
     // check if supermodule/chamber is installed
     if( !fGeom->GetSMstatus(sm) ||
@@ -756,23 +786,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;
     }
 
@@ -780,24 +810,25 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     Double_t loc[] = {AliTRDgeometry::AnodePos()- driftLength, 0., 0.};
     Double_t glb[] = {0., 0., 0.};
     matrix->LocalToMaster(loc, glb);
+    AliDebug(3, Form("Propagate to det[%3d] x_anode[%7.2f] (%f %f)", det, glb[0]+driftLength, glb[1], glb[2]));
 
     // Propagate to the radial distance of the current layer
     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;
     }
@@ -813,7 +844,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;
       }
@@ -824,44 +855,39 @@ 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;
     }
-    // mark track as entering the FIDUCIAL volume of TRD
-    if(kStoreIn){
-      t.SetTrackIn(); 
-      kStoreIn = kFALSE;
-    }
 
     ptrTracklet  = tracklets[ily];
     if(!ptrTracklet){ // BUILD TRACKLET
       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;
       }      
@@ -875,12 +901,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));
@@ -894,7 +924,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));
@@ -915,26 +945,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;
     }
@@ -943,7 +973,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);
@@ -957,11 +987,36 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       }
       AliDebug(4, Form("Failed Chi2[%f]", chi2));
       continue; 
+    }
+     if(fkReconstructor->IsDebugStreaming()){
+      Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+      //      TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
+      AliExternalTrackParam param0(t);
+      AliExternalTrackParam param1(t);
+      param1.Update(p, cov);
+      TVectorD vcov(3,cov);
+      TVectorD vpar(3,p);
+      (*cstreamer) << "FollowBackProlongationInfo"
+               << "EventNumber="       << eventNumber
+               << "chi2="<<chi2
+               << "iplane="<<ily
+               << "vcov.="<<&vcov
+               << "vpar.="<<&vpar
+               << "tracklet.="      << ptrTracklet
+               << "param0.="           << &param0
+               << "param1.="           << &param1
+               << "\n";
+    }
+
+    // mark track as entering the FIDUCIAL volume of TRD
+    if(kStoreIn){
+      t.SetTrackIn();
+      kStoreIn = kFALSE;
     }
     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);
@@ -983,7 +1038,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);
@@ -1155,7 +1213,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);
 
@@ -1315,13 +1376,13 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
   covarPolZ.Invert();
 
   // Update the tracklets
-  Double_t x1, dy, dz;
+  Double_t dy, dz;
   Double_t cov[15];
   memset(cov, 0, sizeof(Double_t) * 15);
   for(Int_t iLayer = 0; iLayer < AliTRDtrackerV1::kNPlanes; iLayer++) {
 
     x  = tracklets[iLayer].GetX0();
-    x1 = x - xref;
+//    x1 = x - xref;
     y  = 0;
     z  = 0;
     dy = 0;
@@ -1803,13 +1864,14 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
   //
 
   // Current track X-position
-  Double_t xpos = t.GetX();
+  Double_t xpos = t.GetX()/*,
+           mass = t.GetMass()*/;
 
   // Direction: inward or outward
   Double_t dir  = (xpos < xToGo) ? 1.0 : -1.0;
 
   while (((xToGo - xpos) * dir) > AliTRDReconstructor::GetEpsilon()) {
-
+//    printf("to go %f\n", (xToGo - xpos) * dir);
     Double_t xyz0[3];
     Double_t xyz1[3];
     Double_t param[7];
@@ -1837,10 +1899,45 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
     // Calculate the mean material budget between start and
     // end point of this prolongation step
     if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) return 0;
-
+    
     // 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
+    Double_t dEdx(0.),
+             bg(TMath::Abs(t.GetP()/mass));
+    if(AliLog::GetDebugLevel("TRD", "AliTRDtrackerV1")>=3){
+      const char *pn[] = {"rho", "x/X0", "<A>", "<Z>", "L", "<Z/A>", "Nb"};
+      printf("D-AliTRDtrackerV1::PropagateTo(): x[%6.2f] bg[%6.2f]\n", xpos, bg);
+      printf("     param :: %s[%e] %s[%e] %s[%e] %s[%e] %s[%e] %s[%e] %s[%e]\n"
+          , pn[0], param[0]
+          , pn[1], param[1]
+          , pn[2], param[2]
+          , pn[3], param[3]
+          , pn[4], param[4]
+          , pn[5], param[5]
+          , pn[6], param[6]);
+    }  
+    switch(fgBB){
+    case kSolid:
+      dEdx = AliExternalTrackParam::BetheBlochSolid(bg);
+      break;
+    case kGas:
+      dEdx = AliExternalTrackParam::BetheBlochGas(bg);
+      break;
+    case kGeant:
+      { // mean exitation energy (GeV)
+        Double_t mee = ((param[3] < 13.) ? (12. * param[3] + 7.) : (9.76 * param[3] + 58.8 * TMath::Power(param[3],-0.19))) * 1.e-9;
+        Double_t mZA = param[5]>1.e-5?param[5]:(param[3]/param[2]);
+        if(AliLog::GetDebugLevel("TRD", "AliTRDtrackerV1")>=3) printf("D-AliTRDtrackerV1::PropagateTo(): Mee[%e] <Z/A>[%e]\n", mee, mZA);
+        // protect against failed calculation of rho in MeanMaterialBudget()
+        dEdx = AliExternalTrackParam::BetheBlochGeant(bg, param[0]>1.e-6?param[0]:2.33, 0.2, 3., mee, mZA);
+      }
+      break;
+    }
+    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;
 
@@ -1853,7 +1950,6 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
 
 }
 
-
 //_____________________________________________________________________________
 Bool_t AliTRDtrackerV1::ReadClusters(TTree *clusterTree)
 {
@@ -1872,12 +1968,14 @@ Bool_t AliTRDtrackerV1::ReadClusters(TTree *clusterTree)
     return kFALSE;
   }
   branch->SetAddress(&clusterArray); 
-  
+
   if(!fClusters){ 
     Float_t nclusters =  fkRecoParam->GetNClusters();
     if(fkReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector;
     fClusters = new TClonesArray("AliTRDcluster", Int_t(nclusters));
     fClusters->SetOwner(kTRUE);
+    SetClustersOwner();
+    AliInfo(Form("Tracker owning clusters @ %p", (void*)fClusters));
   }
   
   // Loop through all entries in the tree
@@ -1911,27 +2009,26 @@ Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree)
   
   fkRecoParam = fkReconstructor->GetRecoParam(); // load reco param for this event
 
-  if(!fkReconstructor->IsWritingClusters()){ 
-    fClusters = AliTRDReconstructor::GetClusters();
+//  if(!fkReconstructor->IsWritingClusters()) AliInfo(Form("IsWritingClusters[%c]", fkReconstructor->IsWritingClusters()?'y':'n'));
+  if(!(fClusters = AliTRDReconstructor::GetClusters())){
+    AliWarning("Clusters unavailable from TRD reconstructor. Trying reading from tree ...");
   } else {
     if(!ReadClusters(cTree)) {
-      AliError("Problem with reading the clusters !");
+      AliError("Reading clusters from tree failed.");
       return 1;
     }
   }
-  SetClustersOwner();
 
   if(!fClusters || !fClusters->GetEntriesFast()){ 
     AliInfo("No TRD clusters");
     return 1;
-  }
+  } else AliInfo(Form("Using :: clusters[%d] onl.tracklets[%d] onl.tracks[%d]",
+    fClusters?fClusters->GetEntriesFast():0,
+    AliTRDReconstructor::GetTracklets()?AliTRDReconstructor::GetTracklets()->GetEntriesFast():0,
+    AliTRDReconstructor::GetTracks()?AliTRDReconstructor::GetTracks()->GetEntriesFast():0));
 
-  //Int_t nin = 
   BuildTrackingContainers();  
 
-  //Int_t ncl  = fClusters->GetEntriesFast();
-  //AliInfo(Form("Clusters %d [%6.2f %% in the active volume]", ncl, 100.*float(nin)/ncl));
-
   return 0;
 }
 
@@ -1945,17 +2042,14 @@ Int_t AliTRDtrackerV1::LoadClusters(TClonesArray * const clusters)
   if(!clusters || !clusters->GetEntriesFast()){ 
     AliInfo("No TRD clusters");
     return 1;
-  }
+  } else AliInfo(Form("Using :: external.clusters[%d]", clusters->GetEntriesFast()));
+
 
   fClusters = clusters;
-  SetClustersOwner();
 
   fkRecoParam = fkReconstructor->GetRecoParam(); // load reco param for this event
   BuildTrackingContainers();  
 
-  //Int_t ncl  = fClusters->GetEntriesFast();
-  //AliInfo(Form("Clusters %d [%6.2f %% in the active volume]", ncl, 100.*float(nin)/ncl));
-
   return 0;
 }
 
@@ -2004,14 +2098,15 @@ void AliTRDtrackerV1::UnloadClusters()
     fTracklets->Delete();
     if(HasRemoveContainers()){delete fTracklets; fTracklets = NULL;}
   }
-  if(fClusters){ 
-    if(IsClustersOwner()) fClusters->Delete();
-    
-    // save clusters array in the reconstructor for further use.
-    if(!fkReconstructor->IsWritingClusters()){
-      AliTRDReconstructor::SetClusters(fClusters);
-      SetClustersOwner(kFALSE);
-    } else AliTRDReconstructor::SetClusters(NULL);
+  if(fClusters && IsClustersOwner()){
+    AliInfo(Form("tracker[%p] clearing %d own clusters @ %p", (void*)this, fClusters->GetEntries(), (void*)fClusters));
+    fClusters->Delete();
+//     
+//     // save clusters array in the reconstructor for further use.
+//     if(!fkReconstructor->IsWritingClusters()){
+//       AliTRDReconstructor::SetClusters(fClusters);
+//       SetClustersOwner(kFALSE);
+//     } else AliTRDReconstructor::SetClusters(NULL);
   }
 
   for (int i = 0; i < AliTRDgeometry::kNsector; i++) fTrSec[i].Clear();
@@ -2066,7 +2161,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
@@ -2587,7 +2682,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 * c
   Double_t chi2[4];
 
   // this should be data member of AliTRDtrack TODO
-  Double_t seedQuality[kMaxTracksStack];
+//  Double_t seedQuality[kMaxTracksStack];
   
   // unpack control parameters
   Int_t config  = ipar[0];
@@ -2811,8 +2906,8 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 * c
         }
         //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
       
-        // book preliminary results
-        seedQuality[ntracks] = like;
+        // book preliminry results
+        //seedQuality[ntracks] = like;
         fSeedLayer[ntracks]  = config;/*sLayer;*/
       
         // attach clusters to the extrapolation seeds
@@ -3115,7 +3210,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++) {
@@ -4042,7 +4139,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