]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDtrackerV1.cxx
Change THnSparseD to THnF in truncated mean related calibration code to reduce the...
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerV1.cxx
index 6afdb846085661776859166e4fc839bec6df7c3e..7f166c88a568d057275f8afcd0cb0d6630e8286a 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)
@@ -350,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
@@ -369,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
@@ -727,14 +726,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 +744,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 +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;
     }
 
@@ -780,24 +778,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 +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;
       }
@@ -824,44 +823,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 +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));
@@ -894,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));
@@ -915,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;
     }
@@ -943,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);
@@ -958,10 +956,15 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       AliDebug(4, Form("Failed Chi2[%f]", chi2));
       continue; 
     }
+    // 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 +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);
@@ -1155,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);
 
@@ -1803,13 +1812,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 +1847,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(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 +1898,6 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
 
 }
 
-
 //_____________________________________________________________________________
 Bool_t AliTRDtrackerV1::ReadClusters(TTree *clusterTree)
 {
@@ -2066,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
@@ -3115,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++) {
@@ -4042,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