trigger names changed (Markus)
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDresolution.cxx
index 085705d..8fb2acb 100644 (file)
@@ -62,6 +62,7 @@
 #include <TLegend.h>
 #include <TGraphErrors.h>
 #include <TGraphAsymmErrors.h>
+#include <TLinearFitter.h>
 #include <TMath.h>
 #include <TMatrixT.h>
 #include <TVectorT.h>
@@ -73,6 +74,7 @@
 #include "AliLog.h"
 #include "AliESDtrack.h"
 #include "AliMathBase.h"
+#include "AliTrackPointArray.h"
 
 #include "AliTRDresolution.h"
 #include "AliTRDgeometry.h"
@@ -105,6 +107,12 @@ Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
     ,"TRDout2MC"
     ,"TRD2MC"
 };
+Char_t const * AliTRDresolution::fgParticle[11]={
+  " p bar", " K -", " #pi -", " #mu -", " e -",
+  " No PID",
+  " e +", " #mu +", " #pi +", " K +", " p",
+};
+
 // Configure segmentation for y resolution/residuals
 Int_t const AliTRDresolution::fgkNresYsegm[3] = {
   AliTRDgeometry::kNsector
@@ -122,6 +130,7 @@ AliTRDresolution::AliTRDresolution()
   ,fIdxPlot(0)
   ,fIdxFrame(0)
   ,fPtThreshold(1.)
+  ,fDyRange(1.5)
   ,fDBPDG(NULL)
   ,fGraphS(NULL)
   ,fGraphM(NULL)
@@ -144,6 +153,7 @@ AliTRDresolution::AliTRDresolution(char* name)
   ,fIdxPlot(0)
   ,fIdxFrame(0)
   ,fPtThreshold(1.)
+  ,fDyRange(1.5)
   ,fDBPDG(NULL)
   ,fGraphS(NULL)
   ,fGraphM(NULL)
@@ -193,7 +203,9 @@ void AliTRDresolution::UserCreateOutputObjects()
 //________________________________________________________
 void AliTRDresolution::InitExchangeContainers()
 {
-  fCl = new TObjArray();
+// Init containers for subsequent tasks (AliTRDclusterResolution)
+
+  fCl = new TObjArray(200);
   fCl->SetOwner(kTRUE);
   fMCcl = new TObjArray();
   fMCcl->SetOwner(kTRUE);
@@ -216,13 +228,11 @@ void AliTRDresolution::UserExec(Option_t *opt)
 
   fCl->Delete();
   fMCcl->Delete();
-/*  fTrklt->Delete();
-  fMCtrklt->Delete();*/
   AliTRDrecoTask::UserExec(opt);
 }
 
 //________________________________________________________
-Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt)
+Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt) const
 {
 // Helper function to calculate pulls in the yz plane 
 // using proper tilt rotation
@@ -271,14 +281,15 @@ TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
     if(!fTracklet->IsOK()) continue;
     Float_t x0 = fTracklet->GetX0();
-    Float_t dq, dl;
+    Float_t dqdl, dl;
     for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
       if(!(c = fTracklet->GetClusters(itb))){ 
         if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
       }
-      dq = fTracklet->GetdQdl(itb, &dl);
+      dqdl = fTracklet->GetdQdl(itb, &dl);
+      if(dqdl<1.e-5) continue;
       dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
-      (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
+      (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dqdl);
     }
 
 //     if(!HasMCdata()) continue;
@@ -312,9 +323,64 @@ TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 
   Int_t sgm[3];
   Double_t covR[7], cov[3], dy[2], dz[2];
-  Float_t pt, x0, y0, z0, dydx, dzdx;
+  Float_t pt, x0, y0, z0, dydx, dzdx, tilt(0.);
   const AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
   AliTRDseedV1 *fTracklet(NULL);  TObjArray *clInfoArr(NULL);
+  // do LINEAR track refit if asked by the user
+  // it is the user responsibility to check if B=0T
+  Float_t param[10];  memset(param, 0, 10*sizeof(Float_t));
+  Int_t np(0), nrc(0); AliTrackPoint clusters[300];
+  if(HasTrackRefit()){
+    Bool_t kPrimary(kFALSE);
+    for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+      if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
+      if(!fTracklet->IsOK()) continue;
+      x0   = fTracklet->GetX0();
+      tilt = fTracklet->GetTilt();
+      AliTRDcluster *c = NULL;
+      fTracklet->ResetClusterIter(kFALSE);
+      while((c = fTracklet->PrevCluster())){
+        Float_t xyz[3] = {c->GetX(), c->GetY(), c->GetZ()};
+        clusters[np].SetCharge(tilt);
+        clusters[np].SetClusterType(0);
+        clusters[np].SetVolumeID(ily);
+        clusters[np].SetXYZ(xyz);
+        np++;
+      }
+      if(fTracklet->IsRowCross()){
+        Float_t xcross(0.), zcross(0.);
+        if(fTracklet->GetEstimatedCrossPoint(xcross, zcross)){
+          clusters[np].SetCharge(tilt);
+          clusters[np].SetClusterType(1);
+          clusters[np].SetVolumeID(ily);
+          clusters[np].SetXYZ(xcross, 0., zcross);
+          np++;
+          nrc++;
+        }
+      }
+      if(fTracklet->IsPrimary()) kPrimary = kTRUE;
+    }
+    if(kPrimary){
+      clusters[np].SetCharge(tilt);
+      clusters[np].SetClusterType(1);
+      clusters[np].SetVolumeID(-1);
+      clusters[np].SetXYZ(0., 0., 0.);
+      np++;
+    }
+    if(!FitTrack(np, clusters, param)) {
+      AliDebug(1, "Linear track Fit failed.");
+      return NULL;
+    }
+    if(HasTrackSelection()){
+      Bool_t kReject(kFALSE);
+      if(fkTrack->GetNumberOfTracklets() != AliTRDgeometry::kNlayer)  kReject = kTRUE; 
+      if(!kReject && !UseTrack(np, clusters, param)) kReject = kTRUE;
+      if(kReject){
+        AliDebug(1, "Reject track for residuals analysis.");
+        return NULL;
+      }
+    }
+  }
   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
     if(!fTracklet->IsOK()) continue;
@@ -325,13 +391,26 @@ TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
     
     // retrive the track angle with the chamber
-    y0   = fTracklet->GetYref(0);
-    z0   = fTracklet->GetZref(0);
-    dydx = fTracklet->GetYref(1);
-    dzdx = fTracklet->GetZref(1);
+    if(HasTrackRefit()){
+      Float_t par[3];
+      if(!FitTracklet(ily, np, clusters, param, par)) continue;
+      dydx = par[2];//param[3];
+      dzdx = param[4];
+      y0   = par[1] + dydx * (x0 - par[0]);//param[1] + dydx * (x0 - param[0]);
+      z0   = param[2] + dzdx * (x0 - param[0]);
+    } else {
+      y0   = fTracklet->GetYref(0);
+      z0   = fTracklet->GetZref(0);
+      dydx = fTracklet->GetYref(1);
+      dzdx = fTracklet->GetZref(1);
+    }
+    /*printf("RC[%c] Primary[%c]\n"
+           "  Fit : y0[%f] z0[%f] dydx[%f] dzdx[%f]\n"
+           "  Ref:  y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", fTracklet->IsRowCross()?'y':'n', fTracklet->IsPrimary()?'y':'n', y0, z0, dydx, dzdx
+           ,fTracklet->GetYref(0),fTracklet->GetZref(0),fTracklet->GetYref(1),fTracklet->GetZref(1));*/
+    tilt = fTracklet->GetTilt();
     fTracklet->GetCovRef(covR);
-    Double_t tilt(fTracklet->GetTilt())
-            ,t2(tilt*tilt)
+    Double_t t2(tilt*tilt)
             ,corr(1./(1. + t2))
             ,cost(TMath::Sqrt(corr));
     AliTRDcluster *c = NULL;
@@ -398,6 +477,7 @@ TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
     }
   }
   if(clInfoArr) delete clInfoArr;
+
   return (TH3S*)arr->At(0);
 }
 
@@ -487,8 +567,6 @@ TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
         << "\n";
     }
   }
-
-
   return (TH2I*)arr->At(0);
 }
 
@@ -877,14 +955,14 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
   Double_t covR[7]/*, cov[3]*/;
 
   if(DebugLevel()>=3){
-    TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
-    fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
+    TVectorD dX(12), dY(12), dZ(12), vPt(12), dPt(12), cCOV(12*15);
+    fkMC->PropagateKalman(&dX, &dY, &dZ, &vPt, &dPt, &cCOV);
     (*DebugStream()) << "MCkalman"
       << "pdg=" << pdg
       << "dx="  << &dX
       << "dy="  << &dY
       << "dz="  << &dZ
-      << "pt="  << &Pt
+      << "pt="  << &vPt
       << "dpt=" << &dPt
       << "cov=" << &cCOV
       << "\n";
@@ -1809,7 +1887,9 @@ void AliTRDresolution::MakeSummary()
 
   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 
-  if(!HasMCdata()){
+  if(!HasMCdata() ||
+    (!fGraphS->At(kMCcluster) || !fGraphM->At(kMCcluster) ||
+     !fGraphS->At(kMCtracklet) || !fGraphM->At(kMCtracklet))){
     delete cOut;
     return;
   }
@@ -1915,12 +1995,12 @@ void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
     Int_t bmax(h1.GetMaximumBin());
     Int_t jBinMin(1), jBinMax(100);
     for(Int_t ibin(bmax); ibin--;){
-      if(h1.GetBinContent(ibin)==0){
+      if(h1.GetBinContent(ibin)<1.){
         jBinMin=ibin; break;
       }
     }
     for(Int_t ibin(bmax); ibin++;){
-      if(h1.GetBinContent(ibin)==0){
+      if(h1.GetBinContent(ibin)<1.){
         jBinMax=ibin; break;
       }
     }
@@ -1936,6 +2016,8 @@ void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
 //________________________________________________________
 void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
 {
+// Core functionality for MakeSummary function.  
+
   h2->Reset();  
   Double_t x,y;
   TGraphErrors *g(NULL); TAxis *ax(h2->GetXaxis());
@@ -1950,33 +2032,33 @@ void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
 
 
 //________________________________________________________
-Char_t const *fgParticle[11]={
-  " p bar", " K -", " #pi -", " #mu -", " e -",
-  " No PID",
-  " e +", " #mu +", " #pi +", " K +", " p",
-};
-const Color_t fgColorS[11]={
-kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
-kGray,
-kRed, kViolet, kMagenta+1, kOrange-3, kOrange
-};
-const Color_t fgColorM[11]={
-kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
-kBlack,
-kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
-};
-const Marker_t fgMarker[11]={
-30, 30, 26, 25, 24,
-28,
-20, 21, 22, 29, 29
-};
 Bool_t AliTRDresolution::PostProcess()
 {
-  //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
+// Fit, Project, Combine, Extract values from the containers filled during execution
+
+  /*fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));*/
   if (!fContainer) {
     AliError("ERROR: list not available");
     return kFALSE;
   }
+
+  // define general behavior parameters
+  const Color_t fgColorS[11]={
+  kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
+  kGray,
+  kRed, kViolet, kMagenta+1, kOrange-3, kOrange
+  };
+  const Color_t fgColorM[11]={
+  kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
+  kBlack,
+  kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
+  };
+  const Marker_t fgMarker[11]={
+  30, 30, 26, 25, 24,
+  28,
+  20, 21, 22, 29, 29
+  };
+
   TGraph *gm= NULL, *gs= NULL;
   if(!fGraphS && !fGraphM){ 
     TObjArray *aM(NULL), *aS(NULL);
@@ -2202,7 +2284,9 @@ TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool
     Int_t nybins=fgkNresYsegm[fSegmentLevel];
     if(expand) nybins*=2;
     h = new TH3S(hname, htitle, 
-                 48, -.48, .48, 60, -1.5, 1.5, nybins, -0.5, nybins-0.5);
+                 48, -.48, .48,            // phi
+                 60, -fDyRange, fDyRange,  // dy
+                 nybins, -0.5, nybins-0.5);// segment
   } else h->Reset();
   arr->AddAt(h, 0);
   sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
@@ -2294,11 +2378,11 @@ TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
   const Int_t kNpt(14);
   const Int_t kNdpt(150); 
   const Int_t kNspc = 2*AliPID::kSPECIES+1;
-  Float_t Pt=0.1, DPt=-.1, Spc=-5.5;
+  Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
-  for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
-  for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
-  for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
+  for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
+  for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
+  for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
 
   // Pt resolution
   sprintf(hname, "%s_%s_Pt", GetNameId(), name);
@@ -2351,11 +2435,11 @@ TObjArray* AliTRDresolution::Histos()
 
   // binnings for plots containing momentum or pt
   const Int_t kNpt(14), kNphi(48), kNdy(60);
-  Float_t Phi=-.48, Dy=-.3, Pt=0.1;
+  Float_t lPhi=-.48, lDy=-.3, lPt=0.1;
   Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
-  for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
-  for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
-  for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
+  for(Int_t i=0; i<kNphi+1; i++,lPhi+=.02) binsPhi[i]=lPhi;
+  for(Int_t i=0; i<kNdy+1; i++,lDy+=.01) binsDy[i]=lDy;
+  for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
 
   // cluster to track residuals/pulls
   fContainer->AddAt(arr = new TObjArray(2), kCharge);
@@ -2426,6 +2510,53 @@ Bool_t AliTRDresolution::Load(const Char_t *file, const Char_t *dir)
   return kTRUE;
 }
 
+//________________________________________________________
+Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
+{
+// Generic function to process sigma/mean for 2D plot dy(x)
+
+  if(!h2) {
+    if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
+    return kFALSE;
+  }
+  if(!Int_t(h2->GetEntries())){
+    if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
+    return kFALSE;
+  }
+  if(!g || !g[0]|| !g[1]) {
+    if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
+    return kFALSE;
+  }
+  // prepare
+  TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
+  TF1 f("f", "gaus", ay->GetXmin(), ay->GetXmax());
+  Int_t n(0);
+  if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
+  if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
+  TH1D *h(NULL);
+  if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
+
+  // do actual loop
+  for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
+    Double_t x = ax->GetBinCenter(ix);
+    Double_t ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
+    h = h2->ProjectionY("py", ix, ix);
+    if((n=(Int_t)h->GetEntries())<stat){
+      if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
+      continue;
+    }
+    f.SetParameter(1, 0.);
+    f.SetParameter(2, 3.e-2);
+    h->Fit(&f, "QN");
+    g[0]->SetPoint(np, x, f.GetParameter(1));
+    g[0]->SetPointError(np, ex, f.GetParError(1));
+    g[1]->SetPoint(np, x, f.GetParameter(2));
+    g[1]->SetPointError(np, ex, f.GetParError(2));
+    np++;
+  }
+  return kTRUE;
+}
+
 
 //________________________________________________________
 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
@@ -2911,6 +3042,151 @@ Bool_t AliTRDresolution::GetGraphArray(Float_t *bb, ETRDresolutionPlot ip, Int_t
   return kTRUE;
 }
 
+//____________________________________________________________________
+Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
+{
+//
+// Fit track with a staight line using the "np" clusters stored in the array "points".
+// The following particularities are stored in the clusters from points:
+//   1. pad tilt as cluster charge
+//   2. pad row cross or vertex constrain as fake cluster with cluster type 1
+// The parameters of the straight line fit are stored in the array "param" in the following order :
+//     param[0] - x0 reference radial position
+//     param[1] - y0 reference r-phi position @ x0
+//     param[2] - z0 reference z position @ x0
+//     param[3] - slope dy/dx
+//     param[4] - slope dz/dx
+//
+// Attention :
+// Function should be used to refit tracks for B=0T
+//
+
+  if(np<40){
+    if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].", np);
+    return kFALSE;
+  }
+  TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
+
+  Double_t x0(0.);
+  for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
+  x0/=Float_t(np);
+
+  Double_t x, y, z, dx, tilt(0.);
+  for(Int_t ip(0); ip<np; ip++){
+    x = points[ip].GetX(); z = points[ip].GetZ();
+    dx = x - x0;
+    zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
+  }
+  if(zfitter.Eval() != 0) return kFALSE;
+
+  Double_t z0    = zfitter.GetParameter(0);
+  Double_t dzdx  = zfitter.GetParameter(1);
+  for(Int_t ip(0); ip<np; ip++){
+    if(points[ip].GetClusterType()) continue;
+    x    = points[ip].GetX();
+    dx   = x - x0;
+    y    = points[ip].GetY();
+    z    = points[ip].GetZ();
+    tilt = points[ip].GetCharge();
+    y -= tilt*(-dzdx*dx + z - z0);
+    Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
+    yfitter.AddPoint(&dx, y, 1.);
+  }
+  if(yfitter.Eval() != 0) return kFALSE;
+  Double_t y0   = yfitter.GetParameter(0);
+  Double_t dydx = yfitter.GetParameter(1);
+
+  param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
+  if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", x0, y0, z0, dydx, dzdx);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, AliTrackPoint *points, const Float_t param[10], Float_t par[3])
+{
+//
+// Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
+// See function FitTrack for the data stored in the "clusters" array
+
+// The parameters of the straight line fit are stored in the array "param" in the following order :
+//     par[0] - x0 reference radial position
+//     par[1] - y0 reference r-phi position @ x0
+//     par[2] - slope dy/dx
+//
+// Attention :
+// Function should be used to refit tracks for B=0T
+//
+
+  TLinearFitter yfitter(2, "pol1");
+
+  // grep data for tracklet
+  Double_t x0(0.), x[60], y[60], dy[60];
+  Int_t nly(0);
+  for(Int_t ip(0); ip<np; ip++){
+    if(points[ip].GetClusterType()) continue;
+    if(points[ip].GetVolumeID() != ly) continue;
+    Float_t xt(points[ip].GetX())
+           ,yt(param[1] + param[3] * (xt - param[0]));
+    x[nly] = xt;
+    y[nly] = points[ip].GetY();
+    dy[nly]= y[nly]-yt;
+    x0    += xt;
+    nly++;
+  }
+  if(nly<10){
+    if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].", nly);
+    return kFALSE;
+  }
+  // set radial reference for fit
+  x0 /= Float_t(nly);
+
+  // find tracklet core
+  Double_t mean(0.), sig(1.e3);
+  AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
+
+  // simple cluster error parameterization
+  Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
+
+  // fit tracklet core
+  for(Int_t jly(0); jly<nly; jly++){
+    if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
+    Double_t dx(x[jly]-x0);
+    yfitter.AddPoint(&dx, y[jly], 1.);
+  }
+  if(yfitter.Eval() != 0) return kFALSE;
+  par[0] = x0;
+  par[1] = yfitter.GetParameter(0);
+  par[2] = yfitter.GetParameter(1);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDresolution::UseTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
+{
+//
+// Global selection mechanism of tracksbased on cluster to fit residuals
+// The parameters are the same as used ni function FitTrack().
+
+  const Float_t kS(0.6), kM(0.2);
+  TH1S h("h1", "", 100, -5.*kS, 5.*kS);
+  Float_t dy, dz, s, m;
+  for(Int_t ip(0); ip<np; ip++){
+    if(points[ip].GetClusterType()) continue;
+    Float_t x0(points[ip].GetX())
+          ,y0(param[1] + param[3] * (x0 - param[0]))
+          ,z0(param[2] + param[4] * (x0 - param[0]));
+    dy=points[ip].GetY() - y0; h.Fill(dy);
+    dz=points[ip].GetZ() - z0;
+  }
+  TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
+  fg.SetParameter(1, 0.);
+  fg.SetParameter(2, 2.e-2);
+  h.Fit(&fg, "QN");
+  m=fg.GetParameter(1); s=fg.GetParameter(2);
+  if(s>kS || TMath::Abs(m)>kM) return kFALSE;
+  return kTRUE;
+}
+
 //________________________________________________________
 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
 {